© C. Pulsoni et al. 2020

1. 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 are essentially divided 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) (Emsellem et al. 2011; 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, mainly dry minor mergers (Naab et al. 2009; Johansson et al. 2012), 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 and fast rotator (i.e., the core and 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; Buitrago et al. 2017). 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 observed tidal debris from recent accretion events in the halos of many galaxies (e.g. Malin & Carter 1983; Janowiecki et al. 2010; Longobardi et al. 2015a; 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. 3), 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 radii of the observed kinematic transitions to the halo and their dependence on the galaxies’ stellar mass seem to support such an interpretation.

To date, only a few studies of the kinematic properties of stellar halos in simulations are available in the literature. Wu et al. (2014) analysed the kinematics of 42 cosmological zoom simulations of galaxies and found a variety of V/σ profile shapes (rising, flat, or with a maximum), in agreement with observations. However, these early simulations did not reproduce the whole spectrum of properties of observed FRs, especially the fast rotating and extended disks (Emsellem et al. 2011; Pulsoni et al. 2018). Recently, Schulze et al. (2020) using the Magneticum Pathfinder simulations showed that these simulations reproduce the observed kinematic properties of galaxies more closely, and that extended kinematics is a valuable tool for gaining insight into galaxy accretion histories. They also found that the kinematic transition radius is a good estimator of radius of the transition between in-situ and ex-situ dominated regions for a subset of galaxies with decreasing V/σ profiles, especially those that did not undergo major mergers in their evolution.

The goal of this paper is to better understand the structural changes between the centers and stellar halos of ETGs with a large and well-resolved sample of simulated galaxies. We study the stellar halo structure, that is, 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 (Springel et al. 2018; Pillepich et al. 2018a; Naiman et al. 2018; Marinacci et al. 2018; Nelson et al. 2018, 2019a), 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 (Genel et al. 2014; Vogelsberger et al. 2014), using a refined galaxy formation model. For this work we consider two cosmological volumes with side lengths ∼100 Mpc and ∼50 Mpc, which are referred to as TNG100 and TNG50. TNG50 is the highest resolution realization of the IllustrisTNG project (Pillepich et al. 2019; Nelson et al. 2019b) with particle resolution more than 15 times better than TNG100.

The paper is organized as follows. For the comparison of the TNG galaxies properties with observations, we first summarize in Sect. 3 how different ETG surveys select their samples and how physical quantities are measured. Section 4 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 (Sect. 5), we proceed to show the photometric results in Sect. 6 and the kinematic results in Sect. 7. Section 8 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, Sect. 9 summarizes our conclusions.

2. The IllustrisTNG simulations

The IllustrisTNG simulations are a new generation of cosmological magnetohydrodynamical simulations using the moving mesh code AREPO (Springel 2010). Compared to the previous Illustris simulations, they include improvements in the models for chemical enrichment, stellar and black hole feedback, and introduce new physics such as the growth and amplification of seed magnetic fields.

The baryonic physics model contains a new implementation of black hole feedback (Weinberger et al. 2017), as well as updates to the galactic wind feedback, stellar evolution and gas chemical enrichment models (Pillepich et al. 2018b). These modifications, in particular those for the two feedback mechanisms, were required to alleviate some of the tensions between Illustris and observations, such as the large galaxy stellar masses below the knee of the galaxy stellar mass function and the gas fractions within group-mass halos. They in turn also improve on the too large stellar sizes of galaxies and the lack of a strong galaxy color bimodality at intermediate and high galaxy masses in Illustris (Nelson et al. 2015).

The IllustrisTNG fiducial model was chosen by assessing the outcome of many different models against the original Illustris by using additional observables, specifically the halo gas mass fraction and the galaxy half-mass radii, with respect to those used to calibrate the Illustris model against observational findings, such as, the star formation rate density as a function of z, the galaxy stellar mass function at z = 0, the z = 0 black hole mass versus halo mass relation, and the z = 0 stellar-to-halo mass relation.

The new AGN feedback model is responsible for the quenching of galaxies in massive halos and for the production of red and passive galaxies at late times, alleviating the discrepancies with observational data at the massive end of the halo mass function (Weinberger et al. 2017; Nelson et al. 2018; Donnari et al. 2019). The faster and more effective winds in TNG reduce the star formation at all masses and all times, resulting in a suppressed z = 0 galaxy stellar mass function for M_* ≲ 10¹⁰ M_⊙, and smaller galaxy sizes (Pillepich et al. 2018b). Overall the TNG model has been demonstrated to agree satisfactorily with many observational constraints (e.g., Genel et al. 2018; Nelson et al. 2018) and to return a reasonable mix of morphological galaxy types (Rodriguez-Gomez et al. 2019).

In this study we consider two simulation runs, TNG100 and TNG50, which are the two highest resolution realizations of the IllustrisTNG intermediate and small cosmological volumes. TNG100 has a volume and resolution comparable with Illustris, while TNG50 reaches resolutions typical of zoom-in simulations. Table 1 summarizes and compares the characteristic parameters of the two simulations.

The TNG model is calibrated at the resolution of TNG100 and all the TNG runs adopt identical galaxy formation models with parameters that are independent of particle mass and spatial resolution (“strong resolution convergence”, according to Schaye et al. 2015). This imposition results in some of the properties of the simulated galaxies being resolution dependent. As discussed by Pillepich et al. (2018b), this can be primarily explained by the fact that better resolution allows the sampling of higher gas densities, hence more gas mass is eligible for star formation and the star formation rate accelerates. This means that, for example, at progressively better resolution, galaxies tend to have increased stellar masses at fixed halo mass and smaller sizes at fixed stellar mass (see also Pillepich et al. 2019 for a quantification of these effects).

3. 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 (Cappellari et al. 2011), 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. 2015b; 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, and 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 (Cappellari et al. 2011). The Atlas3D 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⁸ − 10¹² M_⊙) and environment (field, group, and clusters). This sample is not morphologically selected, but we use the data from van de Sande et al. (2017) 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¹ (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⁹ − 10¹²; 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 van de Sande et al. (2017) report g − i colors, which we convert to g − r using the transformation equation derived in Appendix 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² catalog (Makarov et al. 2014), and convert to g − r colors using the relations in Appendix A.

Sizes – For the Atlas3D sample we use the effective radii (R_e) values in Table 3 of Cappellari et al. (2011). Those for the MASSIVE galaxies are from Ma et al. (2014, Table 3), where we adopt the NSA measurements, where available, or the 2MASS values corrected using their Eq. (4). The data for MANGA are from Graham et al. (2018). For SAMI we use the data presented in van de Sande et al. (2017), and we circularize the effective semi-major axis by using the reported value for the ellipticity. The half light radii for the ePN.S galaxies are in Table 2 of 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. Section 6.1 discusses the systematic effects in comparing observed effective radii and half-mass radii in simulated galaxies.

Stellar masses – The IllustrisTNG model assumes a Chabrier (2003) initial mass function (IMF). The stellar masses for the SAMI survey in van de Sande et al. (2017) are derived using a color–mass relation, and a Chabrier IMF. For Atlas3D, MASSIVE, 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_Kcorr = 1.07M_K + 1.53, and converted to stellar masses with the formula from van de Sande et al. (2019):

$$\begin{array}{c}\hfill {log}_{10}{M}_{\ast }=10.39-0.46(M_{K_{\mathrm{corr}}}+23),\end{array}$$(1)

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. The missing flux correction takes into account the over-subtraction of the sky background by the 2MASS data reduction pipeline (Schombert & Smith 2012) and the limited 4 R_e aperture of the 2MASS measurement.

For the ePN.S sample we derive stellar masses using integrated luminosities from the most extended photometric profiles available in the literature, 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.

There are several sources of errors in the stellar mass estimates of observed galaxies. The uncertainty in the magnitudes derived from the 2MASS photometry are typically ∼0.25 mag (Scott et al. 2013). The uncertainty in the distances typically translate into an error of 0.1 mag on the absolute magnitudes but can reach up to 0.5 mag (van de Sande et al. 2017). These uncertainties correspond to an error on the stellar mass of typically ∼0.1 dex and up to ∼0.2 dex. In addition the total luminosity, and hence the total stellar mass, can be underestimated if the photometry is not deep enough to measure the faint surface brightness of the stellar halos, especially in massive galaxies with large Sérsic indices or described by multiple Sérsic components. Since the stellar masses of the simulated galaxies are evaluated using the total bound stellar mass (Sect. 5.1), this may cause a systematic difference between observed and simulated stellar masses at the high mass end; see also Sect. 6.1.

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−3 R_e). Ellipticities for the SAMI galaxies are from van de Sande et al. (2017), 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). The measurement errors on the ellipticities are per se very small (O(10⁻³), Kormendy et al. 2009), but the characteristic ellipticities used by different surveys for the same galaxies can differ within a root-mean-square scatter of ∼0.05 (see e.g., Veale et al. 2017 and Fig. 2 from Graham 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, that is, with semi-major axis $R_{\mathrm{e}}\sqrt{1-\epsilon }$, where ε is the ellipticity).

