The stellar halos of ETGs in the IllustrisTNG simulations: the photometric and kinematic diversity of galaxies at large radii

We characterize the photometric and kinematic properties of simulated early-type galaxy (ETG) stellar halos, and compare them to observations. We select a sample of ~1200 ETGs in the TNG100 and TNG50 simulations, spanning a stellar mass range of $10^{10.3}-10^{12}M_{\odot}$ and within the range of (g-r) colour and lambda-ellipticity diagram populated by observed ETGs. We determine photometric parameters, intrinsic shapes, and kinematic observables in their extended stellar halos. We study the variation in kinematics from center to halo and connect it to a change in the intrinsic shape of the galaxies. We find that the simulated galaxy sample reproduces the diversity of kinematic properties observed in ETG halos. Simulated fast rotators (FRs) divide almost evenly in one third having flat lambda profiles and high halo rotational support, a third with gently decreasing profiles, and another third with low halo rotation. Slow rotators (SRs) tend to have increased rotation in the outskirts, with half of them exceeding lambda=0.2. For $M_{*}>10^{11.5}M_{\odot}$ halo rotation is unimportant. A similar variety of properties is found for the stellar halo intrinsic shapes. Rotational support and shape are deeply related: the kinematic transition to lower rotational support is accompanied by a change towards rounder intrinsic shape. Triaxiality in the halos of FRs increases outwards and with stellar mass. Simulated SRs have relatively constant triaxiality profiles. Simulated stellar halos show a large variety of structural properties, with quantitative but no clear qualitative differences between FRs and SRs. At the same stellar mass, stellar halo properties show a gradual transition and significant overlap between the two families, despite the clear bimodality in the central regions. This is in agreement with observations of extended photometry and kinematics. [abridged]


Introduction
The family of early type galaxies (ETGs) encompasses galaxies that have typically ceased their star formation at early times, with red colors and small amounts of cold gas and dust today, and that mainly consist of elliptical and lenticular galaxies (Roberts & Haynes 1994;Kauffmann et al. 2003;Blanton & Moustakas 2009). Ellipticals essentially divide into two classes with distinct physical properties (e.g. Kormendy et al. 2009, and references therein): those with low to intermediate masses and coreless luminosity profiles, that rotate rapidly, are relatively isotropic and oblate-spheroidal, and have high ellipticities and disky-distorted isophotes; and those which are frequently among the most massive galaxies, with cored profiles, mostly non-rotating, anisotropic and triaxial, relatively rounder than than coreless systems, and with boxy-distorted isophotes. Thus the dichotomy in the light distributions of the ellipticals roughly corresponds to different kinematic properties, with coreless disky objects being rotationally supported, and cored boxy galaxies having low rotation (Bender 1987). With the advent of integral field spectroscopy (IFS) the classification of elliptical galaxies has shifted to a kinematics-based division between fast rotators (FR) and slow rotators (SR) Graham et al. 2018). In particular low mass, coreless, FR ellipticals share similar properties with lenticular galaxies, which are are also included in the FR family, while massive cored ellipticals are typically SRs.
The formation of massive ETGs is believed to have occurred in two phases (e.g. Oser et al. 2010). In an initial assembly stage, gas collapses in dark matter halos and forms stars in a brief intense burst which is quickly quenched (e.g. Thomas et al. 2005;Conroy et al. 2014;Peng et al. 2010). Present-day simulations agree in that the progenitors of FR and SR at these high redshifts are indistinguishable (Penoyre et al. 2017;Lagos et al. 2017;Schulze et al. 2018, with Illustris, Eagle, and Magneticum, respectively). At z 1 the accretion-dominated phase overtakes, whereby ETGs grow efficiently in size through a series of merger episodes which enrich the galaxies with accreted (ex-situ) stars. The ΛCDM cosmology predicts that structures form hierarchically, in which more massive systems form through the accretion of less massive objects. This means that more massive galaxies can have accreted fractions larger than 80%, while lower mass galaxies are mostly made of in-situ stars, and the accreted components are mainly deposited in the outskirts (Rodriguez-Gomez et al. 2016;Pillepich et al. 2018a). The slow/fast rotators (i.e. the core/coreless) classes result from different formation pathways characterized by different numbers of mergers, merger mass ratio, timing, and gas fractions (Naab et al. 2014;Penoyre et al. 2017, see also the discussion in Kormendy et al. 2009), although the details still depend on the star formation and AGN feedback models adopted by the numerical models (Naab & Ostriker 2017). In general, the result of a formation history dominated by gas dissipation is most likely a coreless FR, while dry major mergers often result in SRs.
The two-phase formation scenario is supported both by observations of compact red nuggets at z ∼ 2, a factor of 2-4 smaller than present day ellipticals (Daddi et al. 2005;Trujillo et al. 2007;van Dokkum et al. 2008), and by evidence for a subsequent rapid size growth with little or no star formation (e.g. van Dokkum et al. 2010;Damjanov et al. 2011;van der Wel et al. 2014). The merger driven size growth is supported by the observed rate of mergers from pair counts and identified interacting galaxies (Hopkins et al. 2008;Robaina et al. 2010), as well as the observation of tidal debris from recent accretion events in the halos of many galaxies (e.g. Malin & Carter 1983;Janowiecki et al. 2010;Iodice et al. 2017;Mancillas et al. 2019).
A consequence of the two-phase formation is that ETGs are layered structures in which the central regions are the remnants of the stars formed in-situ, while the external stellar halos are principally made of accreted material (Bullock & Johnston 2005;Cooper et al. 2010), even though the details strongly depend on stellar mass (Pillepich et al. 2018a). Because of the different nature of the stellar halos, galaxies are expected to show significant variation of physical properties from central regions to large radii, such as shapes of the light profiles (Huang et al. 2013;D'Souza et al. 2014;Spavone et al. 2017), stellar populations (Pastorello et al. 2014;Zibetti et al. 2020), and kinematics (Coccato et al. 2009;Romanowsky & Fall 2012;Arnold et al. 2014;Foster et al. 2016).
Kinematic measurements in the outer halos of ETGs require alternative kinematic tracers to overcome the limitations from the faint surface brightness in these regions, such as planetary nebulae (PNe) (e.g., the ePN.S survey, Arnaboldi et al. 2017, see Sect. 2), or globular clusters (e.g., the SLUGGS survey, Brodie et al. 2014). Recently Pulsoni et al. (2018) found evidence from the ePN.S survey for a kinematic transition between the central regions and the outskirts of ETGs. Despite the FR/SR dichotomy of their central regions, these ETG halos display a variety of kinematic behaviors. A considerable fraction of the ePN.S FRs show reduced rotational support at large radii, which has been interpreted as the fading of a rotating, disk-like component into a more dispersion dominated spheroid; almost half of the FR sample shows kinematic twists or misalignments at large radii, indicating a variation of their intrinsic shapes, from oblate at the center to triaxial in the halo. SRs instead have increased rotational support at large radii. While a smaller group of FRs stands out for having particularly high V/σ ratio in the halo, most of the ePN.S FRs and SRs have similar V/σ ratio in the halo regions. These results suggest the idea that at large radii the dynamical structure of these galaxies could be much more similar than in their high-density centers: if halos are mainly formed from accreted material, their common origin would explain their similarities. The observed kinematic transition to the halo and its dependence on stellar mass seems to support such an interpretation.
The goal of this paper is to better understand the structural changes between ETG centers and stellar halos with a large and well-resolved sample of simulated galaxies. We study the stellar halo structure, i.e., the rotational support and intrinsic shapes of the simulated galaxies, we compare the results with observations, and we investigate how the radial variations in rotational support relate to changes in the halo shapes. We use the IllustrisTNG simulations Pillepich et al. 2018a;Naiman et al. 2018;Marinacci et al. 2018;Nelson et al. 2018Nelson et al. , 2019b, a suite of magnetohydrodynamical simulations that models the formation and evolution of galaxies within the ΛCDM paradigm. It builds and improves upon the Illustris simulation Vogelsberger et al. 2014), using a refined galaxy formation model. This includes a revised implementation of galactic winds, a new black hole feedback mechanism, an updated stellar evolution and gas chemical enrichment model, magnetohydrodynamics, and improvements in the numerical scheme. The complete description of the galaxy formation model can be found in the two methods papers by Weinberger et al. (2017) and Pillepich et al. (2018b). For this work we consider two cosmological volumes with side lengths ∼ 100 Mpc and ∼ 50 Mpc, which are referred to as TNG100 and TNG50. TNG100 includes 2 × 1820 3 initial resolution elements, the mean baryonic gas mass resolution is 1.4 × 10 6 M , the gravitational softening length for the stars is 0.74 kpc at redshift z = 0. TNG50 is the highest resolution realization of the IllustrisTNG project Nelson et al. 2019a): its particle resolution is more than 15 times better than TNG100, with a baryon mass resolution of 8.5 × 10 4 M and 2 × 2160 3 initial resolution elements. The gravitational softening length for the stars in TNG50 is 0.288 kpc at redshift z = 0.
The paper is organized as follows. For comparing the TNG galaxies properties with observations, we first summarize in Section 2 how different ETG surveys select their samples, and how physical quantities are measured. Section 3 then describes and illustrates our methods to derive photometric and kinematic measurements for the simulated galaxies. After selecting the sample of ETGs from the TNG100 and TNG50 simulations (Section 4), we proceed to show the photometric results in Section 5 and the kinematic results in Section 6. Section 7 relates the variation in the kinematic properties from central regions to halos to the parallel changes in the intrinsic structure of galaxies. In a companion paper we will explore the dependence of these properties on the accretion history of galaxies. Finally Section 8 summarizes our conclusions.

