© The Authors 2026

1. Introduction

The consensus holds that early-type galaxies (ETGs) are generally quiescent systems that have undergone early cessation of star formation, resulting in their current red optical colors and minimal reserves of cold gas and dust (van den Bergh 1999). Morphologically, ETGs appear symmetric and predominantly elliptical, with shapes ranging from near-circular to highly elongated (Sandage 2005). Classified within Hubble’s framework (Hubble 1926, 1936), ETGs comprise two main types, i.e., elliptical and lenticular (S0) galaxies. Elliptical galaxies are smooth triaxial ellipsoids devoid of spiral arms or nebulae. Their surface brightness follows de Vaucouleurs’ law (de Vaucouleurs 1948), featuring concentrated central peaks and sharp outward declines (Kormendy 1977). Lenticular galaxies retain disk-like structures but lack spiral arms, often exhibiting central bars or rings (e.g., Sandage et al. 1970; van den Bergh 1976; Graham & Guzmán 2003).

Integral field spectroscopy (IFS) enables the classification of dust-unaffected ETGs into “fast rotators (FRs)” and “slow rotators (SRs)” based on projected kinematic properties (Emsellem et al. 2004, 2007, 2011). This division originated from qualitative observations in the SAURON survey of 48 ETGs (Emsellem et al. 2004), which revealed a dichotomy: galaxies exhibited either regular rotation patterns or complex or irregular stellar velocity maps. They are defined as regular rotators (RR) or non-regular rotators (NRR). Quantification of kinematic regularity is achieved through harmonic term analysis (e.g., the higher-order term of the Fourier decomposition of velocity fields $\overline{k_5}$ (Krajnović et al. 2006, 2008); see details in Section 4), distinguishing RRs from NRRs (Krajnović et al. 2011; Cappellari 2016). The correspondence between RR and FR as well as NRR and SR led to the adoption of the spin parameter λ_R (Emsellem et al. 2011, 2007), with SRs defined by ${\lambda }_R<0.31\sqrt{\epsilon }$ within the effective radius R_e, where ε is the ellipticity. Fast rotators display high λ_R and ε, resembling disk galaxies, while SRs feature weak rotation (λ_R ≲ 0.15). Integral-field spectroscopy surveys such as SAURON (de Zeeuw et al. 2002), ATLAS^3D (Cappellari et al. 2011a), CALIFA (Sánchez et al. 2012), SAMI (Croom et al. 2012; Bryant et al. 2015), MASSIVE (Ma et al. 2014), and MaNGA (Bundy et al. 2015) provide spatially resolved kinematic maps, facilitating studies of ordered versus disordered motion via λ_R–ε diagnostics. However, the robustness of this classification has not been systematically examined in numerical simulations. And the correlation between $\overline{k_5}$ and other kinematic parameters, such as the relative importance of cylindrical rotation energy (Sales et al. 2010) widely used in simulations, is poorly quantified.

Recent cosmological hydrodynamical simulations, such as EAGLE (Schaye et al. 2015; Crain et al. 2015), IllustrisTNG (TNG hereafter, e.g., Nelson et al. 2018), Horizon-AGN (Dubois et al. 2014, 2016), New-Horizon (Dubois et al. 2021), Magneticum (Dolag et al. 2025), and SIMBA (Davé et al. 2019), have made remarkable progress in successfully reproducing realistic visual morphologies and other properties. This remarkable progress has enabled statistical studies of galaxies in a fully self-consistent cosmological environment. On the one hand, a systematic comparison between IFS results and kinematic properties in simulations is crucial for constraining the potential shortage of simulations. On the other hand, numerical simulations can be used to examine the kinematic models and classification of SRs and FRs. It must be pointed out that, in addition to reproducing fairly well the morphological properties observed, modern cosmological hydrodynamic simulations – including the TNG suite – demonstrate a certain degree of reliability in reproducing the observed kinematic properties. Santucci et al. (2024) first established that the EAGLE simulation successfully models the stellar mass fraction variation across orbital families as functions of stellar mass and spin parameter at z = 0, consistent with SAMI observations. This alignment extends to the correlation between specific angular momentum j_* and stellar mass M_*, where Lagos et al. (2017) reported reasonable agreement between EAGLE results and ATLAS^3D data. Similarly, analyses of the Magneticum Pathfinder simulations reveal highly congruent λ_R(R_e) and ε distributions for ETGs compared to the ATLAS^3D, CALIFA, SAMI, and SLUGGS surveys (Schulze et al. 2018). Further validation by Schulze et al. (2020) confirms that Magneticum’s stellar spin parameter profiles closely match observations from SLUGGS (Arnold et al. 2014; Forbes et al. 2017) and ePN.S (Pulsoni et al. 2018). Building on CALIFA’s kinematic "Hubble sequence" derived from stellar orbit-circularity distributions, Xu et al. (2019) analyzed z = 0 TNG galaxies to explore connections between stellar orbital compositions and fundamental galaxy properties. They found TNG100 broadly reproduces observed orbital-component fractions and their stellar-mass dependencies. Similarly, Zhang et al. (2025) further showed that TNG galaxies well match the luminosity fractions of cold, warm, and hot orbital components of CALIFA galaxies extracted via orbit-superposition Schwarzschild modeling (e.g., Schwarzschild 1979; Valluri et al. 2004; van den Bosch et al. 2008; Zhu et al. 2018b,a,c). Moreover, IllustrisTNG reproduces galaxies with unusual structural and kinematic properties, such as low surface brightness galaxies (Zhu et al. 2018d), and jellyfish galaxies (Yun et al. 2019). Therefore, to a certain extent, TNG simulations can be employed for further research on the kinematic parameters and properties of galaxies.

Structural mass ratios derived from kinematic decomposition and cylindrical rotational properties can characterize galaxy dynamical properties and their underlying formation histories. However, observational measurement of spins, κ_rot, and precise structural mass ratios are challenging. For simulations, kinematic phase-space analysis enables robust decomposition: the auto-GMM method, an unsupervised algorithm developed by Du et al. (2019), identifies structures such as cold disks, warm disks, bulges, and halos in TNG galaxies by clustering stellar particles in a 3D kinematic phase space of normalized azimuthal angular momentum (circularity), non-azimuthal angular momentum, and binding energy. Galaxies with minimal stellar halos exhibit tight scaling relations among j_*, M_*, size, and metallicity owing to their quiescent merger histories (Du et al. 2022, 2024; Ma et al. 2024). Despite the auto-GMM’s success in isolating intrinsic structures formed by distinct mechanisms – such as separating stellar halos from bulges (Du et al. 2020, 2021) – the kinematic decomposition cannot be directly applied to observational datasets. Consequently, linking kinematically derived structures as well as galactic spin to IFS observables is critical. Moreover, while κ_rot = 0.5 or a spheroid mass ratio of 0.5 serves as a common threshold for classifying elliptical versus disk galaxies (e.g., Zhao et al. 2020; Fall & Rodriguez-Gomez 2023; Pérez-Montaño et al. 2024; Kataria & Vivek 2024), its relationship to kinematically distinct FRs and SRs in IFS surveys remains unclear.

In this paper, we perform a comparative analysis of galaxy kinematics using TNG simulations and IFS data from the MaNGA and ATLAS^3D surveys, with a focus on classifying ETGs into FRs and SRs, as well as interpreting the definitions of FRs and SRs. The structure of this paper is organized as follows. Section 2 describes the ETGs sample selection and data extraction from IFS observations. Section 3 describes mock observations and sample selection of TNG simulation galaxies. Section 4 demonstrates the kinematic parameters quantifying rotation properties of galaxies. Section 5 investigates the relation between different kinematic parameters. Section 6 examines the bimodality and the classification of SRs and FRs in kinematics. Section 7 compares the radial distribution of kinematic properties between TNG simulations and observations. Our conclusions are discussed and summarized in Sections 8 and 9, respectively.

2. Sample selection and data extraction from IFS observations

2.1. Data extraction of color, morphology, and rotation

This study compares the kinematic properties of simulated TNG galaxies with IFS measurements from the ATLAS^3D and MaNGA surveys for ETGs. ATLAS^3D1 is a complete, multiwavelength, volume-limited survey of nearby ETGs, and its final sample encompasses 260 nearby ETGs. The Sloan Digital Sky Survey-IV (SDSS-IV: York et al. 2000; Gunn et al. 2006; Blanton et al. 2017) Mapping Nearby Galaxies at Apache Point Observatory (MaNGA², Bundy et al. 2015; Yan et al. 2016) is an IFS survey designed to obtain spectral measurements across the surfaces of approximately 10 000 nearby galaxies. The data of color, morphology, and stellar mass of our sample are extracted as follows.

Color g − r − The g − r colors for MaNGA galaxies and most ATLAS^3D galaxies were obtained from the NASA-Sloan Atlas³ (NSA) catalog (Blanton & Roweis 2007; Blanton et al. 2011). For the subset of ATLAS^3D galaxies not included in the NSA catalog, we used B − V colors corrected for Galactic extinction from the Hyperleda⁴ catalog (Makarov et al. 2014). These B − V colors were then converted to g − r colors using the relation provided in Pulsoni et al. (2020).

Size R_e − The effective radius R_e of the ATLAS^3D sample was adopted from the Table 1 in Cappellari et al. (2013b). The data of MaNGA is from Zhu et al. (2023). MaNGA provides spatially resolved spectra, covering radial ranges out to 1.5R_e for the Primary+ sample (about two-thirds of the total sample) and out to 2.5R_e for the secondary sample (about one-third of the total sample).

Stellar mass M_* − The stellar masses for the ATLAS^3D and MaNGA sample were derived using a color–mass relation assuming the Chabrier (2003) initial mass function. For the ATLAS^3D dataset, total absolute K-band luminosities (L_K) were adopted directly from Cappellari et al. (2011a). For the MaNGA dataset, L_K values originate from Graham et al. (2018). Both datasets use luminosities sourced from the 2MASS extended source catalog (Jarrett et al. 2003), with corrections applied for Galactic extinction. These luminosities were further calibrated to address sky-background over-subtraction artifacts using L_Kcorr = 1.07L_K + 1.53 as established by Scott et al. (2013). Finally, stellar masses were calculated via log₁₀(M_*/M_⊙) = 10.39 − 0.46(L_Kcorr + 23) (van de Sande et al. 2019), where M_* is consistent with that derived by the Jeans anisotropic model (JAM) dynamical methodology (Cappellari et al. 2013a; Cappellari 2013).