The uncertainties on the measured λ_e for the Atlas3D galaxies are generally small, Δλ_e ≃ 0.01 (Emsellem et al. 2011). Similar errors apply for the MASSIVE sample, Δλ_e ≲ 0.01 (Veale et al. 2017). SAMI and MANGA instead target objects at larger distances with lower apparent sizes and spatial resolution. For these galaxies the measurement uncertainties are combined with seeing effects, which generally tend to systematically decrease λ_e. In the SAMI galaxies, for a typical seeing of 2 arcsec, van de Sande et al. (2017) find that measurement errors (Δλ_e ∼ 0.01) and seeing effects cancel out for galaxies with λ_e <  0.2, while for λ_e >  0.2 seeing is the dominant effect and causes a median decrease in λ_e of 0.05. For the MANGA regular rotators in the cleaned sample, Graham et al. (2018) estimate mean Δλ_e = [0.005, −0.041] and median errors Δλ_e = [0.004, −0.027].

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 the procedure described in Sect. 4.4 directly to the velocity fields from Emsellem et al. (2004) and Cappellari et al. (2011), giving a median error on V/σ of the order of 0.03.

For the ePN.S galaxies the procedure is applied to the PN velocity fields, whereas for the central regions we use the V_rot and σ from kinemetry analysis on IFS data from Krajnović et al. (2008, 2011); Foster et al. (2016), when available. In the other cases we use V_rot and σ from major axis slits (see references in the ePN.S paper). For the ePN.S galaxies the measurement uncertainties on the V/σ profiles are dominated by the statistical error on the PN velocity fields. The median Δ(V/σ) = 0.08.

4. 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_ij:

$$\begin{array}{c}\hfill I_{\mathit{ij}}=\frac{{\sum }_n{M}_n{x}_{n,i}{x}_{n,j}}{{\sum }_n{M}_n},\end{array}$$(2)

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, for example, 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. 6.1.

4.1. Intrinsic shapes

The three-dimensional intrinsic shapes of the galaxies are evaluated by diagonalizing the inertia tensor I_ij 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 8 R_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 ${\widehat{e}}_j$ (with j = a, b, c) of the inertia tensor.

We also use the triaxiality parameter,

$$\begin{array}{c}\hfill T(r)=\frac{1-p{(r)}^2}{1-q{(r)}^2},\end{array}$$(3)

to quantify the intrinsic shape.