Observed parameters of ETGs
In this paper we compare the kinematic results for the central regions of the simulated TNG galaxies with IFS measurements from the surveys Atlas3D , MANGA (Bundy et al. 2015), SAMI (Croom et al. 2012) and MASSIVE (Ma et al. 2014).
Kinematics measurement at large radii are notably difficult to obtain for ETGs, and therefore discrete kinematic tracers such as planetary nebulae (PNe) and globular clusters (GCs) are typically used to overcome the limitations of absorption line spectroscopy, which is restricted to the central 1-2 R e . PNe are established probes of the stellar kinematics in ETG halos (Hui et al. 1995;Arnaboldi et al. 1996;Méndez et al. 2001;Coccato et al. 2009;Cortesi et al. 2013), out to very large radii (Longobardi et al. 2015;Hartke et al. 2018). Since they are drawn from the main stellar population, their kinematics traces the bulk of the host-galaxy stars, and are directly comparable to integrated light measurements. The relation between GCs and the underlying galaxy stellar population is less straightforward (Forbes & Remus 2018). In general GCs do not necessarily follow the surface brightness distribution and kinematics of the stars (e.g. Brodie & Strader 2006;Coccato et al. 2013;Veljanoski et al. 2014), although there is growing evidence for red, metal-rich GCs to be tracers of the host galaxy properties (Fahrion et al. 2020;Dolfi et al. 2020). Therefore we here compare the kinematics of the simulated galaxies and their stellar halos at large radii with PN kinematic results from the ePN.S early-type galaxy survey (Arnaboldi et al. 2017, Arnaboldi et. al., in prep.).
Below we describe the sample properties for the different surveys, and we give details and sources of the measured quantities used though out this paper.
Sample properties -The Atlas3D survey selected ETGs from a volume-limited sampe of galaxies, with distance within 42 Mpc, and sky declination δ such that |δ − 29 • | < 35 • ), brighter than M K < −21.5 mag. From this parent sample ETGs were morphologically selected as all the galaxies without visible spiral structure. This morphological selection is broadly similar to a selection of the red sequence . The At-las3D ETG sample contains 68 Es and 192 S0s. The SAMI survey (Croom et al. 2012) selected a volume and magnitude limited sample of galaxies in the redshift range 0.004 < z < 0.095, covering a broad range in galaxy stellar mass (M * = 10 8 − 10 12 M ) and environment (field, group, and clusters). This sample is not morphologically selected, but we use the data from  where the quality cuts and the imposed threshold on the velocity dispersion σ > 70 km/s bias the sample towards the ETGs (82%). The galaxies of the MANGA survey (Bundy et al. 2015) are selected from the NASA-Sloan Atlas 1 (NSA) catalog (which is based on the Sloan Digital Sky Survey (SDSS) Data Release 8, Aihara et al. 2011) at low redshift (0.01 < z < 0.15), to follow a flat distribution in stellar mass in the range M * = 10 9 − 10 12 ; in this paper we will compare only with MANGA's galaxies classified as ellipticals or lenticulars as in Graham et al. (2018). The MASSIVE survey (Ma et al. 2014) targets all the most massive ETGs (M * 10 11.5 M ) within a distance of 108 Mpc. Finally, the ePN.S sample of ETGs is magnitude limited M K −23, and includes objects with different structural parameters. This ensures the sample to be a representative group of nearby ETGs. The ePN.S kinematic results (Pulsoni et al. 2018) combine PN kinematics in the halos with literature absorption line data for the central regions.
Colors -The MANGA galaxies, and most of the Atlas3D and MASSIVE objects have measured g − r colors in the NSA catalog. For the SAMI galaxies  report g−i colors, which we convert to g − r using the transformation equation derived in App. A. For all of the ePN.S sample, and some of the Atlas3D and MASSIVE galaxies that are not in the NSA catalog, we use B − V colors corrected for galactic extinction from the Hyperleda 2 catalog (Makarov et al. 2014), and convert to g − r colors using the relations in App. A.
Sizes -For the Atlas3D sample we use the effective radii (R e ) values in Table 3 of . Those for the MAS-SIVE galaxies are from Ma et al. (2014, Graham et al. (2018). For SAMI we use the data presented in , and we circularise the effective semi-major axis by using the reported value for the ellip-ticity. The half light radii for the ePN.S galaxies are in Table 2 from Pulsoni et al. (2018). These are effective semi-major axis distances measured from the most extended photometric profiles available from the literature, extrapolated to very large radii with a Sérsic fit. The ellipticity assumed is in their Table 1.
Stellar masses -The IllustrisTNG model assumes a Chabrier (2003) initial mass function (IMF). The stellar masses for the SAMI survey in  are derived using a color-mass relation, and a Chabrier IMF. For Atlas3D, MAS-SIVE, and MANGA we use the total absolute K-band luminosity M K from the same tables referenced above, which are derived from the 2MASS extended source catalog (Jarrett et al. 2003), and already corrected for galactic extinction. The luminosities M K are then corrected for missing flux as in Scott et al. (2013), M K corr = 1.07M K + 1.53, and converted in stellar masses with the formula from van de Sande et al. (2019) log 10 M * = 10.39 − 0.46(M K corr + 23), which uses the stellar population model-based mass-to-light ratio from Cappellari et al. (2013), their [log(M/L) Salp ], converted to a Chabrier IMF. For the ePN.S sample we derive stellar masses using the integrated luminosities from extended photometric profiles, extrapolated to infinity with a Sérsic fit (references in Pulsoni et al. 2018). We convert the integrated values to stellar masses by using the non-dereddened relations between colors and mass-to-light ratios for ellipticals and S0 galaxies from García-Benito et al. (2019), which assume a Chabrier IMF.
Ellipticities -For the Atlas3D galaxies we use the ellipticity ε measurements within 1 R e reported in Table B1 of Emsellem et al. (2011). 17 out of 260 Atlas3D objects have obvious bar components: for these cases the ellipticity is measured at larger radii (typically 2.5 -3R e ). Ellipticities for the SAMI galaxies are from , and are average ellipticities of the galaxies within 1 R e . MANGA's ellipticities from Graham et al. (2018) are also measurements within the 1 R e isophote, while for the MASSIVE sample Veale et al. (2017) uses ellipticities from NSA where available, and from 2MASS otherwise, which are globally fitted values. The ellipticity profiles for the ePN.S galaxies are referenced in Pulsoni et al. (2018).
Angular momentum parameters λ e -The parameter λ e is derived in the different surveys using different integration areas. While Emsellem et al. 2011 (Atlas3D) and Veale et al. 2017 (MASSIVE) use circular apertures of radius R e , van de Sande et al. 2017 (SAMI) prefer elliptical apertures with semi-major axis R e , and Graham et al. 2018 (MANGA) integrate over the half-light ellipse (an ellipse covering the same area as a circle with radius R e , i.e. with semi-major axis R e √ 1 − ε, where ε is the ellipticity ).
V/σ profiles -The V/σ profiles for the Atlas3D and the ePN.S galaxies are derived from the ratio of the rotation velocity V rot and the azimuthally averaged velocity dispersion in elliptical radial bins. For the Atlas3D galaxies we apply directly the procedure described in Sect. 3.4 to the velocity fields from Emsellem et al. (2004) and Cappellari et al. (2011). In the case of the ePN.S galaxies the procedure is applied to the PN velocity fields, while for the central regions we use the V rot and σ from kinemetry analysis on IFS data from Krajnović et al. (2008); Krajnović et al. (2011);Foster et al. (2016), when available. In the other cases we use V rot and σ from major axis slits (references in the ePN.S paper).

Methods: IllustrisTNG photometry and kinematics
In this section we describe the method for measuring photometry and kinematics in the IllustrisTNG galaxies. For each simulated galaxy we define a coordinate system (x, y, z) aligned with the axes of the simulation box, and centered at the position of the most bound particle in the galaxy. Galaxies are observed both edge-on and along a random fixed line-of-sight (LOS) direction. The edge-on projection is obtained by rotating the particles according to the principal axes of the moment of the inertia tensor I i j where the sum is performed over the 50% most bound stellar particles; x n,i is their coordinates, M n their mass. The random LOS direction is arbitrarily chosen to be the z axis of the simulation box. In this work we will indicate with the lowercase letters x i , v i , and r i the 3D coordinates, velocities and radii, and we reserve capital letters for the corresponding 2D quantities projected on the sky. The coordinate r indicates the intrinsic semi-major axis distance, while R indicates the projected semi-major axis distance.
For any projection, we rotate the galaxies so that the X axis corresponds to the projected major axis, and the Y axis to the projected minor axis. This is done by evaluating the inertia tensor in Eq.2 using the 2D projected coordinates, and summing over the 50% most bound particles. We choose to weight quantities by the mass and not by luminosity, as the former are not affected by uncertainties from stellar population modeling and attenuation effects e.g. from dust. The difference between mass weighted and luminosity weighted quantities, such as in the K band, is generally small for old stellar populations (e.g. Forbes et al. 2008). Radial profiles are shown in units of effective radii R e , which are evaluated as described in Sect. 5.1.

Intrinsic shapes
The three-dimensional intrinsic shapes of the galaxies are evaluated by diagonalizing the inertia tensor I i j in Eq. (2), summed over stellar particles enclosed in elliptical shells. This definition of I i, j without any weight factors is shown by Zemp et al. (2011) to be the least biased method for measuring the local intrinsic shape of a distribution of particles, and we refer to their work for a detailed description of the procedure.
In brief, the galaxies are divided in spherical shells of radii r and r + ∆r. In each shell we calculate the tensor I i, j : the square root of the ratio of its eigenvalues give the axis ratios p and q (with p ≥ q) of the principal axes, the eigenvectors their directions. The spherical shell is subsequently deformed to a homeoid of semi-axes a = r, b = pa and c = qa. We repeat the procedure iteratively until the homeoid is adjusted to the iso-density surface, and the fractional difference between two iteration steps in both axis ratios is smaller than 1%. The values of p and q as functions of the principal major axis length r give the intrinsic shape profiles of the galaxies. We require a minimum number of 1000 particles in each shell as suggested by Zemp et al. (2011), which assures small errors from particle statistics, and, at the same time, the possibility of measuring intrinsic shape profiles out to at least 8R e for ∼ 96% of the selected TNG galaxies. The directions of the principal axes of the galaxies as functions of the galactocentric distance r are given by the eigenvectorsê j (with j = a, b, c) of the inertia tensor.  Fig. 1. Photometric measurements. Top: intrinsic shape of an example galaxy from TNG100 as a function of the intrinsic major axis distance r/R e ; axis ratios and triaxiality parameter are shown in separate panels. This simulated galaxy is oblate with q ∼ 0.4 and p ∼ 0.95 in the central 3R e , and becomes near-prolate (T > 0.7) in the outer halo. In this galaxy 9r soft correspond to 0.48R e . Bottom: ellipticity and photometric position angle profiles for two example galaxies from TNG100 as a function of the projected major axis distance R/R e . The quantities derived from the inertia tensor are shown with solid symbols, those from the mock images with open symbols.
We also use the triaxiality parameter to quantify the intrinsic shape. In App.B we find that shape measurements at 1R e are affected by the resolution of gravitational forces only for the lowest mass galaxies, for which the absolute error on p and q is ∼ 0.1. At r ∼ 9r soft , i.e. r ∼ 3.5R e for the lowest mass galaxies, and r ∼ 1.16R e for M * = 10 11 M , these resolution effects are negligible, and the error on the shape measurements is then due to particle noise and is ∼ 0.02. This uncertainty translates into an error of ∆T = 0.2 on the T parameter for typical values of the axis ratios in fast rotator ETGs (i.e. p = 0.9 and q = 0.5). As discussed in App.B, we consider the triaxiality profiles reliable starting from r = 9r soft ; at smaller radii, where ∆T is larger, we quantify shapes using p and q which are better defined.
In the paper we will consider halos as near-oblate when T ≤ 0.3, and near-prolate when T > 0.7. Halos with intermediate values of T parameter are designated as triaxial. Figure 1 (top panel) shows the principal axis ratios q(r) and p(r) as a function of the major axis distance r for one example TNG galaxy, normalized by the R e of the edge-on projection. The galaxy shown in the example is close to oblate in the central regions, with q(1R e ) = 0.43 and p(1R e ) = 0.95 (T < 0.3). At large radii the galaxy becomes close to prolate with q(10R e ) = 0.76, p(10R e ) = 0.83, and triaxiality parameter T > 0.7. For the galaxy shown 9r soft = 0.48R e .

Ellipticity and photometric position angle profiles
Mass weighted photometry is derived by diagonalizing the 2D inertia tensor (Eq. (2)) using the projected coordinates for a given LOS. We use an iterative procedure similar to the one described in Sect. 3.1 for the 3D intrinsic shape. The square root of the ratio of the two eigenvalues of I i j gives the projected flattening, hence the projected ellipticity ε(R); the components of the eigenvectors define the photometric position angle PA phot (R). The zero point of the PA phot (R) is chosen to be the X axis of the galaxies.
As an independent check on the results, we derived ε(R) and PA phot (R) also from fitting ellipses to mock images of the galaxies, and obtained very similar results. The bottom panels of Fig.  1 shows the ε(R) and PA phot (R) profiles obtained from the inertia tensor (solid symbols) and from the images (open symbols) for two example galaxies. The galaxy TNG100-511175, shown with blue symbols, is the same as the one shown in the top panels of Fig. 1: the increased axis ratio q(r) at r ≥ 3R e is reflected in a decreased projected ellipticity. The example also shows that at low ellipticities the uncertainty on the measured PA phot (R) becomes larger, as is well known. We quantified that our method allows us to measure reliably position angles down to ellipticities ε = 0.1, where the error ∆PA phot (R) from particle noise is ∼ 6 • . Below 0.1 ∆PA phot (R) increases exponentially when ε decreases towards 0.

Central kinematics
For each TNG galaxy we build projected mean velocity and dispersion fields for two projections (edge-on and random LOS). We use a resolution of 0.2 kpc, which corresponds to 2 arcsecs for a galaxy observed at 20 Mpc, comparable to present day IFS surveys (e.g. Law et al. 2016, for MANGA). The stellar particles are binned on a regular spatial grid centered on the galaxy and 8R e wide.
The binned data are then combined into Voronoi bins as described in Cappellari & Copin (2003), so that each bin contains at least 100 stellar particles. In each i-th bin we calculate the projected mean velocity and the mean velocity dispersion as the weighted averages: where the index n runs over the particles in the bin, and N i is the number of particles in the i-th bin. The top panel of Fig. 2 shows the result for one example galaxy; the middle and bottom panels show the halo kinematics and the derived kinematic parameters as described in the next section. The example illustrates that in the central regions, where the density of particles is highest, the velocity field is sampled at the highest resolution. At larger radii the Voronoi bins combine the data in progressively larger bins in order to reach the required minimum number of particles. The systemic velocity of the galaxy is derived by fitting a harmonic expansion as in Pulsoni et al. (2018, their section 4) to the central regions (i.e. at R ≤ 2R e ) of the projected velocity field. The fitted constant term is then subtracted from the velocity fields V i .
From the velocity fields we calculate the angular momentum parameter λ e following the definition of Emsellem et al. (2011) where the weighting with the flux is substituted here with a weighting with the mass M i of each Voronoi bin of index i, M i = n M n,i , and R circ,i is the circular radius of the i-th bin. The cumulative λ e is derived by summing over all the Voronoi bins contained inside an elliptical aperture of semi-major axis R e and flattening given by the ellipticity ε(1R e ). By comparison, the differential λ(R) is summed in elliptical shells.

Halo kinematics
The mean velocity and velocity dispersion fields at large radii are derived using the adaptive smoothing kernel technique (Coccato et al. 2009), used by Pulsoni et al. (2018) to derive halo velocity fields from the discrete velocities of planetary nebulae in the ePN.S survey. For the simulated galaxies, the discrete velocities of the particles at R > 2R e are smoothed with a fully adaptive kernel (A = 1, B = 0), and their stellar masses are included in the weighting. We verified that the kinematic measurements from the adaptively smoothed and the Voronoi binned velocity fields return consistent values in the regions of spatial overlap. The bottom panel of Fig. 2 shows the rotation velocity V rot , kinematic position angle PA kin , and velocity dispersion σ profiles derived from the Voronoi binned velocity fields (in orange), and from the smoothed velocity fields (in blue). V rot and PA kin are derived from fitting a harmonic expansion as in Pulsoni et al. (2018), and σ is azimuthally averaged in elliptical annuli whose flattening follows the ellipticity profile of the galaxies. The zero point of PA kin is defined to be the X axis of the galaxies, consistently with the zero point of PA phot . Error bars on the σ(R) profiles are derived from the standard deviation of the σ values inside each annulus. The values obtained with the smoothed velocity fields are very well consistent with those from the Voronoi binned velocity fields. We also evaluated differential λ profiles using Eq. (5), where the summation is performed over the Voronoi bins and the particles, each weighted by their mass, in elliptical annuli.
The example galaxy shown in Fig. 2 has a massive disk (q 0.4, see Fig. 1, top panels) embedded in a spheroidal halo with high T (T 0.7). The variation in intrinsic shape from nearoblate in the center to strongly triaxial at large radii is accompanied by a modest photometric twist (Fig. 1, bottom panels), and a much larger kinematic twist (Fig. 2) which follows the rotation along the projected minor axis visible in the top panel. At the same radii the rotation velocity V rot is observed to drop, together with the local λ parameter.

Selection in color and mass
The purpose of this paper is to study the stellar halos of a volume-and stellar mass-limited sample of simulated ETGs, and compare with observations. Nelson et al. (2018) verified that TNG100 reproduces well the (g − r) color of 10 10 < M * /M < 10 12.5 galaxies at z = 0, by comparing with the observed distribution from SDSS (Strateva et al. 2001). They also showed that redder galaxies have lower star formation rates, gas fractions, gas metallicities, and older stellar populations, and that they correspond to earlier morphological types (their Figure 13).  Fig. 3. Selection of galaxies in color and stellar mass. Top: the g − r color -stellar mass diagram of the simulated galaxies in TNG50 and TNG100. Red sequence galaxies are selected above the tilted dashed line, in the mass range 10 10.3 < M * /M < 10 12 . Bottom: ETGs from recent IFS surveys as indicated. Most of the observed ETGs in this mass range fall in the same red sequence region.
Thus we extract our sample of ETGs from the TNG50 and TNG100 snapshots at z = 0 in the color-stellar mass diagram, isolating galaxies in the red sequence. To obtain a sample of galaxies in the same area occupied by the Atlas3D and the ePN.S samples (see Sect. 2), we choose (g − r) ≥ 0.05 log 10 (M * /M ) + 0.1mag (6) We do not include any dust extinction model in the calculation of the simulated colors in order to avoid the contamination from dust-reddened late type galaxies. Even in this case, this sample of simulated galaxies unavoidably contains some red disks, while in Atlas3D some of the disks have been removed (see Sect. 2). We limited the sample stellar mass range to 10 10.3 ≤ M * ≤ 10 12 M . This choice assures that the TNG100 galaxies are resolved by at least 2×10 4 stellar particles. By comparison, the minimum number of stellar particles in the selected TNG50 galaxies is 36×10 4 . In addition we impose that the galaxies' effective radius (see Sect. 5.1) R e ≥ 2r soft , to guarantee that the region at r = R e is well resolved for all simulated galaxies. For TNG100 r soft = 0.74 kpc at z = 0, which excludes 38 galaxies at the low mass end (see Fig. 7). In TNG50 all the galaxies have R e > 2 × r soft , where r soft = 0.288 kpc. These criteria select a sample of 2250 galaxies in TNG100 and 168 galaxies in TNG50. Figure 3 shows the color-stellar mass diagram for the simulated galaxies from TNG100 and TNG50, and for observed . Stellar mass function of the TNG galaxies compared with those of the Atlas3D and ePN.S surveys. Top: stellar mass function of the TNG galaxies selected in color and stellar mass, and whose effective radii are well-resolved, as described in Sect. 4.1. The stellar mass functions of both TNG50 and TNG100 agree well with Atlas3D. Bottom: stellar mass function of the final sample of ETGs, selected in color, stellar mass, and intrinsic shape as described in Sect. 4.1 and Sect. 4.2. The removal of centrally elongated objects mostly changes the low-mass part of the TNG50 mass function, while that of TNG100 is still in good agreement with Atlas3D. galaxies from several IFS surveys. Our selection criteria are highlighted with dashed lines. Most of the observed ETGs, including the SAMI galaxies and the MANGA ellipticals and lenticulars are in the selected region of the diagram.
The histograms in the top panel of Fig. 4 show the stellar mass functions for the color-mass-selected samples. The bottom panel instead shows the stellar mass functions of the final samples as defined by adding constraints from the lambda-ellipticity diagram in Sect. 4.2. The red and hatched histograms show the Atlas3D and ePN.S samples, respectively. Here we consider the Atlas3D sample properties to validate our selection criteria, as this survey is especially targeted to study a volume-limited sample of ETGs. The ePN.S sample, which will be used to compare with properties at large radii, is also shown, and it contains on average higher mass galaxies. Both TNG50 and TNG100 are in reasonable agreement with Atlas3D. We remark here that a more generous color selection including bluer galaxies would produce a too large number of high ellipticity galaxies especially in TNG50.
In the following, whenever we compare simulated and observed galaxy samples, we will apply to the observed galaxies the same color and stellar mass selection criteria that we used for the TNG sample. 4.2. Selection of ETGs in the λ-ellipticity diagram: fast and slow rotators Figure 5 shows the λ e -ε(1R e ) diagram for the simulated ETGs in three stellar mass bins, and compares with observed ETG samples. The top row features the diagram for the TNG50 (crosses) and TNG100 (circles) galaxies selected as described in Sect. 4.1, and projected along a random LOS. The middle row shows again the TNG50 and TNG100 galaxies after the additional selection discussed in this section. The bottom row shows the similar diagram for the observed ETG samples (selected in various ways as described in Sect. 2), in the same color and stellar mass region as defined in Sect. 4.1. Here we also include for comparison the spiral and irregular galaxies from the MANGA sample (marked as LTGs).
We observe that a significant fraction of the TNG galaxies shown in the top row populate a region to the right of the λ eε(1R e ) diagram where there are no observed counterparts, i.e. below the magenta line and with ε(1R e ) > 0.5. By color coding the galaxies according to their intrinsic axis ratios at r ∼ 1R e , we find that these galaxies have elongated, triaxial shapes. These systems occur at all values of λ e , i.e. some rotate as rapidly as the MANGA disk galaxies, but others do not show any rotation It is possible that some of the rapidly rotating elongated systems are barred galaxies. Rosas-Guevara et al. (2020) showed that within a dynamically selected sample of disk galaxies the TNG100 simulation produces barred systems in fractions consistent with observational results. The majority of these systems, all characterized (per definition) by high rotation, are quenched and hence will overlap with the colour range of our sample of red galaxies. Some barred galaxies are also expected to be present among the observed ETG samples. For example, in the Atlas3D sample 7% of the galaxies show a clear bar component. For these objects the ε-values shown in Fig. 5 were measured at larger radii, to avoid the influence of the bar on the estimate of ε (Emsellem et al. 2011). However, if their actual ε(1R e ) values were used and placed these objects in the region of the λ e -ε(1R e ) populated by the centrally elongated (at r ∼ 1R e ) TNG galaxies, their fraction would not be large enough to explain the abundance of simulated galaxies in the same region, and none of these have λ e < 0.2. Therefore the presence of a large fraction of centrally elongated galaxies with high ellipticity and intermediate to low λ e in the TNG sample cannot be explained as a simple sample selection bias (note also that resolution effects on the intrinsic shapes at 1R e are at most of the order of 0.1, for the low mass galaxies; see App.B).
In App.C we discuss the properties of these galaxies further, and suggest that they are likely a class of galaxies that are produced by the simulation but are not present in nature. These galaxies occupy a particular mass range that depends on resolution and they are the reddest systems for their mass. We found no similar concentration of elongated systems among the red galaxies in the Illustris simulation, indicating that the galaxy formation model is involved in the production of these centrally elongated galaxies in TNG. The elongated components typically extend up to 3 R e and are embedded in near-oblate spheroids with a wide range of flattening q, with lower median value in TNG50 (q ∼ 0.3) than in TNG100 (q ∼ 0.45), indicating a relation to disk building and bar instability. However, some of these do not contain a disk component (Fig.C.1), and they populate a wide range of rotation (λ e ) approximately uniformly all the way from edge-on λ e = 0.7 to no rotation (Fig.C.2). Therefore we suggest that the centrally elongated galaxies in TNG may be systems that   5. λ e -ε(1R e ) diagrams for TNG and observed galaxy samples, in stellar mass bins. Top row: the sample of TNG50 and TNG100 ETGs, selected by color and mass as in Fig. 3, observed along a random LOS. The galaxies are color coded according to their intermediate to major axis ratios p at 1Re. Middle row: the λ-ellipticity diagram for the TNG galaxies after the additional selection p(1R e ) ≥ 0.6. Bottom row: λ-ellipticity plots for observed ETG samples as shown in the legend, including also late-type galaxies from the MANGA survey for comparison. Their λ e and ε(1R e ) are as described in Sect. 2. The solid black line is the threshold separating fast and slow rotators as defined in Emsellem et al. (2011): The magenta line shows the semi-empirical prediction for edge-on axisymmetric rotators with anisotropy parameter δ = 0.70ε intr from Cappellari et al. (2007). After the additional removal of centrally elongated systems, the colour and mass selected TNG ETG samples give a good representation of the observed ETG samples in the λ e -ε(1R e )-plane, with the exception of the high-λ e (> 0.7) S0 disks specific to the MANGA survey.
were in the process of forming a disk, whose evolution has been interrupted or derailed by rapid dynamical instability, star formation, and feedback in the simulations, in the particular mass range in which they occur.
For these reasons we exclude the centrally elongated objects from our sample of galaxies. We do this by performing a selection in intrinsic shape, and reject all galaxies with intermediate to major axis ratio p < 0.6 at r ∼ 1R e . This choice is motivated by the fact that the intrinsic shape distribution of real galaxies is known (Weijmans et al. 2014;Foster et al. 2017;Li et al. 2018;Ene et al. 2018) although with large uncertainties (Bassett & Foster 2019), and galaxies with p < 0.6 are rare, even among the the slow rotators. By applying this selection criterion we obtain our final sample of simulated ETG galaxies, 1114 objects in TNG100 and 80 in TNG50.
The middle row of Fig. 5 shows that the final selected sample of ETGs lies in the region of the diagram populated by the observed galaxies. The fraction of simulated galaxies in the region of avoidance (i.e. above the black and below the magenta lines) is in agreement with the λ e − ε observations. The loca- tion of the simulated galaxies in the plane follows closely the Atlas3D, SAMI, and MASSIVE galaxies. In the MANGA sample there is a large fraction S0 galaxies with λ e > 0.7 that are not present in the other surveys, and are likely due to differences in the data analysis, possibly to the beam corrections applied by Graham et al. (2018) on the MANGA data (see discussion in Falcón-Barroso et al. 2019). Figure 6 demonstrates that the distribution of stellar halo properties which we are interested in, i.e. λ and triaxiality parameter, are not affected by the sample selection based on p(1R e ). The distributions do not systematically depend on the intrinsic shape of the central regions of the galaxies. As discussed in a companion paper, the properties of the galaxies at large radii are mainly set by their accretion history and not by the details of the star formation in the central regions of galaxies.
The bottom panel of Fig. 4 shows the stellar mass function for the final sample of ETGs, compared with observations. The stellar mass function of the TNG100 ETGs is still similar to Atlas3D. For TNG50 the additional selection has excluded a large fraction of red galaxies in the stellar mass range ∼ 10 10.3 − 10 11 M . This results in a stellar mass function skewed towards high masses (and so more similar to ePN.S).
Henceforth we classify galaxies as slow rotators (SRs) and fast rotators (FRs), using the dividing line introduced by Emsellem et al. (2011), shown in Fig. 5 with the black line: galaxies above this threshold are FRs, and galaxies below are SRs. To reduce the effects of inclination, we choose to classify the simulated galaxies using the values of λ e and ellipticity for their edge-on projection (shown in Fig. C.2).

Summary of the sample selection criteria
The sample of ETG galaxies used in the remainder of this paper is extracted from the TNG50 and TNG100 simulations by selecting galaxies in the stellar mass range 10 10.3 ≤ M * ≤ 10 12 M and with red (g − r) color as in Eq. (6) (see Fig. 3); excluding a small number of objects with R e < 2r soft , to assure sufficient resolution at 1R e ; finally, removing a class of centrally elongated, triaxial galaxies with p(1R e ) < 0.6, which are systems not present in the observed ETG samples, that probably became barunstable and quenched during the process of (central) disk formation.
The selected sample has a distribution of λ e − ε(1R e ) similar to observed ETGs (Fig. 5) and halo properties that are unbiased by the selection in intrinsic shape (Fig. 6). The mass functions of the selected TNG50 and TNG100 ETG samples are given in Fig. 4, and the mass-size relations are shown in Fig. 7.

Photometric properties of the TNG ETG samples
In this section we study the photometric properties of the selected sample of TNG galaxies and how they vary with radius. Section 5.1 discusses the measured galaxy sizes and how our definition of effective radii compares with effective radii inferred from ETG photometry. Section 5.2 compares the distribution of projected ellipticities at 1R e with that from ETG surveys and validates our sample selection. Section 5.3 studies the TNG ellipticity profiles out to the stellar halo, Sect. 5.4 explores the intrinsic shape distribution of stellar halos and its dependence on stellar mass, and Sect. 5.5 investigates the dependence of galaxy triaxiality on radius and stellar mass in the simulated samples. Finally 10.50 10.75 11.00 11.25 11.50 11.75 12.00 10 1   8. Ellipticity distribution and profiles of the selected TNG ETGs. a) Distribution of the random LOS ε(1R e ) compared with Atlas3D and ePN.S as indicated in the legend. Despite some differences between the Atlas3D and the TNG SRs, the good agreement of FR distributions shows that the selected sample of TNG galaxies contains a mixture of disk and spheroidal galaxies consistent with Atlas3D. b) Ellipticity profiles for the ePN.S galaxies (top) and 10 randomly selected example galaxies from TNG100 and 10 from TNG50, projected along a random LOS (bottom). The comparison between TNG and ePN.S galaxies highlights the variety of projected ellipticity profiles in both samples. c) Distribution of the edge-on projected ellipticity values at different radii predicted by TNG. SRs have a rather constant ellipticity distribution with radius. The FRs have a large variety of shapes at large radii as shown by the broadening of the distribution, while the shift of the peak to lower ellipticities shows the tendency of most galaxies to become rounder in their halos.
Sect. 5.6 tests the ability of photometric twist measurements to establish the underlying triaxiality in the TNG galaxies.