Ellipticity ε − For ATLAS^3D galaxies, ellipticities within 1R_e are from Emsellem et al. (2011) Table B1. Among the 260 galaxies analyzed, 17 with prominent bars have ellipticity measurements specifically at 2.5–3R_e to avoid biased inner regions (Emsellem et al. 2011; Pulsoni et al. 2020). The ellipticities of MaNGA galaxies were obtained from Zhu et al. (2023), which are uniformly measured at 1R_e.

Spin spectroscopic parameter λ_R − The spin spectroscopic parameter λ_R was calculated across different surveys using slightly varying integration areas. In the ATLAS^3D survey, Emsellem et al. (2011) used circular apertures with a radius of 1R_e, whereas the MaNGA survey integrates over a half-light ellipse (Zhu et al. 2023). This ellipse covers an area equivalent to a circle of radius 1R_e, with a semimajor axis of $R_{\mathrm{e}}/\sqrt{1-\epsilon }$. Meanwhile, we adopted elliptical apertures for all λ_R calculations in our study.

V_rot/σ profiles − The V_rot/σ profiles for both the ATLAS^3D and MaNGA galaxies were calculated as the ratio of the rotation velocity V_rot(R) (see Equation (3)) and the azimuthally averaged velocity dispersion σ(R) at a given in elliptical radial bins (Pulsoni et al. 2020). For the ATLAS^3D sample, kinematic maps were adopted from Emsellem et al. (2004) and Cappellari et al. (2011a). For the MaNGA galaxies, publicly available stellar kinematics from the SDSS DR17⁵ (Abdurro’uf et al. 2022) were used. We then applied KINEMETRY⁶ to fit the velocity and velocity dispersion fields separately, deriving radial profiles of V_rot and σ for these galaxies.

Inclination angle i − Inclination angles of ATLAS^3D and MaNGA galaxies were adopted from Cappellari et al. (2013b) and Zhu et al. (2023), respectively.

2.2. Sample selection of early-type galaxies

The ATLAS^3D survey selected a sample of bright ETGs with high completeness within a distance of < 42 Mpc (z ≲ 0.01) and a sky declination δ satisfying |δ − 29° | < 35°. Additionally, these galaxies required an absolute K-band magnitude brighter than −21.5 mag (M_* ≳ 6 × 10⁹ M_⊙), and away from the dusty region near the Galaxy equatorial plane with latitude |b|< 15°. Early-type galaxies were identified morphologically by excluding galaxies with visible spiral structures, resulting in a sample similar to galaxies on the red sequence (Cappellari et al. 2011a). The ATLAS^3D ETG sample comprises 68 elliptical galaxies and 192 lenticular galaxies. The final sample comprises 260 nearby ETGs, with their kinematic data acquired using the SAURON IFS instrument (Bacon et al. 2001) mounted on the William Herschel Telescope. The SAURON IFS has a field of view of approximately 33″ × 41″ and is equipped with a square microlens array of 0.94″ × 0.94″. For the ATLAS^3D galaxies, Cappellari et al. (2011a) employed a uniform analysis pipeline to analyze the spectral data cubes, adopting a minimum signal-to-noise (S/N) threshold of 40 for adaptive binning (Cappellari & Copin 2003). Stellar kinematics were extracted using the ppxf⁷ software (Cappellari & Emsellem 2004) in conjunction with the MILES⁸ stellar template library (Sánchez-Blázquez et al. 2006).

The MaNGA survey targets galaxies within a redshift range of 0.01 < z < 0.15 and a stellar mass range of 10⁹ M_⊙ to 10¹² M_⊙, achieving a spatial resolution of 2.5″. For the MaNGA sample, stellar kinematics were derived from spectral data cubes using the Data Analysis Pipeline (DAP; Belfiore et al. 2019; Westfall et al. 2019) with a subset of the MILES stellar library (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011). Prior to kinematic extraction, the spectra were Voronoi-binned (Cappellari & Copin 2003) to a target S/N of 10 to ensure reliable stellar velocity dispersion measurements. The final SDSS DR17 includes over 10 000 galaxy data cubes after the removal of non-galaxy targets. We only selected data cubes with good quality⁹ (drp3qual = 1), as determined by MaNGA’s Data Reduction Pipeline (DRP; Law et al. 2016; Yan et al. 2016), and imposed additional criteria: redshift z < 0.05 (the sample is considered to be approximately complete) and stellar masses between 10¹⁰ M_⊙ and 10^11.5 M_⊙. This yields a final sample of 3100 MaNGA galaxies comparable with the TNG sample.

We selected ETGs using the criterion (g − r)≥0.05log₁₀(M_*/M_⊙)+0.1 mag (Pulsoni et al. 2020) and limited the sample stellar mass range to 10¹⁰ ≤ M_* ≤ 10^11.5 M_⊙, yielding 193 ATLAS^3D (red points) ETGs, shown in the upper panel of Figure 1. Relying solely on color for classification proved inadequate for MaNGA galaxies (blue points), as dust reddening in spiral galaxies can significantly alter their observed colors. To address this limitation, we supplemented g − r color cuts with morphological classifications from the MaNGA Deep Learning catalog¹⁰, which utilizes a convolutional neural network trained on data from Nair & Abraham (2010). Therefore, to select ETGs from the MaNGA data, two criteria must be met simultaneously: first, these galaxies should have relatively red colors; second, they should be classified as ETGs according to the Deep Learning catalog. Adopting the methodology of Domínguez Sánchez et al. (2022), we selected ETGs using the same criteria. This approach identifies 599 elliptical and 334 lenticular galaxies.

Fig. 1. g − r color–stellar mass diagram. Top: observational data (from MaNGA and ATLAS3D). Bottom: simulated data (from TNG50 and TNG100). Both panels encompass all samples within the mass range of 1010 − 1011.5 M⊙. Red sequence ETGs are selected above the dashed line.

3. Sample selection and mock observations of IllustrisTNG simulation galaxies

The TNG simulations use a Lambda cold dark matter (ΛCDM) setup with the cosmological parameters from the Planck survey (Planck Collaboration I 2016), and it exploits all the advantages of the unstructured moving-mesh hydrodynamic method Arepo¹¹ (Springel 2010) but improves the numerical methods, the subgrid physical model, and the recipe for galaxy feedback both from stars and Active Galactic Nucleus (AGNs). In particular, TNG is equipped with a novel dual mode (thermal and kinetic) AGN feedback that shapes and regulates the stellar component within massive systems, maintaining a realistic gas fraction (Weinberger et al. 2017). Also, the feedback model from galactic winds has been improved to better represent low- and intermediate-mass galaxies (Pillepich et al. 2018b).

The TNG suite comprises three runs using different simulation volumes and resolutions, namely TNG50, TNG100, and TNG300 (Nelson et al. 2018; Naiman et al. 2018; Marinacci et al. 2018; Pillepich et al. 2018a, 2019; Springel et al. 2018). Each run begins at z = 127 using the Zeldovich approximation and evolves down to z = 0. In this study, we used TNG50-1 (hereafter, simply TNG50) and TNG100-1 (hereafter, simply TNG100) runs. TNG50 provides a high number of galaxies for statistical analyses and a “zoom-in”-like resolution. TNG50 includes 2 × 2160³ initial resolution elements in a ∼50 comoving Mpc box, with the baryon mass resolution of 8.5 × 10⁴ M_⊙ and a gravitational softening radius r_soft for stars of about 0.3 kpc at z = 0. Dark matter (DM) was resolved with particles of mass 4.5 × 10⁵ M_⊙. Meanwhile, the minimum gas softening reaches 74 comoving parsecs. TNG100 includes approximately 2 × 1820³ resolution elements in a ∼110 comoving Mpc box. The DM and baryonic mass resolutions of TNG100 are m_DM = 7.5 × 10⁶ M_⊙ and m_b = 1.4 × 10⁶ M_⊙, respectively (Springel & Hernquist 2003). The softening length employed for TNG100 for both the DM and stellar components is ϵ = 0.74 kpc, while an adaptive gas gravitational softening is used, with a minimum ϵ_gas, min = 0.185 kpc. TNG50 thus has roughly 15 times better mass resolution, and 2.5 times better spatial resolution, than TNG100 (Nelson et al. 2019).

The identification of galaxies in the simulations was performed using the friends-of-friends (FoF) group finding algorithm (Davis et al. 1985) and the SUBFIND algorithm (Springel et al. 2001; Dolag et al. 2009). The central galaxy (subhalo) is the first (most massive) subhalo of each FoF group, while the other galaxies within the FoF group are its satellites. All components, including gas, stars, DM, and black holes, gravitationally bound to a single galaxy are associated with their host subhalo. Each galaxy was positioned at the minimum of its respective gravitational potential well. This study utilized galaxies from snapshot #99 (corresponding to z = 0) in both TNG50 and TNG100.

3.1. Sample selection of early-type galaxies

Nelson et al. (2018) showed that TNG simulations reproduce realistic galaxy colors matching SDSS observations (Strateva et al. 2001). Redder galaxies in TNG exhibit typical ETG-like properties: suppressed star formation, low gas fractions, higher metallicities, and older stellar populations. Firstly, we selected all galaxies with stellar masses within the range of 10^10.0 ≤ M_* ≤ 10^11.5 M_⊙, including both central and satellite galaxies. Then we identified ETGs using the color-mass diagram (Figure 1, lower panel), applying the same color selection criteria (dashed lines) and the same stellar mass range as for the observations. While TNG color measurements ignore dust effects, we validated our ETG selection using a deep learning classifier (Huertas-Company et al. 2019) trained on Nair & Abraham (2010) data¹². This morphological analysis confirmed that 85% of color-selected TNG100 galaxies are genuine ETGs (1693 E/S0 systems), with the remaining 15% being spiral galaxies. Tests show this contamination negligibly impacts our results, justifying our color-only selection approach for TNG samples.

