The stellar halos of ETGs in the IllustrisTNG simulations: the photometric and kinematic diversity of galaxies at large radii
C. Pulsoni, O. Gerhard, M. Arnaboldi, A. Pillepich, D. Nelson, L. Hernquist, V. Springel
AAstronomy & Astrophysics manuscript no. TNGkinematics c (cid:13)
ESO 2020June 19, 2020
The stellar halos of ETGs in the IllustrisTNG simulations: thephotometric and kinematic diversity of galaxies at large radii
C. Pulsoni , , O. Gerhard , M. Arnaboldi , A. Pillepich , D. Nelson , L. Hernquist , and V. Springel Max-Planck-Institut für extraterrestrische Physik, Giessenbachstraße, 85748 Garching, Germany Excellence Cluster Universe, Boltzmannstraße 2, 85748 Garching, Germany European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748 Garching, Germany Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA17 June 2020
ABSTRACT
Context.
Early-type galaxies (ETGs) are found to follow a wide variety of merger and accretion histories in cosmological simulations.
Aims.
We characterize the photometric and kinematic properties of simulated ETG stellar halos, and compare them to observations.
Methods.
We select a sample of 1114 ETGs in the TNG100 simulation, and 80 in the higher-resolution TNG50. These ETGs span astellar mass range of 10 . − M (cid:12) and are selected within the range of g − r colour and λ -ellipticity diagram populated by observedETGs. We determine photometric parameters, intrinsic shapes, and kinematic observables in their extended stellar halos. We comparethe results with central IFU kinematics and ePN.S planetary nebula velocity fields at large radii, study the variation in kinematics fromcenter to halo, and connect it to a change in the intrinsic shape of the galaxies. Results.
We find that the simulated galaxy sample reproduces the diversity of kinematic properties observed in ETG halos. Simulatedfast rotators (FRs) divide almost evenly in one third having flat λ profiles and high halo rotational support, a third with gentlydecreasing profiles, and another third with low halo rotation. However, the peak of rotation occurs at larger R than in observed ETGsamples. Slow rotators (SRs) tend to have increased rotation in the outskirts, with half of them exceeding λ = .
2. For M ∗ > . M (cid:12) halo rotation is unimportant. A similar variety of properties is found for the stellar halo intrinsic shapes. Rotational support and shapeare deeply related: the kinematic transition to lower rotational support is accompanied by a change towards rounder intrinsic shape.Triaxiality in the halos of FRs increases outwards and with stellar mass. Simulated SRs have relatively constant triaxiality profiles. Conclusions.
Simulated stellar halos show a large variety of structural properties, with quantitative but no clear qualitative di ff erencesbetween FRs and SRs. At the same stellar mass, stellar halo properties show a more gradual transition and significant overlap betweenthe two families, despite the clear bimodality in the central regions. This is in agreement with observations of extended photometryand kinematics. Key words.
Galaxies: elliptical and lenticular, cD – Galaxies: halos – Galaxies: kinematics and dynamics – Galaxies: photometry –Galaxies: structure
1. Introduction
The family of early type galaxies (ETGs) encompasses galax-ies that have typically ceased their star formation at early times,with red colors and small amounts of cold gas and dust to-day, and that mainly consist of elliptical and lenticular galax-ies (Roberts & Haynes 1994; Kau ff mann et al. 2003; Blan-ton & Moustakas 2009). Ellipticals essentially divide into twoclasses with distinct physical properties (e.g. Kormendy et al.2009, and references therein): those with low to intermediatemasses and coreless luminosity profiles, that rotate rapidly, arerelatively isotropic and oblate-spheroidal, and have high ellip-ticities and disky-distorted isophotes; and those which are fre-quently among the most massive galaxies, with cored profiles,mostly non-rotating, anisotropic and triaxial, relatively rounderthan than coreless systems, and with boxy-distorted isophotes.Thus the dichotomy in the light distributions of the ellipticalsroughly corresponds to di ff erent kinematic properties, with core-less disky objects being rotationally supported, and cored boxygalaxies having low rotation (Bender 1987). With the advent ofintegral field spectroscopy (IFS) the classification of elliptical galaxies has shifted to a kinematics-based division between fastrotators (FR) and slow rotators (SR) (Emsellem et al. 2011; Gra-ham et al. 2018). In particular low mass, coreless, FR ellipticalsshare similar properties with lenticular galaxies, which are arealso included in the FR family, while massive cored ellipticalsare typically SRs.The formation of massive ETGs is believed to have occurredin two phases (e.g. Oser et al. 2010). In an initial assemblystage, gas collapses in dark matter halos and forms stars in abrief intense burst which is quickly quenched (e.g. Thomas et al.2005; Conroy et al. 2014; Peng et al. 2010). Present-day sim-ulations agree in that the progenitors of FR and SR at thesehigh redshifts are indistinguishable (Penoyre et al. 2017; Lagoset al. 2017; Schulze et al. 2018, with Illustris, Eagle, and Mag-neticum, respectively). At z (cid:46) ffi ciently in size through a se-ries of merger episodes, mainly dry minor mergers (Naab et al.2009; Johansson et al. 2012), which enrich the galaxies with ac-creted (ex-situ) stars. The Λ CDM cosmology predicts that struc-tures form hierarchically, in which more massive systems formthrough the accretion of less massive objects. This means that
Article number, page 1 of 25 a r X i v : . [ a s t r o - ph . GA ] J un & A proofs: manuscript no. TNGkinematics more massive galaxies can have accreted fractions larger than80%, while lower mass galaxies are mostly made of in-situ stars,and the accreted components are mainly deposited in the out-skirts (Rodriguez-Gomez et al. 2016; Pillepich et al. 2018a). Theslow / fast rotators (i.e. the core / coreless) classes result from dif-ferent formation pathways characterized by di ff erent numbers ofmergers, merger mass ratio, timing, and gas fractions (Naab et al.2014; Penoyre et al. 2017, see also the discussion in Kormendyet al. 2009), although the details still depend on the star forma-tion and AGN feedback models adopted by the numerical mod-els (Naab & Ostriker 2017). In general, the result of a formationhistory dominated by gas dissipation is most likely a corelessFR, while dry major mergers often result in SRs.The two-phase formation scenario is supported both by ob-servations of compact red nuggets at z ∼
2, a factor of 2-4smaller than present day ellipticals (Daddi et al. 2005; Trujilloet al. 2007; van Dokkum et al. 2008), and by evidence for a sub-sequent rapid size growth with little or no star formation (e.g.van Dokkum et al. 2010; Damjanov et al. 2011; van der Welet al. 2014; Buitrago et al. 2017). The merger driven size growthis supported by the observed rate of mergers from pair countsand identified interacting galaxies (Hopkins et al. 2008; Robainaet al. 2010), as well as the observed tidal debris from recent ac-cretion events in the halos of many galaxies (e.g. Malin & Carter1983; Janowiecki et al. 2010; Longobardi et al. 2015b; Iodiceet al. 2017; Mancillas et al. 2019).A consequence of the two-phase formation is that ETGs arelayered structures in which the central regions are the remnantsof the stars formed in-situ, while the external stellar halos areprincipally made of accreted material (Bullock & Johnston 2005;Cooper et al. 2010), even though the details strongly depend onstellar mass (Pillepich et al. 2018a). Because of the di ff erent na-ture of the stellar halos, galaxies are expected to show signifi-cant variation of physical properties from central regions to largeradii, 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 (Coc-cato et al. 2009; Romanowsky & Fall 2012; Arnold et al. 2014;Foster et al. 2016).Kinematic measurements in the outer halos of ETGs requirealternative kinematic tracers to overcome the limitations fromthe faint surface brightness in these regions, such as planetarynebulae (PNe) (e.g., the ePN.S survey, Arnaboldi et al. 2017, seeSect. 3), or globular clusters (e.g., the SLUGGS survey, Brodieet al. 2014). Recently Pulsoni et al. (2018) found evidence fromthe ePN.S survey for a kinematic transition between the centralregions and the outskirts of ETGs. Despite the FR / SR dichotomyof their central regions, these ETG halos display a variety ofkinematic behaviors. A considerable fraction of the ePN.S FRsshow reduced rotational support at large radii, which has beeninterpreted as the fading of a rotating, disk-like component intoa more dispersion dominated spheroid; almost half of the FRsample shows kinematic twists or misalignments at large radii,indicating a variation of their intrinsic shapes, from oblate at thecenter to triaxial in the halo. SRs instead have increased rota-tional support at large radii. While a smaller group of FRs standsout for having particularly high V /σ ratio in the halo, most of theePN.S FRs and SRs have similar V /σ ratio in the halo regions.These results suggest the idea that at large radii the dynamicalstructure of these galaxies could be much more similar than intheir high-density centers: if halos are mainly formed from ac-creted material, their common origin would explain their sim-ilarities. The radii of the observed kinematic transitions to the halo and their dependence on the galaxies’ stellar mass seem tosupport such an interpretation.Up to date only a few studies of the kinematic propertiesof stellar halos in simulations are available in the literature. Wuet al. (2014) analysed the kinematics of 42 cosmological zoomsimulations of galaxies and found a variety of V /σ profile shapes(rising, flat, or with a maximum), in agreement with observa-tions. However these early simulations did not reproduce thewhole spectrum of properties of observed FRs, especially the fastrotating and extended disks (Emsellem et al. 2011; Pulsoni et al.2018). Recently, Schulze et al. (2020) using the MagneticumPathfinder simulations showed that these simulations reproducethe observed kinematic properties of galaxies more closely, andthat extended kinematics is a valuable tool for gaining insightinto galaxy accretion histories. They also found that the kine-matic transition radius is a good estimator of radius of the tran-sition between in-situ and ex-situ dominated regions for a subsetof galaxies with decreasing V /σ profiles, especially those thatdid not undergo major mergers in their evolution.The goal of this paper is to better understand the structuralchanges between the centers and stellar halos of ETGs with alarge and well-resolved sample of simulated galaxies. We studythe stellar halo structure, i.e., the rotational support and intrin-sic shapes of the simulated galaxies, we compare the resultswith observations, and we investigate how the radial variationsin rotational support relate to changes in the halo shapes. Weuse the IllustrisTNG simulations (Springel et al. 2018; Pillepichet al. 2018a; Naiman et al. 2018; Marinacci et al. 2018; Nelsonet al. 2018, 2019b), a suite of magnetohydrodynamical simula-tions that models the formation and evolution of galaxies withinthe Λ CDM paradigm. It builds and improves upon the Illustrissimulation (Genel et al. 2014; Vogelsberger et al. 2014), using arefined galaxy formation model. For this work we consider twocosmological volumes with side lengths ∼
100 Mpc and ∼ ff erent ETG surveys select their samples, and howphysical quantities are measured. Section 4 then describes and il-lustrates our methods to derive photometric and kinematic mea-surements for the simulated galaxies. After selecting the sampleof ETGs from the TNG100 and TNG50 simulations (Section 5),we proceed to show the photometric results in Section 6 and thekinematic results in Section 7. Section 8 relates the variation inthe kinematic properties from central regions to halos to the par-allel changes in the intrinsic structure of galaxies. In a compan-ion paper we will explore the dependence of these properties onthe accretion history of galaxies. Finally Section 9 summarizesour conclusions.
2. The IllustrisTNG simulations
The IllustrisTNG simulations are a new generation of cosmo-logical magnetohydrodynamical simulations using the movingmesh code arepo (Springel 2010). Compared to the previous Il-lustris simulations, they include improvements in the models forchemical enrichment, stellar and black hole feedback, and intro-duce new physics such as the growth and amplification of seedmagnetic fields.The baryonic physics model contains a new implementa-tion of black hole feedback (Weinberger et al. 2017), as well
Article number, page 2 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations
Table 1.
Table of physical and numerical parameters for TNG50 andTNG100. These are the volume of the box, the initial number of par-ticles (gas cells and dark matter particles), the target baryon mass, thedark matter particle mass, and the z = Run Volume N part m baryons m DM r soft name [Mpc ] [10 M (cid:12) ] [10 M (cid:12) ] [pc]TNG50 51 . × .
85 4 . . ×
14 75 738as updates to the galactic wind feedback, stellar evolution andgas chemical enrichment models (Pillepich et al. 2018b). Thesemodifications, in particular those for the two feedback mecha-nisms, were required to alleviate some of the tensions betweenIllustris and observations, such as the large galaxy stellar massesbelow the knee of the galaxy stellar mass function and the gasfractions within group-mass halos. They in turn also improve onthe too large stellar sizes of galaxies and the lack of a stronggalaxy color bimodality at intermediate and high galaxy massesin Illustris (Nelson et al. 2015).The IllustrisTNG fiducial model was chosen by assessing theoutcome of many di ff erent models against the original Illustrisby using additional observables, specifically the halo gas massfraction and the galaxy half-mass radii, with respect to thoseused to calibrate the Illustris model against observational find-ings, i.e. the star formation rate density as a function of z , thegalaxy stellar mass function at z =
0, the z = z = ff ective winds in TNG reduce thestar formation at all masses and all times, resulting in a sup-pressed z = M ∗ (cid:46) M (cid:12) , andsmaller galaxy sizes (Pillepich et al. 2018b). Overall the TNGmodel has been demonstrated to agree satisfactorily with manyobservational constraints (e.g. Genel et al. 2018; Nelson et al.2018) and to return a reasonable mix of morphological galaxytypes (Rodriguez-Gomez et al. 2019).In this study we consider two simulation runs, TNG100 andTNG50, which are the two highest resolution realizations ofthe IllustrisTNG intermediate and small cosmological volumes.TNG100 has a volume and resolution comparable with Illus-tris, while TNG50 reaches resolutions typical of zoom-in simu-lations. Table 1 summarizes and compares the characteristic pa-rameters of the two simulations.The TNG model is calibrated at the resolution of TNG100and all the TNG runs adopt identical galaxy formation modelswith parameters that are independent of particle mass and spatialresolution ("strong resolution convergence" according to Schayeet al. 2015). This imposition results in some of the properties ofthe simulated galaxies being resolution dependent. As discussedby Pillepich et al. (2018b), this can be primarily explained bythe fact that better resolution allows the sampling of higher gasdensities, hence more gas mass is eligible for star formation andthe star formation rate accelerates. This means that, for example,at progressively better resolution, galaxies tend to have increasedstellar masses at fixed halo mass and smaller sizes at fixed stellar mass (see also Pillepich et al. 2019 for a quantification of thesee ff ects).
3. Observed parameters of ETGs
In this paper we compare the kinematic results for the centralregions of the simulated TNG galaxies with IFS measurementsfrom 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 di ffi cultto obtain for ETGs, and therefore discrete kinematic tracers suchas planetary nebulae (PNe) and globular clusters (GCs) are typ-ically used to overcome the limitations of absorption line spec-troscopy, which is restricted to the central 1-2 R e . PNe are estab-lished 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 (Longobardiet al. 2015a; Hartke et al. 2018). Since they are drawn from themain stellar population, their kinematics traces the bulk of thehost-galaxy stars, and are directly comparable to integrated lightmeasurements. The relation between GCs and the underlyinggalaxy stellar population is less straightforward (Forbes & Re-mus 2018). In general GCs do not necessarily follow the surfacebrightness distribution and kinematics of the stars (e.g. Brodie& Strader 2006; Coccato et al. 2013; Veljanoski et al. 2014), al-though there is growing evidence for red, metal-rich GCs to betracers of the host galaxy properties (Fahrion et al. 2020; Dolfiet al. 2020). Therefore we here compare the kinematics of thesimulated galaxies and their stellar halos at large radii with PNkinematic results from the ePN.S early-type galaxy survey (Arn-aboldi et al. 2017, Arnaboldi et. al., in prep.).Below we describe the sample properties for the di ff erentsurveys, and we give details and sources of the measured quan-tities used though out this paper. Sample properties - The Atlas3D survey selected ETGs froma volume-limited sampe of galaxies, with distance within 42Mpc, and sky declination δ such that | δ − ◦ | < ◦ ), brighterthan M K < − . . < z < . M ∗ = − M (cid:12) )and environment (field, group, and clusters). This sample is notmorphologically selected, but we use the data from van de Sandeet al. (2017) where the quality cuts and the imposed thresholdon the velocity dispersion σ >
70 km / s bias the sample to-wards 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 Sur-vey (SDSS) Data Release 8, Aihara et al. 2011) at low redshift(0 . < z < . M ∗ = − ; in this paper we will compareonly with MANGA’s galaxies classified as ellipticals or lenticu-lars as in Graham et al. (2018). The MASSIVE survey (Ma et al.2014) targets all the most massive ETGs ( M ∗ (cid:38) . M (cid:12) ) withina distance of 108 Mpc. Finally, the ePN.S sample of ETGs ismagnitude limited M K (cid:46) −
23, and includes objects with di ff er-ent structural parameters. This ensures the sample to be a rep- http: // & A proofs: manuscript no. TNGkinematics resentative group of nearby ETGs. The ePN.S kinematic results(Pulsoni et al. 2018) combine PN kinematics in the halos withliterature absorption line data for the central regions.
Colors - The MANGA galaxies, and most of the Atlas3D andMASSIVE objects have measured g − r colors in the NSA cata-log. For the SAMI galaxies van de Sande et al. (2017) report g − i colors, which we convert to g − r using the transformation equa-tion derived in App. A. For all of the ePN.S sample, and someof the Atlas3D and MASSIVE galaxies that are not in the NSAcatalog, we use B − V colors corrected for galactic extinctionfrom the Hyperleda catalog (Makarov et al. 2014), and convertto g − r colors using the relations in App. A. Sizes - For the Atlas3D sample we use the e ff ective radii ( R e )values in Table 3 of (Cappellari et al. 2011). Those for the MAS-SIVE galaxies are from Ma et al. (2014, Table 3), where weadopt the NSA measurements, where available, or the 2MASSvalues corrected using their Equation 4. The data for MANGAare from Graham et al. (2018). For SAMI we use the data pre-sented in van de Sande et al. (2017), and we circularise the ef-fective semi-major axis by using the reported value for the el-lipticity. The half light radii for the ePN.S galaxies are in Table2 of Pulsoni et al. (2018). These are e ff ective semi-major axisdistances measured from the most extended photometric profilesavailable from the literature, extrapolated to very large radii witha Sérsic fit. The ellipticity assumed is in their table 1. Section 6.1discusses the systematic e ff ects in comparing observed e ff ectiveradii and half-mass radii in simulated galaxies. Stellar masses - The IllustrisTNG model assumes a Chabrier(2003) initial mass function (IMF). The stellar masses for theSAMI survey in van de Sande et al. (2017) are derived usinga color–mass relation, and a Chabrier IMF. For Atlas3D, MAS-SIVE, and MANGA we use the total absolute K-band luminosity M K from the same tables referenced above, which are derivedfrom the 2MASS extended source catalog (Jarrett et al. 2003),and already corrected for galactic extinction. The luminosities M K are then corrected for missing flux as in Scott et al. (2013), M K corr = . M K + .
53, and converted to stellar masses with theformula from van de Sande et al. (2019)log M ∗ = . − . M K corr + , (1)which uses the stellar population model-based mass-to-light ra-tio from Cappellari et al. (2013), their [log( M / L ) Salp ], convertedto a Chabrier IMF. The missing flux correction takes into accountthe over-subtraction of the sky background by the 2MASS datareduction pipeline (Schombert & Smith 2012) and the limited4 R e aperture of the 2MASS measurement.For the ePN.S sample we derive stellar masses using inte-grated luminosities from the most extended photometric profilesavailable in the literature, extrapolated to infinity with a Sérsic fit(references in Pulsoni et al. 2018). We convert the integrated val-ues to stellar masses by using the non-dereddened relations be-tween colors and mass-to-light ratios for ellipticals and S0 galax-ies from García-Benito et al. (2019), which assume a ChabrierIMF.There are several sources of errors in the stellar mass esti-mates of observed galaxies. The uncertainty in the magnitudesderived from the 2MASS photometry are typically ∼ .
25 mag(Scott et al. 2013). The uncertainty in the distances typicallytranslate into an error of 0 . . ∼ . ∼ . http: // nosity, and hence the total stellar mass, can be underestimatedif the photometry is not deep enough to measure the faint sur-face brightness of the stellar halos, especially in massive galaxieswith large Sérsic indices or described by multiple Sérsic compo-nents. Since the stellar masses of the simulated galaxies are eval-uated using the total bound stellar mass (Section 5.1), this maycause a systematic di ff erence between observed and simulatedstellar masses at the high mass end; see also Section 6.1. Ellipticities - For the Atlas3D galaxies we use the ellipticity ε measurements within 1 R e reported in Table B1 of Emsellemet al. (2011). 17 out of 260 Atlas3D objects have obvious barcomponents: for these cases the ellipticity is measured at largerradii (typically 2.5 - 3 R e ). Ellipticities for the SAMI galaxies arefrom van de Sande et al. (2017), and are average ellipticities ofthe galaxies within 1 R e . MANGA’s ellipticities from Grahamet al. (2018) are also measurements within the 1 R e isophote,while for the MASSIVE sample Veale et al. (2017) uses el-lipticities from NSA where available, and from 2MASS other-wise, which are globally fitted values. The ellipticity profilesfor 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 elliptic-ities used by di ff erent surveys for the same galaxies can di ff erwithin a root-mean-square scatter of ∼ .
05 (see e.g. Veale et al.2017 and figure 2 from Graham et al. 2018).
Angular momentum parameters λ e - The parameter λ e is de-rived in the di ff erent surveys using di ff erent integration areas.While Emsellem et al. (2011, Atlas3D) and Veale et al. (2017,MASSIVE) use circular apertures of radius R e , van de Sandeet al. (2017, SAMI) prefer elliptical apertures with semi-majoraxis R e , and Graham et al. (2018, MANGA) integrate over thehalf-light ellipse (an ellipse covering the same area as a circlewith radius R e , i.e. with semi-major axis R e √ − ε , where ε isthe ellipticity).The uncertainties on the measured λ e for the Atlas3D galax-ies are generally small, ∆ λ e (cid:39) .
01 (Emsellem et al. 2011).Similar errors apply for the MASSIVE sample, ∆ λ e (cid:46) . ff ects, which generally tend to system-atically decrease λ e . In the SAMI galaxies, for a typical seeingof 2 arcsec, van de Sande et al. (2017) find that measurementerrors ( ∆ λ e ∼ .
01) and seeing e ff ects cancel out for galax-ies with λ e < .
2, while for λ e > . ff ect and causes a median decrease in λ e of 0 .
05. For theMANGA regular rotators in the cleaned sample, Graham et al.(2018) estimate mean ∆ λ e = [0 . , − . ∆ λ e = [0 . , − . V /σ profiles - The V /σ profiles for the Atlas3D and theePN.S galaxies are derived from the ratio of the rotation velocity V rot and the azimuthally averaged velocity dispersion in ellipticalradial bins. For the Atlas3D galaxies we apply the procedure de-scribed in Sect. 4.4 directly to the velocity fields from Emsellemet al. (2004) and Cappellari et al. (2011), giving a median erroron V /σ of the order of 0 . V rot and σ from kinemetry analysis on IFS data from Krajnovi´c et al.(2008); Krajnovi´c et al. (2011); Foster et al. (2016), when avail-able. 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 themeasurement uncertainties on the V /σ profiles are dominated Article number, page 4 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations by the statistical error on the PN velocity fields. The median ∆ ( V /σ ) = .
4. Methods: IllustrisTNG photometry andkinematics
In this section we describe the method for measuring photometryand kinematics in the IllustrisTNG galaxies. For each simulatedgalaxy we define a coordinate system ( x , y , z ) aligned with theaxes of the simulation box, and centered at the position of themost bound particle in the galaxy. Galaxies are observed bothedge-on and along a random fixed line-of-sight (LOS) direction.The edge-on projection is obtained by rotating the particles ac-cording to the principal axes of the moment of the inertia tensor I i j I i j = (cid:80) n M n x n , i x n , j (cid:80) n M n , (2)where the sum is performed over the 50% most bound stellar par-ticles; x n , i is their coordinates, M n their mass. The random LOSdirection is arbitrarily chosen to be the z axis of the simulationbox. 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 reservecapital letters for the corresponding 2D quantities projected onthe sky. The coordinate r indicates the intrinsic semi-major axisdistance, while R indicates the projected semi-major axis dis-tance.For any projection, we rotate the galaxies so that the X axiscorresponds to the projected major axis, and the Y axis to theprojected minor axis. This is done by evaluating the inertia ten-sor in Eq.2 using the 2D projected coordinates, and summingover the 50% most bound particles. We choose to weight quan-tities by the mass and not by luminosity, as the former are nota ff ected by uncertainties from stellar population modeling andattenuation e ff ects e.g. from dust. The di ff erence between massweighted and luminosity weighted quantities, such as in the K band, is generally small for old stellar populations (e.g. Forbeset al. 2008). Radial profiles are shown in units of e ff ective radii R e , which are evaluated as described in Sect. 6.1. The three-dimensional intrinsic shapes of the galaxies are eval-uated by diagonalizing the inertia tensor I i j in Eq. (2), summedover stellar particles enclosed in elliptical shells. This definitionof I i , j without any weight factors is shown by Zemp et al. (2011)to be the least biased method for measuring the local intrinsicshape of a distribution of particles, and we refer to their work fora 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 squareroot of the ratio of its eigenvalues give the axis ratios p and q (with p ≥ q ) of the principal axes, the eigenvectors their direc-tions. The spherical shell is subsequently deformed to a homeoidof semi-axes a = r , b = pa and c = qa . We repeat the procedureiteratively until the homeoid is adjusted to the iso-density sur-face, and the fractional di ff erence between two iteration steps inboth axis ratios is smaller than 1%. The values of p and q asfunctions of the principal major axis length r give the intrinsicshape profiles of the galaxies. We require a minimum number of1000 particles in each shell as suggested by Zemp et al. (2011),which assures small errors from particle statistics, and, at thesame time, the possibility of measuring intrinsic shape profiles q ( r ) , p ( r ) TNG100-511175 qp T ( r ) TNG100-511175 [random LOS]TNG100-68 [edge-on] P A p h o t Fig. 1.
Photometric measurements.
Top: intrinsic shape of an examplegalaxy from TNG100 as a function of the intrinsic major axis distance r / R e ; axis ratios and triaxiality parameter are shown in separate panels.This simulated galaxy is oblate with q ∼ . p ∼ .
95 in the cen-tral 3 R e , and becomes near-prolate ( T > .
7) in the outer halo. In thisgalaxy 9 r soft correspond to 0 . R e . Bottom: ellipticity and photomet-ric position angle profiles for two example galaxies from TNG100 as afunction of the projected major axis distance R / R e . The quantities de-rived from the inertia tensor are shown with solid symbols, those fromthe mock images with open symbols. out to at least 8 R e for ∼
96% of the selected TNG galaxies. Thedirections of the principal axes of the galaxies as functions of thegalactocentric distance r are given by the eigenvectors ˆ e j (with j = a , b , c ) of the inertia tensor.We also use the triaxiality parameter T ( r ) = − p ( r ) − q ( r ) (3)to quantify the intrinsic shape.In App.B we find that shape measurements at 1 R e are af-fected by the resolution of gravitational forces only for the low-est mass galaxies, for which the absolute error on p and q is ∼ . r ∼ r soft , i.e. r ∼ . R e for thelowest mass galaxies, and r ∼ . R e for M ∗ = M (cid:12) , theseresolution e ff ects are negligible, and the error on the shape mea-surements is then due to particle noise and is ∼ .
02 in TNG100.This uncertainty translates into an error of ∆ T = . T pa-rameter for typical values of the axis ratios in fast rotator ETGs(i.e. p = . q = . r = r soft ; at Article number, page 5 of 25 & A proofs: manuscript no. TNGkinematics
Table 2.
Absolute uncertainties on the shape measurements in TNG100galaxies: ∆ p , ∆ q , and ∆ T at 2 r soft and 9 r soft (first two rows), which cor-respond to di ff erent multiples of R e for galaxies of di ff erent mass asshown in the last three rows. Triaxiality measurements are consideredunreliable for r < r soft . The radius r out indicates the mean radius of theoutermost shell that contains at least 1000 particles for at least 95% ofthe galaxies within the indicated stellar mass bins. r ∼ r soft r ≥ r soft ∆ p , ∆ q (cid:46) . ∼ . ∆ T − ∼ . r soft r soft r out log M ∗ / M (cid:12) : 10 . − . . R e . R e ∼ R e log M ∗ / M (cid:12) : 10 . − . . R e R e ∼ R e log M ∗ / M (cid:12) > < . R e (cid:46) R e ∼ R e smaller radii, where ∆ T is larger, we quantify shapes using p and q which are better defined. These results for TNG100 aresummarized in Table 2. For the TNG50 galaxies we expect sim-ilar or lower uncertainties.In the paper we will consider halos as near-oblate when T ≤ .
3, and near-prolate when T > .
7. Halos with intermediatevalues 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 exampleTNG galaxy, normalized by the R e of the edge-on projection.The galaxy shown in the example is close to oblate in the centralregions, with q (1 R e ) = .
43 and p (1 R e ) = .
95 ( T < . q (10 R e ) = . p (10 R e ) = .
83, and triaxiality parameter T > .
7. For thegalaxy shown 9 r soft = . R e . Mass weighted photometry is derived by diagonalizing the 2Dinertia tensor (Eq. (2)) using the projected coordinates for a givenLOS. We use an iterative procedure similar to the one describedin Sect. 4.1 for the 3D intrinsic shape. The square root of the ratioof the two eigenvalues of I i j gives the projected flattening, hencethe projected ellipticity ε ( R ); the components of the eigenvectorsdefine the photometric position angle PA phot ( R ). The zero pointof 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 galax-ies, and obtained very similar results. The bottom panels of Fig.1 shows the ε ( R ) and PA phot ( R ) profiles obtained from the inertiatensor (solid symbols) and from the images (open symbols) fortwo example galaxies. The galaxy TNG100-511175, shown withblue symbols, is the same as the one shown in the top panels ofFig. 1: the increased axis ratio q ( r ) at r ≥ R e is reflected in adecreased projected ellipticity. The example also shows that atlow ellipticities the uncertainty on the measured PA phot ( R ) be-comes larger, as is well known. We quantified that our methodallows us to measure reliably position angles down to ellipticities ε = .
1, where the error ∆ PA phot ( R ) from particle noise is ∼ ◦ .Below 0.1 ∆ PA phot ( R ) increases exponentially when ε decreasestowards 0. For each TNG galaxy we build projected mean velocity and dis-persion fields for two projections (edge-on and random LOS).We use a resolution of 0 . Y / R e T N G - [ r a n d o m L O S ]
100 50 0 50 100Velocity [km/s] 60 80 100 120 140Velocity Dispersion [km/s]
10 5 0 5 10X / Re1050510 Y / R e
10 5 0 5 10X / Re T N G - [ r a n d o m L O S ]
100 50 0 50 100Velocity [km/s] 60 80 100 120 140Velocity Dispersion [km/s] V r o t [ k m / s ] TNG100-511175 [random LOS]smoothedVor. binned [ k m / s ] P A k i n [ d e g ] Fig. 2.
Stellar kinematic measurements.
Top : Voronoi binned mean ve-locity fields for an example TNG100 galaxy, with a central bulge anda relatively massive disk. For this galaxy the random LOS projectionalmost coincides with the edge-on projection. The projected major axisis aligned with the X axis. The data points show the projected ( X , Y )positions of the stellar particles, and are colored according to the meanvelocity and velocity dispersion of the corresponding Voronoi bin asshown by the colorbars. Center : Smoothed velocity and velocity disper-sion fields for the same example galaxy above. The central disk structureis embedded in a spheroidal stellar halo.
Bottom : kinematic parametersof the galaxy shown above, derived from the Voronoi binned velocityfields (orange open circles) and from the smoothed velocity fields (bluepoints).Article number, page 6 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations for a galaxy observed at 20 Mpc, comparable to present day IFSsurveys (e.g. Law et al. 2016, for MANGA). The stellar particlesare binned on a regular spatial grid centered on the galaxy and8 R e wide.The binned data are then combined into Voronoi bins as de-scribed in Cappellari & Copin (2003), so that each bin containsat least 100 stellar particles. In each i-th bin we calculate theprojected mean velocity and the mean velocity dispersion as theweighted averages: V i = (cid:80) n M n , i V n , i (cid:80) n M n , i ; σ i = (cid:118)(cid:116) (cid:80) n M n , i (cid:0) V n , i − V i (cid:1) N i N i − (cid:80) j M n , i , (4)where the index n runs over the particles in the bin, and N i is thenumber of particles in the i-th bin. The top panel of Fig. 2 showsthe result for one example galaxy; the middle and bottom panelsshow the halo kinematics and the derived kinematic parametersas described in the next section. The example illustrates that inthe central regions, where the density of particles is highest, thevelocity field is sampled at the highest resolution. At larger radiithe Voronoi bins combine the data in progressively larger bins inorder to reach the required minimum number of particles.The systemic velocity of the galaxy is derived by fitting aharmonic expansion as in Pulsoni et al. (2018, their section 4)to the central regions (i.e. at R ≤ R e ) of the projected velocityfield. The fitted constant term is then subtracted from the velocityfields V i .From the velocity fields we calculate the angular momentumparameter λ e following the definition of Emsellem et al. (2011) λ e = (cid:80) i M i R circ , i | V i | (cid:80) i M i R circ , i (cid:113) V i + σ i , (5)where the weighting with the flux is substituted here with aweighting with the mass M i of each Voronoi bin of index i, M i = (cid:80) n M n , i , and R circ , i is the circular radius of the i-th bin.The cumulative λ e is derived by summing over all the Voronoibins contained inside an elliptical aperture of semi-major axis R e and flattening given by the ellipticity ε (1 R e ). By comparison,the di ff erential λ ( R ) is summed in elliptical shells. As discussedin App. B the angular momentum parameter is not a ff ected byresolution at R (cid:38) R e for the selected sample of galaxies. The mean velocity and velocity dispersion fields at large radii arederived using the adaptive smoothing kernel technique (Coccatoet al. 2009), used by Pulsoni et al. (2018) to derive halo veloc-ity fields from the discrete velocities of planetary nebulae in theePN.S survey. For the simulated galaxies, the discrete velocitiesof the particles at R > R e are smoothed with a fully adaptivekernel ( A = B = V rot , kinematic po-sition angle PA kin , and velocity dispersion σ profiles derivedfrom the Voronoi binned velocity fields (in orange), and fromthe smoothed velocity fields (in blue). V rot and PA kin are derivedfrom fitting a harmonic expansion as in Pulsoni et al. (2018),and σ is azimuthally averaged in elliptical annuli whose flatten-ing follows the ellipticity profile of the galaxies. The zero point Table 3.
Absolute uncertainties on the kinematic parameters λ ( R ) and V /σ ( R ) for di ff erent particle number at the resolution of TNG100. log M ∗ M (cid:12) : 10 . − . . −
11 11 − ∆ λ (2 R e ) ∼ . ∼ . (cid:46) . ∆ ( V /σ )(2 R e ) ∼ . ∼ . ∼ . ∆ λ (8 R e ) ∼ . ∼ . (cid:46) . ∆ ( V /σ )(8 R e ) ∼ . ∼ . ∼ . PA kin is defined to be the X axis of the galaxies, consistentlywith the zero point of PA phot . Error bars on the σ ( R ) profiles arederived from the standard deviation of the σ values inside eachannulus. The values obtained with the smoothed velocity fieldsare very well consistent with those from the Voronoi binned ve-locity fields.We also evaluated di ff erential λ profiles using Eq. (5), wherethe summation is performed over the Voronoi bins and the parti-cles, each weighted by their mass, in elliptical annuli. We es-timated uncertainties on the di ff erential λ ( R ) and on V /σ ( R )in TNG100 by considering a few kinematically representativegalaxies in three stellar mass bins and studied the kinematicparameter distributions derived from 1000 simulations respec-tively, with particle numbers decreased to the typical numbers atdi ff erent multiples of R e . Table 3 lists the standard deviation ofthe 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 (cid:46) .
4, see Fig. 1, top panels) embedded in a spheroidal halowith high T ( T (cid:38) . V rot is observed to drop, togetherwith the local λ parameter.
5. Selection of the sample of ETGs in theIllustrisTNG simulations
The purpose of this paper is to study the stellar halos of avolume- and stellar mass-limited sample of simulated ETGs, andcompare with observations. Nelson et al. (2018) verified thatTNG100 reproduces well the ( g − r ) color of 10 < M ∗ / M (cid:12) < . galaxies at z =
0, by comparing with the observed distri-bution from SDSS (Strateva et al. 2001). They also showed thatredder galaxies have lower star formation rates, gas fractions,gas metallicities, and older stellar populations, and that they cor-respond to earlier morphological types (their Figure 13).Thus we extract our sample of ETGs from the TNG50 andTNG100 snapshots at z = g − r ) ≥ .
05 log ( M ∗ / M (cid:12) ) + . . (6)For M ∗ we use the total bound stellar mass of the galaxies. Wedo not include any dust extinction model in the calculation of thesimulated colors in order to avoid the contamination from dust-reddened late type galaxies. Even in this case, this sample ofsimulated galaxies unavoidably contains some red disks, whilein Atlas3D some of the disks have been removed (see Sect. 3). Article number, page 7 of 25 & A proofs: manuscript no. TNGkinematics C o l o r g - r [ m a g ] TNG100TNG50 C o l o r g - r [ m a g ] MANGA S0MANGA ESAMI A3DMASSIVEePNS
Fig. 3.
Selection of galaxies in color and stellar mass.
Top: the g − r color - stellar mass diagram of the simulated galaxies in TNG50 andTNG100. Red sequence galaxies are selected above the tilted dashedline, in the mass range 10 . < M ∗ / M (cid:12) < . Bottom:
ETGs fromrecent IFS surveys as indicated. Most of the observed ETGs in this massrange fall in the same red sequence region.
We limited the sample stellar mass range to 10 . ≤ M ∗ ≤ M (cid:12) . This choice assures that the TNG100 galaxies are re-solved by at least 2 × stellar particles. By comparison, theminimum number of stellar particles in the selected TNG50galaxies is 36 × .In addition we impose that the galaxies’ e ff ective radius (seeSect. 6.1) R e ≥ r soft , to guarantee that the region at r = R e iswell resolved for all simulated galaxies. For TNG100 r soft = . z =
0, which excludes 38 galaxies at the low mass end(see Fig. 7). In TNG50 all the galaxies have R e > × r soft , where r soft = .
288 kpc. These criteria select a sample of 2250 galaxiesin TNG100 and 168 galaxies in TNG50.Figure 3 shows the color-stellar mass diagram for the sim-ulated galaxies from TNG100 and TNG50, and for observedgalaxies from several IFS surveys. Our selection criteria arehighlighted with dashed lines. Most of the observed ETGs, in-cluding the SAMI galaxies and the MANGA ellipticals andlenticulars are in the selected region of the diagram.The histograms in the top panel of Fig. 4 show the stellarmass functions for the color-mass-selected samples. The bottompanel instead shows the stellar mass functions of the final sam-ples as defined by adding constraints from the lambda-ellipticitydiagram in Sect. 5.2. The red and hatched histograms show theAtlas3D and ePN.S samples, respectively. Here we consider theAtlas3D sample properties to validate our selection criteria, asthis survey is especially targeted to study a volume-limited sam-ple of ETGs. The ePN.S sample, which will be used to comparewith properties at large radii, is also shown, and it contains onaverage higher mass galaxies. Both TNG50 and TNG100 arein reasonable agreement with Atlas3D. We remark here that amore generous color selection including bluer galaxies wouldproduce a too large number of high ellipticity galaxies especiallyin TNG50. pd f A3DePN.STNG50TNG100 pd f A3DePN.STNG50TNG100
Fig. 4.
Stellar mass function of the TNG galaxies compared with thoseof the Atlas3D and ePN.S surveys.
Top : stellar mass function of theTNG galaxies selected in color and stellar mass, and whose e ff ectiveradii are well-resolved, as described in Sect. 5.1. The stellar mass func-tions of both TNG50 and TNG100 agree well with Atlas3D. Bottom :stellar mass function of the final sample of ETGs, selected in color, stel-lar mass, and intrinsic shape as described in Sect. 5.1 and Sect. 5.2. Theremoval of centrally elongated objects mostly changes the low-masspart of the TNG50 mass function, while that of TNG100 is still in goodagreement with Atlas3D.
In the following, whenever we compare simulated and ob-served galaxy samples, we will apply to the observed galaxiesthe same color and stellar mass selection criteria that we usedfor the TNG sample. λ -ellipticity diagram: fast and slow rotators Figure 5 shows the λ e - ε (1 R e ) diagram for the simulated ETGs inthree stellar mass bins, and compares with observed ETG sam-ples. 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 againthe TNG50 and TNG100 galaxies after the additional selectiondiscussed in this section. The bottom row shows the similar dia-gram for the observed ETG samples (selected in various ways asdescribed in Sect. 3), in the same color and stellar mass regionas defined in Sect. 5.1. Here we also include for comparison thespiral and irregular galaxies from the MANGA sample (markedas LTGs).We observe that a significant fraction of the TNG galaxiesshown 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. Article number, page 8 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations e [ R a n d o m L O S ] p(1Re) M M * < 10 M M M * < 10 M TNG100TNG50 M * M0.2 0.4 0.6 0.80.00.10.20.30.40.50.60.70.80.9 e [ R a n d o m L O S ] TNG100 TNG50 TNG100 TNG50 TNG100 TNG50 e [ R a n d o m L O S ] MANGA EMANGA S0 MANGA LTGsSAMI A3D
MANGA EMANGA S0 MANGA LTGsSAMI A3DMASSIVE
MANGA EMANGA S0 MANGA LTGsSAMI A3DMASSIVE
Fig. 5. λ e - ε (1 R e ) diagrams for TNG and observed galaxy samples, in stellar mass bins. Top row : the sample of TNG50 and TNG100 ETGs,selected by color and mass as in Fig. 3, observed along a random LOS. The galaxies are color coded according to their intermediate to major axisratios p at 1Re. Middle row : the λ -ellipticity diagram for the TNG galaxies after the additional selection p (1 R e ) ≥ . Bottom row : λ -ellipticityplots for observed ETG samples as shown in the legend, including also late-type galaxies from the MANGA survey for comparison. Their λ e and ε (1 R e ) 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): λ e = . √ ε . The magenta line shows the semi-empirical prediction for edge-on axisymmetric rotators with anisotropy parameter δ = . ε intr from Cappellari et al. (2007). After the additional removal of centrally elongated systems, the colour and mass selected TNG ETG samples givea good representation of the observed ETG samples in the λ e - ε (1 R e )-plane, with the exception of the high- λ e ( > .
7) S0 disks specific to theMANGA survey. below the magenta line and with ε (1 R e ) > .
5. By color codingthe galaxies according to their intrinsic axis ratios at r ∼ R e , wefind that these galaxies have elongated, triaxial shapes. Thesesystems occur at all values of λ e , i.e. some rotate as rapidly asthe MANGA disk galaxies, but others do not show any rotation(Fig. C.1).It is possible that some of the rapidly rotating elongated sys-tems are barred galaxies. Rosas-Guevara et al. (2020) showedthat within a dynamically selected sample of disk galaxies theTNG100 simulation produces barred systems in fractions con- sistent with observational results. The majority of these systems,all characterized (per definition) by high rotation, are quenchedand hence will overlap with the colour range of our sample of redgalaxies. Some barred galaxies are also expected to be presentamong the observed ETG samples. For example, in the Atlas3Dsample 7% of the galaxies show a clear bar component. For theseobjects the ε -values shown in Fig. 5 were measured at largerradii, to avoid the influence of the bar on the estimate of ε (Em-sellem et al. 2011). However, if their actual ε (1 R e ) values wereused and placed these objects in the region of the λ e - ε (1 R e ) pop- Article number, page 9 of 25 & A proofs: manuscript no. TNGkinematics ulated by the centrally elongated (at r ∼ R e ) TNG galaxies,their fraction would not be large enough to explain the abun-dance of simulated galaxies in the same region, and none of thesehave λ e < .
2. Therefore the presence of a large fraction of cen-trally elongated galaxies with high ellipticity and intermediateto low λ e in the TNG sample cannot be explained as a simplesample selection bias (note also that resolution e ff ects on the in-trinsic shapes at 1 R e are at most of the order of 0.1, for the lowmass galaxies, see App.B).In App.C we discuss the properties of these galaxies fur-ther, and suggest that they are likely a class of galaxies thatare produced by the simulation but are not present in nature.These galaxies occupy a particular mass range that depends onresolution and they are the reddest systems for their mass. Wefound no similar concentration of elongated systems among thered galaxies in the Illustris simulation, and the λ − ε diagramsfor simulated galaxies in Magneticum (Schulze et al. 2018) andEAGLE (Walo-Martín et al. 2020) do not contain many objectswith large ellipticities and intermediate to low λ e . This indi-cates that the new galaxy formation model in TNG is involved inthe occurrence of these centrally elongated galaxies. The elon-gated components typically extend up to 3 R e and are embed-ded in near-oblate spheroids with a wide range of flattening q ,with lower median value in TNG50 ( q ∼ .
3) than in TNG100( q ∼ . λ e ) ap-proximately uniformly all the way from edge-on λ e = . p < . r ∼ R e . This choice is motivatedby the fact that the intrinsic shape distribution of real galaxies isknown (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 < . λ e − ε observations. The loca-tion of the simulated galaxies in the plane follows closely theAtlas3D, SAMI, and MASSIVE galaxies. In the MANGA sam-ple there is a large fraction S0 galaxies with λ e > . ff erencesin the data analysis, possibly to the beam corrections applied byGraham et al. (2018) on the MANGA data (see discussion inFalcón-Barroso et al. 2019).Figure 6 demonstrates that the distribution of stellar haloproperties which we are interested in, i.e. λ and triaxiality param-eter, are not a ff ected by the sample selection based on p (1 R e ).The distributions do not systematically depend on the intrinsicshape of the central regions of the galaxies. As discussed in acompanion paper, the properties of the galaxies at large radii are pd f T N G e )0.00.51.01.52.02.5 pd f e )01234 T N G Fig. 6.
Distribution of λ parameter and triaxiality in the stellar halos(i.e. at R ≥ R e ) for di ff erent intrinsic shapes in the central regions,parametrized by the intermediate to major axis ratio p (1 R e ). TNG100and TNG50 galaxies are shown separately, in the top and bottom panelsrespectively. These distributions are not systematically dependent on p (1 R e ), and thus the stellar halo spin and triaxiality remain unbiasedafter implementing a sample selection based on p (1 R e ). mainly set by their accretion history and not by the details of thestar formation in the central regions of galaxies.The bottom panel of Fig. 4 shows the stellar mass func-tion for the final sample of ETGs, compared with observations.The stellar mass function of the TNG100 ETGs is still sim-ilar to Atlas3D. For TNG50 the additional selection has ex-cluded a large fraction of red galaxies in the stellar mass range ∼ . − M (cid:12) . This results in a stellar mass function skewedtowards high masses (and so more similar to ePN.S).Henceforth we classify galaxies as slow rotators (SRs) andfast rotators (FRs), using the dividing line introduced by Em-sellem et al. (2011), λ e = . √ ε e , (7)shown in Fig. 5 with the black line: galaxies above this thresholdare FRs, and galaxies below are SRs. To reduce the e ff ects of in-clination, we choose to classify the simulated galaxies using thevalues of λ e and ellipticity for their edge-on projection (shownin Fig. C.2). The sample of ETG galaxies used in the remainder of this paperis extracted from the TNG50 and TNG100 simulations by – selecting galaxies in the stellar mass range 10 . ≤ M ∗ ≤ M (cid:12) and with red ( g − r ) color as in Eq. (6) (see Fig. 3); – excluding a small number of objects with R e < r soft , to as-sure su ffi cient resolution at 1 R e ; – finally, removing a class of centrally elongated, triaxialgalaxies with p (1 R e ) < .
6, which are systems not present
Article number, page 10 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations in the observed ETG samples, that probably became bar-unstable and quenched during the process of (central) diskformation.The selected sample has a distribution of λ e − ε (1 R e ) similarto observed ETGs (Fig. 5) and halo properties that are unbiasedby the selection in intrinsic shape (Fig. 6). The mass functionsof the selected TNG50 and TNG100 ETG samples are given inFig. 4, and the mass-size relations are shown in Fig. 7.
6. Photometric properties of the TNG ETG samples
In this section we study the photometric properties of the se-lected sample of TNG galaxies and how they vary with radius.Section 6.1 discusses the measured galaxy sizes and how ourdefinition of e ff ective radii compares with e ff ective radii inferredfrom ETG photometry. Section 6.2 compares the distribution ofprojected ellipticities at 1 R e with that from ETG surveys and val-idates our sample selection. Section 6.3 studies the TNG elliptic-ity profiles out to the stellar halo, Sect. 6.4 explores the intrinsicshape distribution of stellar halos and its dependence on stellarmass, and Sect. 6.5 investigates the dependence of galaxy triaxi-ality on radius and stellar mass in the simulated samples. FinallySect. 6.6 tests the ability of photometric twist measurements toestablish the underlying triaxiality in the TNG galaxies. We first discuss the adopted measurement of the e ff ective radiusfor the simulated galaxies, which we will use in the paper asgalactocentric distance unit.The e ff ective radius R e is derived for each projection (edge-on or random LOS) of the galaxies by using cumulative mass C i r c u l a r i z e d R e [ k p c ] soft TNG100 soft TNG50TNG100 TNG50 C i r c u l a r i z e d R e [ k p c ] MANGA S0MANGA EA3D MASSIVESAMIePN.S
Fig. 7.
Circularized e ff ective radii, i.e. R e √ − ε , of the selected ETGsamples in the TNG100 and TNG50 simulations as a function of stel-lar mass, and comparison with observations (bottom panel). The blackcurves show the median profile of the distribution of galaxy R e from allsurveys and the 10th and 90th percentiles, in both panels. The dashedlines in the top panel indicate the resolution of the two simulations. profiles in elliptical apertures: R e is the major axis radius of theaperture that contains half of the total bound stellar mass.Figure 7 shows the circularized R e as a function of M ∗ forthe final samples of ETGs, and compares it to the distributionof observed e ff ective radii from the di ff erent surveys. The R e inTNG100 are larger than most of the observed R e at M ∗ (cid:38) . ,but they are in reasonable agreement with the ePN.S measure-ments. TNG50 produces smaller galaxies compared to TNG100and to observations at intermediate stellar masses (Pillepich et al.2019). This is purely a resolution e ff ect as discussed in Sect. 2.On the other hand, comparisons to observed R e strongly dependon the operational definitions of galaxy sizes, as discussed byGenel et al. (2018).Observers measure R e by integrating light profiles fitted tothe bright central regions to large radii. This definition of R e tends to underestimate the size (and at the same time the to-tal stellar mass) of the galaxies if the photometric data are notdeep enough to sample the light distribution in the halos, es-pecially in massive galaxies with high Sérsic indices. Pulsoniet al. (2018) determined R e of the ePN.S galaxies from the mostextended photometric profiles available in the literature, usinga 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 fac-tor of ∼ M ∗ > M (cid:12) .However it does not take into account the possibility of an extrahalo component / intra-group or intra-cluster light (ICL) at largeradii. For the simulated galaxies, defining the stellar content ofthe galaxies as all the bound stellar particles identified by the subfind algorithm, automatically includes also ICL stars in themost massive halos, thus overestimating both R e and M ∗ .To quantify these e ff ects requires separating a galaxy fromthe surrounding ICL. A kinematic separation of the ICL similarto Longobardi et al. (2015a) is beyond the scope of this paper.However, Kluge et al. (2020) recently found that if the ICL com-ponent in bright cluster galaxies is identified as the outer com-ponent of a double Sérsic fit, the radius at which it starts dom-inating is ∼
100 kpc with a very large scatter (5 to 400 kpc intheir figure 16). We evaluated the di ff erences in R e and M ∗ thatwe would obtain if instead of using the whole bound stellar masswe limit the galaxy to the mass within 100 kpc. We find that inTNG100 galaxies with M ∗ < . the e ff ects are negligible; in10 . M (cid:12) < M ∗ < M (cid:12) objects the di ff erences in the derived R e and stellar masses are within 10% and 5% respectively, whilebetween 10 M (cid:12) < M ∗ < . M (cid:12) they are within 30% and15%. At higher masses the di ff erences in R e can be larger than50% and those in M ∗ larger than 25%, with a very large scatter.These e ff ects are half as pronounced in TNG50. A size-stellarmass diagram analogous to Fig. 7 using the 100 kpc aperture in-stead of the total bound mass shows an improved agreement withthe observed R e , but TNG100 galaxies with M ∗ > . are stilllarger on average. This may indicate that TNG100 predicts toolarge sizes for high mass galaxies (see also Genel et al. 2018).Because of the somewhat arbitrary choice of the 100 kpclimit on one hand, and the uncertainties in the observed R e distri-bution on the other (from di ff erences in sample selection, qualityof the photometric data, definition of total stellar light, the mass-to-light ratio to obtain total stellar masses), we choose here todefine R e for the simulated galaxies as the half mass radius ofthe total bound stellar mass, and consider the above uncertain-ties in the discussion of the results where relevant. Article number, page 11 of 25 & A proofs: manuscript no. TNGkinematics pd f Slow rotators
A3DePN.S-1ReTNG50TNG100
A3DePN.S-1ReTNG50TNG100 a) pd f Slow rotators 0.00 0.25 0.50 0.75 [edge on]Fast rotators
TNG-1ReTNG-3ReTNG-6ReTNG-8ReTNG-10Re c) E lli p t i c i t y ePN.S ePN.S Slow rotators ePN.S
Fast rotators0 5 10R / Re0.00.20.40.60.8 E lli p t i c i t y [ r a n d o m L O S ] TNG50 + TNG100
TNG50 + TNG100 b) Fig. 8.
Ellipticity distribution and profiles of the selected TNG ETGs. a) Distribution of the random LOS ε (1 R e ) compared with Atlas3D andePN.S as indicated in the legend. Despite some di ff erences between the Atlas3D and the TNG SRs, the good agreement of FR distributions showsthat the selected sample of TNG galaxies contains a mixture of disk and spheroidal galaxies consistent with Atlas3D. b) Ellipticity profiles for theePN.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 theedge-on projected ellipticity values at di ff erent radii predicted by TNG. SRs have a rather constant ellipticity distribution with radius. The FRshave 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 showsthe tendency of most galaxies to become rounder in their halos. Figure 8a shows the distributions of the ellipticities measuredat 1 R e for the final sample of selected ETGs, compared withAtlas3D and ePN.S. In the top panels are the SRs, and in thebottom the FRs.The TNG50 and TNG100 simulations predict a significantfraction of SR galaxies with ε (1 R e ) > .
4, while the observedSRs are relatively rounder. This is a common feature of currentsimulations (Naab et al. 2014; Schulze et al. 2018), and its originis still to be understood. In the case of the FR class, the ellipticitydistributions are rather flat-topped, and in good agreement withAtlas3D. By comparison the ePN.S sample contains on averagerounder (and also more massive, see the bottom panel of Fig. 4)galaxies, and hence a lower number of disk galaxies: none of theePN.S FRs has ellipticity higher than 0.7.Overall Fig. 8a shows that the selected sample of ETGs con-tains a mixture of galaxy types consistent with observations, witha similar balance between disks and spheroids.
The ellipticity profiles of the TNG galaxies are compared overan extended radial range with those of the ePN.S galaxies in Fig.8b. There we show profiles for randomly selected sub-samples ofthe simulated galaxies. Figure 8c instead shows the distributionof ellipticities at di ff erent radii for the fast and the slow rotatorsseparately.The observed profiles for the ePN.S SRs generally mildlyincrease with radius, reaching ε ∼ . 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 di ff erent radii.Most of the simulated FRs have decreasing ellipticity profileswith radius, while a fraction have high ellipticity also at largeradii, as also shown by the ePN.S galaxies. Thus, Fig. 8c showsthat 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 sup-ports the idea of a change in structure of these galaxies at largeradii. The large range of flattening in the stellar halos indicatesa variety in the stellar halo properties. By comparison, the SRsshow only small structural variations. In this section we quantify the distribution of the galaxy intrinsicshapes at large radii. Here we refer to stellar halo as the outer re-gions of the galaxies, where the physical properties are markedlydi ff erent from those of the central regions. While this region maybegin at di ff erent radii in each ETG, we will see in the next sec-tion that the median triaxiality profiles for our sample reach con-stant values beyond ∼ R e . Hence we measure the stellar halointrinsic shape distributions by deriving the intrinsic axes ratiosin a shell around 8 R e , 1.5 R e thick, which is the maximal radius atwhich also the lowest mass systems contain enough particles toreliably measure intrinsic shapes, see Sect.4.1. Di ff erent choicesof the shell thickness, or slightly di ff erent choices of the radius(for example 7 instead of 8 R e ) at which we measure p and q deliver similar results. Article number, page 12 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations p(8Re) q ( R e ) O b l a t e T r i a x i a l P r o l a t e Spherical
TNG100TNG50 0.00.20.40.60.81.0 T ( R e ) pd f - F R s q(8Re) pd f - S R s q q Y Yq q Y Y 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 R e . TNG50 and TNG100 galax-ies are shown with di ff erent symbols as in the legend. Bottom : Intrinsicshape distribution for the halos ( r ∼ R e ) of the FRs (top panels) andSRs (bottom panels) in mass intervals, shown with di ff erent colors asindicated in the figure. The fitted functions are shown with solid lines,and the fitted parameters are reported in the legend. The vertical dashedlines show the comparison with the photometric model used in Pulsoniet al. (2018) to reproduce the observed photometric twists and averageellipticities for the ePN.S survey. Most low-mass TNG galaxies havenear-oblate stellar halos (top), changing towards more triaxial shapeswith increasing stellar mass (bottom panels). 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 alarge variety of possible shapes, from very flat near-oblate with q ∼ .
3, to prolate with q ∼ p ∼ .
5. The majority of (low-mass) galaxies appear to have near-oblate stellar halos, with alarge scatter in minor to major axis ratio q . The bottom panelsof Fig. 9 show the intrinsic shape distributions for the fast andslow rotators separately. The distributions of the minor to major axis ratio q resemble Gaussians, and fitted as such the FRs havemean µ q ∼ . σ q ∼ .
16 in all stellar mass bins.The SRs have µ q ∼ . σ q ∼ .
15, with a tendency for thehighest mass galaxies to be flatter.The distribution of the intermediate to major axis ratio p canbe approximated by a log-normal distribution in Y = ln(1 − p ).The shape of this distribution shows a dependence on stel-lar mass: at higher stellar masses µ Y increases, together withthe width of the distribution. This means that at higher stellarmasses, in both the FR and SR classes, the fraction of near-oblategalaxies with p ∼ µ q = . µ p = . σ of the mean values obtainedfrom the distribution of simulated FRs. We can study how the intrinsic shapes of galaxies vary as a func-tion of radius by looking at their triaxiality profiles. We recallfrom Sect. 4.1 that because of the error due to resolution e ff ectsin the central regions, T profiles are considered well-definedonly beyond r = r soft , where their error ∆ T ∼ . r > . R e for the lowest mass objects, and for r > . R e for galaxies with M ∗ ∼ M (cid:12) .The left panels of Fig. 10 show the median triaxiality profilesfor FRs and SRs. These median profiles were built by binning thegalaxies according to their values of the triaxiality parameter Tat 8 R e , and are plotted against the intrinsic major axis distance r .The scale radius R e that normalizes the three dimensional radiusis the 2D projected e ff ective radius for the edge-on projection ofeach galaxy. The right panels show the median of the distributionof the triaxiality parameter measured at 8 R e as a function of thestellar mass.FRs are characterized by increasing T profiles, which tendto plateau at r > R e where the TNG galaxies show a broadrange 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 afunction of stellar mass. This is visible in the right hand side ofFig. 10, for FRs and SRs separately. At lower masses the TNGgalaxies have preferentially near-oblate shapes, with T (cid:46) . M ∗ > M (cid:12) there is a non-negligible fraction ofgalaxies with prolate-triaxial halos, even among the FRs. Wenote a systematic di ff erence between TNG50 and TNG100 inthe triaxiality of the SR stellar halos. In TNG50 the SRs tend tobe much more oblate, indicating some higher degree of dissipa-tion involved in their evolution compared to the SRs in TNG100.The statistical significance of this di ff erence is marginal sinceTNG50 contains only 14 SRs. Isophotal twists in photometry are generally considered to besignatures of triaxiality. This is because the projection on the skyof coaxial triaxial ellipsoids with varying axis ratios approximat-ing the constant luminosity / mass surfaces of ETGs can result intwisting isophotes (e.g. Benacchio & Galletta 1980). However,the e ff ects of triaxiality on the PA phot profiles are model depen- Article number, page 13 of 25 & A proofs: manuscript no. TNGkinematics T ( r ) oblatetriaxialprolate M )0.00.20.40.60.81.0 T ( R e ) oblatetriaxialprolate TNG100TNG50 T ( r ) oblatetriaxialprolate M )0.00.20.40.60.81.0 T ( R e ) oblatetriaxialprolate TNG100TNG50
Fig. 10. Left panels : Median triaxiality profiles for the fast rotators (top) and the slow rotators (bottom). The numbers on the left indicate thepercentage of FRs or SRs populating each profile. The vertical dotted lines show r = r soft for TNG100 galaxies with M ∗ ∼ . M (cid:12) , the solidlines show r = r soft for TNG100 galaxies with M ∗ ∼ M (cid:12) . At radii larger than these the profiles are not a ff ected by resolution issues. Rightpanels : The median of the stellar halo triaxiality distribution measured at 8 R e as a function of stellar mass. The shaded areas are contoured bythe 25th and 75th percentiles of the distribution. TNG100 and TNG50 are shown separately, FRs are on top, SRs on the bottom. On averagethe triaxiality parameter increases with radius and with stellar mass. The FRs have increasing T ( r ) profiles, while SRs have either constant ordecreasing profiles. dent, that is they depend on axis ratio, on how much the axisratios changes with radius, as well as on the viewing angles.Figure 11 shows the distribution of maximum photometrictwist, i.e the maximum variation of PA phot ( R ), measured between1 and 8 R e , as a function of the halo triaxiality at 8 R e , i.e. wherethe triaxiality profiles have reached a constant value. Each sym-bol in the diagram is color coded by the median projected ellip-ticity between 1 and 8 R e . Galaxies with ellipticity lower than0.1 have the photometric position angle poorly determined, andare shown with smaller symbols.We observe that the amplitude of the photometric twists isonly 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 smalltwists irrespective of T (8 R e ). This is explained by the fact thatlarge twists can be measured for viewing angles close enough toface-on (Pulsoni et al. 2018, that is lower ellipticities in Fig. 11),at which even small values of T can produce large twists. FromFig. 11 we conclude that the amplitude of the photometric twists is a poor indicator for galaxy triaxiality, and that very small pho-tometric twists are intrinsically compatible with triaxial shapes. Article number, page 14 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations M a x p h o t o m e t r i c t w i s t [ d e g ] oblate triaxial prolate FRSR
Fig. 11.
Maximum measured photometric twist between 1 and 8 R e asa function of the halo triaxiality T (8 R e ) for the TNG50 and TNG100galaxies projected along a random LOS. FRs and SRs are shown withdi ff erent symbols as in the legend. Smaller symbols are used to repre-sent galaxies with ellipticity < .
1, for which the errors on PA phot arelarge. The color-bar indicates the median ellipticity measured between 1and 8 R e . The gray solid line shows the median of the photometric twistdistribution as a function of T (8 R e ); the shaded area encloses the 25thand 75th percentiles. The amplitude of the photometric twists in TNGgalaxies is only weakly dependent on the triaxiality parameter.
7. The kinematics properties
In this section we study how the kinematic properties of the TNGgalaxies vary with radius. In Sect. 7.1 we derive median di ff er-ential λ ( R ) profiles to quantify the variety of kinematic behav-iors in the outskirts of FRs and SRs. Sections 7.2 and 7.3 com-pare the shapes of the V rot /σ profiles of the simulated ETGs withthe observed galaxies in the Atlas3D and ePN.S surveys. FinallySect. 7.4 uses the simulated galaxies to assess kinematic mis-alignments and twists as signatures of triaxial shapes in galaxies. The top panels of Fig. 12 show the median di ff erential λ profilesfor FRs and SRs in their edge-on projection. Galaxies are binnedtogether according to the shape of their profiles. We achieve thisby binning the FRs according to their values of λ at R = R e and at R = R e . For the SRs the shape of the λ ( R ) profiles isgenerally a monotonic function of R : in this case we binned theprofiles according to their λ (7 R e ).Most of the FRs reach their maximum λ ( R ) around 3 R e ; only7% of the galaxies increase λ ( R ) between 3 and 10 R e . FRs di-vide almost evenly among a third (34%) that have flat λ profileswith λ (8 R e ) > .
6, a third (40%) with gently decreasing profilesand 0 . > λ (10 R e ) ≥ .
6, and another third (26%) with very lowrotation in the halo ( λ (10 R e ) ≤ . λ (7 R e ) < .
2) and a half with increased rotationalsupport at large radii compared to the central regions. We ob-serve that a small fraction of the SRs (5% of the TNG100 SRsand 15% of the TNG50 SRs) reach very high values of λ at largeradii ( λ (7 R e ) > . λ profiles to observed SRs like NGC 3608 (Pul-soni et al. 2018). The others are galaxies with a clear extendeddisk structure characterized by rapid rotation and low velocitydispersion, but whose central kinematics is dominated by a non-rotating bulge. There are no observed counterparts for the latterin both the ePN.S (33 galaxies) and the SLUGGS surveys (Fosteret al. 2016, 25 galaxies, of which 18 are in common with ePN.S).A larger sample of galaxies with extended kinematic data wouldbe needed to confirm or rule out these objects.For the SRs we note a mismatch in the amount of halo rota-tional support between TNG50 and TNG100, analogous to thedi ff erence observed for the halo triaxiality (Sect. 6.5): on aver-age TNG50 SR halos rotate faster and have more oblate shapes.These di ff erences in halo properties might be due to the depen-dence of the galaxy formation model on the resolution of thesimulations, although the number of SRs in TNG50 (only 14galaxies) is too small to draw strong conclusions.For both FRs and SRs, and in both TNG50 and TNG100, thestellar halo rotational support depends weakly on stellar massup to M ∗ ∼ . M (cid:12) (Fig. 12). However, at high stellar massesthe fraction of galaxies with significant rotation in the halo de-creases, so that at M ∗ > . M (cid:12) most of the galaxies have nonrotating halos. The broad range of possible λ ( R ) profile shapesin Fig. 12 shows that the IllustrisTNG simulations encompass, ifnot exceed, the observed variety of halo rotational support foundin the ePN.S survey, of which one of the key results was the largekinematic diversity of stellar halos. Figure 13 shows the median V rot /σ profiles of the TNG100 andTNG50 galaxies compared with median profiles from Atlas3Dand individual galaxy profiles from the ePN.S sample for 0 ≤ R ≤ R e . Here we normalize the radii by the circularized R e , i.e. R e √ − (cid:15) (1 R e ), for an appropriate comparison with the Atlas3D R e . The profiles of the simulated SRs are similar to the observedprofiles. The FRs instead show a di ff erence in how quickly V rot /σ rises with radius: observed FR galaxies have on aver-age more steeply rising V rot /σ profiles than the simulated ETGs.Very few TNG FRs reach V rot /σ ≥ R e compared to theAtlas3D galaxies, and almost none exceeds V rot σ (1 R e ) = . ff erent shapes of the V rot /σ profiles in observationsand simulations cannot be explained with resolution e ff ects, asin both TNG50 and TNG100 the V rot /σ profiles tend to peak ata median radius of R ∼ R e . By comparison, the ePN.S FRs tendto peak at smaller fractions of R e , at a median 1.3 R e .This di ff erence between the observed and simulated FRs isnot a consequence of the selection functions of the samples. TheTNG galaxies are selected according to color, mass and inter-mediate to major axis ratio p . The selection in p removes cen-trally elongated galaxies, most of which have intermediate tolow λ (1 R e ) (Fig. 5). The Atlas3D sample, selected as describedin Sect. 3, has some of the disk galaxies removed which as inMANGA (Fig. 5) will be mostly located at high λ (1 R e ). We re- Article number, page 15 of 25 & A proofs: manuscript no. TNGkinematics
R / Re ( R ) log10(Stellar mass / M ) ( R e ) TNG100TNG50
R / Re ( R ) log10(Stellar mass / M ) ( R e ) TNG100TNG50
Fig. 12. Left : Median λ profiles for fast rotators (upper left) and slow rotators (lower left panel) in their edge on projection. The median profilesare built by binning the galaxies according to the shape of their λ profiles; see text. The shaded areas are bounded by the 25th and 75th percentilesof the distributions. The numbers on the left indicate the percentage of FRs or SRs populating each profile. Right : Distribution of λ (8 R e ) for thefull 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. call that in the comparison we consistently matched the colorselection and mass range of the Atlas3D sample to the selectioncriteria adopted for the TNG galaxies. Thus the TNG sampleshould in principle contain a larger number of late type galaxieswith strong disks, i.e. high V rot /σ , by comparison with Atlas3D.Figure 13 shows instead that the TNG100 and TNG50 ETG sam-ples lack galaxies with high rotational support at 1 R e .In Sect. 6.1 we discussed that the e ff ective radii used tonormalize the radial scales in TNG would be expected to besystematically slightly overestimated compared to the observed R e since they are defined using the total bound stellar masswhich is often inaccessible in observations. If taken into account,this e ff ect would increase the gap between simulated and ob-served samples. On the other hand, since the mass-size relationis roughly reproduced in the simulations (Fig. 7), the di ff erentsteepness of the V rot /σ profiles in observations and simulationsimplies a di ff erent distribution of the angular momentum as afunction of radius in the simulated galaxies. This could be dueto a too e ffi cient condensation of the gas into stars that does notallow the gas to dissipate and collapse to small enough radii. Figure 14 compares the relation between V rot /σ in the centralregions and V rot /σ in the stellar halos for the simulated andobserved galaxies. The observed galaxies are from the ePN.Ssurvey (Pulsoni et al. 2018, their figure 9) and are reported inthe top left panel. Their measurements use absorption line dataat 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 cen-tral and right top panels show V rot /σ (6 R e ) versus V rot /σ (1 R e )for TNG100 and TNG50 ETGs, respectively. Galaxies close tothe dashed 1:1 lines show similar rotational support in the cen-tral regions and in the outskirts; galaxies below the lines haveincreased rotational support in their stellar halos; galaxies abovethe lines instead have reduced rotation at large radii.The position of the observed SRs below the one-to-one lineis 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 ) (cid:38)
1, some of which are actually extended disks withprominent bulges at the center. These represent a small fraction
Article number, page 16 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations V r o t / A3DePN.S
Slow rotators Fast rotators
R / (Circularized Re) V r o t / [ r a n d o m L O S ] TNG100TNG50 0 1 2 3
R / (Circularized Re)
Fig. 13.
Median V rot /σ profiles for slow rotators (left) and fast rotators(right panels). Top : observations from Atlas3D and ePN.S. For the for-mer we show median profiles and shaded areas bounded by the 20thand the 80th percentiles of the distribution. For the ePN.S galaxies weshow individual profiles.
Bottom : median profiles and 20%-80% distri-butions for the TNG50 and TNG100 ETGs, projected along a randomLOS. Vertical and horizontal dashed lines are for guiding the eyes incomparing top with bottom panels. This shows that the average V rot /σ profiles of the simulated ETGs are shallower than the observed profiles. of the SR family, and do not have observed counterparts in theePN.S and SLUGGS surveys.For the TNG FRs we find a di ff erent distribution when com-paring rotational support at the same radii: most of these galax-ies fill the bottom half of the diagram, with very few reaching V rot /σ (1 R e ) > .
5, and a large fraction having significant V rot /σ at 6 R e . This di ff erence 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 topeak at larger radii than the observed ETGs, there are almost noobjects that reach V rot /σ > . R e .In the bottom panels of Fig. 14 we show the comparisonat adjusted radii: V rot /σ (3 R e ) (i.e. at the median radius wherethe TNG V rot /σ profiles tend to reach the maximum) versus V rot /σ (8 R e ) (i.e. where the decreasing rotation profiles finallyreach their minimum, see also Fig. 12), and find better agree-ment. Now the FRs spread in the region of the diagram abovethe one-to-one line as they do for the ePN.S observations inthe top left panel of Fig. 14, with a sub-population of objectsshowing V / rot σ (8 R e ) > V /σ (3 R e ). We note that FRs with near-prolate stellar halos occupy mostly the lower left corner withlow V rot /σ (3 R e ) and low V rot /σ (8 R e ). Galaxies with triaxial ha-los tend to distribute on the left side of the diagram, galaxies withoblate halos tend to have higher V rot /σ (8 R e ). It is interesting that,aside for the di ff erences in the central regions at R (cid:46) R e , theobserved galaxies with and without signatures for triaxial stel-lar halos distribute similarly as the simulated triaxial and near-oblate stellar halos, respectively, with the triaxial halos havingon average lower V rot /σ (8 R e ).From Fig. 14 we conclude that even though the quantita-tive details between simulated and observed ETGs galaxies donot agree, the simulations do reproduce the observed kinematictransitions between central regions and outskirts, as well as thevariety of halo kinematic classes. Kinematic signatures of galaxy triaxiality can be found from ob-servations of minor axis rotation, kinematic twists, and misalign-ment 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 therandom LOS projected TNG galaxies. The solid lines show theunweighted root-mean-square deviation from zero of the datapoints as a function of the ellipticity for the FR and the SRclasses separately, computed by mirroring the data points aroundzero (i.e. using [ Ψ , − Ψ ]). The dashed lines show the same quan-tities for observed galaxies in Atlas3D and MANGA (Krajnovi´cet al. 2011; Graham et al. 2018). Simulated and observed galax-ies show very similar trends with ellipticity, with very flat galax-ies being strongly aligned, and kinematic misalignment increas-ing with rounder shapes. Simulated FRs are found to be muchmore aligned than SRs, in agreement with observations.Data points in Fig. 15 are color coded according to the in-termediate to major axis ratio p at 1 R e . At these radii the errorson the axis ratios are at most 0.1 for the low mass systems (seeApp.B), so p (1 R e ) measurements are generally well defined. Wefind that many TNG galaxies, both SRs and FRs show a highdegree of alignment ( | Ψ | < ◦ ) even though they are far frombeing oblate (i.e. with p < . Ψ . Al-though we can not draw conclusions on the shape distributionof the real galaxies, Fig. 15 implies that near-alignment of thekinematic and photometric position angles does not exclude tri-axiality for the IllustrisTNG FR galaxies even at 1 R e .Simulated FRs and SRs do show kinematic position anglevariations as a function of radius. An example is the objectshown in Figs. 1 and 2, a galaxy with central disk embeddedin a triaxial stellar halo, which shows a variation in the directionof rotation corresponding to the sudden change in the triaxialityprofile around r ∼ R e .The top panel of Fig. 16 shows the distribution of the max-imum kinematic twists, i.e. the maximum variation of PA kin ,measured between 1 and 8 R e for all the TNG galaxies, as afunction of the stellar halo triaxiality measured at 8 R e , i.e. wherethe median T profiles are constant with radius (see Fig. 10).We find that the large majority of TNG galaxies have smallkinematic twists compared to typical measurement errors of PA kin from discrete tracers (the median error for the ePN.Sgalaxies is ∼ ◦ ). Aside for a small group of near-oblate galax-ies with counter-rotating disks (twist ∼ ◦ ), the main depen-dence of the amplitude of the twist with the triaxiality parameteris such that galaxies with high T are more likely to have largekinematic twists, as shown by the solid gray line in Fig. 10 rep-resenting the median twist as a function of T . This highlights theimportance of kinematics as a fundamental tool to investigatethe intrinsic structure of galaxies besides photometry alone (seeSect. 6.6). On the other hand Fig. 10 points out that not all thetriaxial and prolate galaxies in IllustrisTNG display kinematictwists. This suggests that also in observations a galaxy’s triaxi-ality may not be easily revealed by kinematic misalignments.IFS kinematics show that the central regions of FR have PA kin and PA phot aligned within ∼
10 degrees, while SRs aregenerally misaligned (Krajnovi´c et al. 2011; Fogarty et al. 2015;Ene et al. 2018; Graham et al. 2018). This di ff erence in the mis-alignment distribution was interpreted as signature of a di ff er- Article number, page 17 of 25 & A proofs: manuscript no. TNGkinematics V / ( R e ) NGC0584NGC0821 NGC1023NGC1316NGC1344NGC1399 NGC2768NGC2974NGC3115NGC3377NGC3379 NGC3384NGC3489NGC3608NGC3923NGC4278 NGC4339NGC4365NGC4374NGC4472NGC4473 NGC4494NGC4552 NGC4564NGC4594NGC4636NGC4649NGC4697 NGC4742NGC5128NGC5846 NGC5866 NGC7457 ePN.S Galaxies
FRsFRs - triaxial haloFRs - small number PNeSRs
SRT(8Re)
T(8Re)
T(8Re) > 0.7
SRT(8Re)
T(8Re)
T(8Re) > 0.7 V / ( R e ) TNG 100 - random LOS projection
T(8Re)
T(8Re)
T(8Re) > 0.7
T(8Re)
T(8Re)
T(8Re) > 0.7
Fig. 14. V /σ ratio in the central regions compared with the V /σ at large radii. Top : V /σ (1 R e ) versus V /σ (6 R e ) for the ePN.S galaxies (left),TNG100 (middle), and TNG50 (right). Bottom : V /σ (3 R e ) versus V /σ (8 R e ) for TNG100 (left) and TNG50 (right). Di ff erent colors in the ePN.Ssample mark the SRs (in red), the FRs with kinematic signatures for triaxial halos (i.e. with kinematic twists or misaligments, in green), and FRswith PA kin aligned with PA phot (in blue). The gray points stand for ePN.S FRs for which the measured V rot /σ ( < R e > ) is a lower limit. For thesimulated galaxies the di ff erent colors mark SRs, and FRs with di ff erent halo T (8 R e ): FRs that have near-oblate halos (in blue); FRs with triaxialhalos (in green); near-prolate FRs (in orange). The gray symbols show simulated galaxies with too few particles at 8 R e for measuring T . The sizeof the data points is proportional to the stellar mass. The gray dashed lines show the 1:1 relation. The TNG simulations reproduce the diversityof 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. ent intrinsic shape distribution for the two classes, with the FRsbeing consistent with having oblate shapes, and the SRs beingmoderately triaxial. This interpretation is not supported by theresults in Fig. 15.Even though the central regions of FRs have been found tohave PA kin well aligned with PA phot , kinematic twists are ob-served in the halos. By extending the kinematic study at largerradii using planetary nebulae, Pulsoni et al. (2018) found thatkinematic twists are relatively frequent in the ePN.S sample ofFRs ( ∼ ff erent sample selection between simulations andePN.S, with the former potentially containing a larger fraction ofdisk galaxies even after the matching of the mass functions. Fig-ures 4 and 8a in fact show that the ePN.S sample is on averagemore massive and contains rounder galaxies than TNG100. Article number, page 18 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations ( R e ) [ d e g ] obsTNGFRsSRs Fig. 15.
Distribution of the misalignments
Ψ = PA kin − PA phot at 1 R e as a function of the ellipticity for both TNG100 and TNG50 galaxies.Data points are color coded according to their intermediate to majoraxis ratio p (1 R e ); open symbols are SRs, filled symbols are FRs. Thesolid lines show the scatter of the misalignments as a function of the el-lipticity for the TNG galaxies, calculated on the mirrored ( Ψ , − Ψ ) data.The dashed lines show the scatter for the MANGA and Atlas3D galax-ies. Red lines are for the SRs, blue lines for the FRs. The triaxialityof the IllustrisTNG FRs is consistent with the near-alignment of theirkinematic and photometric position angles. M a x k i n e m a t i c t w i s t [ d e g ] oblate triaxial prolate pd f TNG100+TNG50TNG matched mass func.TNG matched w/ error 35ePN.S
Fig. 16. Top : maximum kinematic twist measured between 1 and 8 R e asa function of the stellar halo triaxiality T (8 R e ). FRs are shown in blue,SRs in red; both TNG50 and TNG100 samples are shown. The gray lineand shaded band show the median of the twist distribution as a functionof T (8 R e ) and the 25th and 75th percentiles. Even though the major-ity of the TNG galaxies have small kinematic twists, high T galaxiesare more like to display a kinematic twist. Bottom : Distribution of themaximum kinematic twist for the TNG50 and TNG100 samples (greenstep histogram) compared with the ePN.S sample (hatched histogram).The filled green histogram shows the distribution of the TNG galaxiesafter their mass function is matched to that of the ePN.S sample; thegray histogram is the mass-matched TNG distribution convolved with ameasurement error and is consistent with the data.Article number, page 19 of 25 & A proofs: manuscript no. TNGkinematics
8. Stellar halo angular momentum and shape
In Sect. 7.3 we observed that the rotational support in the stellarhalo correlates with the intrinsic shape: stellar halos with high V rot /σ are likely near-oblate, while halos with lower V rot /σ canhave larger triaxiality. In this section we connect the variation ofrotational support in the stellar halos to variations of their intrin-sic shape.The top panel of Fig. 17 shows the di ff erential λ ( R ) (mea-sured for the edge on projection), the axis ratio q ( r ), and thetriaxiality parameter profile T ( r ) for an example galaxy. It has adecreasing λ ( R ) profile in the halo, while q ( r ) and T ( r ) increase:for this object the decreased rotational support marks a varia-tion in the intrinsic shape of the galaxy from relatively flat andnear-oblate at R ∼ R e , where the rotation is highest, to triaxialspheroidal in the outskirts. In observations the drop in rotationof the ePN.S FRs is often related to a decrease in ellipticity. Thisled to the idea that central fast rotating regions of FRs are em-bedded in a more dispersion dominated spheroidal stellar halo.We verify this conclusion in Fig. 17. Here the outskirts ofgalaxies are divided into 6 ellipsoidal shells of major axis r be-tween 3 . R e . This radial range is motivated by the require-ment that r is large enough so that the errors on the T parameterare small for galaxies of all masses, but also such that most of thegalaxies (96%) have enough particles at large radii (8 R e ) in or-der to measure their intrinsic shapes. In each ellipsoidal shell ofmajor axis r and width ∆ r we measure the axis ratio q and the tri-axiality 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 shellsand all TNG galaxies in four galaxy stellar mass bins, as a func-tion of the minor to major axis ratio q ( r ) and of the triaxialityparameter T ( r ).We find indeed a relation between λ and shape. High ro-tational support is related to flattened (i.e. low q ) near-oblate(i.e. low T ) shapes. Where λ decreases q grows towards morespheroidal shapes, although the scatter in possible stellar haloshapes is large ( σ q ∼ . λ , but lower T is com-patible with all values of λ . This means that the decrease in therotational support of the TNG FRs follows a change in the struc-ture of the galaxies, which become more spheroidal in the out-skirts. These outer spheroidal components can span all valuesof triaxiality. By comparison to the FRs, the SRs show smallervariations in both intrinsic shapes (Fig. 10) and λ .Figure 17 confirms the dependence of the stellar halo intrin-sic shape and rotational support on stellar mass already describedin Sects. 6.4, 6.5, and 7.1. Lower mass galaxies host more rota-tionally supported stellar halos with near-oblate shapes. At pro-gressively higher masses the fraction of dispersion dominatedstellar halos increases as well as the fraction of halos with hightriaxiality.It is interesting to note that there is no clear separation in Fig-ure 17 between the stellar halo properties of the FRs and SRs, butrather a continuity of properties among the two classes, with thelow T -high λ extreme dominated by the FRs, the high T -low λ limit by the SRs, and the relative importance of the two popu-lations gradually changing with stellar mass. This result impliesthat there is no qualitative di ff erence between the structure of thegalaxies at large radii despite the bimodality of the FR / SR classi-fication of the centers. The IllustrisTNG galaxies agree with theePN.S observations in that FRs and SRs tend be more similar atlarge radii, especially at intermediate to high masses. ( R ) , q ( r ) (R) [edge on]q(r) TNG50-461785
T(r)
TNG50-461785 oblatetriaxialprolate
Fig. 17. Top : λ ( R ), axis ratio q ( r ) and triaxiality T ( r ) profiles for anexample galaxy. Bottom : λ for the edge-on projection versus minor tomajor axis (left panels) and triaxiality parameter (right panels) in stel-lar mass bins. Each data point is a local measurement within a shell ofmajor axis R = r : each TNG galaxy is represented by ∼ . − R e . FRs and SRs are shown with di ff erentcolors as in the legend. The mass bins are labelled on the right mar-gins. Lower rotational support λ in the stellar halos is related to roundershapes, with a wide range of triaxiality which depends on stellar mass.FRs and SRs exhibit a continuity in stellar halo properties rather than abimodality. The similarity of the overall shapes of the λ ( R ) (or V rot /σ ( R ))profiles and of the ε ( R ) profiles between ePN.S and TNG galax-ies described in in Sects. 6.3 and 7.3 suggests that intrinsic shapevariations similar to those found in the TNG galaxies also existin real galaxies. Article number, page 20 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations
9. Summary and conclusions
In this paper we have analysed the kinematic and photometricproperties of ∼ λ e − ε diagram (Fig. 5). There we excluded simulated objects that donot match the observed properties in the central regions of fastrotators (FRs) and slow rotators (SRs), and that appear as highlycentrally elongated red galaxies. We verified that this does nota ff ect our results on the stellar halo properties of the simulatedgalaxies. The resulting ETG sample has mass-, size-, and centralkinematics distributions consistent with observations.For the selected sample we determined mean velocity fields,kinematic, and photometric profiles, and studied the intrinsicshapes of the simulated galaxies from the central regions intotheir outskirts, up to 15 R e . The purpose of the paper is to studythe kinematic properties of stellar halos and connect them tovariations in the structural properties of their galaxies from thecentral regions to the halos. Our conclusions are as follows:1) The di ff erential λ profiles and the triaxiality profiles (Figs.10 and 12) successfully reproduce the diversity of kinematic andphotometric properties of stellar halos observed in the ePN.S sur-vey. – We find that simulated FRs divide almost evenly among athird with flat λ profiles, a third with gently decreasing pro-files, and another third with very low rotation in the halo.Half of the SRs do not show any rotation in the halo, whilethe 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 lowmass fast rotators) have stellar halos consistent with oblateshapes. – Both halo triaxiality and rotational support are found to de-pend on stellar mass, with higher mass galaxies being moretriaxial and more dispersion dominated at large radii.2) Halo intrinsic shape and rotational support are stronglyrelated (Fig. 17): – High λ is related to flattened oblate shapes. – Where λ decreases with radius, galaxies tend to becomerounder, but with a wide range of triaxiality.3) The FR class in TNG shows the largest variety in stellarhalo properties and the largest variations with radius in both in-trinsic shapes and rotational support. Among these galaxies wecan find rotationally supported stellar halos with flattened oblateshapes, as well as FRs that have central rotating disk-like struc-tures embedded in more spheroidal components. For a subset ofthese the stellar halos can reach high triaxiality values. SRs, bycomparison, display milder changes in structure with radius.4) The TNG FRs exhibit shallower V rot /σ ( R ) profiles thanthe Atlas3D and ePN.S galaxies (Fig. 13). Both TNG50 andTNG100 FRs tend to reach a peak in rotation at a median ra-dius of ∼ R e , compared to ∼ − . R e for the Atlas3D andePN.S galaxies. This result implies a more extended distributionof the angular momentum with radius in the TNG galaxies thanobserved.5) However, even though the V rot /σ ( R ) profiles do not agreequantitatively between simulated and observed ETGs galaxies,the simulations do reproduce the observed kinematic transitions between central regions and outskirts. The similarity betweenthe scaled shapes of the V rot /σ ( R ), and the ε ( R ) profiles betweenePN.S and TNG galaxies (Figs. 8b and 14) suggests that also theobserved variations in the kinematics between central regionsand halos trace changes in the intrinsic structure of the galaxies.6) We find that most of the triaxial TNG galaxies displaymodest 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 bycurrently available data.7) By comparing the distributions of the minor-to-major axisratios q , triaxiality parameters T , and angular momentum pa-rameters λ (Fig. 17), we find that lower rotational support in thestellar halos is related to rounder shapes, with a wide range oftriaxiality which depends on stellar mass. In this there is no qual-itative di ff erence between the FRs and SRs. Rather, despite thebimodality of the central regions, the two classes show a continu-ity 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 di ff erent sides of the distri-bution gradually changes with stellar mass. This is in agreementwith ePN.S observations of ETG halos.In a companion paper we will investigate the dependence ofthe stellar halo parameters on the accretion history of galaxies,and explore the relation between stellar and dark matter haloproperties. Acknowledgements.
We thank the anonymous referee for helpful comments that improved the clar-ity of the paper. C.P. is extremely grateful to F. Hofmann for his support. Thisresearch has made use of the NASA / IPAC Extragalactic Database (NED).
References
Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29Arnaboldi, M., Freeman, K. C., Mendez, R. H., et al. 1996, ApJ, 472, 145Arnaboldi, M., Pulsoni, C., Gerhard, O., & PN. S Consortium. 2017, in IAUSymposium, Vol. 323, Planetary Nebulae: Multi-Wavelength Probes of Stellarand Galactic Evolution, ed. X. Liu, L. Stanghellini, & A. Karakas, 279–283Arnold, J. A., Romanowsky, A. J., Brodie, J. P., et al. 2014, ApJ, 791, 80Bassett, R. & Foster, C. 2019, MNRAS, 487, 2354Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57Benacchio, L. & Galletta, G. 1980, MNRAS, 193, 885Bender, R. 1987, Mitteilungen der Astronomischen Gesellschaft Hamburg, 70,226Binney, J. 1985, MNRAS, 212, 767Blanton, M. R. & Moustakas, J. 2009, ARA&A, 47, 159Brodie, J. P., Romanowsky, A. J., Strader, J., et al. 2014, ApJ, 796, 52Brodie, J. P. & Strader, J. 2006, ARA&A, 44, 193Buitrago, F., Trujillo, I., Curtis-Lake, E., et al. 2017, MNRAS, 466, 4888Bullock, J. S. & Johnston, K. V. 2005, ApJ, 635, 931Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ, 798, 7Cappellari, M. & Copin, Y. 2003, MNRAS, 342, 345Cappellari, M., Emsellem, E., Bacon, R., et al. 2007, MNRAS, 379, 418Cappellari, M., Emsellem, E., Krajnovi´c, D., et al. 2011, MNRAS, 413, 813Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2013, MNRAS, 432, 1862Chabrier, G. 2003, PASP, 115, 763Chua, K. T. E., Pillepich, A., Vogelsberger, M., & Hernquist, L. 2019, MNRAS,484, 476Coccato, L., Arnaboldi, M., & Gerhard, O. 2013, MNRAS, 436, 1322Coccato, L., Gerhard, O., Arnaboldi, M., et al. 2009, MNRAS, 394, 1249Conroy, C., Graves, G. J., & van Dokkum, P. G. 2014, ApJ, 780, 33Cooper, A. P., Cole, S., Frenk, C. S., et al. 2010, MNRAS, 406, 744Cortesi, A., Arnaboldi, M., Coccato, L., et al. 2013, A&A, 549, A115Croom, S. M., Lawrence, J. S., Bland-Hawthorn, J., et al. 2012, MNRAS, 421,872Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680Damjanov, I., Abraham, R. G., Glazebrook, K., et al. 2011, ApJ, 739, L44Dolfi, A., Forbes, D. A., Couch, W. J., et al. 2020, MNRAS, 495, 1321Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817
Article number, page 21 of 25 & A proofs: manuscript no. TNGkinematics
D’Souza, R., Kau ff man, G., Wang, J., & Vegetti, S. 2014, MNRAS, 443, 1433Emsellem, E., Cappellari, M., Krajnovi´c, D., et al. 2011, MNRAS, 414, 888Emsellem, E., Cappellari, M., Peletier, R. F., et al. 2004, MNRAS, 352, 721Ene, I., Ma, C.-P., Veale, M., et al. 2018, MNRAS, 479, 2810Fahrion, K., Lyubenova, M., Hilker, M., et al. 2020, A&A, 637, A26Falcón-Barroso, J., van de Ven, G., Lyubenova, M., et al. 2019, A&A, 632, A59Fogarty, L. M. R., Scott, N., Owers, M. S., et al. 2015, MNRAS, 454, 2050Forbes, D. A., Lasky, P., Graham, A. W., & Spitler, L. 2008, MNRAS, 389, 1924Forbes, D. A. & Remus, R.-S. 2018, MNRAS, 479, 4760Foster, C., Pastorello, N., Roediger, J., et al. 2016, MNRAS, 457, 147Foster, C., van de Sande, J., D’Eugenio, F., et al. 2017, MNRAS, 472, 966Franx, M., Illingworth, G., & de Zeeuw, T. 1991, ApJ, 383, 112García-Benito, R., González Delgado, R. M., Pérez, E., et al. 2019, A&A, 621,A120Genel, S., Nelson, D., Pillepich, A., et al. 2018, MNRAS, 474, 3976Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175Graham, M. T., Cappellari, M., Li, H., et al. 2018, MNRAS, 477, 4711Hartke, J., Arnaboldi, M., Gerhard, O., et al. 2018, A&A, 616, A123Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356Huang, S., Ho, L. C., Peng, C. Y., Li, Z.-Y., & Barth, A. J. 2013, ApJ, 768, L28Hui, X., Ford, H. C., Freeman, K. C., & Dopita, M. A. . 1995, ApJ, 449, 592Iodice, E., Spavone, M., Capaccioli, M., et al. 2017, ApJ, 839, 21Janowiecki, S., Mihos, J. C., Harding, P., et al. 2010, ApJ, 715, 972Jarrett, T., Chester, T., Cutri, R., Schneider, S., & Huchra, J. 2003, \ aj, 125, 525Jester, S., Schneider, D. P., Richards, G. T., et al. 2005, AJ, 130, 873Johansson, P. H., Naab, T., & Ostriker, J. P. 2012, ApJ, 754, 115Kau ff mann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 54Kluge, M., Neureiter, B., Ri ff eser, A., et al. 2020, ApJS, 247, 43Kormendy, J., Fisher, D. B., Cornell, M. E., & Bender, R. 2009, ApJS, 182, 216Krajnovi´c, D., Bacon, R., Cappellari, M., et al. 2008, MNRAS, 390, 93Krajnovi´c, D., Emsellem, E., Cappellari, M., et al. 2011, MNRAS, 414,2923–2949Lagos, C. d. P., Theuns, T., Stevens, A. R. H., et al. 2017, MNRAS, 464, 3850Law, D. R., Cherinka, B., Yan, R., et al. 2016, AJ, 152, 83Li, H., Mao, S., Cappellari, M., et al. 2018, ApJ, 863, L19Longobardi, A., Arnaboldi, M., Gerhard, O., & Hanuschik, R. 2015a, A&A, 579,A135Longobardi, A., Arnaboldi, M., Gerhard, O., & Mihos, J. C. 2015b, A&A, 579,L3Ma, C.-P., Greene, J. E., McConnell, N., et al. 2014, ApJ, 795, 158Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., & Vauglin, I. 2014, A&A,570, A13Malin, D. F. & Carter, D. 1983, ApJ, 274, 534Mancillas, B., Duc, P.-A., Combes, F., et al. 2019, A&A, 632, A122Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113Méndez, R. H., Ri ff eser, A., Kudritzki, R.-P., et al. 2001, ApJ, 563, 135Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJ, 699, L178Naab, T., Oser, L., Emsellem, E., et al. 2014, MNRAS, 444, 3357Naab, T. & Ostriker, J. P. 2017, ARA&A, 55, 59Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206Nelson, D., Pillepich, A., Genel, S., et al. 2015, Astronomy and Computing, 13,12Nelson, D., Pillepich, A., Springel, V., et al. 2019a, MNRAS, 490, 3234Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624Nelson, D., Springel, V., Pillepich, A., et al. 2019b, Computational Astrophysicsand Cosmology, 6, 2Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ,725, 2312Pastorello, N., Forbes, D. A., Foster, C., et al. 2014, MNRAS, 442, 1003Peng, Y.-j., Lilly, S. J., Kovaˇc, K., et al. 2010, ApJ, 721, 193Penoyre, Z., Moster, B. P., Sijacki, D., & Genel, S. 2017, MNRAS, 468, 3883Pillepich, A., Nelson, D., Hernquist, L., et al. 2018a, MNRAS, 475, 648Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196Pillepich, A., Springel, V., Nelson, D., et al. 2018b, MNRAS, 473, 4077Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14Pulsoni, C., Gerhard, O., Arnaboldi, M., et al. 2018, A&A, 618, A94Robaina, A. R., Bell, E. F., van der Wel, A., et al. 2010, ApJ, 719, 844Roberts, M. S. & Haynes, M. P. 1994, ARA&A, 32, 115Rodriguez-Gomez, V., Pillepich, A., Sales, L. V., et al. 2016, MNRAS, 458,2371–2390Rodriguez-Gomez, V., Snyder, G. F., Lotz, J. M., et al. 2019, MNRAS, 483, 4140Romanowsky, A. J. & Fall, S. M. 2012, ApJS, 203, 17Rosas-Guevara, Y., Bonoli, S., Dotti, M., et al. 2020, MNRAS, 491, 2547Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521Schombert, J. & Smith, A. K. 2012, PASA, 29, 174Schulze, F., Remus, R.-S., Dolag, K., et al. 2020, MNRAS, 493, 3778Schulze, F., Remus, R.-S., Dolag, K., et al. 2018, MNRAS, 480, 4636Scott, N., Graham, A. W., & Schombert, J. 2013, ApJ, 768, 76Spavone, M., Capaccioli, M., Napolitano, N. R., et al. 2017, A&A, 603, A38Springel, V. 2010, MNRAS, 401, 791 Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676Strateva, I., Ivezi´c, Ž., Knapp, G. R., et al. 2001, AJ, 122, 1861Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621,673Trujillo, I., Conselice, C. J., Bundy, K., et al. 2007, MNRAS, 382, 109van de Sande, J., Bland-Hawthorn, J., Fogarty, L. M. R., et al. 2017, ApJ, 835,104van de Sande, J., Lagos, C. D. P., Welker, C., et al. 2019, MNRAS, 484, 869van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28van Dokkum, P. G., Franx, M., Kriek, M., et al. 2008, ApJ, 677, L5van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018Veale, M., Ma, C.-P., Greene, J. E., et al. 2017, MNRAS, 471, 1428Veljanoski, J., Mackey, A. D., Ferguson, A. M. N., et al. 2014, MNRAS, 442,2929Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518Walo-Martín, D., Falcón-Barroso, J., Dalla Vecchia, C., Pérez, I., & Negri, A.2020, MNRAS, 494, 5652Weijmans, A.-M., de Zeeuw, P. T., Emsellem, E., et al. 2014, MNRAS, 444, 3340Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291Wu, X., Gerhard, O., Naab, T., et al. 2014, MNRAS, 438, 2701Zemp, M., Gnedin, O. Y., Gnedin, N. Y., & Kravtsov, A. V. 2011, ApJS, 197, 30Zibetti, S., Gallazzi, A. R., Hirschmann, M., et al. 2020, MNRAS, 491, 3562 Article number, page 22 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations
Appendix A: Color transformation equations C o l o r g - r [ m a g ] g-r = 0.62(g-i)+0.06 TNG50TNG100MANGAA3DMASSIVE C o l o r g - r [ m a g ] g-r = 0.79(B-V)+0.04 J+05TNG50TNG100A3DMASSIVE
Fig. A.1.
Relation between colors for the observed galaxies, as shownin the legend. The simulated galaxies with M ∗ ≥ M (cid:12) are also shownbut without any modeling for dust attenuation, which probably explainswhy they tend to follow a di ff erent color relation than the observations.The solid lines show the linear fit to the observations, while the dashedline shows the color transformation equation from Jester et al. (2005). In Sect. 5.1 we selected a sample of ETGs from TNG50 andTNG100 by using a cut in g − r color and stellar mass, and com-pare with observations. Most of the observed galaxies have g − r colors in the NSA catalog. For a small fraction of the obser-vations, only B − V or g − i colors are available. Jester et al.(2005) published transformation equations between the SDSSand UBVRcIc magnitudes valid for stars and quasars, whichprovide reasonable results for galaxies. We note, however, thatthese transformation equations do not readily apply to the ob-served galaxies. In Fig. A.1 we show the relation between colorsin 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 colorsare from Hyperleda. We find that the observed ETGs follow therelations: g − r = . g − i ) − . g − r = . B − V ) − . . (A.1)Hence we use these derived transformation equations to trans-form the measured B − V or g − i colors in g − r . Appendix B: Accuracy on the measured intrinsicshapes and angular momentum parameter λ e Physical properties measured in simulated galaxies are a ff ectedby resolution e ff ects coming from the discrete particle nature ofthese systems. In a collisionless dark matter-only (DMO) simu-lation the resolution of the gravitational potential depends on thesoftening length and the particle mass resolution (i.e. the numberof particles), which particularly a ff ect measurements in the cen-tral regions of the simulated systems (e.g. Power et al. 2003). Ina full physics simulation instead, like TNG100 and TNG50, theresolution has additional e ff ects on the baryon physics, which re-quire the models to be calibrated on observations (Schaye et al.2015; Pillepich et al. 2018b). In this section we are not concernedwith the latter e ff ects, which also have an impact on the galaxyproperties. Here we aim at quantifying the e ff ects from the reso-lution of the gravitational potential on shape and spin measure-ments at 1 R e and beyond.Chua et al. (2019) analyzed the convergence of intrinsicshape profiles in the Illustris-DMO simulation, with the proce-dure recommended by Zemp et al. (2011), also used in this pa-per and outlined in Sect. 4.1. From their Figure 1 we find that the shape profiles are converged within 0.1 in both p and q al-ready at r ∼ r soft . Since TNG100 has similar particle resolu-tion as Illustris-DMO (and smaller r soft ), we can apply a simi-lar convergence criterion in TNG100. We verified that the twosimulations have comparable particle numbers at r ∼ r soft , i.e.where full convergence is achieved according to the prescriptionof Chua et al. (2019). At smaller radii the full physics TNG100simulation has more particles than the DMO simulation since thebaryons are subjected to dissipation; hence we expect similar, ifnot better, convergence in TNG100.Twice the softening radius in both TNG50 and TNG100corresponds to a radius (cid:46) R e for all the selected galaxies (seeFig. 7, and note that in TNG100 we excluded the galaxies with R e < r soft from the sample). Since R e increases with stellarmass, the more massive systems are better resolved. For exam-ple, in TNG100 the median e ff ective radius at M ∗ = M (cid:12) is1 R e ∼ . r soft , where the resolution e ff ects on both p and q are (cid:46) .
02. Thus the absolute error on the axes ratios p and q mea-sured at r = R e is mass dependent, and it is at most 0.1 for thelow mass systems. At r (cid:38) r soft , which corresponds to r (cid:38) . R e at the low mass end, and to 1 . R e for M ∗ = M (cid:12) , the soft-ening e ff ects are negligible compared to the errors coming fromparticle statistics. For a minimum required number of 1000 parti-cles per ellipsoidal shell we found that these errors are generally (cid:46) .
02 on both p and q . Hence at r (cid:38) r soft the uncertainty onthe axes ratios is of the order of 0.02.Throughout the paper we quantify the intrinsic shapes ofgalaxies also with the triaxiality parameter. Because of its defi-nition (Eq. (3)), the error ∆ T is shape dependent, as well as massdependent as discussed above. In this work we consider reliable T measurements performed at r ≥ r soft , where the uncertaintieson the axes ratios are ∼ .
02, corresponding to ∆ T ∼ . p = . q = . ∆ T grows large for the low mass galaxies, we quantify theintrinsic shapes using only the better determined p and q (e.g. inFig. 15).The angular momentum parameter λ e is evaluated by inte-grating over all the particles within R ≤ R e , hence, by defini-tion, it is derived more reliably than the flattening. This is appar-ent in the convergence study of Lagos et al. (2017) on the EA-GLE simulations, which have particle mass resolution and grav-itational softening length very similar to TNG100. The studyshows that the angular momentum within 1 R e is well convergedalready at stellar masses above M ∗ = . M (cid:12) , i.e. within galax-ies resolved by about 2000 particles. The selected galaxies in theTNG100 sample are resolved with more than 2 × particles,hence the measured λ e is independent of the softening and parti-cle mass resolution.In conclusion we find that both angular momentum andshapes are not (or only marginally in case of the shape) a ff ectedby the resolution of the gravitational potential and by the particlenumber at r (cid:38) 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 excludingan excess of centrally elongated galaxies not present in the ob-served samples. The inconsistency is revealed in the λ e − ε (1 R e )diagram of Fig. 5, in which the centrally elongated galaxies dis-tribute in a region at high ellipticity where there are few observedcounterparts. In App.B we showed that the resolution e ff ects at1 R e on the intrinsic shapes are at most of the order of 0.1: theseare not enough to explain the di ff erences with observations in the Article number, page 23 of 25 & A proofs: manuscript no. TNGkinematics
X / Re Z / R e TNG100-17212
X / Re V / X / Re Z / R e TNG100-17
X / Re V / Fig. C.1.
Stellar mass distribution (left) and V /σ fields (right) of thecentral 3 Re for two examples of the centrally elongated galaxiesprojected edge-on. In the top panel the elongated central component( p (1 R e ) = . q (1 R e ) = .
33, and λ e = .
54) is embedded in adisk structure ( q (5 R e ) = . p (1 R e ) = . q (1 R e ) = .
30, and λ e = .
10) is embedded in aspheroid ( q (5 R e ) = . e [ E dg e O n ] TNG100TNG50 I n t e r m e d i a t e / M a j o r A x i s Fig. C.2. λ e -ellipticity (1 R e ) diagram for the TNG galaxies selected incolor and stellar mass as in Sect. 5.1, and - di ff erently from Fig. 5 - pro-jected edge-on. The symbols are colored according to the intermediateto major axis ratio measured at 1 R e . The black and magenta lines are asin Fig. 5. λ e -ellipticity diagram, nor the extreme values measured for theflattening p (1 R e ).Figure C.1 shows example velocity fields for two of thesecentrally elongated objects (i.e. with p (1 R e ) < .
6) projectededge-on. The colors show the line-of-sight mean velocity di- p ( r ) - T N G pd f Fig. C.3.
Top: Median p ( r ) profiles for the galaxies with p (1 Re ) < . q (5 R e ) in the centrally elongated galaxies of TNG100(left) and TNG50 (right). Dashed and dotted vertical lines show themean and the median of the distributions, respectively. The galaxieswith p (1 Re ) < . ∼ − r ∼ R e ) they are mostly near-oblate ( p ∼ .
9) with medianflattening q ∼ . − .
4, depending on the simulation. vided by the velocity dispersion at the position of each parti-cle. The two galaxies shown have similar flattening in the cen-tral regions p (1 R e ) ∼ . q (1 R e ∼ . p (1 R e ) < . ff erent degrees of edge-on rotationat 1 R e and some of them do not rotate at all (see also bottompanel of Fig. C.1), which is at odds with what is typically seenin barred galaxies (e.g. Falcón-Barroso et al. 2019).The top panels of Fig. C.3 show the median p ( r ) profiles forthe centrally elongated galaxies. The profiles are built by bin-ning 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 fractionof centrally elongated TNG100 or TNG50 galaxies that popu-late the median profile. The elongated regions can extend up toa few Re, typically ∼ R e . Outside these radii, galaxies are near-oblate, with median p (5 R e ) ∼ .
9, and a large scatter on the flat-tening q (5 R e ) (bottom panels in Fig. C.3). The TNG100 galaxieshave rather spheroidal shapes with q (5 R e ) ∼ .
45; the TNG50galaxies tend to be flatter with q (5 R e ) ∼ .
3, and a peak in thedistribution at ∼ . p (1 R e ) is shown by the color of the data points. The centrallyelongated systems occur prominently among the redder simu-lated galaxies, and in specific stellar mass ranges. In TNG100 Article number, page 24 of 25. Pulsoni et al.: The stellar halos of ETGs in the IllustrisTNG simulations C o l o r ( g - r ) [ m a g ] TNG10010.0 10.5 11.0 11.5 12.0log10 Stellar mass / M0.00.20.40.60.8 C o l o r ( g - r ) [ m a g ] TNG500.2 0.4 0.6 0.8p (1Re)
Fig. C.4.
Color (g-r) as a function of the stellar mass. The top panelshows TNG100, the bottom panel TNG50. The dashed lines illustratethe selection of the sample of ETGs with masses 10 . ≤ M ∗ / M (cid:12) ≤ , and red color (Sect. 5.1). The colors indicate the axis ratios p (1 R e ). they are numerous in the interval 10 . M (cid:12) (cid:46) M ∗ (cid:46) . M (cid:12) ,while in TNG50 they have 10 . M (cid:12) (cid:46) M ∗ (cid:46) . M (cid:12) .The fact that these centrally elongated galaxies are prefer-entially produced in a particular mass range, and that most ofthem are old, red systems, points towards these objects beinga class of galaxies that are produced by the simulation but arenot present in nature, probably related to the way the simulatedgalaxies accrete gas, form stars, and quench during a rapid dy-namical evolution in this specific mass range. The di ff erencebetween TNG100 and TNG50 in the stellar mass range wherethe centrally elongated systems are produced could be explainedby the higher star formation e ffi ciency in the higher resolutionsimulation (Pillepich et al. 2018b). This hypothesis is strength-ened by the absence of a similar concentration of centrally elon-gated systems among the red galaxies of intermediate masses inthe Illustris simulation: hence the presence of this population ofgalaxies is due to the galaxy formation model. We note that inTNG100 this mass range approximately coincides with the kneein the stellar mass-halo mass relation (e.g. Behroozi et al. 2013,and reference therein), where accretion could be particularly ef-ficient.As Fig. C.3 showed, these centrally elongated systems areembedded in an near-oblate component with various flattening q , and q tends to be smaller in TNG50. This di ff erence in flat-tening is likely due to the better resolution of TNG50, whichallows to resolve thinner disks (Pillepich et al. 2019, their ap-pendix B and C). From App. B above the error on q due to thespatial resolution is of the order of 0.1 at ∼ R e for the lowestmass systems, so that the measured thickness of a thin disk withradius ∼ R e and q ∼ . λ e values seenin Fig. C.2.This suggests that these centrally elongated galaxies may besystems that were in the process of forming a disk, whose evo-lution was interrupted or derailed by rapid dynamical instability,star formation, and feedback in the simulations, in the partic-ular mass range in which they occur. Fig. 6 in the main textdemonstrates that the distributions of angular momentum andtriaxiality in the surrounding stellar halos are not a ff ected by thecentral elongation p (1 R e ) of the simulated galaxies. Thus theseover-elongated systems can simply be excluded from the sampleselected for our main analysis.) of the simulated galaxies. Thus theseover-elongated systems can simply be excluded from the sampleselected for our main analysis.