Sizes of the TNG galaxies
We first discuss the adopted measurement of the effective radius for the simulated galaxies, which we will use in the paper as galactocentric distance unit.
The effective radius R e is derived for each projection (edgeon or random LOS) of the galaxies by using cumulative mass profiles in elliptical apertures: R e is the major axis radius of the aperture that contains half of the total bound stellar mass. Figure 7 shows the circularized R e as a function of M * for the final samples of ETGs, and compares it to the distribution of observed effective radii from the different surveys. The R e in TNG100 are larger than most of the observed R e at M * 10 10.75 , but they are in reasonable agreement with the ePN.S measurements. The higher resolution TNG50 produces smaller galaxies compared to TNG100 and to observations at intermediate stellar masses; see also Pillepich et al. (2019). On the other hand, comparisons to observed R e strongly depend on the operational definitions of galaxy sizes, as discussed by Genel et al. (2018).
Observers measure R e by integrating light profiles fitted to the bright central regions to large radii. This definition of R e tends to underestimate the size (and at the same time the total stellar mass) of the galaxies if the photometric data are not deep enough to sample the light distribution in the halos, especially in massive galaxies with high Sérsic indices. Pulsoni et al. (2018) determined R e of the ePN.S galaxies from the most extended photometric profiles available in the literature, using a Sérsic fit of the outermost regions to integrate to large radii. This approach leads to an average increase of the R e by a factor of ∼ 2 for the most massive objects with M * > 10 11 M . However it does not take into account the possibility of an extra halo component/intra-group or intra-cluster light (ICL) at large radii. For the simulated galaxies, defining the stellar content of the galaxies as all the bound stellar particles identified by the subfind algorithm, automatically includes also ICL stars in the most massive halos, thus overestimating both R e and M * .
To quantify these effects requires separating a galaxy from the surrounding ICL. A kinematic separation of the ICL similar to Longobardi et al. (2015) is beyond the scope of this paper. However, Kluge et al. (2020) recently found that if the ICL component in bright cluster galaxies is identified as the outer component of a double Sérsic fit, the radius at which it starts dominating is ∼ 100 kpc with a very large scatter (5 to 400 kpc in their figure 16). We evaluated the differences in R e and M * that we would obtain if instead of using the whole bound stellar mass we limit the galaxy to the mass within 100 kpc. We find that in TNG100 galaxies with M * < 10 10.5 the effects are negligible; in 10 10.5 M < M * < 10 11 M objects the differences in the derived R e and stellar masses are within 10% and 5% respectively, while between 10 11 M < M * < 10 11.5 M they are within 30% and 15%. At higher masses the differences in R e can be larger than 50% and those in M * larger than 25%, with a very large scatter. These effects are half as pronounced in TNG50. A size-stellar mass diagram analogous to Fig. 7 using the 100 kpc aperture instead of the total bound mass shows an improved agreement with the observed R e , but TNG100 galaxies with M * > 10 10.75 are still larger on average. This may indicate that TNG100 predicts too large sizes for high mass galaxies (see also Genel et al. 2018).
Because of the somewhat arbitrary choice of the 100 kpc limit on one hand, and the uncertainties in the observed R e distribution on the other (from differences in sample selection, quality of the photometric data, definition of total stellar light, the massto-light ratio to obtain total stellar masses), we choose here to define R e for the simulated galaxies as the half mass radius of the total bound stellar mass, and consider the above uncertainties in the discussion of the results where relevant. Figure 8a shows the distributions of the ellipticities measured at 1 R e for the final sample of selected ETGs, compared with Atlas3D and ePN.S. In the top panels are the SRs, and in the bottom the FRs.