In addition, we excluded galaxies where the effective radius R_e fell below twice the gravitational softening length (R_e < 2r_soft) where r_soft is 0.74 kpc and 0.288 kpc in TNG100 and TNG50, respectively. This criterion eliminated five objects (2.0% of the TNG50 parent sample; N = 244) and 245 objects (9.4% of the TNG100 parent sample; N = 2602). We further excluded merger-contaminated galaxies through visual inspection of stellar density maps. After applying these criteria, our final ETG sample consisted of 2342 and 236 galaxies from TNG100 and TNG50, respectively.

The ATLAS^3D parent sample, which includes 260 ETGs as well as 611 spiral galaxies at z ≲ 0.01, has been confirmed as a complete and representative sample of the nearby galaxy population, with an ETG fraction of approximately 30% (260/871) (Cappellari et al. 2011a). In the TNG simulations, our selection criteria yield ETG fractions of 27% (236/873) in TNG50 and 35% (2342/6369) in TNG100. By comparison, an observational sample (MaNGA) at z < 0.05 (Section 2.2) classifies 23% of galaxies as ETGs, a value closely aligned with TNG50. A separate local sample (0.05 < z < 0.1, g < 16 mag) reports a higher ETG fraction of 33% (14 034 galaxies) (Nair & Abraham 2010). TNG100 appears to reproduce a more consistent overall fraction of ETGs, while the discrepancy between TNG100 and TNG50 suggests that TNG50 may underestimate quenching efficiency. Meanwhile, in terms of spatial coverage and environmental representation, TNG50 broadly aligns with the ATLAS^3D survey in sampling low-density cosmological regimes (e.g., Cappellari et al. 2011b), while TNG100 corresponds to MaNGA’s cosmological volume (z < 0.05), which incorporates higher-density environments (e.g., Goddard et al. 2017; Greene et al. 2017). Shown in Figure 2, the stellar mass distributions of TNG50 and TNG100 galaxies show no clear differences attributable to their cosmological environments. And these environmental differences between simulation pairs do not fully account for the observational-simulation discrepancies we found (see details in Section 6.3).

Fig. 2. Mass–size relation. Top: stellar mass distribution function of the chosen sample of ETGs from observational data and the TNG simulations. Bottom: Circularized effective radii Re, of the selected ETG samples in the TNG100 and TNG50 simulations as a function of stellar mass and compared with observations (ATLAS3D and MaNGA). Both TNG50 (solid green line) and TNG100 (solid gray line) galaxies share comparable stellar masses and sizes with the ATLAS3D (dashed red line) and MaNGA (dashed blue line) galaxies. The lines and points represent the median values, with colored regions and error bars indicating the 16th to 84th percentile range of each distribution. To avoid overlap, we shift the error bars slightly for clarity.

3.2. Mock images and data extraction

3.2.1. Mass and size

We define the total stellar mass M_* as the total bound stellar mass of the galaxy. That is, the calculation of the total stellar mass utilizes the complete set of particle data associated with galaxies and subhalos, encompassing all gravitationally bound particles identified through the SUBFIND algorithm. The effective radius R_e is estimated by the circularized effective radius in a similar way to that of observations, where the area of the ellipse that encloses half of M_* is equated to $\pi R_{\mathrm{e}}^2$.

The multi-Gaussian expansion (MGE) fitting is conducted using the Python package MGEFIT¹³ (Cappellari 2002) on the 2D stellar mass map within an area of 60 kpc × 60 kpc. The stellar surface density in the MGE fitting is expressed as

$$\begin{array}{c}\hfill \mathrm{\Sigma }(x\prime ,y\prime )={\displaystyle \sum _{k=1}^N}\frac{{\mathrm{\Sigma }}_k}{2\pi {\sigma }_k^2{q}_k^{\prime }}exp[-\frac{1}{2{\sigma }_k^2}(x{\prime }^2+\frac{y{\prime }^2}{q_k^{\prime 2}})],\end{array}$$(1)

where Σ_k, σ_k, and $q_k\prime$ represent the total stellar mass, dispersion along the major axis, and axial ratio of the k-th Gaussian component, respectively. MGEFIT outputs the 2D half-mass radius R_e and semimajor axis $R_{\mathrm{e}}^{\mathrm{maj}}$ based on the 2D stellar mass distribution in the MGE formalism. It should be noted that the MGE approach does not extrapolate the light or density of a galaxy to infinite radii when determining R_e. Beyond three times the dispersion of the largest MGE Gaussian component, the model’s flux effectively drops to zero. As a result, this method yields a slightly smaller R_e compared to other techniques (de Vaucouleurs 1948; Burstein et al. 1987; de Vaucouleurs et al. 1991). The R_e and $R_e^{\mathrm{maj}}$ values used in this paper are thus scaled by a factor of 1.35, as suggested by Cappellari et al. (2013b), which subtracts potential systematic errors due to different methods.

We compare the R_e and M_* distributions of TNG galaxies with ATLAS^3D and MaNGA samples in Figure 2. The TNG galaxies (dashed profiles) agree with observed (solid profiles) R_e and M_* values well. Those with M_* ≥ 10^10.8 M_⊙ may be slightly higher (∼2 kpc), but this has minimal impact on our results. Moreover, we also examined the measurement using r−band luminosity images. Our analysis reveals that R_e measured by r-band fitting yields a slightly higher value by ≲1 kpc than those of mass-based measurements, propagating to a marginal increase in Δλ_R(R_e)≲0.03. Though quantifiable, these differences negligibly impact our kinematic conclusions.

3.2.2. Kinematic maps and their Fourier expansion models

We computed the line-of-sight (LOS) mean velocity and velocity dispersion fields for TNG galaxies from a randomly oriented perspective. All images are binned into 0.2 kpc intervals, equivalent to 2 arc seconds for a galaxy at a distance of 20 Mpc. This creates mock kinematic observations comparable to the ATLAS^3D survey (Cappellari et al. 2011a). Examples of fast (left panels) and slow (right panels) rotators are shown in Figure 3. Stellar particles are grouped into a regular grid centered on the galaxy, spanning a width of 40 kpc. Voronoi binning¹⁴ (Cappellari & Copin 2003) is then applied to merge bins with an insufficient number of stellar particles. Due to the Poissonian noise in particle counts, the target S/N relates to the number of particles per bin as $\mathrm{S}/\mathrm{N}=\sqrt{N_{\mathrm{p}}}$. To achieve robust kinematics, we set S/N ≥ 10, requiring at least 100 stellar particles per bin. The mean velocity and velocity dispersion of the i-th bin are calculated as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \overline{V_i}=\frac{{\sum }_n{M}_{n,i}{V}_{n,i}}{{\sum }_n{M}_{n,i}},{\sigma }_i=\sqrt{\frac{{\displaystyle \sum _n}{M}_{n,i}{(V_{n,i}-\overline{V_i})}^2}{\frac{N_i}{N_i-1}{\displaystyle \sum _n}{M}_{n,i}}},\end{array}\end{array}$$(2)

Fig. 3. Velocity maps and models of examples of a FR (left) and a SR (right), with raw data (top), a reconstructed velocity field employing a pure cosine law along ellipses (middle), and residuals depicting deviations (bottom). Our visualization emphasizes structural features through isodensity ellipses spaced at regular Re intervals from Re to 3Re. The major axis orientations – with solid lines representing photometrically derived morphological axes and dashed lines indicating kinematically determined axes – are both measured at Re.

where the index n runs over the stellar particles with mass M_n, i within ith Voronoi bin, N_i is the particle number of this bin.

We further applied the KINEMETRY package for modeling the velocity maps through a truncated Fourier expansion along elliptical contours (Krajnović et al. 2006, 2008; Foster et al. 2016):

$$\begin{array}{c}\hfill V(\psi )=V_0+{\displaystyle \sum _{n=1}^k}[A_{n}sin(n\psi )+B_{n}cos(n\psi )],\end{array}$$(3)

where ψ denotes the eccentric anomaly, V₀ represents the systemic velocity, and A_n, B_n are the Fourier coefficients. In the case of an ideal rotating disk, symmetry requirements force all coefficients except B₁ to vanish, simplifying the velocity profile to V(ψ) = V₀ + B₁cos(ψ). The cosine term captures the characteristic signature of regular rotation B₁ = V_rot (the second row of Figure 3), where V_rot denotes the velocity amplitude, while nonzero higher-order coefficients reveal departures from axisymmetry in the velocity field (the third row of Figure 3).

4. Parameters quantifying rotation properties of galaxies

The relative importance of rotation is a key factor that determines the morphology, assembly history, and internal dynamics of galaxies. Several dimensionless parameters are commonly employed to characterize the rotational properties of galaxies, such as λ_R (the classical spin parameter), κ_rot (cylindrical rotation energy), and structural mass ratios derived from kinematic decomposition. Meanwhile, the higher-order term $\overline{k_5}$ from the Fourier decomposition of velocity fields, which is used to describe kinematic anomalies, has been used to classify fast and SRs in ETGs (Emsellem et al. 2004, 2011), despite lacking a well-defined physical connection to the relative importance of rotation. $\overline{k_5}$ is defined as the mass-weighted average of the fifth-order kinematic Fourier coefficient $k_5=\sqrt{A_5^2+B_5^2}/V_{\mathrm{rot}}$ within the effective radius R_e, where V_rot represents the velocity amplitude, A₅ and B₅ are the fifth sine and cosine terms, respectively, in the truncated Fourier expansion for modeling the velocity field along elliptical contours (see Equation (3)). The relationships among these parameters and their consistency with observational data are not yet fully understood.

This study systematically estimates and compares these parameters through: (1) rigorous evaluation of their diagnostic capabilities, (2) quantitative analysis of their inter-parameter correlations, and (3) direct comparison with observational results. This comprehensive approach seeks to establish a universal framework for characterizing galactic rotation and classifying of slow and FRs of galaxies in both observations and simulations.

4.1. The λ_R − ε diagram

The spin parameter λ_R is a dimensionless parameter that quantifies the rotation of galaxies as defined by Emsellem et al. (2011)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\lambda }_R\equiv \frac{⟨R\mid \overline{V}⟩}{⟨R\sqrt{{\overline{V}}^2+{\sigma }^2}⟩}=\frac{{\sum }_i{M}_i{R}_i|\overline{V_i}|}{{\sum }_i{M}_i{R}_i\sqrt{{\overline{V_i}}^2+{\sigma }_i^2}},\end{array}\end{array}$$(4)