In Appendix B we find that shape measurements at 1 R_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 the resolution on TNG100. At r ∼ 9 r_soft, which is r ∼ 3.5 R_e for the lowest mass galaxies, and r ∼ 1.1 R_e for M_* = 10¹¹ M_⊙, these resolution effects are negligible, and the error on the shape measurements is then due to particle noise and is ∼0.02 in TNG100. 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 Appendix B, we consider the triaxiality profiles reliable starting from r = 9 r_soft; at smaller radii, where ΔT is larger, we quantify shapes using p and q which are better defined. These results for TNG100 are summarized in Table 2. For the TNG50 galaxies, we expect similar or lower uncertainties.

In the paper, we 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(1 R_e) = 0.43 and p(1 R_e) = 0.95 (T <  0.3). At large radii the galaxy becomes close to prolate with q(10 R_e) = 0.76, p(10 R_e) = 0.83, and triaxiality parameter T >  0.7. For the galaxy shown, 9 r_soft = 0.48 R_e.

Fig. 1. Photometric measurements. Top: intrinsic shape of an example galaxy from TNG100 as a function of the intrinsic major axis distance r/Re; 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 3 Re, and becomes near-prolate (T >  0.7) in the outer halo. In this galaxy 9 rsoft correspond to 0.48 Re. Bottom: ellipticity and photometric position angle profiles for two example galaxies from TNG100 as a function of the projected major axis distance R/Re. The quantities derived from the inertia tensor are shown with solid symbols, those from the mock images with open symbols.

4.2. 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. 4.1 for the 3D intrinsic shape. The square root of the ratio of the two eigenvalues of I_ij 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 ≥ 3 R_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.

4.3. 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 arcsec 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 8 R_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 ith bin we calculate the projected mean velocity and the mean velocity dispersion as the weighted averages:

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

where the index n runs over the particles in the bin, and N_i is the number of particles in the ith 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.

Fig. 2. Stellar kinematic measurements. Top: Voronoi binned mean velocity fields for an example TNG100 galaxy, with a central bulge and a relatively massive disk. For this galaxy the random LOS projection almost coincides with the edge-on projection. The projected major axis is aligned with the X axis. The data points show the projected (X, Y) positions of the stellar particles, and are colored according to the mean velocity and velocity dispersion of the corresponding Voronoi bin as shown by the colorbars. Center: smoothed velocity and velocity dispersion fields for the same example galaxy above. The central disk structure is embedded in a spheroidal stellar halo. Bottom: kinematic parameters of the galaxy shown above, derived from the Voronoi binned velocity fields (orange open circles) and from the smoothed velocity fields (blue points).

The systemic velocity of the galaxy is derived by fitting a harmonic expansion as in Pulsoni et al. (2018, their Sect. 4) to the central regions (i.e., at R ≤ 2 R_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)

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{e}}=\frac{{\sum }_i{M}_i{R}_{\mathrm{circ},i}|V_i|}{{\sum }_i{M}_i{R}_{\mathrm{circ},i}\sqrt{V_i^2+{\sigma }_i^2}},\end{array}$$(5)

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 = ∑_nM_n, i, and R_circ, i is the circular radius of the ith bin. The cumulative λ_endixe 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 ε(1 R_e). By comparison, the differential λ(R) is summed in elliptical shells. As discussed in Appendix B the angular momentum parameter is not affected by resolution at R ≳ 1 R_e for the selected sample of galaxies.

4.4. 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 >  2 R_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. We estimated uncertainties on the differential λ(R) and on V/σ(R) in TNG100 by considering a few kinematically representative galaxies in three stellar mass bins and studied the kinematic parameter distributions derived from 1000 simulations respectively, with particle numbers decreased to the typical numbers at different multiples of R_e. Table 3 lists the standard deviation of the distributions for typical numbers of particles at 2 R_e and 8 R_e.

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 near-oblate 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.

5. Selection of the sample of ETGs in the IllustrisTNG simulations

5.1. 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¹⁰ <  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 Fig. 13).

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. 3), we chose

$$\begin{array}{c}\hfill (g-r)\ge 0.05{log}_{10}(M_{\ast }/M_{\odot })+0.1\mathrm{mag}.\end{array}$$(6)

For M_* we use the total bound stellar mass of the galaxies. 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. 3).

We limited the sample stellar mass range to 10^10.3 ≤ M_* ≤ 10¹² M_⊙. This choice assures that the TNG100 galaxies are resolved by at least 2 × 10⁴ stellar particles. By comparison, the minimum number of stellar particles in the selected TNG50 galaxies is 36 × 10⁴.

In addition, we impose that the galaxies’ effective radius (see Sect. 6.1) R_e ≥ 2 r_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 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.

Fig. 3. Selection of galaxies in color and stellar mass. Top: 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 1010.3 <  M*/M⊙ <  1012. Bottom: ETGs from recent IFS surveys as indicated. Most of the observed ETGs in this mass range fall in the same red sequence region.

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. 5.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 that a more generous color selection, including bluer galaxies, would produce a too large number of high ellipticity galaxies, especially in TNG50.

Fig. 4. 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. 5.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 Sects. 5.1 and 5.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.

In the following, whenever we compare simulated and observed galaxy samples, we apply to the observed galaxies the same color and stellar mass selection criteria that we used for the TNG sample.

5.2. Selection of ETGs in the λ-ellipticity diagram: fast and slow rotators

Figure 5 shows the λ_e–ε(1 R_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. 5.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. 3), in the same color and stellar mass region as defined in Sect. 5.1. Here we also include for comparison the spiral and irregular galaxies from the MANGA sample (marked as LTGs).

Fig. 5. λe–ε(1 Re) diagrams for TNG and observed galaxy samples, in stellar mass bins. Top row: λ-ellipticity diagram for 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 1 Re. Middle row: λ-ellipticity diagram for the TNG galaxies after the additional selection p(1 Re)≥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 ε(1 Re) are as described in Sect. 3. The solid black line is the threshold separating fast and slow rotators as defined in Emsellem et al. (2011): ${\lambda }_{\mathrm{e}}=0.31\sqrt{\epsilon }$. 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–ε(1 Re)-plane, with the exception of the high-λe (> 0.7) S0 disks specific to the MANGA survey.