Ellipticity distribution in the central regions
The TNG50 and TNG100 simulations predict a significant fraction of SR galaxies with ε(1R e ) > 0.4, while the observed SRs are relatively rounder. This is a common feature of current simulations (Naab et al. 2014;Schulze et al. 2018), and its origin is still to be understood. In the case of the FR class, the ellipticity distributions are rather flat-topped, and in good agreement with Atlas3D. By comparison the ePN.S sample contains on average rounder (and also more massive, see the bottom panel of Fig. 4) galaxies, and hence a lower number of disk galaxies: none of the ePN.S FRs has ellipticity higher than 0.7.
Overall Fig. 8a shows that the selected sample of ETGs contains a mixture of galaxy types consistent with observations, with a similar balance between disks and spheroids.

Ellipticity profiles
The ellipticity profiles of the TNG galaxies are compared over an extended radial range with those of the ePN.S galaxies in Fig.  8b. There we show profiles for randomly selected sub-samples of the simulated galaxies. Figure 8c instead shows the distribution of ellipticities at different radii for the fast and the slow rotators separately.
The observed profiles for the ePN.S SRs generally mildly increase with radius, reaching ε ∼ 0.3 at 4R e . By comparison, the simulated SRs have more nearly constant ellipticity profiles. This can also be seen in the histograms of Fig. 8c where the ε distribution is almost unvaried between different radii.
Most of the simulated FRs have decreasing ellipticity profiles with radius, while a fraction have high ellipticity also at large radii, as also shown by the ePN.S galaxies. Thus, Fig. 8c shows that at larger radii the FR ellipticity distribution peaks at smaller ε and, at the same time, it broadens.
The decrease in ellipticity of the majority of the FRs supports the idea of a change in structure of these galaxies at large radii. The large range of flattening in the stellar halos indicates a variety in the stellar halo properties. By comparison, the SRs show only small structural variations.