where the weighting with the flux is substituted here with a weighting with the mass contained within the ith Voronoi bin, M_i. The value R_i is the distance of the ith bin to the galaxy center. The values λ_R(R_e) and λ_R(3R_e) are computed by summing over all Voronoi bins within ellipses of one and three times R_e, with semimajor axis $R_{\mathrm{e}}^{\mathrm{maj}}$ (Graham et al. 2018)¹⁵.

Galaxies with higher rotational velocities typically exhibit larger spin parameters (λ_R) and lower ellipticities (ε), as illustrated in the λ_R − ε diagrams of Figures 4 and 5 for both TNG simulations and IFS observations.

Fig. 4. λR(Re)−ε(3Re) diagram color-coded by κrot(all), fspheroid(all), fhalo(all), and fbulge(all) for TNG50 (upper) and TNG100 (lower). Here, κrot(all), fspheroid(all), fhalo(all), and fbulge(all) represent the significance of cylindrical rotation, the spheroidal mass fraction, the stellar halo mass fraction, and the bulge mass fraction, respectively (see Sections 4.3 and 4.4). The magenta curve shows the λR − ε relation for edge-on, axisymmetric galaxies derived by Cappellari et al. (2007). The dotted black curves illustrate this relation across varying i, with steps of Δi = 10°. Additionally, the thin dashed black lines represent the theoretical distribution for the different rotations spaced at intervals of ΔλR, intr = 0.1. The bold black curve represents our proposed new threshold for distinguishing between SRs and FRs. This threshold is derived from the tensor virial theorem for oblate galaxies with intrinsic parameters λR, intr(Re) = 0.4 and εintr = 0.525, assuming an anisotropy of δ = 0.7εintr = 0.367. For comparison, the thin solid black line represents the standard empirical threshold defined by ${\lambda }_R(R_{\mathrm{e}})=0.31\times \sqrt{\epsilon }$ proposed in Emsellem et al. (2011, Eq. (3)). We present the proportion of SRs classified by our new threshold alongside the standard old threshold at the upper-left corner. We smoothed the color data using the locally weighted regression method LOESS (Cleveland & Devlin 1988; Cappellari et al. 2013a) with a smoothing factor of 0.1. The horizontal dashed blue line represents the empirical threshold defined by λR(Re) = 0.1 according to Emsellem et al. (2007); the black polygon denotes the empirical threshold set by λR(Re) < 0.08 + ε/4  and  ε < 0.4 as per Cappellari (2016); and the green polygon signifies the empirical threshold given by ${\lambda }_R(R_{\mathrm{e}})<{\lambda }_{R\mathrm{start}}+\epsilon /4\mathrm{and}\epsilon <0.35+\frac{{\lambda }_{R\mathrm{start}}}{1.538}$ (with λRstart = 0.12 when the IFS data quality is comparable to that of SAMI) based on van de Sande et al. (2021).

Our analysis confirms a substantial galaxy population in the high-ellipticity region (ε ≳ 0.4) with intermediate-to-low λ_R(R_e) values–occupying the right portion of the λ_R(R_e)–ε(R_e) plane. Pulsoni et al. (2020) argued that such galaxies (∼50% of total ETG samples) should be excluded as nonphysical via applying an axis-ratio cut at R_e as suggested by observations (Weijmans et al. 2014; Foster et al. 2017; Li et al. 2018; Ene et al. 2018). Excluding such galaxies may induce severe selection bias. While ellipticity (ε) is commonly used to quantify inclination of galaxies, we find that these specific systems are predominantly barred galaxies, whose central ε(R_e) are enhanced by the bar structure. There is no robust justification for excluding such galaxies, as the presence of elongated bar structures is unlikely to significantly alter overall galactic kinematics. To mitigate inclination errors from this central triaxiality of bars and the associated bias of excluding such galaxies, we adopted ε measured at 3R_e. Seventeen of the 260 galaxies from ATLAS^3D were measured at 2.5R_e − 3R_e, and the majority of observational data was measured at R_e. This issue should not significantly affect our conclusion for fast-rotating cases. Emsellem et al. (2011) introduced an empirical threshold ${\lambda }_R(R_{\mathrm{e}})=0.31\sqrt{\epsilon }$ (thin black curve in Figure 5) to distinguish FRs from SRs. This criterion is initially defined by the fifth-order kinematic Fourier coefficient $\overline{k_5}$, which quantifies non-axisymmetric rotational components. For instance, the bottom-right panel of Figure 3 showcases a SR with elevated $\overline{k_5}$ values (see details in Section 4). This classical criterion yields about 12% SRs in TNG simulations, which is consistent with the result of both ATLAS^3D and MaNGA.

Fig. 5. λR(Re)−ε(Re) diagram of the ATLAS3D(left panel) and MaNGA (right panel) samples. The proportions of SRs and FRs classified using our new threshold, compared to the standard old definition, are provided at the top. The definitions of all the lines are the same as those in Figure 4.

To assess the uncertainty in kinematic measurements at different galactocentric radii, we computed the deviation between λ_R values measured within R_e and 3R_e, as shown in Figure 6. The differences between morphologically derived position angles (PA_morph) and kinematically derived ones (PA_kin) are generally minor (i.e., < 20 degrees), as PA_morph and PA_phot can be derived by applying KINEMETRY to perform elliptical fitting on simulated galaxy images. In contrast, λ_R(R_e) values are systematically lower than λ_R(3R_e), which is expected due to the kinematic suppression often caused by central slow-rotating components such as bulges. This radial dependence highlights the importance of well-defined apertures in kinematic analyses, particularly for galaxies with structurally complex components. For consistency, we adopted 3R_e measurements to represent accurate galactic properties, while R_e measurements were retained for direct comparison with observational studies.

Fig. 6. Distributions of the differences between various galactic parameters measured within (at) 3Re and Re isodensity ellipses. From left to right: λR(3Re)−λR(Re), ε(3Re)−ε(Re), PAmorph(3Re)−PAmorph(Re), PAkin(3Re)−PAkin(Re), and i(3Re)−i(Re). Here, λR(3Re), ε(3Re), PAmorph(3Re), and PAkin(3Re), i(3Re) represent the respective parameters measured within (at) 3Re, while λR(Re), ε(Re), PAmorph(Re), PAkin(Re), i(Re) denote the same parameters but measured within (at) Re isophoto. The red and blue histograms signify the galaxies in TNG50 and TNG100, respectively.

In addition, the inclinations for the observational data were derived from JAM modeling, which is quite time-consuming. Therefore, we used the simple correlation $i=\arccos \sqrt{1-\epsilon }$ that allowed us to apply our method to more IFS data where dynamical modeling was not available. We computed the deviation between i values measured within R_e and 3R_e, as shown in the right-most panel of Figure 6. The differences between i(3R_e) and i(R_e) are generally minor.

4.2. Approximation of the intrinsic spin parameter λ_R, intr

The location of galaxies within the λ_R − ε diagram is influenced not only by their rotational properties but also by their inclination angles. To correct for this projection effect, we estimated an intrinsic spin parameter, λ_R, intr, which approximates the spin value for an edge-on configuration (i = 90°), under the assumption that galaxies are axisymmetric systems. Following the method of Binney & Tremaine (1987), we derived the edge-on (V/σ) as

$$\begin{array}{c}\hfill {\left(\frac{V}{\sigma }\right)}_{\mathrm{e}}={\left(\frac{V}{\sigma }\right)}_{\mathrm{obs}}\frac{\sqrt{1-\delta {cos}^{2}i}}{sini},\end{array}$$(5)

where the observed ratio ${(\frac{V}{\sigma })}_{\mathrm{obs}}=\sqrt{\frac{{\displaystyle \sum _j}{M}_j{\overline{V_j}}^2}{{\displaystyle \sum _j}{M}_j{\sigma }_j^2}}$ is computed along the line of sight for a given inclination angle i. The value δ denotes the velocity anisotropy randomly sampled from a uniform distribution in the empirically motivated range [0, 0.7ε_intr], as suggested by Cappellari et al. (2007). We derived the edge-on ellipticity by ${\epsilon }_{\mathrm{intr}}=1-\sqrt{1+\epsilon (\epsilon -2)/{\sin }^{2}i}$. We then approximated the intrinsic spin parameter λ_R, intr by the empirical formula from Emsellem et al. (2007, 2011):

$$\begin{array}{c}\hfill {\lambda }_{R,\mathrm{intr}}\approx \frac{a{(\frac{V}{\sigma })}_{\mathrm{e}}}{\sqrt{1+a^2{(\frac{V}{\sigma })}_{\mathrm{e}}^2}},\end{array}$$(6)

where a = 1.1. λ_R, intr roughly corresponds to the values of λ_R projected onto the magenta curve in Figures 4 and 5.

To assess how accurately this method recovers λ_R, intr, we calculated the deviation Δλ_R, intr between the intrinsic spin parameter and approximated its real values in edge-on views, specifically for the TNG galaxies, as shown in Figure 7. The distributions for TNG50 (red histograms) and TNG100 (blue histograms) reveal no significant differences in λ_R, intr measurements within R_e (left panel) versus within 3R_e (right panel). The results demonstrate that the method reliably recovers the intrinsic spin parameter, with deviations falling within acceptable limits.

Fig. 7. Distribution of the deviation ΔλR, intr within different regions between the intrinsic spin parameter and its edge-on projection. Left panel: Distribution of the deviation ΔλR, intr(Re). Right panel: Distribution of the deviation ΔλR, intr(3Re). The red histogram signifies the galaxy sample derived from TNG50, whereas the blue histogram represents the galaxy sample from TNG100. The vertical dashed lines show the median of the intrinsic spin spectroscopic parameter estimation error. The median and the 16th and 84th percentiles are listed at the top right of each subplot.

4.3. Mass fractions of kinematically derived structures using auto-GMM: spheroid f_spheroid and stellar halo f_halo