We observe that a significant fraction of the TNG galaxies shown in the top row populate a region to the right of the λ_e–ε(1 R_e) diagram where there are no observed counterparts, i.e. below the magenta line and with ε(1 R_e) > 0.5. By color coding the galaxies according to their intrinsic axis ratios at r ∼ 1 R_e, we find that these galaxies have elongated, triaxial shapes. These systems occur at all values of λ_e, that is, some rotate as rapidly as the MANGA disk galaxies, but others do not show any rotation (Fig. C.1).

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 ε(1 R_e) values were used and placed these objects in the region of the λ_e–ε(1 R_e) populated by the centrally elongated (at r ∼ 1 R_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 1 R_e are at most of the order of 0.1, for the low mass galaxies, see Appendix B).

In Appendix 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, and the λ − ε diagrams for simulated galaxies in Magneticum (Schulze et al. 2018) and EAGLE (Walo-Martín et al. 2020) do not contain many objects with large ellipticities and intermediate to low λ_e. This indicates that the new galaxy formation model in TNG is involved in the occurrence of these centrally elongated galaxies. 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 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 ∼ 1 R_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 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 location 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, that is, λ and triaxiality parameter, are not affected by the sample selection based on p(1 R_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.

Fig. 6. Distribution of λ parameter and triaxiality in the stellar halos (i.e. at R ≥ 4 Re) for different intrinsic shapes in the central regions, parametrized by the intermediate to major axis ratio p(1 Re). TNG100 and TNG50 galaxies are shown separately, in the top and bottom panels respectively. These distributions are not systematically dependent on p(1 Re), and thus the stellar halo spin and triaxiality remain unbiased after implementing a sample selection based on p(1 Re).

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¹¹ 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),

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{e}}=0.31\sqrt{{\epsilon }_{\mathrm{e}}},\end{array}$$(7)

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).

5.3. 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¹² M_⊙ and with red (g − r) color as in Eq. (6) (see Fig. 3);

• excluding a small number of objects with R_e <  2 r_soft, to assure sufficient resolution at 1 R_e;

• finally, removing a class of centrally elongated, triaxial galaxies with p(1 R_e) < 0.6, which are systems not present in the observed ETG samples that probably became bar-unstable and quenched during the process of (central) disk formation.

The selected sample has a distribution of λ_e − ε(1 R_e) that is 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.

Fig. 7. Circularized effective radii, i.e. $R_{\mathrm{e}}\sqrt{1-\epsilon }$, of the selected ETG samples in the TNG100 and TNG50 simulations as a function of stellar mass, and comparison with observations (bottom panel). The black curves show the median profile of the distribution of galaxy Re from all surveys and the 10th and 90th percentiles, in both panels. The dashed lines in the top panel indicate the resolution of the two simulations.

6. 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 6.1 discusses the measured galaxy sizes and how our definition of effective radii compares with effective radii inferred from ETG photometry. Section 6.2 compares the distribution of projected ellipticities at 1 R_e with that from ETG surveys and validates our sample selection. Section 6.3 studies the TNG ellipticity profiles out to the stellar halo, Sect. 6.4 explores the intrinsic shape distribution of stellar halos and its dependence on stellar mass, and Sect. 6.5 investigates the dependence of galaxy triaxiality on radius and stellar mass in the simulated samples. Finally, Sect. 6.6 tests the ability of photometric twist measurements to establish the underlying triaxiality in the TNG galaxies.

6.1. 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 (edge-on 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. TNG50 produces smaller galaxies compared to TNG100 and to observations at intermediate stellar masses (Pillepich et al. 2019). This is purely a resolution effect, as discussed in Sect. 2. 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¹¹ 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. (2015b) 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 Fig. 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¹¹ M_⊙ objects the differences in the derived R_e and stellar masses are within 10% and 5% respectively, while between 10¹¹ 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 the 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 mass-to-light ratio to obtain total stellar masses), here we 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.

6.2. Ellipticity distribution in the central regions

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.

Fig. 8. Ellipticity distribution and profiles of the selected TNG ETGs. a: distribution of the random LOS ε(1 Re) 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.