Intrinsic shape distribution of the stellar halos
In this section we quantify the distribution of the galaxy intrinsic shapes at large radii. Here we refer to stellar halo as the outer regions of the galaxies, where the physical properties are markedly different from those of the central regions. While this region may begin at different radii in each ETG, we will see in the next section that the median triaxiality profiles for our sample reach con- r = 9r soft for TNG100 galaxies with M * ∼ 10 10.3 M , the solid lines show r = 9r soft for TNG100 galaxies with M * ∼ 10 11 M . At radii larger than these the profiles are not affected by resolution issues. Right panels: The median of the stellar halo triaxiality distribution measured at 8R e as a function of stellar mass. The shaded areas are contoured by the 25th and 75th percentiles of the distribution. TNG100 and TNG50 are shown separately, FRs are on top, SRs on the bottom. On average the triaxiality parameter increases with radius and with stellar mass. The FRs have increasing T (r) profiles, while SRs have either constant or decreasing profiles. 7 instead of 8R e ) at which we measure p and q deliver similar results.
The top panel of Fig. 9 shows the minor to major axis ratio q as a function of the intermediate to major axis ratio p: we find a large variety of possible shapes, from very flat near-oblate with q ∼ 0.3, to prolate with q ∼ p ∼ 0.5. The majority of (lowmass) galaxies appear to have near-oblate stellar halos, with a large scatter in minor to major axis ratio q. The bottom panels of Fig. 9 show the intrinsic shape distributions for the fast and slow rotators separately. The distributions of the minor to major axis ratio q resemble Gaussians, and fitted as such the FRs have mean µ q ∼ 0.5 and dispersion σ q ∼ 0.16 in all stellar mass bins. The SRs have µ q ∼ 0.6 and σ q ∼ 0.15, with a tendency for the highest mass galaxies to be flatter.
The distribution of the intermediate to major axis ratio p can be approximated by a log-normal distribution in Y = ln(1 − p).
The shape of this distribution shows a dependence on stellar mass: at higher stellar masses µ Y increases, together with the width of the distribution. This means that at higher stellar masses, in both the FR and SR classes, the fraction of near-oblate galaxies with p ∼ 1 decreases.
The vertical dashed lines in Fig. 9 shows the comparison with the photometric model used by Pulsoni et al. (2018) to reproduce the distribution of maximum photometric twists versus mean ellipticity of the observed FRs. Their model parameters µ q = 0.6 and µ p = 0.9 are within 1σ of the mean values obtained from the distribution of simulated FRs.

Triaxiality profiles
We can study how the intrinsic shapes of galaxies vary as a function of radius by looking at their triaxiality profiles. We recall from Sect. 3.1 that because of the error due to resolution effects in the central regions, T profiles are considered well-defined only beyond r = 9r soft , where their error ∆T ∼ 0.2 for typical FR axis ratios (App.B). Thus we show profiles only for r > 3.5R e for the lowest mass objects, and for r > 1.16R e for galaxies with M * ∼ 10 11 M .
The left panels of Fig. 10 show the median triaxiality profiles for FRs and SRs. These median profiles were built by binning the galaxies according to their values of the triaxiality parameter T at 8R e , and are plotted against the intrinsic major axis distance r. The scale radius R e that normalizes the three dimensional radius is the 2D projected effective radius for the edge-on projection of each galaxy. The right panels show the median of the distribution of the triaxiality parameter measured at 8R e as a function of the stellar mass.
FRs are characterized by increasing T profiles, which tend to plateau at r > 5R e where the TNG galaxies show a broad range of intrinsic shapes despite all having near-oblate centers. SRs tend to have flatter profiles.
We find that the stellar halo intrinsic shape distribution is a function of stellar mass. This is visible in the right hand side of Fig. 10, for FRs and SRs separately. At lower masses the TNG galaxies have preferentially near-oblate shapes, with T 0.2, but at larger masses the median triaxiality parameter increases, so that at M * > 10 11 M there is a non-negligible fraction of galaxies with prolate-triaxial halos, even among the FRs. We note a systematic difference between TNG50 and TNG100 in the triaxiality of the SR stellar halos. In TNG50 the SRs tend to be much more oblate, indicating some higher degree of dissipation involved in their evolution compared to the SRs in TNG100. The statistical significance of this difference is marginal since TNG50 contains only 14 SRs.

Photometric twists and triaxiality in TNG ETGs
Isophotal twists in photometry are generally considered to be signatures of triaxiality. This is because the projection on the sky of coaxial triaxial ellipsoids with varying axis ratios approximating the constant luminosity/mass surfaces of ETGs can result in twisting isophotes (e.g. Benacchio & Galletta 1980). However, the effects of triaxiality on the PA phot profiles are model dependent, that is they depend on axis ratio, on how much the axis ratios changes with radius, as well as on the viewing angles. Figure 11 shows the distribution of maximum photometric twist, i.e the maximum variation of PA phot (R), measured between 1 and 8 R e , as a function of the halo triaxiality at 8 R e , i.e. where the triaxiality profiles have reached a constant value. Each symbol in the diagram is color coded by the median projected ellipticity between 1 and 8 R e . Galaxies with ellipticity lower than 0.1 have the photometric position angle poorly determined, and are shown with smaller symbols.
We observe that the amplitude of the photometric twists is only weakly dependent on the triaxiality. Near-oblate and nearprolate galaxies are slightly less likely to have constant PA phot than triaxial galaxies, but the majority of the galaxies have small twists irrespective of T (8R e ). This is explained by the fact that large twists can be measured for viewing angles close enough to face-on (Pulsoni et al. 2018, that is lower ellipticities in Fig. 11), at which even small values of T can produce large twists. From Fig. 11 we conclude that the amplitude of the photometric twists is a poor indicator for galaxy triaxiality, and that very small photometric twists are intrinsically compatible with triaxial shapes.  Fig. 11. Maximum measured photometric twist between 1 and 8 R e as a function of the halo triaxiality T (8R e ) for the TNG50 and TNG100 galaxies projected along a random LOS. FRs and SRs are shown with different symbols as in the legend. Smaller symbols are used to represent galaxies with ellipticity < 0.1, for which the errors on PA phot are large. The color-bar indicates the median ellipticity measured between 1 and 8R e . The gray solid line shows the median of the photometric twist distribution as a function of T (8R e ); the shaded area encloses the 25th and 75th percentiles. The amplitude of the photometric twists in TNG galaxies is only weakly dependent on the triaxiality parameter.

The kinematics properties
In this section we study how the kinematic properties of the TNG galaxies vary with radius. In Sect. 6.1 we derive median differential λ(R) profiles to quantify the variety of kinematic behaviors in the outskirts of FRs and SRs. Sections 6.2 and 6.3 compare the shapes of the V rot /σ profiles of the simulated ETGs with the observed galaxies in the Atlas3D and ePN.S surveys. Finally Sect. 6.4 uses the simulated galaxies to assess kinematic misalignments and twists as signatures of triaxial shapes in galaxies.