Recent studies (Du et al. 2019, 2020, 2021) introduced auto-GMM, an automated unsupervised method for accurately decomposing galaxy kinematic structures in simulation data. This approach analyzes stellar particles in a 3D kinematic phase space comprising the circularity parameter ϵ = j_z/j_c(e) (Abadi et al. 2003), the non-azimuthal angular momentum ratio j_p/j_c(e), and the normalized binding energy e/|e|_max (Doménech-Moral et al. 2012). The value j_z denotes the specific angular momentum component aligned with the galaxy’s minor axis (z-axis) in cylindrical rotation. This quantity is normalized by j_c(e), the maximum possible circular angular momentum corresponding to the particle’s binding energy e. The parameter j_z/j_c(e) thus quantifies rotation aligned with the galaxy’s net angular momentum, whereas j_p/j_c(e) measures misaligned rotation, and e/|e|_max indicates how tightly bound each stellar particle is.

When applied to TNG simulations, this method successfully identified distinct structural components, including cold disk, warm disk, bulge, and stellar halo. The total spheroidal mass fraction f_spheroid was obtained by combining the contributions from both the bulge (f_bulge) and the stellar halo (f_halo), with data taken from Du et al. (2020) and Du et al. (2021)¹⁶. The measurements are provided at three different scales: f(R_e) for stars within 1R_e, f(3R_e) within 3R_e, and f(all) for all stars in each galaxy. f_halo defined in this way correlates closely with the strength of merger history (Du et al. 2021; Proctor et al. 2024). One of the key advantages of auto-GMM is its allowance for weakly rotating spheroids, differing from other studies (Zana et al. 2022; Cristiani et al. 2024; Liang et al. 2025), which assume zero rotation in bulges and stellar halos–a likely unrealistic assumption.

The spheroidal mass fraction f_spheroid serves as a quantitative measure for determining the relative dominance of spheroidal components in galaxies. Based on this metric, we can classify galaxies with f_spheroid < 0.5 as FRs, which means disk structures dominate them, while those with f_spheroid > 0.5 are identified as SRs.

4.4. Assessing the significance of cylindrical rotation energy κ_rot

The significance of cylindrical rotation was quantified by the parameter ${\kappa }_{\mathrm{rot}}=\frac{K_{\mathrm{rot}}}{K}=\frac{1}{K}{\displaystyle \sum \left[\frac{1}{2}{m}_i{\left(\frac{j_{z,i}}{R_i}\right)}^2\right]}$, as defined in (Sales et al. 2010), where K and m represent the total kinetic energy and mass of each stellar particle, respectively. This widely used metric provides a clear physical interpretation of ordered rotation in simulated galaxies and is computationally straightforward to evaluate.

However, several important considerations should be noted. First, while commonly adopted, the threshold of κ_rot = 0.5 for distinguishing between disk galaxies (corresponding to FRs) and elliptical galaxies (SRs) lacks strong theoretical justification (Zhao et al. 2020; Roshan et al. 2021; Lu et al. 2025). Second, for a completely dispersion-dominated system, theoretical considerations predict κ_rot = 1/3. Finally, the physical interpretation remains ambiguous for both disk-dominated galaxies (κ_rot > 0.5) and elliptical galaxies (κ_rot < 0.5). These limitations highlight the need for careful interpretation when applying κ_rot criteria to different galactic systems.

5. Approximating κ_rot and f_spheroid from scaling relations with λ_R, intr

The intrinsic spin parameter λ_R, intr, together with the λ − ε diagram, offers a relatively accessible observational metric. A higher λ_R, intr value signifies that a galaxy’s kinematic structure is dominated by ordered rotation, while a lower value indicates a greater influence of random motion. Consequently, Figure 4 reveals a clear trend: both f_spheroid and κ_rot decrease as λ_R, intr declines. In this section, we quantitatively examine the scaling relations connecting λ_R, intr, κ_rot, and f_spheroid.

We begin by investigating the relationship between λ_R, intr and κ_rot. The fundamental difference in their velocity dependence – linear for λ_R, intr versus quadratic for κ_rot – motivates a quadratic fitting function of the form:

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{rot}}=A{\lambda }_{R,\mathrm{intr}}^2+B{\lambda }_{R,\mathrm{intr}}+C\end{array}$$(7)

where A, B, and C represent the coefficients to be determined through regression analysis. This functional form naturally accommodates their distinct velocity scalings while enabling quantitative comparisons of their kinematic relationship. The first row of Figure 8 reveals a robust correlation between λ_R, intr and κ_rot (see the fitting result at the bottom-right corners). This relationship holds consistently across both TNG50 (red) and TNG100 (blue) simulations, regardless of the measurement aperture - whether considering stars within R_e, 3R_e, or the entire system. Notably, as demonstrated in the right-most panel, λ_R, intr(R_e) alone can be used to approximate the global κ_rot(all) with an acceptable uncertainty of ∼0.05, despite being measured only within the effective radius (R_e). Furthermore, the criterion κ_rot > 0.5 has become a standard threshold for identifying disk-dominated systems in numerical simulations, approximately equivalent to λ_R, intr(R_e)≳0.4.

Fig. 8. Distribution and correlation of κrot vs. λR, intr (upper) and fspheroid vs. λR, intr (lower) measured within different regions from left to right. The solid red and blue curves show the best-fit results for the TNG50 and TNG100 datasets, respectively, while the dashed curves represent their 2σ confidence interval. The corresponding fitting coefficients and σ values are listed in the bottom-right corner of each panel. In the lower panels, we additionally provide the Pearson correlation coefficient (rp). The rp ∼ −1 value indicates a strong anticorrelation between λR, intr and fspheroid.

Our analysis further reveals a linear scaling relation between λ_R, intr and f_spheroid, as shown in the lower panels of Figure 8. This correlation robustly validates λ_R, intr as an effective diagnostic for quantifying spheroid mass fractions across different measurement apertures; see the fitting result at the bottom-right corners. The observed anticorrelation follows f_spheroid ≈ 1 − λ_R, intr (or equivalently, f_disk ≈ λ_R, intr), indicating that λ_R, intr directly traces disk dominance with a typical uncertainty of σ ∼ 0.1. This establishes λ_R, intr as a powerful observational proxy for connecting IFS kinematics to the intrinsic structural decomposition of galaxies. Even when limited to one effective radius (λ_R, intr(R_e)) – the typical spatial coverage of IFS observations – λ_R, intr retains diagnostic power for characterizing global galaxy rotation.

The parameter λ_R thus can serve as a proxy for estimating λ_R, intr, κ_rot, and f_spheroid. Empirically, λ_R, intr(R_e) provides a reasonable approximation of the global rotational properties, specifically κ_rot(all) and f_spheroid(all). Moreover, it is worth highlighting that the λ_R, intr − f_spheroids relationship is primarily determined by f_halo, shown in the third column of Figure 4. Bulges only contribute to scatter without altering the fundamental scaling relation, as shown in the fourth column. In the fourth column, we further leverage the empirical correlations between λ_R, intr(R_e) and κ_rot(all), f_spheroid(all) derived from TNG50 simulations to infer approximate values of κ_rot(all) and f_spheroid(all) based on λ_R, intr(R_e)¹⁷ for observed galaxies. Details on a subset of MaNGA galaxies, specifically nine randomly selected sample galaxies, are listed in Table 1. The complete content of MaNGA, ATLAS^3D, TNG50, and TNG100 data is publicly available.

6. Kinematic classification of galaxies into fast (disk) and slow (elliptical) rotating galaxies

In this study, we equate disk galaxies with FRs and elliptical galaxies with SRs. Regarding this, it should be noted that morphological classifications (disk vs. elliptical galaxies) and kinematic classifications (FRs vs. SRs) describe fundamentally distinct galaxy properties. Disks and ellipticals denote visual structures observed photometrically (e.g., spiral arms vs. smooth ellipsoids; Hubble 1926, 1936; Sandage 2005), while FRs and SRs categorize galaxies based on dynamical properties quantified by spin parameter λ_R (Emsellem et al. 2007; Cappellari 2016). Historically, these classifications strongly correlated–ellipticals were presumed pressure-supported SRs, while disks were rotationally dominated FRs–driving early unification in literature. However, IFS surveys such as ATLAS^3D, SAMI, CALIFA, and MaNGA have revealed significant exceptions: part of morphological ellipticals exhibit FR kinematics (Emsellem et al. 2011; Fogarty et al. 2015), and a minority of disks show SR-like suppression of rotation (Falcón-Barroso et al. 2015; Graham et al. 2018). Thus, equating disks and ellipticals directly with FRs and SRs overlooks dynamically complex populations and risks oversimplification.

Therefore, this work introduces unified classification thresholds that reconcile galaxy morphology and kinematics (see Section 6.2). Our criteria reclassify approximately 4% of previously identified disk galaxies as ellipticals and approximately 16% of FRs as SRs. With the associated publicly released catalog (Table 1), we establish a consistent disk and FR as well as elliptical and SR mapping applicable to both observational surveys and cosmological simulations.

6.1. Failure of the empirical classification of slow/fast rotators: The high-order Fourier kinematic components $\overline{k_5}$ and the empirical λ_R(R_e)−ε criterion

To evaluate the reliability of the standard FRs and SRs classification in previous studies (e.g., Krajnović et al. 2011; Emsellem et al. 2011), we also calculated the $\overline{k_5}$ (see Section 4). However, we find that $\overline{k_5}$ exhibits no significant correlation with either κ_rot or f_spheroid, suggesting it poorly traces kinematic disk structure or global rotation (Figure 9), where it is assumed that galaxies with f_spheroid < 0.5 or κ_rot > 0.5 are classified as FRs (disk-dominated), while those with f_spheroid > 0.5 or κ_rot < 0.5 are classified as SRs (spheroid-dominated). This implies that both $\overline{k_5}$ and the standard λ_R − ε criterion (thin black profiles in Figures 4 and 5) inadequately distinguish SRs from FRs, though they successfully identify extreme slow-rotating cases (Lu et al. 2023).

Fig. 9. κrot(all), fspheroid(all) of the selected ETGs samples in TNG simulation, plotted as a function of $\overline{k_5}(R_{\mathrm{e}})$. No clear correlation appears between κrot(all) and fspheroid(all) of the selected ETGs samples in the TNG simulations. The red and blue dots represent the TNG50 and TNG100 samples, respectively.