The TNG50 and TNG100 simulations predict a significant fraction of SR galaxies with ε(1 R_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 yet 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 have 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.

6.3. 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 4 R_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.

6.4. 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 constant values beyond ∼5 R_e. Hence we measure the stellar halo intrinsic shape distributions by deriving the intrinsic axes ratios in a shell around 8 R_e, 1.5 R_e thick, which is the maximal radius at which also the lowest mass systems contain enough particles to reliably measure intrinsic shapes, see Sect. 4.1. Different choices of the shell thickness, or slightly different choices of the radius (for example 7 instead of 8 R_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 (low-mass) 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.

Fig. 9. Intrinsic shape distribution of ETG stellar halos in TNG. Top: minor to major axis ratio q versus intermediate to major axis ratio p coloured by triaxiality as measured at 8 Re. TNG50 and TNG100 galaxies are shown with different symbols as in the legend. Bottom: intrinsic shape distribution for the halos (r ∼ 8 Re) of the FRs (top panels) and SRs (bottom panels) in mass intervals, shown with different colors as indicated in the figure. The fitted functions are shown with solid lines, and the fitted parameters are reported in the legend. The vertical dashed lines show the comparison with the photometric model used in Pulsoni et al. (2018) to reproduce the observed photometric twists and average ellipticities for the ePN.S survey. Most low-mass TNG galaxies have near-oblate stellar halos (top), changing towards more triaxial shapes with increasing stellar mass (bottom panels).

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.

6.5. 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. 4.1 that because of the error due to resolution effects in the central regions, T profiles are considered well-defined only beyond r = 9 r_soft, where their error ΔT ∼ 0.2 for typical FR axis ratios (Appendix B). Thus we show profiles only for r >  3.5 R_e for the lowest mass objects, and for r >  1.16 R_e for galaxies with M * ∼10¹¹ 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 8 R_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 8 R_e as a function of the stellar mass.

Fig. 10. Left panels: median triaxiality profiles for the fast rotators (top) and the slow rotators (bottom). The numbers on the left indicate the percentage of FRs or SRs populating each profile. The vertical dotted lines show r = 9 rsoft for TNG100 galaxies with M* ∼ 1010.3 M⊙, the solid lines show r = 9 rsoft for TNG100 galaxies with M* ∼ 1011 M⊙. At radii larger than these the profiles are not affected by resolution issues. Right panels: median of the stellar halo triaxiality distribution measured at 8 Re 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.

FRs are characterized by increasing T profiles, which tend to plateau at r >  5 R_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¹¹ 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.

6.6. 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, that is, 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, 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.

Fig. 11. Maximum measured photometric twist between 1 and 8 Re as a function of the halo triaxiality T(8 Re) 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 PAphot are large. The color-bar indicates the median ellipticity measured between 1 and 8 Re. The gray solid line shows the median of the photometric twist distribution as a function of T(8 Re); 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.

We observe that the amplitude of the photometric twists is only weakly dependent on the triaxiality. Near-oblate and near-prolate 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(8 R_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.

7. The kinematics properties

In this section, we study how the kinematic properties of the TNG galaxies vary with radius. In Sect. 7.1, we derive median differential λ(R) profiles to quantify the variety of kinematic behaviors in the outskirts of FRs and SRs. Sections 7.2 and 7.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. 7.4 uses the simulated galaxies to assess kinematic misalignments and twists as signatures of triaxial shapes in galaxies.

7.1. 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 = 4 R_e and at R = 8 R_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 λ(7 R_e).

Fig. 12. Left: median λ profiles for fast rotators (upper left) and slow rotators (lower left panel) in their edge on projection. The median profiles are built by binning the galaxies according to the shape of their λ profiles; see text. The shaded areas are bounded by the 25th and 75th percentiles of the distributions. The numbers on the left indicate the percentage of FRs or SRs populating each profile. Right: distribution of λ(8 Re) for the full sample, as a function of stellar mass. FRs and SRs are shown separately, as well as TNG50 and TNG100. FR galaxies show a large variety of λ profiles, whereas SRs have either constant or increased λ in the stellar halo. Overall the halo rotational support decreases at high stellar masses.

Most of the FRs reach their maximum λ(R) around 3 R_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 λ(8 R_e) > 0.6, a third (40%) with gently decreasing profiles and 0.3 >  λ(10 R_e)≥0.6, and another third (26%) with very low rotation in the halo (λ(10 R_e)≤0.3).

The SRs essentially divide between a half (53%) with non-rotating halos (λ(7 R_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 (λ(7 R_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 non-rotating 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. 6.5): on average TNG50 SR halos rotate faster and have more oblate shapes. These differences in halo properties might be due to the dependence 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.

7.2. Simulated versus observed rotation profiles – central regions

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 ≤ 4 R_e. Here we normalize the radii by the circularized R_e, i.e. $R_{\mathrm{e}}\sqrt{1-ϵ(1R_{\mathrm{e}})}$, for an appropriate comparison with the Atlas3D R_e.

Fig. 13. Median Vrot/σ profiles for slow rotators (left) and fast rotators (right panels). Top: observations from Atlas3D and ePN.S. For the former we show median profiles and shaded areas bounded by the 20th and the 80th percentiles of the distribution. For the ePN.S galaxies we show individual profiles. Bottom: median profiles and 20%−80% distributions for the TNG50 and TNG100 ETGs, projected along a random LOS. Vertical and horizontal dashed lines are for guiding the eyes in comparing top with bottom panels. This shows that the average Vrot/σ profiles of the simulated ETGs are shallower than the observed profiles.

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 1 R_e compared to the Atlas3D galaxies, and almost none exceeds V_rotσ(1 R_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 ∼ 3 R_e. By comparison, the ePN.S FRs tend to peak at smaller fractions of R_e, at a median 1.3 R_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 λ(1 R_e) (Fig. 5). The Atlas3D sample, selected as described in Sect. 3, has some of the disk galaxies removed which as in MANGA (Fig. 5) will be mostly located at high λ(1 R_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 1 R_e.

In Sect. 6.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.

7.3. Simulated versus observed profiles – outskirts

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 Fig. 9) and are reported in the top left panel. Their measurements use absorption line data at 1 R_e and PN data for the halos, which on average cover 6 R_e with a large scatter (minimum 3 R_e, maximum 13 R_e). The central and right top panels show V_rot/σ(6 R_e) versus V_rot/σ(1 R_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.

Fig. 14. V/σ ratio in the central regions compared with the V/σ at large radii. Top: V/σ(1 Re) versus V/σ(6 Re) for the ePN.S galaxies (left), TNG100 (middle), and TNG50 (right). Bottom: V/σ(3 Re) versus V/σ(8 Re) 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 PAkin aligned with PAphot (in blue). The gray points stand for ePN.S FRs for which the measured Vrot/σ(⟨6 Re⟩) is a lower limit. For the simulated galaxies the different colors mark SRs, and FRs with different halo T(8 Re): 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 8 Re 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.

The position of the observed SRs below the one-to-one line is reproduced by the simulations. As already discussed in Sect. 7.1, there are a few outliers among the simulated SRs with V_rot/σ(6 R_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/σ(1 R_e) > 1.5, and a large fraction having significant V_rot/σ at 6 R_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. 7.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 1 R_e.

In the bottom panels of Fig. 14 we show the comparison at adjusted radii: V_rot/σ(3 R_e) (i.e. at the median radius where the TNG V_rot/σ profiles tend to reach the maximum) versus V_rot/σ(8 R_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σ(8 R_e) > V/σ(3 R_e). We note that FRs with near-prolate stellar halos occupy mostly the lower left corner with low V_rot/σ(3 R_e) and low V_rot/σ(8 R_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/σ(8 R_e). It is interesting that, aside for the differences in the central regions at R ≲ 3 R_e, the observed galaxies with and without signatures for triaxial stellar halos distribute similarly as the simulated triaxial and near-oblate stellar halos, respectively, with the triaxial halos having on average lower V_rot/σ(8 R_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.

7.4. 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(1 R_e)−PA_phot(1 R_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 (Krajnović et al. 2011; Graham et al. 2018). Simulated and observed galaxies show very similar trends with ellipticity, with very flat galaxies 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.

Fig. 15. Distribution of the misalignments Ψ = PAkin − PAphot at 1 Re as a function of the ellipticity for both TNG100 and TNG50 galaxies. Data points are color coded according to their intermediate to major axis ratio p(1 Re); open symbols are SRs, filled symbols are FRs. The solid lines show the scatter of the misalignments as a function of the ellipticity for the TNG galaxies, calculated on the mirrored (Ψ, −Ψ) data. The dashed lines show the scatter for the MANGA and Atlas3D galaxies. Red lines are for the SRs, blue lines for the FRs. The triaxiality of the IllustrisTNG FRs is consistent with the near-alignment of their kinematic and photometric position angles.

Data points in Fig. 15 are color-coded according to the intermediate to major axis ratio p at 1 R_e. At these radii the errors on the axis ratios are at most 0.1 for the low mass systems (see Appendix B), so p(1 R_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 1 R_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 ∼ 5 R_e.

The top panel of Fig. 16 shows the distribution of the maximum kinematic twists, that is, 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 8 R_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. 6.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.

Fig. 16. Top: maximum kinematic twist measured between 1 and 8 Re as a function of the stellar halo triaxiality T(8 Re). 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(8 Re) 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.

IFS kinematics show that the central regions of FR have PA_kin and PA_phot aligned within ∼10°, while SRs are generally misaligned (Krajnović et al. 2011; 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 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°, which is 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.

8. Stellar halo angular momentum and shape

In Sect. 7.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 ∼ 2 R_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.

Fig. 17. Top: λ(R), 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 and 8 Re. 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.

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 and 8 R_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 (8 R_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 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. 6.4, 6.5, and 7.1. Lower mass galaxies host more rotationally 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 Fig. 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. 6.3 and 7.3 suggests that intrinsic shape variations similar to those found in the TNG galaxies also exist in real galaxies.

9. Summary and conclusions

In this paper, we have analysed the kinematic and photometric properties of ∼1200 early type galaxies (ETGs) from the IllustrisTNG 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 15 R_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 ∼3 R_e, compared to ∼1 − 1.3 R_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.

Acknowledgments

We thank the anonymous referee for helpful comments that improved the clarity of the paper. C.P. is extremely grateful to F. Hofmann for his support. This research has made use of the NASA/IPAC Extragalactic Database (NED).

References

Appendix A: Color transformation equations

In Sect. 5.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:

$$\begin{array}{c}\hfill \begin{array}{cc}& g-r=0.79(g-i)-0.04\hfill \\ \hfill & g-r=0.62(B-V)-0.06.\hfill \end{array}\end{array}$$(A.1)

Fig. A.1. Relation between colors for the observed galaxies, as shown in the legend. The simulated galaxies with M* ≥ 1010 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).

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 1 R_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. 4.1. From their Fig. 1, we find that the shape profiles are converged within 0.1 in both p and q already at r ∼ 2 r_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 ∼ 9 r_soft, that is, 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 <  2 r_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¹¹ M_⊙ is 1 R_e ∼ 7.8 r_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 = 1 R_e is mass dependent, and it is at most 0.1 for the low mass systems. At r ≳ 9 r_soft, which corresponds to r ≳ 3.5 R_e at the low mass end, and to 1.16 R_e for M_* = 10¹¹ 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 ≳ 9 r_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 ≥ 9 r_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 ≤ 1 R_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 EAGLE simulations, which have particle mass resolution and gravitational softening length very similar to TNG100. The study shows that the angular momentum within 1 R_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⁴ 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 ≳ 1 R_e, and they are well-determined at larger radii.

Appendix C: Elongated galaxies in IllustrisTNG

In Sect. 5.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 − ε(1 R_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 Appendix B we showed that the resolution effects at 1 R_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 λ_e-ellipticity diagram, nor the extreme values measured for the flattening p(1 R_e).

Figure C.1 shows example velocity fields for two of these centrally elongated objects (i.e. with p(1 R_e) < 0.6) projected edge-on. The colors show the line-of-sight mean velocity divided by the velocity dispersion at the position of each particle. The two galaxies shown have similar flattening in the central regions p(1 R_e)∼0.4 and q(1 R_e ∼ 0.3), but in one case (top panel) the galaxy contains an extended disk, in the other the galaxy is spheroidal in the outskirts (bottom panel). Most of the selected systems with p(1 R_e) < 0.6 have velocity fields with rotation around the intrinsic minor axis (as the example in Fig. C.1). However, Fig. C.2 shows that these galaxies (in blue to yellow colors) can exhibit different degrees of edge-on rotation at 1 R_e and some of them do not rotate at all (see also bottom panel of Fig. C.1), which is at odds with what is typically seen in barred galaxies (e.g. Falcón-Barroso et al. 2019).

Fig. C.1. Stellar mass distribution (left) and V/σ fields (right) of the central 3 Re for two examples of the centrally elongated galaxies projected edge-on. In the top panel the elongated central component (p(1 Re) = 0.42, q(1 Re) = 0.33, and λe = 0.54) is embedded in a disk structure (q(5 Re) = 0.33). In the bottom panel the bar-like structure (p(1 Re) = 0.39, q(1 Re) = 0.30, and λe = 0.10) is embedded in a spheroid (q(5 Re) = 0.61).

Fig. C.2. λe-ellipticity (1 Re) diagram for the TNG galaxies selected in color and stellar mass as in Sect. 5.1, and – differently from Fig. 5 – projected edge-on. The symbols are colored according to the intermediate to major axis ratio measured at 1 Re. The black and magenta lines are as in Fig. 5.

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(1 R_e) and p(6 R_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 ∼3 R_e. Outside these radii, galaxies are near-oblate, with median p(5 R_e)∼0.9, and a large scatter on the flattening q(5 R_e) (bottom panels in Fig. C.3). The TNG100 galaxies have rather spheroidal shapes with q(5 R_e)∼0.45; the TNG50 galaxies tend to be flatter with q(5 R_e)∼0.3, and a peak in the distribution at ∼0.2.

Fig. C.3. Top: median p(r) profiles for the galaxies with p(1 Re) < 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(5 Re) 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(1 Re) < 0.6 have elongated shapes up to ∼3 − 4 Re; outside this regions (r ∼ 5 Re) they are mostly near-oblate (p ∼ 0.9) with median flattening q ∼ 0.3 − 0.4, depending on the simulation.

Finally Fig. C.4 shows the color-stellar mass diagram for the galaxies in TNG100 and TNG50 separately. The axis ratio p(1 R_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 they are numerous in the interval 10^10.4 M_⊙ ≲ M_* ≲ 10^10.9 M_⊙, while in TNG50 they have 10^10.0 M_⊙ ≲ M_* ≲ 10^10.5 M_⊙.

Fig. C.4. Color (g − r) as a function of the stellar mass. Top panel: TNG100, bottom panel: TNG50. The dashed lines illustrate the selection of the sample of ETGs with masses 1010.3 ≤ M*/M⊙ ≤ 1012, and red color (Sect. 5.1). The colors indicate the axis ratios p(1 Re).

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 shows, 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 Appendices B and C). From Appendix B above the error on q due to the spatial resolution is of the order of 0.1 at ∼1 R_e for the lowest mass systems, so that the measured thickness of a thin disk with radius ∼5 R_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. Figure 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(1 R_e) of the simulated galaxies. Thus these over-elongated systems can simply be excluded from the sample selected for our main analysis.

All Tables

All Figures

Fig. 1. Photometric measurements. Top: intrinsic shape of an example galaxy from TNG100 as a function of the intrinsic major axis distance r/Re; 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 3 Re, and becomes near-prolate (T >  0.7) in the outer halo. In this galaxy 9 rsoft correspond to 0.48 Re. Bottom: ellipticity and photometric position angle profiles for two example galaxies from TNG100 as a function of the projected major axis distance R/Re. The quantities derived from the inertia tensor are shown with solid symbols, those from the mock images with open symbols.

In the text

Fig. 2. Stellar kinematic measurements. Top: Voronoi binned mean velocity fields for an example TNG100 galaxy, with a central bulge and a relatively massive disk. For this galaxy the random LOS projection almost coincides with the edge-on projection. The projected major axis is aligned with the X axis. The data points show the projected (X, Y) positions of the stellar particles, and are colored according to the mean velocity and velocity dispersion of the corresponding Voronoi bin as shown by the colorbars. Center: smoothed velocity and velocity dispersion fields for the same example galaxy above. The central disk structure is embedded in a spheroidal stellar halo. Bottom: kinematic parameters of the galaxy shown above, derived from the Voronoi binned velocity fields (orange open circles) and from the smoothed velocity fields (blue points).

In the text

	Fig. 3. Selection of galaxies in color and stellar mass. Top: 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 1010.3 <  M*/M⊙ <  1012. Bottom: ETGs from recent IFS surveys as indicated. Most of the observed ETGs in this mass range fall in the same red sequence region.
In the text

Fig. 4. 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. 5.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 Sects. 5.1 and 5.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.

In the text

Fig. 5. λe–ε(1 Re) diagrams for TNG and observed galaxy samples, in stellar mass bins. Top row: λ-ellipticity diagram for 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 1 Re. Middle row: λ-ellipticity diagram for the TNG galaxies after the additional selection p(1 Re)≥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 ε(1 Re) are as described in Sect. 3. The solid black line is the threshold separating fast and slow rotators as defined in Emsellem et al. (2011): ${\lambda }_{\mathrm{e}}=0.31\sqrt{\epsilon }$. 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–ε(1 Re)-plane, with the exception of the high-λe (> 0.7) S0 disks specific to the MANGA survey.

In the text

Fig. 6. Distribution of λ parameter and triaxiality in the stellar halos (i.e. at R ≥ 4 Re) for different intrinsic shapes in the central regions, parametrized by the intermediate to major axis ratio p(1 Re). TNG100 and TNG50 galaxies are shown separately, in the top and bottom panels respectively. These distributions are not systematically dependent on p(1 Re), and thus the stellar halo spin and triaxiality remain unbiased after implementing a sample selection based on p(1 Re).

In the text

Fig. 7. Circularized effective radii, i.e. $R_{\mathrm{e}}\sqrt{1-\epsilon }$, of the selected ETG samples in the TNG100 and TNG50 simulations as a function of stellar mass, and comparison with observations (bottom panel). The black curves show the median profile of the distribution of galaxy Re from all surveys and the 10th and 90th percentiles, in both panels. The dashed lines in the top panel indicate the resolution of the two simulations.

In the text

Fig. 8. Ellipticity distribution and profiles of the selected TNG ETGs. a: distribution of the random LOS ε(1 Re) 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.

In the text

Fig. 9. Intrinsic shape distribution of ETG stellar halos in TNG. Top: minor to major axis ratio q versus intermediate to major axis ratio p coloured by triaxiality as measured at 8 Re. TNG50 and TNG100 galaxies are shown with different symbols as in the legend. Bottom: intrinsic shape distribution for the halos (r ∼ 8 Re) of the FRs (top panels) and SRs (bottom panels) in mass intervals, shown with different colors as indicated in the figure. The fitted functions are shown with solid lines, and the fitted parameters are reported in the legend. The vertical dashed lines show the comparison with the photometric model used in Pulsoni et al. (2018) to reproduce the observed photometric twists and average ellipticities for the ePN.S survey. Most low-mass TNG galaxies have near-oblate stellar halos (top), changing towards more triaxial shapes with increasing stellar mass (bottom panels).

In the text

Fig. 10. Left panels: median triaxiality profiles for the fast rotators (top) and the slow rotators (bottom). The numbers on the left indicate the percentage of FRs or SRs populating each profile. The vertical dotted lines show r = 9 rsoft for TNG100 galaxies with M* ∼ 1010.3 M⊙, the solid lines show r = 9 rsoft for TNG100 galaxies with M* ∼ 1011 M⊙. At radii larger than these the profiles are not affected by resolution issues. Right panels: median of the stellar halo triaxiality distribution measured at 8 Re 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.

In the text

Fig. 11. Maximum measured photometric twist between 1 and 8 Re as a function of the halo triaxiality T(8 Re) 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 PAphot are large. The color-bar indicates the median ellipticity measured between 1 and 8 Re. The gray solid line shows the median of the photometric twist distribution as a function of T(8 Re); 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.

In the text

Fig. 12. Left: median λ profiles for fast rotators (upper left) and slow rotators (lower left panel) in their edge on projection. The median profiles are built by binning the galaxies according to the shape of their λ profiles; see text. The shaded areas are bounded by the 25th and 75th percentiles of the distributions. The numbers on the left indicate the percentage of FRs or SRs populating each profile. Right: distribution of λ(8 Re) for the full sample, as a function of stellar mass. FRs and SRs are shown separately, as well as TNG50 and TNG100. FR galaxies show a large variety of λ profiles, whereas SRs have either constant or increased λ in the stellar halo. Overall the halo rotational support decreases at high stellar masses.

In the text

Fig. 13. Median Vrot/σ profiles for slow rotators (left) and fast rotators (right panels). Top: observations from Atlas3D and ePN.S. For the former we show median profiles and shaded areas bounded by the 20th and the 80th percentiles of the distribution. For the ePN.S galaxies we show individual profiles. Bottom: median profiles and 20%−80% distributions for the TNG50 and TNG100 ETGs, projected along a random LOS. Vertical and horizontal dashed lines are for guiding the eyes in comparing top with bottom panels. This shows that the average Vrot/σ profiles of the simulated ETGs are shallower than the observed profiles.

In the text

Fig. 14. V/σ ratio in the central regions compared with the V/σ at large radii. Top: V/σ(1 Re) versus V/σ(6 Re) for the ePN.S galaxies (left), TNG100 (middle), and TNG50 (right). Bottom: V/σ(3 Re) versus V/σ(8 Re) 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 PAkin aligned with PAphot (in blue). The gray points stand for ePN.S FRs for which the measured Vrot/σ(⟨6 Re⟩) is a lower limit. For the simulated galaxies the different colors mark SRs, and FRs with different halo T(8 Re): 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 8 Re 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.

In the text

Fig. 15. Distribution of the misalignments Ψ = PAkin − PAphot at 1 Re as a function of the ellipticity for both TNG100 and TNG50 galaxies. Data points are color coded according to their intermediate to major axis ratio p(1 Re); open symbols are SRs, filled symbols are FRs. The solid lines show the scatter of the misalignments as a function of the ellipticity for the TNG galaxies, calculated on the mirrored (Ψ, −Ψ) data. The dashed lines show the scatter for the MANGA and Atlas3D galaxies. Red lines are for the SRs, blue lines for the FRs. The triaxiality of the IllustrisTNG FRs is consistent with the near-alignment of their kinematic and photometric position angles.

In the text

Fig. 16. Top: maximum kinematic twist measured between 1 and 8 Re as a function of the stellar halo triaxiality T(8 Re). 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(8 Re) 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.

In the text

Fig. 17. Top: λ(R), 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 and 8 Re. 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.

In the text

Fig. A.1. Relation between colors for the observed galaxies, as shown in the legend. The simulated galaxies with M* ≥ 1010 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 the text

Fig. C.1. Stellar mass distribution (left) and V/σ fields (right) of the central 3 Re for two examples of the centrally elongated galaxies projected edge-on. In the top panel the elongated central component (p(1 Re) = 0.42, q(1 Re) = 0.33, and λe = 0.54) is embedded in a disk structure (q(5 Re) = 0.33). In the bottom panel the bar-like structure (p(1 Re) = 0.39, q(1 Re) = 0.30, and λe = 0.10) is embedded in a spheroid (q(5 Re) = 0.61).

In the text

	Fig. C.2. λe-ellipticity (1 Re) diagram for the TNG galaxies selected in color and stellar mass as in Sect. 5.1, and – differently from Fig. 5 – projected edge-on. The symbols are colored according to the intermediate to major axis ratio measured at 1 Re. The black and magenta lines are as in Fig. 5.
In the text

Fig. C.3. Top: median p(r) profiles for the galaxies with p(1 Re) < 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(5 Re) 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(1 Re) < 0.6 have elongated shapes up to ∼3 − 4 Re; outside this regions (r ∼ 5 Re) they are mostly near-oblate (p ∼ 0.9) with median flattening q ∼ 0.3 − 0.4, depending on the simulation.

In the text

	Fig. C.4. Color (g − r) as a function of the stellar mass. Top panel: TNG100, bottom panel: TNG50. The dashed lines illustrate the selection of the sample of ETGs with masses 1010.3 ≤ M*/M⊙ ≤ 1012, and red color (Sect. 5.1). The colors indicate the axis ratios p(1 Re).
In the text