Lambda profiles
The top panels of Fig. 12 show the median differential λ profiles for FRs and SRs in their edge-on projection. Galaxies are binned together according to the shape of their profiles. We achieve this by binning the FRs according to their values of λ at R = 4R e and at R = 8R e . For the SRs the shape of the λ(R) profiles is generally a monotonic function of R: in this case we binned the profiles according to their λ(7R e ).
Most of the FRs reach their maximum λ(R) around 3R e ; only 7% of the galaxies increase λ(R) between 3 and 10 R e . FRs divide almost evenly among a third (34%) that have flat λ profiles with λ(8R e ) > 0.6, a third (40%) with gently decreasing profiles and 0.3 > λ(10R e ) ≥ 0.6, and another third (26%) with very low rotation in the halo (λ(10R e ) ≤ 0.3).  The SRs essentially divide between a half (53%) with nonrotating halos (λ(7R e ) < 0.2) and a half with increased rotational support at large radii compared to the central regions. We observe that a small fraction of the SRs (5% of the TNG100 SRs and 15% of the TNG50 SRs) reach very high values of λ at large radii (λ(7R e ) > 0.6). The majority of these galaxies are genuine slow rotators with strongly rotating halos similar in terms of velocity fields and λ profiles to observed SRs like NGC 3608 (Pulsoni et al. 2018). The others are galaxies with a clear extended disk structure characterized by rapid rotation and low velocity dispersion, but whose central kinematics is dominated by a nonrotating bulge. There are no observed counterparts for the latter in both the ePN.S (33 galaxies) and the SLUGGS surveys (Foster et al. 2016, 25 galaxies, of which 18 are in common with ePN.S). A larger sample of galaxies with extended kinematic data would be needed to confirm or rule out these objects.
For the SRs we note a mismatch in the amount of halo rotational support between TNG50 and TNG100, analogous to the difference observed for the halo triaxiality (Sect. 5.5): on average TNG50 SR halos rotate faster and have more oblate shapes. These differences in halo properties might be due to the depen-dence of the galaxy formation model on the resolution of the simulations, although the number of SRs in TNG50 (only 14 galaxies) is too small to draw strong conclusions.
For both FRs and SRs, and in both TNG50 and TNG100, the stellar halo rotational support depends weakly on stellar mass up to M * ∼ 10 11.3 M (Fig. 12). However, at high stellar masses the fraction of galaxies with significant rotation in the halo decreases, so that at M * > 10 11.5 M most of the galaxies have non rotating halos. The broad range of possible λ(R) profile shapes in Fig. 12 shows that the IllustrisTNG simulations encompass, if not exceed, the observed variety of halo rotational support found in the ePN.S survey, of which one of the key results was the large kinematic diversity of stellar halos. Figure 13 shows the median V rot /σ profiles of the TNG100 and TNG50 galaxies compared with median profiles from Atlas3D and individual galaxy profiles from the ePN.S sample for 0 ≤ R ≤ 4R e . Here we normalize the radii by the circularized R e , i.e. R e √ 1 − (1R e ), for an appropriate comparison with the Atlas3D R e .

Simulated versus observed rotation profiles -central regions
The profiles of the simulated SRs are similar to the observed profiles. The FRs instead show a difference in how quickly V rot /σ rises with radius: observed FR galaxies have on average more steeply rising V rot /σ profiles than the simulated ETGs. Very few TNG FRs reach V rot /σ ≥ 1 within 1R e compared to the Atlas3D galaxies, and almost none exceeds V rot σ(1R e ) = 1.5 in either TNG100 or TNG50.
The different shapes of the V rot /σ profiles in observations and simulations cannot be explained with resolution effects, as in both TNG50 and TNG100 the V rot /σ profiles tend to peak at a median radius of R ∼ 3R e . By comparison, the ePN.S FRs tend to peak at smaller fractions of R e , at a median 1.3R e .
This difference between the observed and simulated FRs is not a consequence of the selection functions of the samples. The TNG galaxies are selected according to color, mass and intermediate to major axis ratio p. The selection in p removes centrally elongated galaxies, most of which have intermediate to low λ(1R e ) (Fig. 5). The Atlas3D sample, selected as described in Sect. 2, has some of the disk galaxies removed which as in MANGA (Fig. 5) will be mostly located at high λ(1R e ). We recall that in the comparison we consistently matched the color selection and mass range of the Atlas3D sample to the selection criteria adopted for the TNG galaxies. Thus the TNG sample should in principle contain a larger number of late type galaxies with strong disks, i.e. high V rot /σ, by comparison with Atlas3D. Figure 13 shows instead that the TNG100 and TNG50 ETG samples lack galaxies with high rotational support at 1R e .
In Sect. 5.1 we discussed that the effective radii used to normalize the radial scales in TNG would be expected to be systematically slightly overestimated compared to the observed R e since they are defined using the total bound stellar mass which is often inaccessible in observations. If taken into account, this effect would increase the gap between simulated and observed samples. On the other hand, since the mass-size relation is roughly reproduced in the simulations (Fig. 7), the different steepness of the V rot /σ profiles in observations and simulations implies a different distribution of the angular momentum as a function of radius in the simulated galaxies. This could be due to a too efficient condensation of the gas into stars that does not allow the gas to dissipate and collapse to small enough radii. Figure 14 compares the relation between V rot /σ in the central regions and V rot /σ in the stellar halos for the simulated and observed galaxies. The observed galaxies are from the ePN.S survey (Pulsoni et al. 2018, their figure 9) and are reported in the top left panel. Their measurements use absorption line data at 1R e and PN data for the halos, which on average cover 6R e with a large scatter (minimum 3R e , maximum 13R e ). The central and right top panels show V rot /σ(6R e ) versus V rot /σ(1R e ) for TNG100 and TNG50 ETGs, respectively. Galaxies close to the dashed 1:1 lines show similar rotational support in the central regions and in the outskirts; galaxies below the lines have increased rotational support in their stellar halos; galaxies above the lines instead have reduced rotation at large radii.

Simulated versus observed profiles -outskirts
The position of the observed SRs below the one-to-one line is reproduced by the simulations. As already discussed in Sect. 6.1 there are a few outliers among the simulated SRs with V rot /σ(6R e ) 1, some of which are actually extended disks with prominent bulges at the center. These represent a small fraction of the SR family, and do not have observed counterparts in the ePN.S and SLUGGS surveys.
For the TNG FRs we find a different distribution when comparing rotational support at the same radii: most of these galaxies fill the bottom half of the diagram, with very few reaching V rot /σ(1R e ) > 1.5, and a large fraction having significant V rot /σ at 6R e . This difference is explained by the shallower V rot /σ(R) profiles of the simulated galaxies compared to the observed FRs, as discussed in Sect. 6.2. Since the simulated galaxies tend to peak at larger radii than the observed ETGs, there are almost no objects that reach V rot /σ > 1.5 already at 1R e .
In the bottom panels of Fig. 14 we show the comparison at adjusted radii: V rot /σ(3R e ) (i.e. at the median radius where the TNG V rot /σ profiles tend to reach the maximum) versus V rot /σ(8R e ) (i.e. where the decreasing rotation profiles finally reach their minimum, see also Fig. 12), and find better agreement. Now the FRs spread in the region of the diagram above the one-to-one line as they do for the ePN.S observations in the top left panel of Fig. 14, with a sub-population of objects showing V/ rot σ(8R e ) > V/σ(3R e ). We note that FRs with nearprolate stellar halos occupy mostly the lower left corner with low V rot /σ(3R e ) and low V rot /σ(8R e ). Galaxies with triaxial halos tend to distribute on the left side of the diagram, galaxies with oblate halos tend to have higher V rot /σ(8R e ). It is interesting that, aside for the differences in the central regions at R 3R e , the observed galaxies with and without signatures for triaxial stellar halos distribute similarly as the simulated triaxial and nearoblate stellar halos, respectively, with the triaxial halos having on average lower V rot /σ(8R e ).
From Fig. 14 we conclude that even though the quantitative details between simulated and observed ETGs galaxies do not agree, the simulations do reproduce the observed kinematic transitions between central regions and outskirts, as well as the variety of halo kinematic classes. TNG 50 -random LOS projection 14. V/σ ratio in the central regions compared with the V/σ at large radii. Top: V/σ(1R e ) versus V/σ(6R e ) for the ePN.S galaxies (left), TNG100 (middle), and TNG50 (right). Bottom: V/σ(3R e ) versus V/σ(8R e ) for TNG100 (left) and TNG50 (right). Different colors in the ePN.S sample mark the SRs (in red), the FRs with kinematic signatures for triaxial halos (i.e. with kinematic twists or misaligments, in green), and FRs with PA kin aligned with PA phot (in blue). The gray points stand for ePN.S FRs for which the measured V rot /σ(< 6R e >) is a lower limit. For the simulated galaxies the different colors mark SRs, and FRs with different halo T (8R e ): FRs that have near-oblate halos (in blue); FRs with triaxial halos (in green); near-prolate FRs (in orange). The gray symbols show simulated galaxies with too few particles at 8R e for measuring T . The size of the data points is proportional to the stellar mass. The gray dashed lines show the 1:1 relation. The TNG simulations reproduce the diversity of observed ETG halo kinematics and they echo the observed kinematic transitions between the central regions and outskirts of the FR galaxies, albeit at larger radii.