Moreover, SRs are typically identified using other empirically derived thresholds in the λ_R(R_e)−ε plane. Common selection criteria include: 1) Emsellem et al. (2007)’s constant threshold λ_R(R_e) < 0.1 (dashed blue line in Figure 4); (2) Emsellem et al. (2011)’s scaling relation ${\lambda }_R(R_{\mathrm{e}})<0.31\sqrt{\epsilon }$ (black curve); 3) Cappellari (2016)’s polygon criterion requiring λ_R(R_e) < 0.08 + ε/4 with ε < 0.4 (black polygon); and 4) van de Sande et al. (2021)’s adaptable threshold λ_R(R_e) < λ_Rstart + ε/4 where ε < 0.35 + λ_Rstart/1.538 (green polygon), adopting λ_Rstart = 0.12 for SAMI-like data quality. Given the functional similarity of these thresholds–particularly their convergence near ${\lambda }_R(R_{\mathrm{e}})=0.31\sqrt{\epsilon }$–we designate the Emsellem et al. (2011) relation as our reference criterion for distinguishing classical FRs from SRs. A revised threshold, better quantifying galactic rotation, is needed to reliably classify SRs in observations and enable robust comparisons with simulations.

6.2. Reclassification of fast and slow rotators through the kinematic bimodality

As established in the Section 5, λ_R, intr exhibits strong correlations with both f_spheroid and κ_rot, confirming its utility as a fundamental parameter for characterizing galactic rotation and classifying slow versus FRs.

Figure 10 shows a pronounced bimodality in kinematic properties, with well-separated peaks for fast and SRs. Both the ATLAS^3D and MaNGA surveys exhibit this dichotomy:

• Fast rotators (FRs) / disk galaxies are characterized by high intrinsic spin (λ_R, intr(R_e)∼0.7), thus strong rotational support (κ_rot ∼ 0.7) and low spheroid fractions (f_spheroid ∼ 0.3).

• Slow rotators (SRs) / elliptical galaxies have low spin (λ_R, intr(R_e)∼0.2), weak rotation (κ_rot ∼ 0.3), and dominant spheroid components (f_spheroid ∼ 0.8).

The vertical dashed black lines in Figure 10 mark the new thresholds we suggest for distinguishing slow and FRs, giving λ_R, intr ∼0.4, κ_rot ∼ 0.5, and f_spheroid ∼ 0.6. These thresholds roughly correspond to the minimum probability of the bimodal distributions. The threshold λ_R, intr ∼0.4 is exactly same to that suggested by Wang et al. (2024). This bimodal distribution persists across all mass ranges, revealing a significant population of rapidly rotating ETGs. Notably, the κ_rot ∼ 0.5 threshold – commonly used in simulations to select disk galaxies – proves observationally robust for identifying FRs and disk galaxies. We maintain the conventional designations of SRs and FRs but introduce refined kinematic thresholds to better distinguish between these populations. This update ensures consistency with established nomenclature while improving the physical accuracy of the classification scheme.

Fig. 10. Probability density distributions of λR, intr(Re), λR(Re), κrot(all), fspheroid(all), and fhalo(all). The solid green, solid gray, dashed blue, and dashed red profiles represent ETG samples of TNG50, TNG100, MaNGA, and ATLAS3D, respectively. From left to right: Three stellar mass ranges of 1010 − 1010.4 M⊙, 1010.4 − 1010.8 M⊙, and 1010.8 − 1011.5 M⊙. In the top right corner of each upper panel, we give the proportion and absolute count of galaxies within each mass range relative to the total ETG sample. From top to bottom, these vertical dashed black lines represent λR, intr(Re) = 0.4, λR(Re) = 0.3, κrot(all) = 0.5, fspheroid(all) = 0.6, and $f_{\mathrm{halo}}(\mathrm{all})=0.4$, respectively.

The distribution of spin parameter (λ_R) is shown in the second row of Figure 10. As λ_R constitutes the projected counterpart of intrinsic λ_R, intr, it follows that λ_R ∼ λ_R, intr sin i when assuming the galaxy as an isotropic thin disk (see Equation (6)). This systematic projection effect depresses observed λ_R values below their intrinsic values. For FRs, the peak shifts from an intrinsic value of λ_R, intr ∼ 0.65 to a lower observed value near λ_R ∼ 0.55. The threshold decreases to 0.3. It further weakens the bimodality in the distribution and reduces the diagnostic power of λ_R alone, making it less accurate than λ_R, intr for distinguishing between fast and SRs. Furthermore, van de Sande et al. (2021) selected ETGs in the SAMI Galaxy Survey using visual morphology, reporting a similarly bimodal distribution of λ_Re values. As shown in their Figure 8, these bimodality peaks occur near λ_R(R_e)∼0.5 and λ_R(R_e)∼0.1 within a comparable stellar mass range to our sample – a result consistent with our findings from both MaNGA and ATLAS^3D.

Figure 10 clearly shows a distinct population divergence between FRs classified according to our new threshold (thick solid black profiles in Figure 4) and those identified using the conventional criterion (thin solid black profiles in Figure 4). Under the previous classification scheme, only galaxies with minimal disk components and weak rotational support–specifically elliptical systems meeting the conditions κ_rot < 0.4, f_spheroid > 0.8, and f_halo > 0.6 – were categorized as SRs. The TNG simulations exhibit approximately 10% more SRs compared to observational samples. Nevertheless, the overall fraction of SRs in TNG remains consistent with those derived from both the ATLAS^3D and MaNGA surveys. Moreover, the proposed threshold offers a universal criterion for classifying galaxies of all types – including both early-type and late-type galaxies – based on their rotational properties.

6.3. Kinematic bimodality deficiency in TNG: Systematic under-rotation of simulated ETGs

While the TNG simulations reproduce a similar overall proportion of SRs to observations – as detailed in Section 6.2 – they exhibit a less distinct kinematic bimodality. This discrepancy arises from an overabundance of galaxies with intermediate rotational support, as shown in Figure 10. Both TNG50 and TNG100 show a significant deficit of FRs, leading to a systematic shift toward lower intrinsic spin parameter (λ_R, intr) values throughout the simulated population. This deficiency of fast-rotating ETGs in TNG simulations results in systematically different kinematic properties compared to observations: simulated ETGs exhibit weaker rotation (lower κ_rot), higher spheroid fractions (f_spheroid), and higher stellar halo fractions (f_halo), as shown in Figure 10. This result manifests two primary shortcomings: (1) a near-total absence of rapidly rotating systems with λ_R, intr > 0.6, and (2) a displacement of the characteristic λ_R, intr peak from the observed value of ∼0.7 to below 0.5. Correspondingly, the simulations contain an overabundance of slow-rotating ETGs and yield a higher fraction of ETGs with significant spheroidal components (f_spheroid > 0.5), many of which represent disk galaxies embedded within massive stellar halos (He et al. 2025).

It is worth highlighting that the ATLAS^3D survey contains approximately 16% fewer massive galaxies compared to MaNGA (27% vs. 43%), as evidenced by the sparse population in the top-right corner of Figure 10 and the stellar mass distribution shown in Figure 2 (top panels). This disparity likely reflects MaNGA’s selection bias toward more massive systems, but this issue has no clear effect on the kinematic bimodality, as shown in Figure 10. TNG100 and TNG50 reproduces the ATLAS^3D mass distribution reasonably well, though TNG100 slightly overproduces galaxies in the 10^10.4 − 10^10.8 M_⊙ range and underrepresents the less massive ETGs. The less pronounced kinematic bimodality is thus unlikely to be caused by the discrepancy of the proportion of ETGs within each mass range relative to the total ETG sample.

This systematic discrepancy suggests that TNG simulations produce galaxies with higher spheroid and stellar halo mass fractions than those observed. It is still unclear what reasons cause this difference. It may stem from an insufficient numerical resolution to accurately capture structural properties and their dynamical evolution, particularly in central regions, such as bars, gas inflows, and other nuclear fast-rotating structures. Furthermore, the treatment approaches of the interstellar medium and feedback subgrid models also have an impact on the kinematics in especially the inner regions of galaxies. Active galactic nuclei and stellar feedback shape the star formation history and even kinematic features via regulating gas cooling and quenching (Weinberger et al. 2018; Donnari et al. 2021). Although TNG simulations successfully reproduce many galactic properties, accurately matching observed kinematic distributions remains challenging. It is also important to note that observational IFU measurements are typically limited to central regions (within approximately R_e), which cannot fully represent the overall rotational characteristics of galaxies and introduce significant uncertainties. Consequently, the apparent bimodality observed in kinematic distributions may also be subject to reliability concerns. Moreover, as mentioned in Section 3.1, TNG50 appears to produce a smaller proportion of ETGs (27%) compared to observations (33%). This result also indicates that disk galaxies in TNG simulations exhibit lower quenching probabilities into early-type morphologies compared to observations. Consequently, in the real universe, systems with significant disk components appear more susceptible to quenching processes that transform them into ETGs–a transition that current simulations may not fully capture.

6.4. A novel method for approximating stellar halo mass fractions using λ_R, intr

The ubiquitous presence of massive stellar halos fundamentally challenges conventional galaxy classification and structural decomposition methods (Gadotti & Sánchez-Janssen 2012; Du et al. 2020; He et al. 2025). These halos primarily form through merger-driven processes – particularly major mergers and strong tidal interactions – which efficiently disrupt the disks of both progenitor and satellite galaxies (Du et al. 2021). In halo-embedded disk systems, also called Sombrero-like galaxies, where low-mass disks are embedded within massive stellar halos, traditional morphological analyses often misclassify halo mass as part of the disk structure (Du et al. 2020; He et al. 2025), leading to systematic underestimation. While classical bulges have served as historical proxies for merger activity (e.g., Toomre 1977; Aguerri et al. 2001), TNG simulations demonstrate their inadequacy in capturing the full scope of merger histories (Du et al. 2021), necessitating more robust diagnostics.