Relation of kinematic misalignments and twists with triaxiality in TNG ETGs
Kinematic signatures of galaxy triaxiality can be found from observations of minor axis rotation, kinematic twists, and misalignment of the kinematic position angle PA kin with the photometric PA phot (see e.g. Binney 1985;Franx et al. 1991). Figure 15 shows the distribution of misalignments Ψ = PA kin (1R e ) − PA phot (1R e ) as a function of the ellipticity for the random LOS projected TNG galaxies. The solid lines show the unweighted root-mean-square deviation from zero of the data points as a function of the ellipticity for the FR and the SR classes separately, computed by mirroring the data points around zero (i.e. using [Ψ, −Ψ]). The dashed lines show the same quantities for observed galaxies in Atlas3D and MANGA Graham et al. 2018). Simulated and observed galaxies show very similar trends with ellipticity, with very flat galax-ies being strongly aligned, and kinematic misalignment increasing with rounder shapes. Simulated FRs are found to be much more aligned than SRs, in agreement with observations. Data points in Fig. 15 are color coded according to the intermediate to major axis ratio p at 1R e . At these radii the errors on the axis ratios are at most 0.1 for the low mass systems (see App.B), so p(1R e ) measurements are generally well defined. We find that many TNG galaxies, both SRs and FRs show a high degree of alignment (|Ψ| < 10 • ) even though they are far from being oblate (i.e. with p < 0.9). This result agrees with the analysis of Bassett & Foster (2019) of Illustris galaxies, among which they found many triaxial and prolate objects with small Ψ. Although we can not draw conclusions on the shape distribution of the real galaxies, Fig. 15 implies that near-alignment of the kinematic and photometric position angles does not exclude triaxiality for the IllustrisTNG FR galaxies even at 1R e . Simulated FRs and SRs do show kinematic position angle variations as a function of radius. An example is the object shown in Figs. 1 and 2, a galaxy with central disk embedded in a triaxial stellar halo, which shows a variation in the direction of rotation corresponding to the sudden change in the triaxiality profile around r ∼ 5R e .
The top panel of Fig. 16 shows the distribution of the maximum kinematic twists, i.e. the maximum variation of PA kin , measured between 1 and 8 R e for all the TNG galaxies, as a function of the stellar halo triaxiality measured at 8R e , i.e. where the median T profiles are constant with radius (see Fig. 10).
We find that the large majority of TNG galaxies have small kinematic twists compared to typical measurement errors of PA kin from discrete tracers (the median error for the ePN.S galaxies is ∼ 35 • ). Aside for a small group of near-oblate galaxies with counter-rotating disks (twist ∼ 180 • ), the main dependence of the amplitude of the twist with the triaxiality parameter is such that galaxies with high T are more likely to have large kinematic twists, as shown by the solid gray line in Fig. 10 representing the median twist as a function of T . This highlights the importance of kinematics as a fundamental tool to investigate the intrinsic structure of galaxies besides photometry alone (see Sect. 5.6). On the other hand Fig. 10 points out that not all the triaxial and prolate galaxies in IllustrisTNG display kinematic twists. This suggests that also in observations a galaxy's triaxiality may not be easily revealed by kinematic misalignments.
IFS kinematics show that the central regions of FR have PA kin and PA phot aligned within ∼ 10 degrees, while SRs are generally misaligned Fogarty et al. 2015;Ene et al. 2018;Graham et al. 2018). This difference in the misalignment distribution was interpreted as signature of a different intrinsic shape distribution for the two classes, with the FRs Top: maximum kinematic twist measured between 1 and 8 R e as a function of the stellar halo triaxiality T (8R e ). FRs are shown in blue, SRs in red; both TNG50 and TNG100 samples are shown. The gray line and shaded band show the median of the twist distribution as a function of T (8R e ) and the 25th and 75th percentiles. Even though the majority of the TNG galaxies have small kinematic twists, high T galaxies are more like to display a kinematic twist. Bottom: Distribution of the maximum kinematic twist for the TNG50 and TNG100 samples (green step histogram) compared with the ePN.S sample (hatched histogram). The filled green histogram shows the distribution of the TNG galaxies after their mass function is matched to that of the ePN.S sample; the gray histogram is the mass-matched TNG distribution convolved with a measurement error and is consistent with the data. being consistent with having oblate shapes, and the SRs being moderately triaxial. This interpretation is not supported by the results in Fig. 15.
Even though the central regions of FRs have been found to have PA kin well aligned with PA phot , kinematic twists are observed in the halos. By extending the kinematic study at larger radii using planetary nebulae, Pulsoni et al. (2018) found that kinematic twists are relatively frequent in the ePN.S sample of FRs (∼ 40%). They concluded that if the central regions of FRs are oblate, kinematic twists would indicate a transition to halos with triaxial shapes.
In the bottom panel of Fig. 16 we compare the distribution of kinematic twists in the IllustrisTNG galaxies with the observed distribution in the ePN.S sample. To do that we match the TNG mass function to that of the ePN.S sample by randomly selecting the appropriate fraction of TNG galaxies in mass bins. The filled green histogram in Fig. 16 shows the distribution of kinematic twists for 100 random realizations of PNS-like samples extracted from the TNG galaxies. We then convolved the resulting distribution with a Gaussian error of 35 degrees, i.e. the median error of the ePN.S measurements (gray histogram). We find that the simulated galaxies show a similar distribution as the ePN.S galaxies, although there is an indication for a lower fraction of galaxies with large kinematic twists in IllustrisTNG. This might be due to a different sample selection between simulations and ePN.S, with the former potentially containing a larger fraction of disk galaxies even after the matching of the mass functions. Figures 4 and 8a in fact show that the ePN.S sample is on average more massive and contains rounder galaxies than TNG100.

Stellar halo angular momentum and shape
In Sect. 6.3 we observed that the rotational support in the stellar halo correlates with the intrinsic shape: stellar halos with high V rot /σ are likely near-oblate, while halos with lower V rot /σ can have larger triaxiality. In this section we connect the variation of rotational support in the stellar halos to variations of their intrinsic shape.
The top panel of Fig. 17 shows the differential λ(R) (measured for the edge on projection), the axis ratio q(r), and the triaxiality parameter profile T (r) for an example galaxy. It has a decreasing λ(R) profile in the halo, while q(r) and T (r) increase: for this object the decreased rotational support marks a variation in the intrinsic shape of the galaxy from relatively flat and near-oblate at R ∼ 2R e , where the rotation is highest, to triaxial spheroidal in the outskirts. In observations the drop in rotation of the ePN.S FRs is often related to a decrease in ellipticity. This led to the idea that central fast rotating regions of FRs are embedded in a more dispersion dominated spheroidal stellar halo.
We verify this conclusion in Fig. 17. Here the outskirts of galaxies are divided into 6 ellipsoidal shells of major axis r between 3.5 to 8R e . This radial range is motivated by the requirement that r is large enough so that the errors on the T parameter are small for galaxies of all masses, but also such that most of the galaxies (96%) have enough particles at large radii (8R e ) in order to measure their intrinsic shapes. In each ellipsoidal shell of major axis r and width ∆r we measure the axis ratio q and the triaxiality parameter. Then we measure the edge on projected λ(R) in an elliptical shell of major axis R = r and width ∆R = ∆r. Figure 17 shows the halo edge-on projected λ(R) for all shells and all TNG galaxies in four galaxy stellar mass bins, as a function of the minor to major axis ratio q(r) and of the triaxiality parameter T (r).
We find indeed a relation between λ and shape. High rotational support is related to flattened (i.e. low q) near-oblate (i.e. low T ) shapes. Where λ decreases q grows towards more spheroidal shapes, although the scatter in possible stellar halo shapes is large (σ q ∼ 0.15). Stellar halos of both FRs and SRs with high triaxiality generally have low λ, but lower T is compatible with all values of λ. This means that the decrease in the rotational support of the TNG FRs follows a change in the structure of the galaxies, which become more spheroidal in the outskirts. These outer spheroidal components can span all values , axis ratio q(r) and triaxiality T (r) profiles for an example galaxy. Bottom: λ for the edge-on projection versus minor to major axis (left panels) and triaxiality parameter (right panels) in stellar mass bins. Each data point is a local measurement within a shell of major axis R = r: each TNG galaxy is represented by ∼ 6 data points measured between 3.5 − 8R e . FRs and SRs are shown with different colors as in the legend. The mass bins are labelled on the right margins. Lower rotational support λ in the stellar halos is related to rounder shapes, with a wide range of triaxiality which depends on stellar mass. FRs and SRs exhibit a continuity in stellar halo properties rather than a bimodality. of triaxiality. By comparison to the FRs, the SRs show smaller variations in both intrinsic shapes (Fig. 10) and λ. Figure 17 confirms the dependence of the stellar halo intrinsic shape and rotational support on stellar mass already described in Sects. 5.4, 5.5, and 6.1. Lower mass galaxies host more rota-tionally supported stellar halos with near-oblate shapes. At progressively higher masses the fraction of dispersion dominated stellar halos increases as well as the fraction of halos with high triaxiality.
It is interesting to note that there is no clear separation in Figure 17 between the stellar halo properties of the FRs and SRs, but rather a continuity of properties among the two classes, with the low T -high λ extreme dominated by the FRs, the high T -low λ limit by the SRs, and the relative importance of the two populations gradually changing with stellar mass. This result implies that there is no qualitative difference between the structure of the galaxies at large radii despite the bimodality of the FR/SR classification of the centers. The IllustrisTNG galaxies agree with the ePN.S observations in that FRs and SRs tend be more similar at large radii, especially at intermediate to high masses.
The similarity of the overall shapes of the λ(R) (or V rot /σ(R)) profiles and of the ε(R) profiles between ePN.S and TNG galaxies described in in Sects. 5.3 and 6.3 suggests that intrinsic shape variations similar to those found in the TNG galaxies also exist in real galaxies.