Stellar halos and suppressed rotation serve as sensitive tracers of galactic merger events, as mergers dynamically convert ordered rotation into random motion, although some disk will survive from mergers with special orbital configuration (Zeng et al. 2021). This process simultaneously builds up stellar halos and erases disk kinematics. The first row of Figure 11 reveals a moderate anticorrelation between λ_R, intr and the stellar halo fraction f_halo derived from auto-GMM analysis, though with a substantial scatter of σ ∼ 0.15.

Fig. 11. Correlations of λR, intr versus fhalo and f′halo versus fhalo. Top row: Relationships between λR, intr and fhalo measured within Re (first column), 3Re (second column), and for all stars (third column). Fourth column: Relation between λR, intr(Re) and fhalo(all). The solid red and blue lines, along with their corresponding dashed lines, represent the linear fits and 2σ confidence interval for the TNG50 and TNG100 datasets, respectively. Fit coefficients, σ values, and Pearson correlation coefficients (rp) are provided in the top-right corner of each panel. Second row: Similar layout, but illustrating the relationship between f′halo and fhalo. Here, the solid black and green lines, with accompanying dashed lines, indicate the linear fits and 2σ confidence interval for the TNG50 and TNG100 datasets, respectively.

We therefore introduce a novel methodology that leverages the λ_R, intr–f_spheroids relation established in Section 5 to estimate stellar halo masses, despite the considerable associated uncertainty. Our approach computes f′_halo = f_spheroids − f′_bulge, where bulges are operationally defined as bright central concentrations according to classical morphological criteria (Hubble 1936; Buta 2013), without distinguishing between classical and pseudo-bulge subtypes. The bulge mass fraction f′_bulge is approximated by integrating the stellar surface density profile over the central 0–1 kpc region in face-on projections, while f_spheroids is approximated using the scaling relation with λ_R, intr. Meanwhile, we examine whether f′_bulge deviates further from f_bulge if f′_bulge is approximated by integrating the stellar surface density profile within the central 0–2 kpc or 0–3 kpc regions.

We anticipate improved accuracy with future IFS measurements extending to larger radii, as suggested by the second and third columns of Figure 11. It is noteworthy that currently, tracer populations such as planetary nebulae can, to a certain extent, already facilitate the extension of IFS measurements to larger radii (e.g., Pulsoni et al. 2018). For f′_halo we obtain a somewhat acceptable agreement with dynamically derived f_halo values, exhibiting a scatter of approximately σ ∼ 0.15. It is worth mentioning that this method is unable to measure stellar halos with f′_halo ≲ 0.3, where the stellar halo is too small. This method enables a somewhat statistical analysis of ETG merger histories with massive stellar halos using IFS measurements (e.g., Huang et al. 2013; Deason et al. 2016; D’Souza & Bell 2018; Monachesi et al. 2019; Du et al. 2021; Zhu et al. 2022), overcoming limitations inherent in purely morphological stellar halo identification.

7. Radial distribution of V_rot/σ in the inner regions of selected fast rotators

Fast-rotating ETGs, characterized by λ_R, intr > 0.4 as defined in Section 6.2, serve as ideal tracers of minimally disturbed galactic evolution (van Dokkum et al. 2015; Cappellari 2016; Wang et al. 2020; Zhu et al. 2024). Their kinematic radial profiles, relatively unaffected by major mergers or strong tidal interactions, provide critical benchmarks for studying intrinsic galaxy evolution in both simulations and observations. We thus compare the radial distribution of V_rot/σ in Figure 12. To systematically account for the effects of stellar mass and inclination, we categorize our galaxy sample into three mass bins from left to right. Each mass bin is further divided into two inclination subgroups, i.e., i = 30° −60° and i = 60° −90°, from top to bottom. This comprehensive sampling strategy enables us to isolate fundamental kinematic trends while controlling for projection effects and mass-dependent behaviors.

Fig. 12. Median Vrot/σ profiles for FRs selected using our new threshold (κrot(all) > 0.5), projected along a random LOS. Top row: Galaxy samples that fall within mass ranges [1010.0, 1010.4] M⊙, with galaxy samples presented at different inclinations: 30° −60° and 60° −90°, respectively, from left to right. The solid green lines, solid gray lines, dashed red lines, and dashed blue lines, along with their corresponding shaded regions in the figure, represent the data from TNG50, TNG100, ATLAS3D, and MaNGA, respectively. The datasets from the simulation and the observation showcase median profiles, with shaded regions indicating the 16th and 84th percentiles of the distribution. Vertical and horizontal dashed lines are included as visual aids for comparing simulations and observations. The profiles of the simulated ETG FRs closely resemble those observed. Middle and bottom rows: same as the top row, but with galaxy samples that fall within mass ranges [1010.4, 1010.8] M⊙ and [1010.8, 1011.5] M⊙.

We utilized KINEMETRY to separately fit the galaxy’s velocity field V and velocity dispersion field σ, to obtain the radial distributions of V_rot and σ for each galaxy. Figure 12 shows that the kinematic properties of fast-rotating galaxies in TNG50 (green) match reasonably well with ATLAS^3D galaxies. TNG100 galaxies exhibit slightly smaller V_rot/σ values in their central regions, likely due to enhanced numerical heating effects stemming from insufficient resolution, particularly in lower-mass systems. The MaNGA survey similarly shows somewhat depressed V_rot/σ ratios compared to ATLAS^3D, consistent with expectations given its poorer spatial resolution. In particular, the observed sharp rise in V_rot/σ within the central regions (R < 0.5R_e) of some massive ATLAS^3D galaxies (bottom-left panel) results from fast-rotating structures. Moreover, it is worth highlighting that V_rot/σ keeps rising toward R > R_e, IFS kinematic measurements in R < R_e clearly cannot characterize the full dynamical state of FRs. The potential existence of massive stellar halos around disks (Gadotti & Sánchez-Janssen 2012; He et al. 2025), for example the Sombrero galaxy, may significantly change the overall λ_R and V_rot/σ at large radius.

Our results indicate that the TNG simulations produce a substantially larger population of galaxies with intermediate rotational support, as detailed in Section 6.2. Although TNG50 reproduces the kinematic properties of fast-rotating ETGs – consistent with previous studies (Xu et al. 2019; Zhang et al. 2025) – it is important to note that those earlier studies did not exclude the potential effect of SRs in their samples, a methodological choice that may obscure kinematic distinctions via inducing large scatters. It may diminish the visibility of potential discrepancies. Furthermore, differences in the inclination angles among galaxy samples can also influence the analysis. This overabundance of intermediate-rotators in TNG could explain the shallower V_rot/σ(R) radial profiles reported in Pulsoni et al. (2020).

8. Discussions: evaluating kinematic analysis in advanced cosmological simulations

Similar analyses have suggested consistency between TNG simulations and observations. However, when incorporating careful prevalence estimates and inclination corrections, certain discrepancies emerge. The random orientation of galaxies may introduce nonnegligible scatter in λ_R − ε diagrams, potentially accentuating apparent simulation-observation agreement – as exemplified by the close but uncorrected alignment shown in Figures 4 and 5. The differences in probability density profiles are discernible (Figure 10, second row). A typical inclination angle i varies between 30 and 60 degrees, which reduces the spin parameter we measured along the line of sight by a factor of sin(i)≈0.5 − 0.8. The λ_R values of FRs, generally within the range of 0.3–0.7, can thus be significantly underestimated. This geometric effect further reduces the potential kinematic discrepancies between observations and simulations. The λ_R − ε diagram thus alone cannot conclusively establish kinematic consistency between datasets.

Our results highlight the importance of refined analyses when evaluating simulation-observation consistency for not only TNG simulations. Using the Magneticum Pathfinder simulation, Schulze et al. (2018) reported general agreement in λ_R(R_e)−ε distributions between simulated and observed ETGs. Schulze et al. (2020) further extended this approach to ∼500 galaxies (2 × 10¹⁰ M_⊙ < M_* < 1.7 × 10¹² M_⊙), noting reasonable matches in two-point gradient relationships. Similar concordance was observed in EAGLE simulations by Walo-Martín et al. (2020) and van de Sande et al. (2019) (M > 5 × 10⁹ M_⊙,  z = 0), though the latter study noted systematically lower values in Horizon-AGN and Magneticum. Our findings suggest that such comparisons may benefit from explicit treatment of projection effects to ensure conclusive interpretations.

Kinematic properties serve as a fossil record of a galaxy’s formation history. Consequently, a first-order consistency between simulations and observations may be insufficient; moreover, the nonnegligible mixture of SRs and FRs in the previous classification schemes is also likely to mitigate potential discrepancy. In this study, we thus performed detailed statistical comparisons of their projection-corrected distributions across different mass ranges within a unified quantitative framework of both simulations and observations.

9. Summary and conclusions

This study presents a comparative analysis of galaxy kinematics between the TNG simulations (TNG50 and TNG100) and observational data from the MaNGA and ATLAS^3D surveys. We revisited the classification of ETGs into FRs (disk-dominated systems) and SRs (dispersion-supported ellipticals) using intrinsic dynamical parameters – including the intrinsic spin parameter λ_R, intr, rotational energy fraction κ_rot, and spheroid mass fraction f_spheroid. Our results show that widely adopted criteria, such as the empirical ${\lambda }_R(R_{\mathrm{e}})=0.31\sqrt{\epsilon }$ relation and the Fourier harmonic coefficient $\overline{k_5}$, are inadequate for robust FR and SR separation due to their weak correlation with fundamental dynamical properties. We instead proposed revised thresholds (λ_R, intr ∼0.4, κ_rot ∼ 0.5, f_spheroid ∼ 0.6) derived from bimodal kinematic distributions, offering a physically motivated universal distinction between rotation- and dispersion-dominated galaxies.

We also examined inter-parameter correlations among commonly used kinematic quantities and applied scaling relations derived from TNG simulations to estimate λ_R, intr, κ_rot, and f_spheroid in observed MaNGA and ATLAS^3D galaxies. Furthermore, we introduced a novel method for inferring stellar halo mass fractions from integral-field spectra (IFS) kinematic measurements, though current estimates remain subject to significant uncertainties.

Our comparison reveals a notable deficiency in the TNG simulations: they exhibit less distinct kinematic bimodality than observations, producing fewer FRs and showing systematically lower λ_R, intr values. An overabundance of galaxies with intermediate rotational support and elevated spheroid or halo fractions drives this discrepancy. Although the global fractions of SRs and FRs are roughly consistent between simulations (∼38% SRs) and observations (∼30%), the weakened bimodality in TNG reflects a failure to fully capture the dynamical diversity of real ETGs. These findings highlight the need for improved numerical resolution and more sophisticated subgrid physics in future simulations, as well as more extensive observational constraints, to advance our understanding of galaxy formation and evolution.

Data availability

Table 1 is only available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/707/A137

Acknowledgments

The authors gratefully acknowledge insightful discussions with Luis C. Ho, Lars Hernquist, Bitao Wang, and Ling Zhu, and also extend our heartfelt gratitude to the referee for their insightful comments, which have significantly improved the quality of the paper. This work is supported by the National Natural Science Foundation of China under grant No. 12573010, the National Key R&D Program of China (No. XY-2025-1459), the Fundamental Research Funds for the Central Universities (No. 20720230015), the Science Fund for Creative Research Groups of the National Natural Science Foundation of China (No. 12221003), and the China Manned Space Program (No. CMS-CSST-2025-A10). The TNG50/TNG100 simulation used in this work, one of the flagship runs of the IllustrisTNG project, has been run on the HazelHen Cray XC40-system at the High Performance Computing Center Stuttgart as part of project GCS-ILLU of the Gauss Centers for Supercomputing (GCS). This work is also strongly supported by the Computing Center in Xi’an, China. The authors also acknowledge the assistance of DeepSeek’s AI language model in improving the language and clarity of this manuscript.

References

All Tables

All Figures

	Fig. 1. g − r color–stellar mass diagram. Top: observational data (from MaNGA and ATLAS3D). Bottom: simulated data (from TNG50 and TNG100). Both panels encompass all samples within the mass range of 1010 − 1011.5 M⊙. Red sequence ETGs are selected above the dashed line.
In the text

Fig. 2. Mass–size relation. Top: stellar mass distribution function of the chosen sample of ETGs from observational data and the TNG simulations. Bottom: Circularized effective radii Re, of the selected ETG samples in the TNG100 and TNG50 simulations as a function of stellar mass and compared with observations (ATLAS3D and MaNGA). Both TNG50 (solid green line) and TNG100 (solid gray line) galaxies share comparable stellar masses and sizes with the ATLAS3D (dashed red line) and MaNGA (dashed blue line) galaxies. The lines and points represent the median values, with colored regions and error bars indicating the 16th to 84th percentile range of each distribution. To avoid overlap, we shift the error bars slightly for clarity.

In the text

Fig. 3. Velocity maps and models of examples of a FR (left) and a SR (right), with raw data (top), a reconstructed velocity field employing a pure cosine law along ellipses (middle), and residuals depicting deviations (bottom). Our visualization emphasizes structural features through isodensity ellipses spaced at regular Re intervals from Re to 3Re. The major axis orientations – with solid lines representing photometrically derived morphological axes and dashed lines indicating kinematically determined axes – are both measured at Re.

In the text

Fig. 4. λR(Re)−ε(3Re) diagram color-coded by κrot(all), fspheroid(all), fhalo(all), and fbulge(all) for TNG50 (upper) and TNG100 (lower). Here, κrot(all), fspheroid(all), fhalo(all), and fbulge(all) represent the significance of cylindrical rotation, the spheroidal mass fraction, the stellar halo mass fraction, and the bulge mass fraction, respectively (see Sections 4.3 and 4.4). The magenta curve shows the λR − ε relation for edge-on, axisymmetric galaxies derived by Cappellari et al. (2007). The dotted black curves illustrate this relation across varying i, with steps of Δi = 10°. Additionally, the thin dashed black lines represent the theoretical distribution for the different rotations spaced at intervals of ΔλR, intr = 0.1. The bold black curve represents our proposed new threshold for distinguishing between SRs and FRs. This threshold is derived from the tensor virial theorem for oblate galaxies with intrinsic parameters λR, intr(Re) = 0.4 and εintr = 0.525, assuming an anisotropy of δ = 0.7εintr = 0.367. For comparison, the thin solid black line represents the standard empirical threshold defined by ${\lambda }_R(R_{\mathrm{e}})=0.31\times \sqrt{\epsilon }$ proposed in Emsellem et al. (2011, Eq. (3)). We present the proportion of SRs classified by our new threshold alongside the standard old threshold at the upper-left corner. We smoothed the color data using the locally weighted regression method LOESS (Cleveland & Devlin 1988; Cappellari et al. 2013a) with a smoothing factor of 0.1. The horizontal dashed blue line represents the empirical threshold defined by λR(Re) = 0.1 according to Emsellem et al. (2007); the black polygon denotes the empirical threshold set by λR(Re) < 0.08 + ε/4  and  ε < 0.4 as per Cappellari (2016); and the green polygon signifies the empirical threshold given by ${\lambda }_R(R_{\mathrm{e}})<{\lambda }_{R\mathrm{start}}+\epsilon /4\mathrm{and}\epsilon <0.35+\frac{{\lambda }_{R\mathrm{start}}}{1.538}$ (with λRstart = 0.12 when the IFS data quality is comparable to that of SAMI) based on van de Sande et al. (2021).

In the text

	Fig. 5. λR(Re)−ε(Re) diagram of the ATLAS3D(left panel) and MaNGA (right panel) samples. The proportions of SRs and FRs classified using our new threshold, compared to the standard old definition, are provided at the top. The definitions of all the lines are the same as those in Figure 4.
In the text

Fig. 6. Distributions of the differences between various galactic parameters measured within (at) 3Re and Re isodensity ellipses. From left to right: λR(3Re)−λR(Re), ε(3Re)−ε(Re), PAmorph(3Re)−PAmorph(Re), PAkin(3Re)−PAkin(Re), and i(3Re)−i(Re). Here, λR(3Re), ε(3Re), PAmorph(3Re), and PAkin(3Re), i(3Re) represent the respective parameters measured within (at) 3Re, while λR(Re), ε(Re), PAmorph(Re), PAkin(Re), i(Re) denote the same parameters but measured within (at) Re isophoto. The red and blue histograms signify the galaxies in TNG50 and TNG100, respectively.

In the text

Fig. 7. Distribution of the deviation ΔλR, intr within different regions between the intrinsic spin parameter and its edge-on projection. Left panel: Distribution of the deviation ΔλR, intr(Re). Right panel: Distribution of the deviation ΔλR, intr(3Re). The red histogram signifies the galaxy sample derived from TNG50, whereas the blue histogram represents the galaxy sample from TNG100. The vertical dashed lines show the median of the intrinsic spin spectroscopic parameter estimation error. The median and the 16th and 84th percentiles are listed at the top right of each subplot.

In the text

Fig. 8. Distribution and correlation of κrot vs. λR, intr (upper) and fspheroid vs. λR, intr (lower) measured within different regions from left to right. The solid red and blue curves show the best-fit results for the TNG50 and TNG100 datasets, respectively, while the dashed curves represent their 2σ confidence interval. The corresponding fitting coefficients and σ values are listed in the bottom-right corner of each panel. In the lower panels, we additionally provide the Pearson correlation coefficient (rp). The rp ∼ −1 value indicates a strong anticorrelation between λR, intr and fspheroid.

In the text

	Fig. 9. κrot(all), fspheroid(all) of the selected ETGs samples in TNG simulation, plotted as a function of $\overline{k_5}(R_{\mathrm{e}})$. No clear correlation appears between κrot(all) and fspheroid(all) of the selected ETGs samples in the TNG simulations. The red and blue dots represent the TNG50 and TNG100 samples, respectively.
In the text

Fig. 10. Probability density distributions of λR, intr(Re), λR(Re), κrot(all), fspheroid(all), and fhalo(all). The solid green, solid gray, dashed blue, and dashed red profiles represent ETG samples of TNG50, TNG100, MaNGA, and ATLAS3D, respectively. From left to right: Three stellar mass ranges of 1010 − 1010.4 M⊙, 1010.4 − 1010.8 M⊙, and 1010.8 − 1011.5 M⊙. In the top right corner of each upper panel, we give the proportion and absolute count of galaxies within each mass range relative to the total ETG sample. From top to bottom, these vertical dashed black lines represent λR, intr(Re) = 0.4, λR(Re) = 0.3, κrot(all) = 0.5, fspheroid(all) = 0.6, and $f_{\mathrm{halo}}(\mathrm{all})=0.4$, respectively.

In the text

Fig. 11. Correlations of λR, intr versus fhalo and f′halo versus fhalo. Top row: Relationships between λR, intr and fhalo measured within Re (first column), 3Re (second column), and for all stars (third column). Fourth column: Relation between λR, intr(Re) and fhalo(all). The solid red and blue lines, along with their corresponding dashed lines, represent the linear fits and 2σ confidence interval for the TNG50 and TNG100 datasets, respectively. Fit coefficients, σ values, and Pearson correlation coefficients (rp) are provided in the top-right corner of each panel. Second row: Similar layout, but illustrating the relationship between f′halo and fhalo. Here, the solid black and green lines, with accompanying dashed lines, indicate the linear fits and 2σ confidence interval for the TNG50 and TNG100 datasets, respectively.

In the text

Fig. 12. Median Vrot/σ profiles for FRs selected using our new threshold (κrot(all) > 0.5), projected along a random LOS. Top row: Galaxy samples that fall within mass ranges [1010.0, 1010.4] M⊙, with galaxy samples presented at different inclinations: 30° −60° and 60° −90°, respectively, from left to right. The solid green lines, solid gray lines, dashed red lines, and dashed blue lines, along with their corresponding shaded regions in the figure, represent the data from TNG50, TNG100, ATLAS3D, and MaNGA, respectively. The datasets from the simulation and the observation showcase median profiles, with shaded regions indicating the 16th and 84th percentiles of the distribution. Vertical and horizontal dashed lines are included as visual aids for comparing simulations and observations. The profiles of the simulated ETG FRs closely resemble those observed. Middle and bottom rows: same as the top row, but with galaxy samples that fall within mass ranges [1010.4, 1010.8] M⊙ and [1010.8, 1011.5] M⊙.

In the text