Summary and conclusions
In this paper we have analysed the kinematic and photometric properties of ∼ 1200 early type galaxies (ETGs) from the Illus-trisTNG cosmological simulations TNG100 and TNG50, with a focus on their stellar halos. The sample of simulated ETGs was selected in stellar mass and in color (Fig. 3) and in the λ e − ε diagram (Fig. 5). There we excluded simulated objects that do not match the observed properties in the central regions of fast rotators (FRs) and slow rotators (SRs), and that appear as highly centrally elongated red galaxies. We verified that this does not affect our results on the stellar halo properties of the simulated galaxies. The resulting ETG sample has mass-, size-, and central kinematics distributions consistent with observations.
For the selected sample we determined mean velocity fields, kinematic, and photometric profiles, and studied the intrinsic shapes of the simulated galaxies from the central regions into their outskirts, up to 15R e . The purpose of the paper is to study the kinematic properties of stellar halos and connect them to variations in the structural properties of their galaxies from the central regions to the halos. Our conclusions are as follows: 1) The differential λ profiles and the triaxiality profiles (Figs. 10 and 12) successfully reproduce the diversity of kinematic and photometric properties of stellar halos observed in the ePN.S survey.
-We find that simulated FRs divide almost evenly among a third with flat λ profiles, a third with gently decreasing profiles, and another third with very low rotation in the halo. Half of the SRs do not show any rotation in the halo, while the other half has increased rotational support at large radii. -FRs generally tend to show increased triaxiality with radius, although the majority (partially driven by the numerous low mass fast rotators) have stellar halos consistent with oblate shapes. -Both halo triaxiality and rotational support are found to depend on stellar mass, with higher mass galaxies being more triaxial and more dispersion dominated at large radii.
2) Halo intrinsic shape and rotational support are strongly related (Fig. 17): -High λ is related to flattened oblate shapes.
-Where λ decreases with radius, galaxies tend to become rounder, but with a wide range of triaxiality.
3) The FR class in TNG shows the largest variety in stellar halo properties and the largest variations with radius in both intrinsic shapes and rotational support. Among these galaxies we can find rotationally supported stellar halos with flattened oblate shapes, as well as FRs that have central rotating disk-like structures embedded in more spheroidal components. For a subset of these the stellar halos can reach high triaxiality values. SRs, by comparison, display milder changes in structure with radius.
4) The TNG FRs exhibit shallower V rot /σ(R) profiles than the Atlas3D and ePN.S galaxies (Fig. 13). Both TNG50 and TNG100 FRs tend to reach a peak in rotation at a median radius of ∼ 3R e , compared to ∼ 1 − 1.3R e for the Atlas3D and ePN.S galaxies. This result implies a more extended distribution of the angular momentum with radius in the TNG galaxies than observed. 5) However, even though the V rot /σ(R) profiles do not agree quantitatively between simulated and observed ETGs galaxies, the simulations do reproduce the observed kinematic transitions between central regions and outskirts. The similarity between the scaled shapes of the V rot /σ(R), and the ε(R) profiles between ePN.S and TNG galaxies (Figs. 8b and 14) suggests that also the observed variations in the kinematics between central regions and halos trace changes in the intrinsic structure of the galaxies. 6) We find that most of the triaxial TNG galaxies display modest photometric twists that only weakly depend on triaxiality (Fig. 11). For these galaxies kinematic twists are larger (Fig. 16), but in many cases they are not large enough to be measured by currently available data. 7) By comparing the distributions of the minor-to-major axis ratios q, triaxiality parameters T , and angular momentum parameters λ (Fig. 17), we find that lower rotational support in the stellar halos is related to rounder shapes, with a wide range of triaxiality which depends on stellar mass. In this there is no qualitative difference between the FRs and SRs. Rather, despite the bimodality of the central regions, the two classes show a continuity of halo properties with the FRs dominating the low T -high λ end of the distribution and the SRs dominating in the high T -low λ extreme. The relative weight of the different sides of the distribution gradually changes with stellar mass. This is in agreement with ePN.S observations of ETG halos.
In a companion paper we will investigate the dependence of the stellar halo parameters on the accretion history of galaxies, and explore the relation between stellar and dark matter halo properties.  Fig. A.1. Relation between colors for the observed galaxies, as shown in the legend. The simulated galaxies with M * ≥ 10 10 M are also shown but without any modeling for dust attenuation, which probably explains why they tend to follow a different color relation than the observations. The solid lines show the linear fit to the observations, while the dashed line shows the color transformation equation from Jester et al. (2005).
In Sect. 4.1 we selected a sample of ETGs from TNG50 and TNG100 by using a cut in g − r color and stellar mass, and compare with observations. Most of the observed galaxies have g − r colors in the NSA catalog. For a small fraction of the observations, only B − V or g − i colors are available. Jester et al. (2005) published transformation equations between the SDSS and UBVRcIc magnitudes valid for stars and quasars, which provide reasonable results for galaxies. We note, however, that these transformation equations do not readily apply to the observed galaxies. In Fig. A.1 we show the relation between colors in galaxies that have measured both g−r and B−V, and g−r and g−i. g−r and g−i colors are from the NSA catalog, B−V colors are from Hyperleda. We find that the observed ETGs follow the relations: Hence we use these derived transformation equations to transform the measured B − V or g − i colors in g − r.
Appendix B: Accuracy on the measured intrinsic shapes and angular momentum parameter λ e Physical properties measured in simulated galaxies are affected by resolution effects coming from the discrete particle nature of these systems. In a collisionless dark matter-only (DMO) simulation the resolution of the gravitational potential depends on the softening length and the particle mass resolution (i.e. the number of particles), which particularly affect measurements in the central regions of the simulated systems (e.g. Power et al. 2003). In a full physics simulation instead, like TNG100 and TNG50, the resolution has additional effects on the baryon physics, which require the models to be calibrated on observations (Schaye et al. 2015;Pillepich et al. 2018b). In this section we are not concerned with the latter effects, which also have an impact on the galaxy properties. Here we aim at quantifying the effects from the resolution of the gravitational potential on shape and spin measurements at 1R e and beyond. Chua et al. (2019) analyzed the convergence of intrinsic shape profiles in the Illustris-DMO simulation, with the procedure recommended by Zemp et al. (2011), also used in this paper and outlined in Sect. 3.1. From their Figure 1 we find that the shape profiles are converged within 0.1 in both p and q already at r ∼ 2r soft . Since TNG100 has similar particle resolution as Illustris-DMO (and smaller r soft ), we can apply a similar convergence criterion in TNG100. We verified that the two simulations have comparable particle numbers at r ∼ 9r soft , i.e. where full convergence is achieved according to the prescription of Chua et al. (2019). At smaller radii the full physics TNG100 simulation has more particles than the DMO simulation since the baryons are subjected to dissipation; hence we expect similar, if not better, convergence in TNG100.
Twice the softening radius in both TNG50 and TNG100 corresponds to a radius R e for all the selected galaxies (see Fig. 7, and note that in TNG100 we excluded the galaxies with R e < 2r soft from the sample). Since R e increases with stellar mass, the more massive systems are better resolved. For example, in TNG100 the median effective radius at M * = 10 11 M is 1R e ∼ 7.8r soft , where the resolution effects on both p and q are 0.02. Thus the absolute error on the axes ratios p and q measured at r = 1R e is mass dependent, and it is at most 0.1 for the low mass systems. At r 9r soft , which corresponds to r 3.5R e at the low mass end, and to 1.16R e for M * = 10 11 M , the softening effects are negligible compared to the errors coming from particle statistics. For a minimum required number of 1000 particles per ellipsoidal shell we found that these errors are generally 0.02 on both p and q. Hence at r 9r soft the uncertainty on the axes ratios is of the order of 0.02.
Throughout the paper we quantify the intrinsic shapes of galaxies also with the triaxiality parameter. Because of its definition (Eq. (3)), the error ∆T is shape dependent, as well as mass dependent as discussed above. In this work we consider reliable T measurements performed at r ≥ 9r soft , where the uncertainties on the axes ratios are ∼ 0.02, corresponding to ∆T ∼ 0.2 for typical FRs axis ratios (i.e. p = 0.9 and q = 0.5). At smaller radii, where ∆T grows large for the low mass galaxies, we quantify the intrinsic shapes using only the better determined p and q (e.g. in Fig. 15).
The angular momentum parameter λ e is evaluated by integrating over all the particles within R ≤ 1R e , hence, by definition, it is derived more reliably than the flattening. This is apparent in the convergence study of Lagos et al. (2017) on the EA-GLE simulations, which have particle mass resolution and gravitational softening length very similar to TNG100. The study shows that the angular momentum within 1R e is well converged already at stellar masses above M * = 10 9.5 M , i.e. within galaxies resolved by about 2000 particles. The selected galaxies in the TNG100 sample are resolved with more than 2 × 10 4 particles, hence the measured λ e is independent of the softening and particle mass resolution.
In conclusion we find that both angular momentum and shapes are not (or only marginally in case of the shape) affected by the resolution of the gravitational potential and by the particle number at r 1R e , and they are well-determined at larger radii.

Appendix C: Elongated galaxies in IllustrisTNG
In Sect. 4.2 we further restrict the sample selection by excluding an excess of centrally elongated galaxies not present in the observed samples. The inconsistency is revealed in the λ e − ε(1R e ) diagram of Fig. 5, in which the centrally elongated galaxies distribute in a region at high ellipticity where there are few observed counterparts. In App.B we showed that the resolution effects at 1R e on the intrinsic shapes are at most of the order of 0.1: these are not enough to explain the differences with observations in the Article number, page 21 of 23 A&A proofs: manuscript no. TNGkinematics  Top: Median p(r) profiles for the galaxies with p(1Re) < 0.6 in TNG100 (left) and TNG50 (right). Only bins with more than 10 galaxies are shown. The percentages indicate the fraction of centrally elongated galaxies that populate each median profile. Bottom: Distribution of the axis ratio q(5R e ) in the centrally elongated galaxies of TNG100 (left) and TNG50 (right). Dashed and dotted vertical lines show the mean and the median of the distributions, respectively. The galaxies with p(1Re) < 0.6 have elongated shapes up to ∼ 3 − 4 Re; outside this regions (r ∼ 5R e ) they are mostly near-oblate (p ∼ 0.9) with median flattening q ∼ 0.3 − 0.4, depending on the simulation. The top panels of Fig. C.3 show the median p(r) profiles for the centrally elongated galaxies. The profiles are built by binning together the galaxies according to the values of p(1R e ) and p(6R e ). The percentages next to each profiles gives the fraction of centrally elongated TNG100 or TNG50 galaxies that populate the median profile. The elongated regions can extend up to a few Re, typically ∼ 3R e . Outside these radii, galaxies are nearoblate, with median p(5R e ) ∼ 0.9, and a large scatter on the flattening q(5R e ) (bottom panels in Fig. C.3). The TNG100 galaxies have rather spheroidal shapes with q(5R e ) ∼ 0.45; the TNG50 galaxies tend to be flatter with q(5R e ) ∼ 0.3, and a peak in the distribution at ∼ 0.2.
Finally Fig. C.4 shows the color-stellar mass diagram for the galaxies in TNG100 and TNG50 separately. The axis ratio p(1R e ) is shown by the color of the data points. The centrally elongated systems occur prominently among the redder simulated galaxies, and in specific stellar mass ranges. In TNG100 The fact that these centrally elongated galaxies are preferentially produced in a particular mass range, and that most of them are old, red systems, points towards these objects being a class of galaxies that are produced by the simulation but are not present in nature, probably related to the way the simulated galaxies accrete gas, form stars, and quench during a rapid dynamical evolution in this specific mass range. The difference between TNG100 and TNG50 in the stellar mass range where the centrally elongated systems are produced could be explained by the higher star formation efficiency in the higher resolution simulation (Pillepich et al. 2018b). This hypothesis is strengthened by the absence of a similar concentration of centrally elongated systems among the red galaxies of intermediate masses in the Illustris simulation: hence the presence of this population of galaxies is due to the galaxy formation model. We note that in TNG100 this mass range approximately coincides with the knee in the stellar mass-halo mass relation (e.g. Behroozi et al. 2013, and reference therein), where accretion could be particularly efficient.
As Fig. C.3 showed, these centrally elongated systems are embedded in an near-oblate component with various flattening q, and q tends to be smaller in TNG50. This difference in flattening is likely due to the better resolution of TNG50, which allows to resolve thinner disks (Pillepich et al. 2019, their appendix B and C). From App. B above the error on q due to the spatial resolution is of the order of 0.1 at ∼ 1R e for the lowest mass systems, so that the measured thickness of a thin disk with radius ∼ 5R e and q ∼ 0.2 may be subject to a similar uncertainty. If the centrally elongated galaxy components in TNG100 and TNG50 actually form within disks, then they could be the result of a bar instability, such that at the particular stellar mass range where the bars are produced, too much cold gas forms stars too quickly, while building massive bar-unstable disks. Then the feedback from the intense star formation would sweep away the remaining gas, quench the star formation, and lead to the formation of preferentially red bar-like inner components. For some of these objects the time-scale for dynamical friction against the surrounding stars and dark matter might be short enough to slow down their rotation, generating the wide range of λ e values seen in Fig. C.2.
This suggests that these centrally elongated galaxies may be systems that were in the process of forming a disk, whose evolution was interrupted or derailed by rapid dynamical instability, star formation, and feedback in the simulations, in the particular mass range in which they occur. Fig. 6 in the main text demonstrates that the distributions of angular momentum and triaxiality in the surrounding stellar halos are not affected by the central elongation p(1R e ) of the simulated galaxies. Thus these over-elongated systems can simply be excluded from the sample selected for our main analysis.