© The Authors 2023

1. Introduction

Galaxies acquire angular momentum primordially by gravitational torques induced by the large-scale tidal fields (Peebles 1969). Then, by approximate conservation of the angular momentum, the collapsing star-forming gas builds up the galaxy rotation (Fall 1979; Fall & Efstathiou 1980). The progenitors of present-day early-type galaxies (ETGs) are mostly fast rotating, disk-like objects (e.g., Penoyre et al. 2017; Lagos 2017). Subsequently, during their evolution, ETGs can undergo several physical processes that lead to a gain or loss in specific Angular Momentum (sAM), that is the angular momentum per unit mass of the stellar component j_*.

Since z ∼ 1, massive ETGs with M_* ≳ 10^10.5 M_⊙ tend, on average, to decrease rotational support (Choi & Yi 2017; Walo-Martín et al. 2020) due to the effect of mergers (e.g., Jesseit et al. 2009; Naab et al. 2014; Lagos 2018). Gas-poor mergers increase the mass and the velocity dispersion of the stars, while destroying ordered rotation (e.g., Jesseit et al. 2009) and redistributing AM mostly to the dark matter halos by dynamical friction (e.g., Barnes 1988). The presence of gas, either in the satellite or in the host galaxy, instead leads to a net spin-up of the merger remnant (e.g., Naab et al. 2014; Penoyre et al. 2017). These different formation paths are thought to be at the base of the evolution of the class of slow rotators (SRs) from the fast rotators (FRs) between z = 1 and 0 (Penoyre et al. 2017; Lagos 2017; Schulze et al. 2018).

Another parameter impacting j_* is the epoch of latest gas accretion. Tidal torque theory predicts that gas infalling at later times has higher sAM (Catelan & Theuns 1996). In addition, simulations also predict a change with time in the main accretion mode of galaxies, from filamentary at z ≳ 1 to gas cooling from a hydrostatic halo, which is more efficient in spinning up galaxies (e.g., Garrison-Kimmel et al. 2018). Therefore, galaxies that form most of their stars early and that do not accrete gas at recent times, for example because they are prevented by active galactic nucleus (AGN) feedback, have systematically lower j_* (Lagos 2017; Rodriguez-Gomez et al. 2022).

The coevolution of mass and AM establishes an empirical relation between the sAM of the stellar component, j_*, and the total stellar mass M_*. Fall (1983) found that galaxies distribute according to a power-law $j_{\ast }\propto M_{\ast }^{\alpha }$ with α ∼ 0.6, also called the Fall relation, with the proportionality constant tightly correlating with the bulge-to-total mass fraction or the Hubble type (Cortese et al. 2016; Fall & Romanowsky 2018). This value of α is remarkably close to the expected α = 2/3 for dark matter halos from tidal torque theory (Peebles 1969; Efstathiou & Jones 1979). Therefore, the observed j_* − M_* relation for disk galaxies has been interpreted to result from the approximate conservation of primordial angular momentum of the stellar component, which is similarly torqued as the dark matter halo (Fall 1983; Romanowsky & Fall 2012).

The j_* − M_* relation for spiral galaxies is now well established for a large range of masses (Posti et al. 2018a; Mancera Piña et al. 2021; Di Teodoro et al. 2023, and references therein), facilitated by the fact that for exponential disks the j_*(≤r) converges rapidly beyond 2R_e. ETGs are found to roughly follow a parallel sequence to the spirals, with approximately five times lower j_* in ellipticals (Fall & Romanowsky 2013) and eight times lower j_* in bulge-only galaxies (Fall & Romanowsky 2018). However, in this case, the measurement of the total j_* is challenging and consequently the j_* − M_* relation for these galaxies is far less explored than for the late types.

The case of massive ETGs is of particular interest since their evolution is dominated by mergers which have a strong effect on j_*. However, as mentioned above, the inclusion of ETGs in the j_* − M_* diagram is challenging. Following the pioneering work of Fall (1983), the only work to date that has attempted such a measurement by integrating velocity profiles of ETGs out to large radii is Romanowsky & Fall (2012). As discussed by these authors, the issue resides in the larger Sérsic indices of ETGs compared to disk galaxies, which imply that a larger fraction of their light, and therefore of their total AM, is distributed in the external regions. Hence, accurate measurements of j_* in ETGs require extended kinematic measurements, out to radii that are inaccessible to stellar absorption-line spectroscopy. These are possible only through alternative kinematic tracers such as planetary nebulae (PNe, e.g., Coccato et al. 2009) or globular clusters (GCs, e.g., Schuberth et al. 2010; Strader et al. 2011).

Extended kinematic studies of ETGs (Foster et al. 2016; Pulsoni et al. 2018; Dolfi et al. 2021) have revealed that these galaxies can display a large variety of kinematic behaviors, including embedded disks, strongly rotating outskirts, twisting velocity fields, and multiple rotating components. The presence of these features in both FRs and SRs suggests that the stellar halos of ETGs are often triaxial. The kinematic diversity in ETG stellar halos emphasizes the importance of an approach based on a two-dimensional kinematic mapping to estimate their total j_*, sufficiently extended to trace the variations of rotation amplitudes and of the rotation direction with radius.

Another complication is the estimate of the projection effects on j_*, because of the three-dimensional geometry of ETGs compared to disk-dominated systems and their complex kinematics at large radii. Romanowsky & Fall (2012) tackled this issue using randomly oriented, simple axisymmetric models with cylindrical velocity fields. These assumptions, however, are not necessarily valid for ETGs and might bias the determination of j_*. For example, the velocity fields of regularly rotating FRs are often characterized by a “spider” morphology, with a rotation amplitude decreasing above and below the projected major axis (e.g., Krajnović et al. 2011). Then assuming a cylindrical morphology with a constant rotation amplitude above and below the major axis systematic overestimates j_*.

This paper is part of the extended Planetary Nebula Spectrograph (ePN.S) survey which uses PNe to sample the kinematics of the stellar halos in ETGs (Arnaboldi et al. 2017). PNe are established probes of the stellar population in ETG halos and are good kinematic traces of the bulk of the host-galaxy stars (e.g., Coccato et al. 2009; Cortesi et al. 2013b). To date, several studies have demonstrated that the PNe spatial distribution follows the surface brightness of the host galaxy and that their kinematics is directly comparable to integrated light measurements (Hui et al. 1995; Arnaboldi et al. 1996; Méndez et al. 2001; Coccato et al. 2009; Cortesi et al. 2013a).

The goal is to use the ePN.S kinematic data out to large radii to measure j_* in 32 ETGs using all of the two-dimensional kinematic information. This increases the sample of ETGs of Romanowsky & Fall (2012) by a factor of four for which the sAM was calculated from similarly extended velocity data¹. To do this, we complemented the PN kinematics with absorption line kinematics from integral field spectroscopy (IFS) in the central regions available in the literature or newly extracted from archive Multi-Unit Spectroscopic Explorer (MUSE) cubes. We corrected for projection effects by using simulated galaxies from the IllustrisTNG cosmological simulation as physically motivated models, which have been found to reproduce the j_* − M_* relation and its dependency on morphology well (Di Teodoro et al. 2023; Rodriguez-Gomez et al. 2022).

Previous j_* determinations for ETGs are based on assuming constant mass-to-light ratios with radius. In this paper, we examine this assumption by considering both blue photometric bands (i.e., B, V, or g) and the infrared emission at 3.6 μm, which is a good proxy for the stellar mass (e.g., Forbes et al. 2017). We also explore the effects of IMF gradients on the distribution of ETGs in the j_* − M_* diagram using results from stellar population studies from the literature.

The paper is structured as follows. Section 2 describes the data used in this work (Sects. 2.1–2.4) and the procedures to reconstruct the two-dimensional velocity fields and images (Sects. 2.5 and 2.6). We explain how we derived differential λ(R) profiles in Sect. 3 and aperture projected j_*(R) profiles in Sect. 4. Section 4.4 contains the main observational result of this paper, the dependence of the projected j_* on the stellar masses of ETGs. In Sect. 5 we compare the distribution of projected j_* of simulated TNG100 ETGs with the ePN.S galaxies (Sect. 5.2) and estimate the correction for the limited radial coverage of the PN to estimate the total projected sAM (Sect. 6). In Sect. 7, the TNG100 ETGs are used to evaluate projection effects on j_*. We explain in Sect. 8 how we derived the total j_* − M_* relation for the ePN.S ETGs and in Sect. 9 we estimate the retained fraction f_j of halo sAM as a function of M_*. Finally in Sect. 10 we draw our conclusions.

2. The data

2.1. The ePN.S survey and the ETG sample

The ePN.S survey aims to investigate the kinematics, the dynamics, the angular momentum, and the mass distribution in the halos of ETGs using PNe as kinematic tracers where the surface brightness is too low for absorption-line spectroscopy. The advantage of using PNe over other tracers is that they sample 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. 2013a), out to very large radii (Longobardi et al. 2015; Hartke et al. 2018, 2022).

The ePN.S survey targets a sample of 32 nearby ETGs with absolute magnitudes −22.38 > M_K > −26.02, distances ≤25 Mpc, and covering a wide range of internal parameters (i.e., luminosity, central velocity dispersion, ellipticity, boxy/diskyness, see Fig. 2 in Arnaboldi et al. 2017). Thus the sample includes a representative subset of nearby bright ETGs. Compared to a magnitude-limited sample of ETGs such as, for example, Atlas3D (Cappellari et al. 2011), the ePN.S galaxies are on average more massive and have lower ellipticities (see Fig. 11 in Pulsoni et al. 2018).

The ePN.S ETGs include 24 fast (FRs) and 9 slow rotators (SRs) according to the classification of Emsellem et al. (2011), such that SRs have λ_e ≤ 0.31ε. In this paper, we follow their definition and refer to FRs as the ensemble of fast rotating ellipticals and S0s, but we also refer separately to the fast rotating ellipticals as E-FRs. The ePN.S sample also contains the two major mergers remnants NGC1316 and NGC5128, which are interesting cases for studying angular momentum transport to the galaxy outskirts by dynamical friction (e.g., Barnes & Efstathiou 1987; Barnes 1988; Navarro & White 1994).

The survey is based on PN observations mostly done with the Planetary Nebula Spectrograph (PN.S) at the William Herschel Telescope in La Palma (Douglas et al. 2002), but also includes two catalogs from Counter Dispersed Imaging with FORS2@VLT, and six further catalogs from the literature (references in Table 1 in Pulsoni et al. 2018), for a total of 32 ETGs. The catalogs contain a total of 8636 PNe, making the ePN.S the largest kinematic survey to-date of extra-galactic PNe in the outer halos of ETGs. The data cover 4, 6, and 8 effective radii (R_e) for, respectively, 85%, 41%, and 17% of the sample, and with median extension of 5.6R_e (see Pulsoni et al. 2018).

The procedure of outlier removal and construction of Bona Fide PNe catalogs is described in Pulsoni et al. (2018). We refer to that paper for a detailed kinematic analysis of the ePN.S sample and to the procedure to derive smoothed velocity and velocity dispersion fields.

2.2. Kinematic data in the central regions

In order to achieve a complete two-dimensional map of the ETG kinematics, we combine the PN smoothed mean velocity and velocity dispersion fields with two-dimensional kinematics maps from IFS for the central regions (R ≲ 1 − 2R_e), where the PN detection is incomplete. A good fraction (24/32) of the ePN.S galaxies is part of the Atlas3D survey (Emsellem et al. 2004; Cappellari et al. 2011), which made available the full two-dimensional velocity fields². In addition to these, NGC3115 has available MUSE IFS data from Guérou et al. (2016).

For NGC0584, NGC1316, NGC1399, and NGC4594 we used reduced data cubes from the ESO science archive³. The analysis was carried out using the GIST pipeline (Bittner et al. 2019). As a first step, the pipeline shifts the spectra to rest-frame, and applies any necessary spatial masks to the data. In most cases we used a target signal-to-noise ratio of 200 − 250 to increase the signal-to-noise ratio of the data using Voronoi binning (Cappellari & Copin 2003), and then derived the stellar kinematics using PPxF (Cappellari & Emsellem 2004) and MILES templates (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011). We restrict our analysis to 4800 − 6000 Å to avoid emission lines that affect the calculation of the velocity dispersion (e.g., Barbosa et al. 2018) and we masked any strong emission lines and residual sky lines. The mean velocity and velocity dispersion profiles for these galaxies are shown in Appendix A. The full two-dimensional velocity fields will be made available in a future publication (Ennis et al., in prep.).

For 2/32 ePN.S ETGs, that is NGC3923 and NGC4742, IFS kinematic data are not available. Major and minor axis slit data from Carter et al. (1998) show that NGC3923 has negligible rotation in the center, therefore most of the contribution to its sAM comes from the outskirts. Therefore for NGC3923 we only use the PN kinematics. For NGC4742 we use the major axis slit data from Davies et al. (1983).

Finally, the central regions of NGC5128 have been observed with several MUSE programs. However, due to the presence of the extended dust lane in the center of the galaxy, it was not possible to derive kinematics from the archival MUSE data covering the central ∼6.5 arcmin. On the other hand, this galaxy has the richest PN catalog of the sample, with 1222 PNe distributed from 2 to more than 50 kpc, that is from ∼0.5 out to almost 20R_e. Because of the excellent coverage of the PN data and because most the AM of this galaxy is distributed at large radii, we do not consider additional kinematic data for this galaxy and derive its sAM from the PN kinematics only.

We also use SLUGGS kinemetry results from Foster et al. (2016), that is rotation velocity, kinematic position angle, and velocity dispersion profiles available for 18/32 ePN.S galaxies. These data are based on kinematic maps from observations using slitlets and extend out to typically ∼3R_e. The kinemetric profiles are used to bridge the radial gap between IFS velocity fields and PN data when necessary, as described in Sect. 2.5. Table 1 summarizes the kinematic data used for each galaxy.

2.3. Photometric data

We use the most radially extended photometric data available in the literature, which are typically in the optical B, V, or g bands (e.g., Caon et al. 1990; Kormendy et al. 2009; Iodice et al. 2017; Spavone et al. 2017, see Table 1). For most of the ePN.S galaxies, extended ellipticity and photometric position angle profiles are also available. For NGC3489, NGC4339, NGC4742, and NGC5128, whose ellipticity and position angle profiles are not available in the literature, we assume constant ellipticity and position angle with radius, equal to the average ⟨ε⟩ and ⟨PA_phot⟩ values listed in Table 1 in Pulsoni et al. (2018). In particular, for NGC5128 we use the ellipticity estimated by Rejkuba et al. (2022).

We also consider light profiles extracted from the infrared (IR) 3.6 μm imaging with the Spitzer Space Telescope and published by Forbes et al. (2017). The Spitzer data are available for 20/32 ePN.S ETGs and typically cover radii out to 100–200 arcsec, depending on the galaxy. Hence, we use the Sérsic fits from Forbes et al. (2017) to extrapolate the light profiles to larger radii, out to the radial extent of the ePN.S data. Ellipticity and position angle profiles at 3.6 μm are unfortunately not provided. Therefore, we assume that the shape of the isophotes at these longer wavelengths are the same as in the bluer bands. Table 1 summarizes and collects the references of the photometric data used in this paper.

2.4. Distances, effective radii, and stellar masses

The distances of the ePN.S galaxies are derived from the surface brightness fluctuation method and are listed and referenced in Table 1 of Pulsoni et al. (2018), together with the adopted effective radii R_e. These are circularized R_e, that is the semi-major axis of the ellipse enclosing half of the galaxy light multiplied by the squared-root of the axis-ratio. The stellar masses M_* are estimated from their total absolute K-band magnitudes M_K obtained from the 2MASS extended source catalog (Jarrett et al. 2003), assuming the distances above. We corrected for the over-subtraction of the sky background by the 2MASS data reduction pipeline (Schombert & Smith 2012) using the formula M_Kcorr = 1.07M_K + 1.53 provided by Scott et al. (2013), and converted to stellar masses following the relation log₁₀M_* = 10.39 − 0.46(M_Kcorr + 23) from van de Sande et al. (2019). This is based on the mass-to-light ratio from the stellar population modeling of Cappellari et al. (2013a), converted to a Chabrier (2003) initial mass function (IMF). These M_* values are in good agreement with the stellar masses derived using a mass-to-light ratio dependent on the (B − V)₀ color as in Fall & Romanowsky (2013, where (B − V)₀ is the extinction corrected color from Hyperleda⁴), and with the stellar masses obtained by Forbes et al. (2017) from the Spitzer 3.6 μm luminosity, but converted from a Kroupa (2001) to a Chabrier IMF by a factor 0.92 (Madau & Dickinson 2014) and adjusted to the adopted distances (see Fig. C.1). The mean variation between M_* values obtained from the three methods for the same galaxies is ∼0.05 dex, which we consider as uncertainty on M_*. This does not include the uncertainty on the distances, which is typically on the order of 0.1 mag (e.g., Blakeslee et al. 2009). We subsequently consider the effect of IMF variations among and within galaxies on M_* in Sect. 4.2.2. Galaxy types are taken from the Hyperleda catalog.

2.5. Reconstructing two-dimensional velocity fields

The IFS and ePN.S mean LOS velocity fields are divided in elliptical annuli with constant ellipticity ⟨ε⟩ and position angle ⟨PA_phot⟩. The velocities in each bin are fitted with the model

$$\begin{array}{cc}& V(a,\varphi )=V_0(a)+V_{\mathrm{rot}}(a)cos[\varphi -P{A}_{\mathrm{kin}}(a)]\hfill \\ \hfill & +V_{\mathrm{s}3}(a)sin[3\varphi -3P{A}_{\mathrm{kin}}(a)]\hfill \\ \hfill & +V_{\mathrm{c}3}(a)cos[3\varphi -3P{A}_{\mathrm{kin}}(a)],\hfill \end{array}$$(1)

where a is the semi-major axis of the bin and ϕ the eccentric anomaly

$$\begin{array}{c}\hfill \varphi =arctan[(y_n/(1-⟨\epsilon ⟩)x_n)];\end{array}$$(2)

for more details, readers can also refer to Pulsoni et al. (2018). The coordinates (x_n, y_n) of the velocity fields are rotated such that x_n is aligned with the photometric major axis given by ⟨PA_phot⟩. The rotation velocity V_rot, the kinematic position angle PA_kin⁵, the amplitudes of the third order harmonics V_s3 and V_c3, and the constant V₀ are free parameters. From V₀(a), we estimate the systemic velocity of each field as the weighted sum

$$\begin{array}{c}\hfill V_{\mathrm{sys}}=\frac{{\sum }_{\mathrm{bins}}{V}_0/{(\mathrm{\Delta }{V}_0)}^2}{{\sum }_{\mathrm{bins}}1/{(\mathrm{\Delta }{V}_0)}^2},\end{array}$$(3)

where ΔV₀ are the errors on V₀(a). For ePN.S fields, ΔV₀ are derived from Monte Carlo simulations as described in Pulsoni et al. (2018); for the IFS fields, ΔV₀ are the errors on the fitted parameters.

The two-dimensional velocity fields are reconstructed by combining together the IFS data in the center and the ePN.S smoothed velocity fields at large radii. This is done by subtracting each velocity field by its systemic velocity V_sys, interpolating them onto a regular grid of pixels of coordinates (x_n, y_n) to create two-dimensional mean velocity V_n(x_n, y_n) and velocity dispersion maps σ_n(x_n, y_n), and “bridging” the radial gap between IFS and PN kinematics using SLUGGS data when available, or a smooth interpolation between the IFS and PN fields.

For the mean velocity fields, V_n(x_n, y_n) are estimated in the radial gap between the IFS and the ePN.S maps using Eq. (1), with the parameters V_rot, PA_kin, V_s3, and V_c3 given by a linear interpolation between the fitted values on the IFS data in the center and the ePN.S data at large radii. If available, we use the V_rot and PA_kin profiles from SLUGGS instead of the interpolated values. For the velocity dispersion fields, σ_n(x_n, y_n) are given by the velocity dispersion values from the IFS maps in the center and from ePN.S in the outskirts. In the radial gap, we assume σ constant in elliptical annuli and use a liner interpolation between the σ(a) profiles from the IFS and ePN.S data, or SLUGGS velocity dispersion profiles. We checked that our linear interpolation between kinematic parameters is consistent with long-slit data whenever available at these intermediate radii. Table 1 summarizes the kinematic data and details the procedure used to reconstruct the velocity fields of each galaxy.

Figure 1 illustrates the reconstruction of the two-dimensional velocity field for an example galaxy. In the top panel, the fitted kinematic parameters on IFS and ePN.S data are plotted as a function of a together with kinemetry of SLUGGS data from Foster et al. (2016). The final reconstructed velocity field V_n(x_n, y_n) is shown in the bottom panel: this is given by the IFS map from Atlas3D at the center, by the ePN.S smoothed velocity field at large radii, and by the reconstructed velocity field using SLUGGS data at intermediate radii.

Fig. 1. Kinematic data of NGC1023. Top: fitted kinematic parameters in elliptical bins as a function of major axis distance a for NGC1023 (see Eq. (1)). Results on the Atlas3D velocity field are shown in red, and those on the ePN.S smoothed velocity field are in green. Blue symbols show rotational velocities and kinematic position angles from kinemetry on SLUGGS data (Foster et al. 2016). These agree within errors with the Vrot and PAkin values measured from Atlas3D. ePN.S values also agree within errors except for the first PN velocity point at 70 arcsec, where the PN systematically underestimate the mean velocity (Pulsoni et al. 2018, Sect. 7.2). In the radial range marked with vertical lines, we therefore use the Vrot(a) and PAkin(a) profiles from SLUGSS, and Vs3(a) and Vc3(a) from the linear interpolation of the values fitted on the Atlas3D and ePN.S velocity fields. Bottom: the two-dimensional rotation velocity field of NGC1023 reconstructed from Atlas3D, SLUGGS, and ePN.S., using Eq. (1). The velocity field is rotated so that the photometric major axis (PAphot = 83.3 deg) is horizontal. The final mean velocity field is given by the Atlas3D data in the center (within the black square), the reconstructed two-dimensional field from SLUGGS at intermediate radii (within the black ellipse), and the ePN.S data at large radii. The PN positions are highlighted by open black circles.

2.6. Reconstructing galaxy images

The spatial distribution of light f(x_n, y_n) in the galaxies is reconstructed from the photometric profiles μ(a), ε(a), and PA_phot(a) available in the literature (references in Pulsoni et al. 2018), where μ(a) is the surface brightness profile along the galaxy semi-major axis a, ε(a) the ellipticity profile, and PA_phot(a) is the photometric position angle profile. This is done by assuming that the galaxy isophotes can be approximated by perfect ellipses with ellipticity and position angle given by ε(a) and PA_phot(a), and to which we assign surface brightness μ(a). Thus, for each galaxy, we create an image f_n(x_n, y_n) = 10^−μ_n/2.5 that represents a map of weights for the pixels of coordinates (x_n, y_n).

In case the photometric profiles are not extended enough in radius to cover the extent of the PN velocity field, we extrapolate them as follows. The surface brightness profile μ(a) is fitted with a Sérsic profile (Caon et al. 1993) which is extrapolated to large radii. In these outer regions, the ellipticity ε(a) and position angle PA_phot(a) are assumed to be constant and equal to the outermost measurement available. Figure 2 shows an example of a reconstructed galaxy image from photometric profiles. The extrapolation of the profiles to large radii is shown with solid lines.

Fig. 2. Photometric data of NGC1023. Top: photometric profiles for NGC1023 extrapolated to large radii. Bottom: the galaxy image of NGC1023 reconstructed on a grid of coordinates from the photometric profiles above.

3. Local λ(R) profiles

The angular momentum-like parameter λ, first introduced by Emsellem et al. (2007), is a commonly used proxy for quantifying the projected rotation field. It is defined as

$$\begin{array}{c}\hfill \lambda (a)=\frac{{\sum }_n{f}_n{R}_n|V_n|}{{\sum }_n{f}_n{R}_n\sqrt{{\sigma }_n^2+V_n^2}},\end{array}$$(4)

where $R_n=\sqrt{x_n^2+{\mathit{y}}_n^2}$ is the circular radius of the n-th pixel of coordinates (x_n, y_n); V_n and σ_n are the mean line-of-sight velocity and velocity dispersion; f_n is the flux.

Galaxies are divided into elliptical radial annuli with constant flattening ⟨ε⟩ and major axis position angle ⟨PA_phot⟩ (reported in Table 1 of Pulsoni et al. 2018). The aperture value of λ_e within 1R_e is used to divide ETGs into FRs and SRs (see Sect. 2.1). The local λ(a) integrated within elliptical annuli of mean semi-major axis a instead quantifies the local rotational support. Here, V_n and σ_n are given by the reconstructed velocity fields (Sect. 2.5), and f_n are the fluxes from the reconstructed images (Sect. 2.6). Prior to this work, Coccato et al. (2009) used PNe to derive the λ(a) profiles of a subsample of the ePN.S ETGs. Their approach, which substitutes the weighting by f_n with the weighting by the completeness-corrected number density of PNe, gives very similar results.

Figure 3 shows the local λ(a) profiles of the FRs (top panel) and the SRs (bottom panel). The λ(a) profiles of the FRs show considerable diversity: after peaking at about 1 − 2R_e, they either stay constant or decline more or less steeply with a. This diversity of λ(a) profiles results in a large range of stellar halo rotational support for this class of galaxies. We also note that, in the ePN.S sample, galaxies with the lowest λ(∼6R_e) are mostly E-FRs, while S0s often contain extended stellar disks or rapidly rotating halos with λ(∼6R_e) reaching values of 0.6–0.8 (solid versus dashed lines in Fig. 3). SRs modestly increase their λ with a and often can exceed λ = 0.2. Hence some of the ePN.S SRs host stellar halos with moderate rotation. The two mergers NGC1316 and NGC5128 are highlighted in Fig. 3 with gray lines. Both galaxies display moderate rotational support with λ ∼ 0.3 − 0.4, increasing mildly with radius.

Fig. 3. Local λ profiles of FRs (top panel) and SRs (bottom panel). S0 galaxies are shown with dashed lines, elliptical FRs with solid lines. The two mergers NGC1316 and NGC5128 are highlighted in gray.

A similar variety of λ(a) profiles has been found in simulated galaxies (Pulsoni et al. 2020; Schulze et al. 2020). Also for the TNG100 ETGs Pulsoni et al. (2020) did not find a dependence of stellar halo λ on M_*, except at M_* > 10^11.5 M_⊙ where strongly rotating outskirts with λ > 0.6 are not present. Similarly, we do not detect any dependence of the stellar halo rotational support on stellar mass but the small size of the ePN.S sample at M_* > 10^11.5 M_⊙ does not allow to draw conclusions on trends at high masses. Figure 3 quantifies the variety of kinematic behaviors of stellar halos explored in Pulsoni et al. (2018) and emphasizes the importance of extended kinematic measurements to quantify the angular momentum content in these systems, as the kinematic properties measured in the central 1 − 2R_e do no necessarily extrapolate to large radii.

4. The projected specific angular momentum of the ePN.S galaxies

In this work, we capitalize on the full two-dimensional kinematic data to determine the total projected sAM j_p. This allows us to account for deviations from cylindrical geometry of the velocity fields, nonaxisymmetric features, and misaligned rotation in the total j_p budget. In a coordinate system centered on the galaxy, with the x-axis aligned with the projected photometric major axis given by ⟨PA_phot⟩ and the z-axis aligned with the line-of-sight, we define a projected sAM vector of a galaxy $\overrightarrow{j_p}$ as⁶:

$$\begin{array}{c}\hfill \overrightarrow{j_p}\equiv \frac{\overrightarrow{J_p}}{M_{\ast }}=\frac{\int \overrightarrow{R}\times \widehat{z}V(x,y)\mathrm{\Sigma }(x,y)\mathrm{d}x\mathrm{d}y}{\int \mathrm{\Sigma }(x,y)\mathrm{d}x\mathrm{d}y},\end{array}$$(5)

where $\widehat{z}$ is the unit vector aligned with the line-of-sight, $\overrightarrow{R}$ is the position vector of coordinates (x, y) in the plane of the sky, V(x, y) is the line-of-sight velocity, and Σ(x, y) the surface stellar mass density of the galaxy. The integral in Eq. (5) can be approximated with a sum over surface elements with coordinates (x_n, y_n), that is the pixels in our kinematic maps, with mean line-of-sight velocity, V_n, and stellar mass, Σ_nΔx_nΔy_n:

$$\begin{array}{c}\hfill \overrightarrow{j_p}\simeq \frac{{\sum }_n{\overrightarrow{R}}_n\times \widehat{z}{V}_n{(M_{\ast }/L_b)}_n{f}_{b,n}\mathrm{\Delta }{x}_n\mathrm{\Delta }{y}_n}{{\sum }_n{(M_{\ast }/L_b)}_n{f}_{b,n}\mathrm{\Delta }{x}_n\mathrm{\Delta }{y}_n},\end{array}$$(6)

where we substituted the stellar mass density per unit area with the flux in a photometric band indexed with b, multiplied by the corresponding mass-to-light ratio: Σ_n = (M_*/L_b)_n f_b, n. The modulus of $\overrightarrow{j_p}$ is

$$\begin{array}{c}\hfill j_p=\sqrt{j_{p,x}^2+j_{p,y}^2}.\end{array}$$(7)

The aperture profile j_p(≤a) is derived using Eq. (6) and summing over elliptical apertures of increasing a. The local j_p(a) is instead derived by summing over pixels within elliptical annuli. Unless otherwise stated, we indicate with j_p the total (aperture) projected sAM, integrated out to the outermost available measurement, at a mean 6R_e (median 5.6R_e) with a range [3, 13]R_e.

If the mass-to-light ratio is a constant quantity within galaxies, the (M_*/L_b) term cancels out in Eq. (6) and the pixels (x_n, y_n) are weighted only by their fluxes. In this case, we define a “light-weighted” j_p, light, calculated by weighting the local j_p(x, y) with light profiles in blue optical bands, similar to previous work on sAM in ETGs (Sect. 4.1).

Estimating the gradients of (M_*/L_b) is difficult because it requires constraining how the stellar population properties change out to large radii. To construct a “stellar-mass-weighted” integrated j_p, mass, we use IR fluxes in the Spitzer 3.6 μm band as proxy for the stellar mass distribution (Sect. 4.2.1).

Finally, a variation of the IMF with radius as indicated by several studies of massive ETGs (see Smith 2020, for a review) would have significant impact on the (M_*/L_b) gradients, although with still large uncertainties. In Sect. 4.2.2, we estimate these effects on j_p using the results of Bernardi et al. (2023). We denote with j_{p, mass + IMF} the mass-weighted total j_p obtained with this model for the IMF-gradients.

4.1. Aperture j_p(≤a) profiles: Weighting by light

We start our analysis with the light-weighted case, following a similar approach as (Romanowsky & Fall 2012). We derive light-weighted j_p(≤a) profiles in apertures using Eq. (6) and a constant (M_*/L_b) throughout the galaxies. We choose light profiles in blue optical bands, which are typically the most radially extended as the contribution from the sky background is lower at these wavelengths: for example, the sky in the I band is 3 mag more luminous than in B (Patat 2003). The blue fluxes are used to determine the f_n(x_n, y_n) maps.

The left panels of Fig. 4 show the resulting j_p(≤a) profiles for the ePN.S FRs and SRs, plotted out to the mean major-axis distance of the outermost kinematic bin (see Pulsoni et al. 2018). The top right panel shows the median j_p(≤a) profiles normalized at j_p(≤6R_e). Galaxies are divided among lenticulars, E-FRs, and SRs. The bottom right panel displays instead the median projected angular momentum profiles J_p(≤a) normalized at J_p(≤6R_e). J_p(≤a)/J_p(≤6R_e) is a cumulative quantity that increases monotonically with radius; therefore the J_p(≤a)/J_p(≤6R_e) profiles are of more immediate interpretation than the

$$\begin{array}{c}\hfill \frac{j_p(\le a)}{j_p(\le 6R_e)}=\frac{J_p(\le a)}{J_p(\le 6R_e)}\frac{M_{\ast }(\le 6R_e)}{M_{\ast }(\le a)}\end{array}$$(8)

Fig. 4. Light-weighted aperture sAM profiles for fast and slow rotators (left panels): ellipticals and S0s are distinguished by solid and dashed lines, respectively. The two mergers, NGC1316 and NGC5128, are shown with dotted gray lines. The right panels show the median jp(≤R)/jp(≤6Re) (top) and the median cumulative angular momentum Jp(≤R)/Jp(≤6Re) profiles (bottom) for galaxies divided into elliptical FRs, S0s, and SRs. The two recent major mergers NGC1316 and NGC5128 are highlighted with gray lines.

profiles, which depend on the relative rate at which both J(≤a) and M_*(≤a) increase with radius.

The aperture j(≤a) of most FRs are monotonically increasing functions that tend to plateau beyond 2R_e. By assuming that the aperture values of j_p and J_p measured at 6R_e are good approximations for the galaxy-integrated quantities, we estimate that the E-FRs reach a median 48% of total j_p within 1R_e and a median 90% within 3R_e. The median j_p profile of the S0s increases more slowly with major axis distance a, reaching 30% of j_p(≤6R_e) at 1R_e and 77% at 3R_e, eventually plateauing beyond 4R_e. Hence for both FR classes the j(≤a) profiles are nearly converged within the radial range of the PN data.

The difference between E-FRs and S0s is also visible, although less marked, in the median total J_p(≤a)/J_p(≤6R_e) profiles and can be explained by the different distributions of rotational support for the two classes in Fig. 3: most elliptical FRs have less rotation in their stellar halos compared to S0s, which often rotate rapidly to large radii. The median J_p(≤a)/J_p(≤6R_e) profiles are shallower compared to the sAM profiles: E-FRs and S0s contain 50% of their AM within 1.7 and 2.2 R_e, respectively, and both reach 80% within 4R_e. The fact that the J_p(≤a)/J_p(≤6R_e) seem to gently flatten at large radii suggests that only a small fraction of the total AM is left in the outskirts. Within our sample, we do not observe significant differences between low and high-mass FRs in both their median j_p(≤a)/j_p(≤6R_e) or J_p(≤a)/J_p(≤6R_e).

Slow rotators exhibit markedly different j_p profiles and radial AM distributions. Their more extended mass distributions and the fact that these galaxies rotate faster at larger radii (see Fig. 3) determine much steeper outer profiles than for FRs. The median j_p(≤a)/j_p(≤6R_e) and J_p(≤a)/J_p(≤6R_e) do not converge within the radial coverage of the ePN.S data, as Fig. 4 clearly shows. Thus a nonnegligible fraction of the total angular momentum of these galaxies is distributed at larger radii. The correction from j_p(≤6R_e) to the total, galaxy-integrated j_p is estimated in Sect. 6 using cosmological simulations.

Figure 5 illustrates the different two-dimensional distributions of the local angular momentum in different classes of galaxies and highlights their different dynamical structure at large radii. It shows the median ratio of the j_p, x component (see Eq. (7)) to the local j_p(a) as a function of a, where the x-axis is aligned with the projected major axis of the galaxy. While the component j_p, y is determined by rotation along the projected major axis, nonzero j_p, x signals the presence of kinematic twists or misalignments contributed by minor axis rotation. The median profiles are derived for the ePN.S SRs, FRs divided between low mass and high mass at M_* = 10^10.6 M_⊙, S0s, and E-FRs.

Fig. 5. Median contribution to the local j from misaligned rotation as a function of major axis distance for different groups of galaxies as in the legend. The median profiles have been smoothed to highlight the radial trends. The profiles for the two mergers NGC1316 and NGC5128 are shown separately.

The contribution of j_p, x to the local j_p(a) through Eq. (7) depends on the rotator class as well as on M_* in the ePN.S sample. In low-mass FRs and in S0s, the contribution from “off axis” rotation to j_p(a) is negligible: j_p(a) comes mostly from major-axis rotation. In high-mass FRs, the contribution of j_p, x increases with a, from 10% at 2R_e to 30% at 6R_e. The ratio ∥j_p, x∥/j_p(a) is even larger for the E-FRs, which increases from 20% at 2R_e to 50% at 6R_e. However, even in this case, the contribution from off-axis rotation to the total AM is small, as most of the total AM is dominated by the central 3R_e (see Fig. 4). On the other hand in SRs, both components j_p, x and j_p, y are equally important to the total AM budget. For these systems extended 2-D kinematic information is essential for measuring their total angular momentum.

The two major mergers show a different distribution of j_p with a in Fig. 4, although in both cases the profiles do not seem to plateau within the radial coverage of the PN data. NGC1316 j_p increases steeply within the central 1R_e and most of its j_p is contributed by the major-axis rotation. The j_p and J_p profiles of NGC5128 rise more slowly with radius, meaning that a larger fraction of its AM is distributed at large radii more similarly to SRs. In this galaxy, a large contribution to j_p comes from minor-axis rotation, as shown in Fig. 5.

4.2. Aperture j_p(≤a) profiles: Weighting by stellar mass

Evaluating the stellar mass associated to the light emitted by the galaxies at each radius would require modeling of the star formation history through spectral analysis. The stellar population mix at each radius determines the age and metallicity distribution which fix the (M_*/L_b) ratio, modulo an assessment of IMF which establishes the overall normalization of the (M_*/L_b) (e.g., Poci et al. 2019). For a constant IMF, stellar population gradients in ETGs imply (M_*/L_b) ratios about 20 − 30% larger in the center than at 0.5R_e (e.g., Domínguez Sánchez et al. 2019; Ge et al. 2021). However, recent studies find that ETGs have significant IMF gradients with radius, with enhanced fractions of low-mass stars in the central regions and standard (Kroupa- or Chabrier-like) IMF beyond ∼0.5R_e (Martín-Navarro et al. 2015; Parikh et al. 2018; La Barbera et al. 2019). The presence of these low-mass stars, which contribute only a few percent to the bolometric light of an old stellar population (see, e.g., Fig. 4 from Conroy 2013), can increase the (M_*/L_b) by a factor of more than 2 in the center (Domínguez Sánchez et al. 2019). Although a direct determination of the stellar mass distribution is beyond the scope of this paper, in this section we aim at evaluating the overall effect of stellar population and IMF gradients on j_p.

4.2.1. Constant IMF

We start by assuming a constant IMF. A good proxy for the stellar mass distribution is the IR-light emission, as the fluxes at these wavelengths are dominated by the emission from the old stars. Hence, compared to bluer wavelengths, they are less sensitive to the emission from younger stars with lower (M_*/L_b) ratios. For this investigation, we considered light profiles from Spitzer 3.6 μm imaging published by Forbes et al. (2017). These data are available for 20/32 ePN.S galaxies and cover their central ∼3R_e. At larger radii, the flux is assumed to continue following the extrapolation of the Sérsic profile fitted to the inner data.

Where they overlap, the 3.6 μm profiles are typically steeper than those in the blue bands. Therefore, in the calculation of j_p weighted with IR fluxes, the central regions become more strongly weighted. This results in aperture j_p profiles with similar shapes but lower amplitudes compared to those weighted with blue-band fluxes. Figure 6 shows the comparison for four example galaxies. The total j_p values (measured integrating within the radial coverage of the PN data) are lower by a mean 13%, independent of stellar mass. Therefore, we use this factor to correct the total j_p from light-weighted to mass-weighted in galaxies that lack Spitzer 3.6 μm profiles.

Fig. 6. Light-weighted and mass-weighted jp(≤a) profiles in four example galaxies, two FRs (top) and two SRs (bottom). For each galaxy, we show the profiles weighted by the blue-band fluxes, by the IR fluxes, by the mass-to-light ratio profile given by the colors, and by the IR fluxes corrected for IMF-driven gradients in the mass-to-light ratio. The error bar shows the error on jp derived from the errors on the colors and the dispersion in the mass-to-light ratio versus color relations.

Another approach is to consider spatially resolved mass-to-light-ratio versus color relations, typically calibrated in the central 1 − 2R_e for large samples of galaxies (García-Benito et al. 2019; Ge et al. 2021), and apply them to extended color profiles for the ePN.S galaxies. Unfortunately, color profiles that cover the radial extent of the kinematic data are available only for 12 ePN.S galaxies (Ho et al. 2011; Watkins et al. 2014; Iodice et al. 2016, 2017; Spavone et al. 2017; Ragusa et al. 2022) to which we applied the relations from García-Benito et al. (2019, which assume a Chabrier IMF) to derive the mass-to-light ratio profiles for the corresponding blue-band fluxes. The resulting j_p profiles also in this case typically show similar shapes and lower amplitudes compared with the blue-light-weighted profiles. The total j_p is consistent with the determinations from the IR-light-weighting within the errors on the colors and on the mass-to-light-ratio versus color relations, as shown in Fig. 6.

4.2.2. With IMF gradients

We estimate the effects of IMF-driven gradients in the (M_*/L_b) ratio using the results of Bernardi et al. (2023), who measure M/L and IMF gradients in spatially resolved spectra of ETGs. In agreement with previous studies, they find that IMF-driven M/L gradients are substantial in the central 1R_e, where the IMF typically goes from standard at R > 0.5 − 1R_e to bottom-heavier in the center, and provide mean mass excess profiles α(a) = M_{*, variable IMF}(a)/M_{*, fixed IMF}(a) separately for E-FRs, S0s, and SRs in bins of luminosity in the r-band and central velocity dispersion σ₀. We convert their values based on Kroupa IMF to a Chabrier IMF by dividing them by a fraction 0.92 (any constant factor is unimportant in the computation of j_p but relevant for correcting the stellar masses, see below).

To associate the mean α(a) profiles from Bernardi et al. (2023) to the ePN.S ETGs, we sort our galaxies in similar bins of luminosity and velocity dispersion. We use M_r magnitudes in the AB system from Cappellari et al. (2013b) for the Atlas3D galaxies, from Brown et al. (2014) for NGC4594, from Buzzo et al. (2021) for NGC3115, from Iodice et al. (2016, 2017) for NGC1316 and NGC1399, and from Sandage & Visvanathan (1978) for NGC1344, NGC3923, NGC4742, and NGC5128. The central σ₀ values quoted in Bernardi et al. (2023) are not corrected for the seeing effects and the velocity dispersion profiles are shown only for R > 1 kpc. Hence, for a fairer comparison with these data, we use as σ₀ the velocity dispersion at a = 1 kpc. The five galaxies NGC3377, NGC3489, NGC4339, NGC4742, and NGC7457, that is the five least massive ETGs in the sample, have too low σ₀ to fall in any of Bernardi et al.’s bins. Therefore, for these objects we do not use an IMF correction. For the other systems, we calculate j_p using Eq. (6) and weighting with the IR fluxes multiplied by the IMF-driven mass excess α(a). The increased mass in the central regions causes an overall decrease in amplitude of the j_p profiles which is mildly mass dependent (see Sect. 4.4). The IMF-corrected j_p profiles for four example galaxies are also shown in Fig. 6.

We also apply IMF corrections to the stellar masses M_* using the mass excess profiles α(a) from Bernardi et al. (2023). From these we derive

$$\begin{array}{c}\hfill M_{\ast },varIMF/M_{\ast },ChabrierIMF=\frac{2\pi \int f(a)\alpha (a)a\mathrm{d}a}{2\pi \int f(a)a\mathrm{d}a},\end{array}$$(9)

where f(a) are the fluxes in the blue-bands. Using the IR fluxes in Eq. (9) gives slightly larger M_{*, var IMF}/M_{*, Chabrier IMF} ratios by ∼2 percent. Finally we correct stellar masses derived in Sect. 2.4 assuming a fixed Chabrier IMF, by multiplying the values M_* by the ratio M_{*, var IMF}/M_{*, Chabrier IMF}.

4.3. Errors on the total projected sAM j_p

The uncertainties on the measured total j_p come from the uncertainties on the stellar mass distribution and on the mean velocities. The first are accounted for in the different determinations of j_p using light-weighting or mass-weighting, which give mean differences on the order of 15%.

We estimate the errors on the measured j_p from the uncertainties on the mean velocities, which are largely dominated by the errors on the PN velocity fields. We use Monte Carlo simulations to evaluate the effect of these errors on j_p. We build simulated PN catalogs by adding to their mean velocity a random value from a Gaussian distribution centered at 0 and with dispersion given by the measurement error and the velocity dispersion added in quadrature (see the discussion in Sect. 3 of Pulsoni et al. 2018). The simulated catalogs are used to produce the simulated mean PN velocity fields. Each simulated PN velocity field is then complemented with a Monte Carlo simulation of the IFS kinematics and of the interpolated velocity field in the center, using their corresponding errors on the velocities. The uncertainties on j_p is then the standard deviation of the distribution of j_p values from the Monte Carlo simulations. These are a median ∼10% for the S0s and a median 20% for the E-FRs and SRs, and similar for light-weighted and mass-weighted values. For the galaxies without IR fluxes available, for which we estimated the mass weighted j_p by reducing the light-weighted j_p using a mean shift (see Sect. 4.2.1), the error on the mass-weighted j_p is taken to be the sum in quadrature of the error on the light-weighted j_p and the standard deviation around the applied mean shift.

4.4. The j_p − M_* diagram

Figure 7 shows the relation between the total projected sAM and the stellar mass of the ePN.S galaxies. We show the blue luminosity-weighted j_p, light, the mass-weighted j_p, mass using IR fluxes, and the mass-weighted j_{p, mass + IMF} corrected for IMF-gradients, as calculated in the previous Sects. In the first two cases, j_p is plotted against M_* derived as described in Sect. 2.4 assuming a constant Chabrier IMF. The j_{p, mass + IMF} values are instead shown against M_{*, var IMF}, that is M_* corrected for IMF gradients (see Sect. 4.2.2). All these quantities are listed in Table 2.

Fig. 7. Projected sAM as a function of stellar mass for the ePN.S galaxies. Lenticulars, fast, and slow rotators ellipticals are shown with different colors and symbols as in the legend. We show with progressively decreasing sizes jp, light vs. M*, jp, mass vs. M*, and jp, mass + IMF vs. M*, var IMF. Arrows connect different measures for the same galaxies. Dashed lines show the power-law fit to jp, light vs. M*; solid lines show the fit to jp, mass vs. M*; dotted lines show the fit to jp, mass + IMF vs. M*, var IMF. Purple lines show the fits to the elliptical galaxies (E-FRs+SRs), while green lines show the fits to the lenticulars.

In all three cases, galaxies are found to have j_p values increasing with stellar mass. The lenticulars show systematically higher j_p compared to the ellipticals (E-FRs and SRs) of similar M_*. The dependence on morphology is in agreement with previous work (e.g., Romanowsky & Fall 2012; Fall & Romanowsky 2013, 2018). The nine SRs in Fig. 7 appear to follow the relation traced by the fast rotating ellipticals for large M_*. Even though for the SRs the measured j_p likely underestimates the total projected j_p as their j_p(≤a) profiles are not converged (see Fig. 4), we estimate in Sect. 6 that the integration of j_p out to 15R_e increases the values for the SRs by only 0.18 dex with respect to the FRs.

Assuming a power-law relation of the form

$$\begin{array}{c}\hfill {log}_{10}{j}_p/j_0=A[{log}_{10}{M}_{\ast }/M_{\odot }-11]\end{array}$$(10)

and performing a separate fit to the lenticular and elliptical galaxies, we find a slope close to 0.6 for most cases, as shown in Table 3. Weighting j_p by the stellar mass does not strongly change the slope A compared to the light-weighting case, but systematically decreases the normalization j₀. Only for the S0 galaxies does the correction for IMF-driven gradients in the mass-to-light ratio introduce a tilt in the slope, mostly driven by the four low-mass S0s for which we did not perform a correction for IMF gradients (see Sect. 4.2.1). The fitted parameters A and log j₀ and their errors are collected in Table 3. The errors are given by the sum in quadrature of the errors on the fit and the standard deviation of the distribution of parameters given the errors on the j_p values. These are obtained by fitting Eq. (10) on Monte Carlo simulations of the j_p values extracted from Gaussian distribution centered on j_p and sigma equal to their errors.

We conclude this section noting that the definition of j_p used in this work is different from the similarly called quantity j_p defined in Romanowsky & Fall (2012). As commented by these authors in their appendix, their j_p does not directly quantify the projection of the total angular momentum j_t on the plane of the sky, but represents an intermediate step in the derivation of j_t from long slit observations along the projected semi-major axis, assuming cylindrical rotation (see their Eq. (3)). Therefore, a direct comparison with these previous (systematically higher) estimates of j_p is not straightforward.

5. The projected sAM of the IllustrisTNG ETGs

The new generation of cosmological hydrodynamical simulations are able to produce a rich variety of galaxy morphologies and to resolve the dynamical and stellar population properties of galaxies. The increased resolution and the inclusion of efficient star formation feedback has proved to be fundamental to overcome the “angular momentum catastrophe” (e.g., Navarro et al. 1995; Sommer-Larsen et al. 1999; Governato et al. 2007) and reproduce realistic galaxies with angular momentum content that matches the observations (Genel et al. 2015; Teklu et al. 2015; Zavala et al. 2016; Lagos 2017, 2018).

IllustrisTNG is a suite of cosmological magneto-hydrodynamical simulations that form and evolve galaxies in a ΛCDM universe. They include prescriptions for star formation and evolution, chemical enrichment of the inter-stellar medium, gas cooling and heating, black hole and supernova feedback (Springel et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Naiman et al. 2018; Marinacci et al. 2018). IllustrisTNG generates a population of galaxies with good mixture of morphological galaxy types (Rodriguez-Gomez et al. 2019) and whose j_* − M_* relation, and its dependence on morphology, is in good agreement with observations (Di Teodoro et al. 2023; Rodriguez-Gomez et al. 2022). Furthermore, the TNG ETGs were demonstrated to reproduce many of the kinematical and morphological properties of observed ETGs out into the stellar halos. They show similar changes in rotational support and flattening with radius and similar fractions of ETGs displaying kinematic twists (Pulsoni et al. 2020). In this work, we use the simulated TNG ETGs as models of the ePN.S galaxies to estimate the fraction of the total angular momentum distributed outside the radial coverage of the ePN.S data and to account for projection effects.

5.1. Sample selection and derivation of physical quantities

The IllustrisTNG simulations adopt a universal Chabrier IMF, consistent with our choice for the ePN.S M_*. We select the sample of simulated z = 0 ETGs from the entire TNG100 volume as in Pulsoni et al. (2020), with stellar masses in the enlarged range 10^10.1 ≤ M_*/M_⊙ ≤ 10¹², red colors (g − r)≥0.05log₁₀(M_*/M_⊙)+0.1 mag, and with effective radii R_e ≥ 2r_soft. Stellar masses and effective radii are derived considering all the stellar particles bound to the galaxies. As in Pulsoni et al. (2020), we perform a further selection in the λ_e − ε_e diagram (where λ_e is integrated within an elliptical aperture of semi-major axis R_e and ε_e is the ellipticity at 1R_e) excluding a fraction of elongated galaxies, whose properties are inconsistent with observations. These criteria selected a sample of 1327 galaxies, of which 1047 are FRs and 280 are SRs. The stellar mass function and ellipticity distribution of the TNG100 ETGs are in reasonable agreement with Atlas3D but, compared to the ePN.S sample, they contain a larger fraction of galaxies high ellipticity and lower masses (Pulsoni et al. 2020).

Each simulated ETG is observed along 100 random line-of-sight directions. For each of these projections, we derive projected j_p(≤a) profiles, ellipticity ε(a) profiles, and rotational support Λ(≤a) profiles. The ellipticity profiles are derived using the 2-D inertia tensor as in Pulsoni et al. (2020). The projected j_p(≤a) is determined by applying Eq. (6) to the discrete velocities of the particles, summed within elliptical apertures and weighted by the particle stellar masses. We find that using the stellar particles instead of the smoothed velocities (as derived in Pulsoni et al. 2020) gives very similar results (Fig. D.1), but we adopt the first approach as it is computationally faster. For the same reason, we define a Λ parameter quantifying the galaxy rotational support similar to the λ parameter in Eq. (4). This is given by the ratio of the projected angular momentum per unit mass and radius and the square-root of the line-of-sight kinetic energy per unit mass:

$$\begin{array}{cc}& \mathrm{\Lambda }=\frac{\sqrt{l_{p,x}^2+l_{p,y}^2}}{\sqrt{K}}\hfill \end{array}$$

where

$$\begin{array}{cc}& {\overrightarrow{l}}_p=\frac{{\sum }_n{\overrightarrow{R}}_n\times \widehat{z}{v}_{z,n}{m}_n}{{\sum }_n|{\overrightarrow{R}}_n|m_n}\mathrm{and}K=\frac{{\sum }_n{m}_n{v}_{z,n}^2}{{\sum }_n{m}_n}.\hfill \end{array}$$(11)

Here $\widehat{z}$ is the (random) line-of-sight direction, x and y are the particle coordinates on the projection plane. ${\overrightarrow{R}}_n$ is the two-dimensional position vector of each stellar particle with mass m_n on the (x, y) plane, and v_z, n is the particle line-of-sight velocity with respect to the center of mass of the galaxy. Here we consider Λ(≤a) calculated within elliptical apertures of semi-major axis a. The advantage of using Λ is that the former can be easily derived from the discrete velocities of the particles and does not necessarily require the intermediate step of producing mean velocity and velocity dispersion fields. And, indeed, applying Eq. (11) to the mean velocity V and velocity dispersion fields σ (where in K the $v_{z,n}^2$ is substituted with the V_rms = V² + σ²) delivers nearly identical results (see Fig. D.2). This allows us to also consistently derive Λ for the ePN.S galaxies from their velocity fields (see Fig. D.3).

5.2. Comparison with the ePN.S ETGs

To compare the projected sAM of simulated ETGs with that of the ePN.S galaxies, we first need to match the selection function of the two samples (see Sects. 2.1 and 5.1). Since j_p changes systematically with morphology as well as stellar mass, we match the ePN.S galaxies with simulated galaxies of similar (projected) types and M_*. That is, for each ePN.S galaxy we select an ensemble of analogs among the 1327 × 100 random projections of TNG100 ETGs such that they belong to the same rotator class (FRs or SRs) and have similar stellar mass M_*, projected ellipticity ε, and rotational support Λ: |M_* − M_*, ePN.S|≤Δ_M*, |ε − ε_ePN.S|≤Δ_ε, |Λ − Λ_ePN.S|≤Δ_Λ. The quantities ε and Λ are measured where Λ is maximum. This choice is justified by the fact that, although the TNG100 ETGs qualitatively reproduce the λ(a) and ε(a) profiles of observed ETGs (Pulsoni et al. 2020), the simulated galaxies have a less steep distribution of the angular momentum with radius in the central 1 − 3R_e (see Appendix B). This is also reflected in the systematically larger radii at which the simulated galaxies reach the peak in rotation compared to observed ePN.S ETGs (Pulsoni et al. 2020). Therefore, the most consistent way to compare observations with simulations is to consider ε and Λ where rotation is maximum, instead of considering values measured at fixed multiples of R_e.

Figure 8 shows the distribution of the differences between the j_p of the ePN.S galaxies and the median of the j_p(≤6R_e) of the simulated ePN.S analogs, divided by the sum in quadrature of the scatter of the distribution and the error on the ePN.S j_p. In this figure Δ_M* = 0.1dex M_⊙, Δ_ε = 0.05, Δ_Λ = 0.05, but the results do not depend on the bin size. The three distributions correspond to the light-weighted ePN.S j_p, light, the mass-weighted ePN.S j_p, mass, and the j_{p, mass + IMF} values corrected for IMF-driven gradients in the mass-to-light ratios. In the case of j_p, light, we use the blue light-weighted Λ_light values (Eq. (11)); in the case of j_p, mass and j_{p, mass + IMF}, we use the mass-weighted Λ_mass through the IR fluxes, which are slightly smaller by a median 2% than the Λ_light (see also Fig. D.3). We do not correct Λ for IMF gradients as Λ is only used to match the ePN.S galaxies to the TNG ETGs, and the TNG simulations adopt a constant IMF. The j_p and Λ values of the TNG100 ETGs are instead always mass weighted.

Fig. 8. Comparison between the projected sAM of each ePN.S galaxy jp,ePN.S and the median of the distribution of jp(≤6Re) of the ePN.S analogs among the TNG100 ETGs, divided by the sum in quadrature of 1σ-scatter of distribution of the TNG analogs and the error on the ePN.S jp,ePN.S. The three histograms correspond to the three jp determinations for the ePN.S galaxies: the light-weighted jp, light, the mass-weighted jp, mass, and the jp, mass + IMF corrected for IMF gradients. The jp values for the TNG100 galaxies are mass-weighted. Vertical dashed lines show the median of the three distributions.

The distribution of the differences is roughly centered at −0.5σ-scatter around the median for j_p, light, indicating that the ePN.S galaxies have slightly lower j_p within ∼6R_e than the TNG100 analogs. Conversely, TNG100 ETGs selected to have M_* and j_p close to the ePN.S galaxies are of an “earlier type”, with slightly lower rotational support Λ and ellipticity ε. For j_{p, mass + IMF}, the median of the distribution shifts to lower values (−1.25σ).

Overall, the TNG100 simulation gives a reasonably good description of the angular momentum content of ETGs within 6R_e. Even though the simulated galaxies have a different distribution of j_p in the central regions compared to observations (see Appendix B), they are converged to the ePN.S values at 6R_e.

6. The contribution of the outskirts (a > 6R_e) to the total j_p

In Sect. 5.2, we showed that the TNG100 galaxies have similar values of j_p to the ePN.S galaxies. Assuming that the simulated galaxies have a similar distribution of j_p at large radii as observed ETGs, we can use the TNG100 ETGs to estimate the contribution to the total, galaxy-integrated, j_p from the regions outside the radial coverage of the ePN.S data, typically beyond a > 6R_e.

Figure 9 shows the median difference between the total, galaxy-integrated j_p and j_p(≤6R_e) as a function of stellar mass in the simulated FRs and SRs. Low-mass galaxies have essentially converged to their total j_p already at 6R_e, especially the FRs. For FRs with M_* > 10^10.5 M_⊙, j_p increases by a median 10% beyond 6R_e, with only the most massive systems with M_* ≳ 10^11.7 M_⊙ increasing j_p by 30 − 40%. SRs instead increase j_p considerably in their outskirts, by ∼40% within 6 and 15R_e, and by ∼60% if we consider all the bound particles.

Fig. 9. Median difference between the total, galaxy-integrated jp and jp ≤ 6Re, divided by jp, as a function of stellar mass in the TNG100 ETGs. Galaxies are divided into FRs (blue lines) and SRs (red lines). Solid thick lines are for jp integrated out to 15Re; the shaded areas report the quartiles of the distribution. The dashed lines are for jp integrated out to the outermost bound particle, while thin solid lines are for apertures of 10 and 12.5Re.

However, including all the bound particles in simulated massive ETGs might lead to an over-estimate of the total j_p, as many of these galaxies are centrals in their group halos and the bound particles as identified by the SUBFIND algorithm also include the intra-group or intra-cluster light (ICL). Separating among galaxies and ICL is a nontrivial task and beyond the scope of this work (see Arnaboldi & Gerhard 2022 for a review). For example, Pillepich et al. (2018) use an operative definition of ICL as all the stellar particles beyond a fixed aperture of 100 kpc, which corresponds to ∼5R_e or less in massive ETGs such as NGC4472 or NGC4365 (e.g., Kormendy et al. 2009). At these radii, the j_p(≤a) profiles of the ePN.S SRs still increase steeply with a. Therefore as a compromise, we consider all the particles within a radius of 15R_e, corresponding to ∼300 kpc for these large galaxies to estimate the total j_p. This is approximately the radius at which j_p profiles of the two SRs of Romanowsky & Fall (2012) converge to their total value, given by a power-law extrapolation of their velocity profiles to infinity. The choice of 15R_e as radial limit of integration of the sAM in the simulated ETGs does not affect the determination of j_p in FRs but it is critical for massive SRs, with differences of the order 30% on j_p if we vary the integration limit to the whole extent of the simulated stellar halo or, say, to 10R_e (see Fig. 9).

7. Deprojecting the galaxy angular momentum

Building on the result that the TNG100 ETGs are reasonable models for the observed ePN.S galaxies, we use the simulated ETGs to estimate the projection effects and determine the total true (three-dimensional) sAM j_t from the measured projected sAM j_p. To correct for projection effects is a nontrivial task, as it requires a full knowledge of the three dimensional rotational velocity field at all radii as well as the three dimensional distribution of stellar mass. A simple way to structure the problem is to define a “deprojection factor” C_i such that

$$\begin{array}{c}\hfill j_t=C_i{j}_p,\end{array}$$(12)

where j_p is the galaxy integrated projected sAM. The factor C_i therefore includes all dependencies on inclination, density, and rotation-velocity profiles.

Romanowsky & Fall (2012) used a similar parametrization as Eq. (12). However, in that case, the quantity indicated with j_p does not quantify the projected angular momentum on the plane of the sky and is by definition different from the j_p derived in this work in Sect. 4 and in Eq. (12). Therefore we can not directly apply the C_i factor derived in Romanowsky & Fall (2012) to our data.

In previous work, C_i was calibrated using theoretical models. By assuming that galaxies are transparent spheroids, with axisymmetric density distributions and cylindrical velocity fields, Romanowsky & Fall (2012) find that C_i depends primarily on the inclination relative to the line-of-sight and little (∼10%) on the detailed shape of the rotation profile, while it shows no dependence on the Sérsic index. With these assumptions on the velocity fields, they estimated an inclination-averaged deprojection factor separately for the elliptical and the lenticular galaxies to take into account the inclination bias in the galaxy classification.

However, the assumptions of axisymmetry and cylindrical velocity fields do not hold for all the ePN.S ETGs. In this sample of galaxies we observe a large variety of velocity fields (Pulsoni et al. 2018). Most FRs display velocity fields with rotation strongly concentrated along the major axis, indicating deviations from the cylindrical geometry. In addition all the ePN.S SRs and 40% of the ePN.S FRs show kinematic signatures of triaxial stellar halos. This diversity of kinematics and intrinsic structure complicates the derivation of C_i from models.

A way forward is offered by the simulated TNG100 ETGs, which are found to qualitatively reproduce the observed ePN.S velocity fields and ellipticity profiles, suggesting a similar intrinsic structure (Pulsoni et al. 2020), and j_p values (Sect. 5.2). Therefore, by assuming that the TNG100 ETGs are reasonable models of the real ETGs, we can use them to directly measure the projection factor C_i connecting j_p to j_t.

7.1. The projection factor as a function of inclination

We derive the deprojection factor C_i using 100 random line-of-sight projections of our sample of 1327 simulated ETGs. For each projection, we derive C_i = j_p/j_t using all the particles within 15R_e (see Sect. 6 for the justification). Figure 10 shows the distribution of C_i values for different inclinations. For inclinations close to edge-on, j_p differs from j_t by a factor of ∼2 with small scatter⁷. At lower inclinations, the mean factor C_i increases as a larger fraction of the angular momentum becomes hidden by projection effects. The width of the C_i distribution becomes correspondingly larger. For near-face-on projections, the projected sAM is (on average) only a small fraction of the total. The C_i distribution in fact stretches from values close to 1, which are for minor-axis rotators, to values above 10. In these cases, j_t can only be recovered with a large uncertainty.

Fig. 10. Distribution of Ci for different inclinations measured on 100 random line-of-sight projections of the TNG100 ETGs, as labeled in the legend. The legend also lists the mean, the median, and the standard deviation of each distribution. The gray histogram show the distribution of Ci for all inclinations.

Figure 10 demonstrates the strong dependency of C_i on inclination and therefore the need of an assessment of the galaxy inclination before any attempt of reconstructing j_t from j_p. In Romanowsky & Fall (2012), the problem is bypassed by considering inclination averaged values for C_i and neglecting possible inclination biases. For the TNG100 ETGs, the median C_i over all the inclinations is 2.3 (mean 3.12, see the gray histogram in Fig. 10).

7.2. The projection factor from observed Λ − ε plane

Galaxy inclinations for ETGs are difficult to estimate. However, the projected shape and rotational support of galaxies are also a function of inclination. Therefore, one can expect a dependence of C_i on Λ and ε, which are directly measurable quantities. Figure 11 shows the variation of C_i with these observables for the entire sample of simulated ETGs, not just the ePN.S analogs, with Λ and ε calculated at the location of maximum Λ (see Sect. 5.2). We divide galaxies into SRs, low, and high-mass FRs (M_* > 10^10.6 M_⊙).

Fig. 11. Distribution of Ci in the Λ − ε plane, where Λ and ε are measured at the location of maximum rotational support for the entire sample of simulated ETGs. Galaxies are divided in low-mass FRs (left), high-mass FRs (center), and SRs (right). The first row of panels shows the median μ(Ci) in each bin. The second row shows the root-mean-square scatter around the mean, σ(Ci). The third row shows the ratio between scatter and median. The location of the ePN.S galaxies in the Λlight − ε plane is also shown with open circles. For each ePN.S ETG, the median Ci value of its “analogs”, selected to have similar Λ, ε, and M*, is used to “de-project” jp (see Sect. 8).

Indeed C_i varies smoothly with the observables Λ and ε. At high Λ and ε the median value of C_i is close to 2 and the scatter in the distribution of C_i is very small, of the order 0.2 or less. This is much smaller than the scatter for edge-on ETGs (Fig. 10), because high Λ and ε single out edge-on galaxies with strong intrinsic rotation and flatter shapes. At decreasing Λ and ε, the values of C_i increase in parallel with decreasing inclinations. There is a difference in the trend of C_i at decreasing Λ and ε for low-mass FRs, high-mass FRs, and SRs. The distribution of C_i reaches high median μ(C_i) and scatter values σ(C_i) for low-mass FRs, up to μ(C_i) of the order 10 when they are observed close to face-on. For the high-mass FRs and SRs the increase in both median and scatter is progressively reduced, indicating a systematic change in the dynamical structure of these galaxies with an increasing contribution to j_t from minor axis rotation (see Fig. 5). For the SRs, μ(C_i) does not exceed 4-5.

To conclude, the TNG100 ETGs reveal that three-dimensional sAM can be well predicted given the projected ellipticity (ε) and velocity field (Λ) of an ETG.

8. The total sAM of the ePN.S galaxies

In this section we derive the total sAM j_t for the ePN.S galaxies from the measured projected j_p. To do this, we need to estimate the increase of sAM at large radii that is “missed” by the spatial coverage of the ePN.S survey and determine the projection factor as defined in Eq. (12) to correct for projection effects. As discussed above, we can consider the TNG100 ETGs as good, physically motivated, models of the ePN.S ETGs and use them to evaluate the corrections on the observed j_p values.

For each ePN.S galaxy, whose PN data extend out to a_max, we derive the total sAM using a median correction factor from the TNG ETGs divided in low-mass FRs, high-mass FRs (threshold mass 10^10.6 M_⊙), and SRs, and with similar Λ and ε:

$$\begin{array}{cc}\hfill j_{t,\mathrm{e}PN.S}& \equiv \mathrm{median}{\left[C_i\frac{j_p(\le 15R_e)}{j_p(\le a_{\max })}\right]}_{\mathrm{TNG}100}{j}_{p,\mathrm{e}PN.S}\hfill \\ \hfill & =\mathrm{median}{\left[\frac{j_t(\le 15R_e)}{j_p(\le a_{\max })}\right]}_{\mathrm{TNG}100}{j}_{p,\mathrm{e}PN.S}\hfill \end{array}$$(13)

since C_i ≡ j_t(≤15R_e)/j_p(≤15R_e). The median C_i adopted for each ePN.S galaxy are shown in Fig. 11. The correction from j_p(≤a_max) to the total j_p is quantified in Fig. 9. When applying Eq. (13) to the blue light-weighted j_p we use correction factors based on the blue light-weighted Λ_light values, while for j_p, mass and j_{p, mass + IMF} we use the mass-weighted Λ_mass. The projected sAM, the stellar masses, and the correction factors are listed in Table 2.

8.1. The j_t − M_* diagram

The left panel of Fig. 12 shows the j_t − M_* relation for the ePN.S galaxies. The ePN.S sample displays the well-known increase of sAM with stellar mass, and the correlation with morphology. Elliptical galaxies have significantly lower average j_t than lenticulars of similar stellar masses.

Fig. 12. jt − M* plane for the ePN.S galaxies (left) and comparison with previous works (right). Left panel: total sAM of the ePN.S galaxies as a function of their stellar mass. For clarity we only show the mass-weighted jt, mass values data-points. Solid lines show the power-law fit to jt, mass vs. M*; dashed lines show the fit to the mass-weighted jt, light vs. M*; dotted lines show the fit to the values corrected for IMF gradients jt, mass + IMF − M*, var IMF using the mean mass excess profiles from Bernardi et al. (2023) as discussed in Sect. 4.2.2. Green lines refer to S0 galaxies, purple lines refer to the ellipticals. Right panel: dashed green and purple lines show the power-law fits to the ePN.S lenticulars and elliptical galaxies, respectively, as in the left panel. Dotted black lines show the jt − M* relations for different bulge fractions β as derived by Fall & Romanowsky (2018). Dashed and solid black lines show the fits of Romanowsky & Fall (2012) and Fall & Romanowsky (2013) to their sample of S0s and ellipticals. In both panels, the magenta line shows the results of Di Teodoro et al. (2023) on their samples of spiral galaxies.

The uncertainties shown as error-bars in Fig. 12 are derived from the errors on the projected j_p (see Sect. 4.3) and the width in the distributions of the correction factor j_t(≤15R_e)/j_p(≤a_max) from the simulated ETGs. The widths of these asymmetric distributions around the medians are estimated using their quartiles.

We fit the power law in Eq. (10) to j_t, light − M_*, j_t, mass − M_*, and j_{t, mass + IMF} − M_{*, var IMF}, separating between the S0s and ellipticals; see Sect. 8.2. The results are reported in Table 4.

For the ellipticals the slope in the mass-weighted case is 0.76 ± 0.23, while for the S0 this is 0.55 ± 0.17. Weighting by light or including the IMF gradients does not strongly impact the value of the slope. Only for the S0s, the slope decreases slightly from 0.54 to 0.45 when including IMF gradients.

The normalization of the power law at 10¹¹ M_⊙ is 2.45 ± 0.10 dex for the ellipticals in the mass-weighted case. This is a factor of two lower than for the S0s, and a factor of 9 lower than spiral galaxies (see also Sect. 8.5). The normalization is systematically higher in the light-weighted case by a factor of 1.2, while it decreases by a factor ∼1.4 when accounting for IMF gradients.

For the ellipticals the scatter is relatively larger and the power law therefore less certain (Table 4). The vertical scatter σ(log j_t) is of order 0.3 (0.37 including IMF variations) and comparable to or slightly larger than the combined scatter of 0.31 dex expected from the distribution of dark matter halo spin (0.23 dex, Macciò et al. 2008), the stellar-mass-halo-mass relation (0.15 dex, Moster et al. 2013), and the median error of 40%j_t in the sAM measurements corresponding to 0.15 dex. The small difference between the observed and expected values likely reflects the different formation and evolution paths characterizing these objects. Table 4 reports the orthogonal scatter ${\sigma }_{\perp }=\sigma (log{j}_t)/\sqrt{1+A^2}$ with respect to the power-law fits to facilitate comparison with previous work.

8.2. Fitting and variations

The results quoted in Table 4 derive from a least square fit to the data-points without weighting them by their uncertainties, as also in previous work on ETGs. This is motivated by the fact that galaxies with lower angular momentum or higher stellar masses are also those with larger formal errors, both from the PN velocity fields and the from width of the distribution of the projection factors. Weighting by the errors would lead the fit to be completely driven by a handful of E-FRs with the smallest errors, biassing log j₀. We tried to overcome this by imposing a minimum value for the uncertainties equal to the scatter of the data-points around the power law. In this case, while the fitted parameters for the S0s are similar to the un-weighted case, for the Es the slope decreases by 10% and the normalization increases by 5% because of still higher weight of lower-mass, faster-rotating ellipticals.

Monte Carlo simulations with similar samples of measurements drawn from a power-law relation, with errors depending on M_* and j_t as in the observed sample, and typical intrinsic scatter were made to test different fit methods. The standard deviation in A and log j₀ for the elliptical galaxy sample were found to be typically ∼0.2 and ∼0.1 dex, respectively, consistent with the errors given in Table 4. Biases in the mean were small for A but can be substantial (up to 0.2 dex) for log j₀. The least biased results in the mean were obtained from unweighted fits. In the discussion above we have therefore quoted the results of the unweighted fit as less biased toward high-sAM galaxies.

The fitted parameters listed in Table 4 are for j_t obtained with correction factors integrated out to 15R_e (see Eq. (13) and Sect. 6). Extending the outer boundary to integrate over all the bound particles of the simulated TNG ETGs leads to larger correction factors for the more massive ellipticals, up to a median 0.06 dex for the most massive SRs. This modest increase in the j_t values at the high-mass end determines a steeper slope of the j_t − M_*, from ∼0.75 to ∼0.85 for the elliptical galaxies, but leaves the normalization unchanged. Hence integrating out to the virial radius, and therefore including the contribution of the ICL in the correction factors, does not strongly impact our conclusions.

8.3. The j_t − M_* diagram for FRs and SRs

Based on IFS studies revealing an ubiquity of rotating components in ETGs (Krajnović et al. 2011; van de Sande et al. 2017; Graham et al. 2018), the classification scheme for these galaxies has shifted from a morphological to a kinematic paradigm, distinguishing between FRs and SRs (Emsellem et al. 2007). Figure 12 (left panel) shows that FRs and SRs are segregated in stellar mass but not in j_t. The j_t values of the SRs are more uncertain than the FRs, because of the larger uncertainties on j_p as well as the larger contribution of j_t distributed beyond the radial coverage of our data (see Sect. 6). The ePN.S FRs, on the other hand, show a larger scatter among galaxies at fixed stellar mass, with j_t values differing by more than an order of magnitude at fixed M_* (see also the values of σ_⊥ for the two families in Table 4). This wide range of j_t is unlikely explainable by projection effects, which should be already accounted for by the dependence of the correction factor C_i on the projected Λ and ε. The differences in the measured j_t is likely intrinsic and driven by differences in the bulge fractions among FRs.

A fit of the power law in Eq. (10) to the FRs and SRs separately yields very similar relations, with nearly identical normalization and slightly steeper slope for the FRs, despite the difference in the stellar mass range probed by the two samples (see Table 4). Our results suggest that there is no fundamental difference in the sAM content of FRs and SRs, but only in the way this is distributed with radius.

8.4. Comparison with previous work

The right panel of Fig. 12 compares our results with the previous determination of the j_t − M_* relation for ETGs in Romanowsky & Fall (2012) and in Fall & Romanowsky (2013), and for galaxies with different bulge fractions as derived in Fall & Romanowsky (2018). These relations are obtained from j_t values weighted with blue photometric bands, as the j_t, light in this work. The stellar masses in Romanowsky & Fall (2012) are derived from K-band photometry assuming a constant mass-to-light ratio M_*/L_K = 1. Fall & Romanowsky (2013, 2018) revise j_t − M_* relation using a color-dependent M_*/L_K, which returns M_* comparable to ours (see Fig. C.1), and which moves the position of the ellipticals slightly upward by 0.05 dex compared to Romanowsky & Fall (2012). These determinations are closely comparable with the j_t, light − M_* relations derived in this work.

Our j_t, light − M_* relation for elliptical galaxies has a systematically lower normalization compared to Romanowsky & Fall (2012) and Fall & Romanowsky (2013) by a factor 1.62 and 1.67 respectively. It is closer but still below the j_t − M_* determination for pure bulges in Fall & Romanowsky (2018, for the limit β = 1, see Fig. 12). The j_t, light − M_* relation for the S0s determined in the present work also has lower normalization and shallower slope than that of Romanowsky & Fall (2012). These differences become more marked with our mass-weighted j_t measurements.

A galaxy-by-galaxy comparison between the j_t, light values derived here and in Romanowsky & Fall (2012) is shown in Fig. C.2 and discussed in Appendix C: the j_t values are systematically higher in Romanowsky & Fall (2012) for the subset of galaxies in common. This is at least partially explained by their assumption of cylindrical velocity fields which would systematically overestimates j_t. An exception is NGC5128, for which the strong minor-axis rotation gives a nonnegligible contribution to j_t that is not accounted for in the major-axis based measure of Romanowsky & Fall (2012). By taking into account the full two-dimensional kinematic information out to large radii and the effects of mass-weighting versus light-weighting, we find that elliptical galaxies have significantly lower j_t than previously estimated.

8.5. Comparison with spiral galaxies

The j_t − M_* relation for spiral galaxies is well established, with different studies returning consistent values for both slope (∼0.6) and normalization (∼3.4 dex at 10¹¹ M_⊙, e.g., Fall 1983; Fall & Romanowsky 2013; Obreschkow & Glazebrook 2014; Posti et al. 2018a; Hardwick et al. 2022). This is because of the rapid convergence of the j_p(≤a) profiles within a few R_e and the better constrained inclination angles compared to ETGs. The pink line in Fig. 12 shows the j_t − M_* relation for spiral galaxies from Di Teodoro et al. (2023), which is based on IR light-profiles and therefore directly comparable with our j_t, mass − M_* relations.

The comparison with the measurements for the ePN.S sample shows that, at fixed stellar mass, earlier morphological types have on average lower j_t, confirming the result of Romanowsky & Fall (2012). Even in the case of SRs, which have their largest sAM in the outskirts, the estimated j_t do not converge to the values for spiral galaxies with similar stellar masses, even when including the contribution from ICL (see previous Section). We find that elliptical galaxies have 9 times (and up to 13 times when including IMF effects) less j_t than spiral galaxies of the same stellar mass.

Another difference between elliptical and spiral galaxies is the larger scatter of the ellipticals in the j_t − M_* diagram. The orthogonal scatter σ_⊥ in log-log space around their best fitting power law is ∼0.24 (and increases up to 0.3 when including IMF variations). This is larger than the scatter typically measured for spirals, which ranges from 0.17–0.22, depending on the sample (Posti et al. 2018a; Hardwick et al. 2022; Di Teodoro et al. 2021, 2023).

8.6. Comparison with simulations

In this section, we compare the distribution of the ePN.S galaxies in the j_t − M_* plane with that of simulated ETGs from TNG100 and from previous works using cosmological hydrodynamical simulations. For consistency, we consider only the mass-weighted measurements with constant Chabrier IMF. Figure 13 shows that overall the ETG samples from current cosmological hydrodynamical simulations give a good representation of the observed distribution of the ePN.S ETGs in the j_t − M_* plane.

Fig. 13. Comparison with simulations. The distribution of ePN.S galaxies on the jt − M* plane (mass-weighted measurements with Chabrier IMF) is compared to the fit on their simulated analogs from TNG100 and to results from other works investigating ETGs in cosmological-hydrodynamical simulations (see text).

We first discuss the comparison of the ePN.S galaxies with their analogs among the TNG100 ETGs. The TNG100 analogs were selected in Sect. 5.2 to have similar stellar masses, projected ellipticity, and rotational support as the ePN.S galaxies. For each ePN.S galaxy, we derive a median “simulated” j_t from its analogs, which is used to obtain a corresponding j_t − M_* relation for simulated ellipticals and lenticulars.

The TNG analog ETGs have overall similar properties as the ePN.S ETGs in the j_t − M_* plane, with the analog S0s having systematically larger j_t than the analog ellipticals and j_t increasing with M_* similarly as observed. Both families of TNG100 analogs have on average larger j_t compared to the ePN.S galaxies, by a factor 1.8 for the ellipticals and 1.9 for the lenticulars, however within ∼1 standard deviation of the distribution as already quantified in Sect. 5.2 and Fig. 8. This is likely related to the different radial distributions of j_p between TNG and observed galaxies, see Fig. B.1. For this reason we used ratios but not absolute j_p values for the TNG galaxies in Eq. (13) for estimating the total j_t of the ePN.S galaxies from their measured j_p. It is unclear whether this offset is entirely due to differences in the physical properties between real and simulated galaxies, or whether it also reflects some bias in the sample selection not accounted for by our definition of analogs. The enlarged but still limited number statistics offered by the ePN.S sample does not allow us to investigate this conclusively.

With a different selection of TNG100 ETGs based either on the kinetic energy associated with ordered rotation, or on a neural network trained morphological classification, Rodriguez-Gomez et al. (2022) found median sAM values tracing roughly those of the ePN.S lenticulars (dotted gray gray lines in Fig. 13), likely indicating a larger average contribution of disk-like components in their sample ETGs compared to the ePN.S ETGs (see Sects. 2.1 and 5.1).

Blue lines in Fig. 13 show the median stellar sAM measured within 5R_e in galaxies from EAGLE by Lagos (2017), selected to have the reddest colors, (u − r) > 2.2, or alternatively the oldest ages in their simulated samples (mass-weighted stellar ages older than 9.5 Gyr). Their measurements follow well the location of the ePN.S ETGs. The smaller aperture used to integrate their sAM should not contribute to important systematic differences at least up to M_* < 10^11.2 M_⊙ where the large majority of ePN.S FRs has essentially converged sAM profiles already at 5R_e (see Fig. 4). The systematic offset with respect to the mean TNG100 values could at least partially be attributed to differences in the selection criteria, as the EAGLE galaxies are on average redder and older, and therefore likely to have lower sAM (Lagos 2017).

Finally, Teklu et al. (2015) analyzed spheroidal galaxies from the Magneticum simulation selected based on their circularity parameter. Their sAM values are in good agreement with the ePN.S ellipticals (mean values shown by the orange line in Fig. 13).

This comparison shows that the current generation of cosmological simulations reproduce well the properties of observed ETGs on the j_t − M_* plane. The vertical offset among lines in Fig. 13 is most likely due to different disk-to-spheroid ratios in the different samples, given the sensitivity of j_t to galaxy type, from sample selections and possibly also different galaxy formation models.

9. The retained fraction of angular momentum

In the previous section, we have derived the total sAM j_t for the ePN.S galaxies which we now consider as our estimate of their total j_*. The empirical j_* − M_* relation of the stellar component in galaxies is often interpreted in conjunction with the properties of their dark matter halos, by defining:

$$\begin{array}{c}\hfill f_j\equiv j_{\ast }/j_h\mathrm{and}{f}_M\equiv M_{\ast }/M_h,\end{array}$$(14)

where j_h and M_h are sAM and mass of the dark halo. f_M is the stellar-mass-halo-mass (SMHM) relation, also referred to as the star formation efficiency, and f_j is the retained fraction of angular momentum.

The connection between the two components is based on the theoretical framework in which the stellar component “inherits” a fraction f_j of the primordial j_h exerted by tidal torques to the collapsing dark matter halos. The tidal torque theory predicts that $j_h\propto M_h^{2/3}$(Peebles 1969; Efstathiou & Jones 1979) which, given Eq. (14), translates into a relation for j_* as a function of M_*:

$$\begin{array}{c}\hfill j_{\ast }\propto f_j{f}_M^{-2/3}{M}_{\ast }^{2/3}.\end{array}$$(15)

Since disk galaxies are observed to follow closely $j_{*}\propto M_{*}^{2/3}$, Eq. (15) leads to the result that the product $f_j{f}_M^{-2/3}~\text{const}$, meaning that galaxies that are “efficient at forming stars are also efficient at retaining angular momentum” (Posti et al. 2018b). Previous results for disks report a retention factor between 0.8 and 1 (Fall & Romanowsky 2013; Posti et al. 2018a; Di Teodoro et al. 2021, 2023; Romeo et al. 2023). For elliptical galaxies, this is more uncertain. In this section, we revisit the estimate of the retention fraction f_j for ETGs given the j_* − M_* derived in this work.

We follow the derivation from Posti et al. (2018b), which assumes a Hubble constant H₀ = 67.7 km s⁻¹ Mpc⁻¹ (Planck Collaboration XIII 2016) and adopts Navarro et al. (1996) dark matter halos. This gives that the stellar j_* can be written as

$$\begin{array}{c}\hfill j_{\ast }=\frac{77.4}{\sqrt{F_E(c)}}\left(\frac{\lambda }{0.035}\right)f_j{f}_M^{-2/3}{\left(\frac{M_{\ast }}{{10}^{10}{M}_{\odot }}\right)}^{2/3}\mathrm{kpc}\mathrm{km}{\mathrm{s}}^{-1}\end{array}$$(16)

where λ is the halo spin parameter (Peebles 1969) and F_E(c) is a dimensionless factor that depends on the concentration parameter c (Mo et al. 1998).

There are several models of the SMHM relation available in the literature. The current consensus is that f_M has a characteristic “bell” shape, so that it peaks around M_* ∼ 10^10.5 M_⊙ and decreases at lower and higher masses by ∼2 order of magnitudes. The decrease at high masses is usually interpreted as signature of AGN feedback reducing the star formation efficiency (Silk & Rees 1998; Croton et al. 2006; Wechsler & Tinker 2018). Here, we consider the model of f_M derived by Moster et al. (2013) from abundance matching, which is a good approximation for ETGs (Posti & Fall 2021), although other models are qualitatively similar (see, e.g., Behroozi et al. 2019, and references therein).

Figure 14 shows the median retention factor f_j and its scatter as a function of the M_* for ellipticals and S0s separately. This is derived as in Posti et al. (2018b) by Monte-Carlo simulating 10⁴ galaxies with a uniform distribution in M_* between 10¹⁰ and 10¹² M_⊙. Based on M_*, we assign to each galaxy a halo mass M_h by sampling the SMHM relation of Moster et al. (2013) with its scatter (0.15 dex), a halo concentration c by sampling its log-normal distribution centered on the power-law dependence on M_h found by Dutton & Macciò (2014), with scatter of 0.11 dex, and a spin parameter λ sampled from a log-normal distribution with mean value λ = 0.035 and scatter of 0.23 dex (Macciò et al. 2008) independent of M_h. Finally, we assign a j_* value assuming a Gaussian distribution of log j_* values with mean value given by the mass-weighted power laws fitted on the ellipticals and the S0s (Table 4) at M_*, and sigma equal to the median uncertainty on the log j_* measurements. The simulated M_*, M_h, c, λ, and j_* are inserted in Eq. (14) solved for f_j, to obtain the distribution of f_j as a function of M_*.

Fig. 14. Ratio of the stellar to halo sAM as a function of stellar mass fj(M*) for ETGs assuming the SMHM relation of Moster et al. (2013). Ellipticals and S0s are shown with thick purple and green lines, respectively. Colored bands show the scatter around the median given by the Monte-Carlo simulations described in the text. The thinner purple lines show the median fj(M*) profiles using fM from Dutton et al. (2010) for ETGs, Behroozi et al. (2013, 2019) for quiescent galaxies as labeled. For comparison, we also show the results of Posti et al. (2019) for spiral galaxies with a magenta line.

At the peak of the SMHM, that is for M_* ∼ 10¹⁰ − 10^10.5 M_⊙, f_j ∼ 0.25 for the S0s and ∼0.08 for the ellipticals. At higher masses, as the star formation efficiency decreases, f_j is reduced by 1.5 order of magnitudes at 10¹² M_⊙. It is worth noting that at high stellar masses the SMHM relations are more poorly constrained and the extrapolations to M_∗ = 10¹² M_⊙ of different models give results that can differ by more than an order of magnitude in f_M because of the uncertainties in the stellar mass function at the high-mass end (see, i.e., Fig. 34 in Behroozi et al. 2019, and their discussion). To illustrate this uncertainty, Fig. 14 also shows the f_j(M_*) relations for the light-weighted j_t, light of the elliptical galaxies using the f_M models from Dutton et al. (2010) and Behroozi et al. (2013, 2019) for quiescent galaxies. This means that the retention factor f_j of ETGs sharply decreases with M_* in tandem with the star formation efficiency; the magnitude of this decrease, however, is highly uncertain. The differences between f_j(M_*) at 10¹⁰ < M_∗ < 10^10.5 M_⊙ instead reflect the differences in the peak of the star formation efficiency between models.

Romanowsky & Fall (2012) also find values f_j ∼ 0.1 for the ellipticals, assuming the SMHM from Dutton et al. (2010). In this case, the higher sAM in their data is compensated by the lower star formation efficiency assumed. Assuming Dutton et al. (2010), we find f_j ∼ 0.04 for low-mass ellipticals and f_j ∼ 0.01 at the high mass end.

It is interesting to compare the f_j(M_*) relation of the ETGs with that of spiral galaxies. The recent work of Posti et al. (2019), Posti & Fall (2021), and Di Teodoro et al. (2023) has shown that, contrary to ETGs, the star formation efficiency of massive spirals increases monotonically with M_* until reaching f_M ∼ f_b above M_∗ = 10¹¹ M_⊙. This implies that the angular momentum retention factor of spirals depends weakly on M_* and it is close to f_j ∼ 0.8 − 1. Hence, the stellar component of these galaxies is consistent with the approximate conservation of the primordial sAM, in stark contrast with the behavior of ETGs which show declining f_j with M_* and f_j ≪ 1.

The difference between ETGs and spirals can be explained by the different formation pathways that characterize these systems. On one side, massive spirals live in low-mass dark matter halos (e.g., Posti & Fall 2021). They follow a gradual evolution, with star formation activity sustained by smooth gas accretion and few minor mergers, and regulated by stellar winds, which promote the formation of high j_* galaxies (e.g., Genel et al. 2015; Zavala et al. 2016; Rodriguez-Gomez et al. 2022). On the other, massive ellipticals inhabit more massive structures. They form most of their stars from low angular momentum gas at early times while their strong AGN feedback inhibits the accretion of high j gas at later times, determining the formation of low j_* galaxies (Genel et al. 2015; Zavala et al. 2016; Lagos 2017; Rodriguez-Gomez et al. 2022). In addition, ETGs undergo a further loss of j_* through their merger activity (Lagos 2018), which redistribute a fraction of their j_* to the dark matter halo by dynamical friction (e.g., Barnes 1988). Since more massive ETGs host more massive supermassive black holes (Kormendy & Ho 2013) and undergo a larger number of massive gas-poor mergers than low-mass ETGs (Rodriguez-Gomez et al. 2016, 2017; Pulsoni et al. 2021), they are likely to “retain” a lower fraction of their primordial j_* compared to lower-mass ETGs, consistent with the declining trend of f_j with M_*.

In addition, the fact that massive ETGs often live in massive environments further complicates the interpretation of our results and the connection between the properties of the stellar and the dark matter component. While the SMHM relation associates halo masses of M_h ∼ 10¹⁵ M_⊙ (approximately the mass of the Coma Cluster) to galaxies of M_∗ ∼ 10^11.5 M_⊙, the j_* measured in this work does certainly not integrate over the whole stellar component associated to the (cluster) dark matter halo. For example, our j_* measurement does not include the contribution from satellite galaxies in the halo, which might contain a large fraction of the total j_*. Indeed, the fraction of the cluster stellar mass locked in the satellite galaxies can range from 20 to 90% (Montes 2022). Therefore, even though the measured j_* converge within a number of effective radii to the total j_* within the galaxies, it may not be representative of the total sAM of the whole stellar component contained in the halo.

10. Conclusions

For this study, we measured local λ(a) and projected sAM profiles in apertures, j_p(≤a), for a sample of 32 nearby ETGs from the ePN.S survey. These galaxies have stellar masses M_* in the range 10^10.1 − 10^11.9 M_⊙ and two-dimensional kinematic data out to a mean radius of six effective radii (6R_e). We used PNe as kinematic tracers of the stellar halos and IFS data for the central regions. This has increased the number of ETGs for which the sAM can be estimated out to large radii by a factor of four over previous work. In addition, the full two-dimensional kinematic information has allowed us to do so without making assumptions about the geometry of the velocity fields and to include the contribution from minor axis rotation.

In order to reconstruct the total sAM j_* from the measured projected j_p, we used simulated ETGs from the IllustrisTNG simulation TNG100. The TNG100 ETGs have velocity fields and projected ellipticity profiles qualitatively similar to the observed ETGs (Pulsoni et al. 2020). Simulated “ePN.S analogs” which were chosen to have a similar ellipticity, rotational support, and stellar mass as the ePN.S ETGs, also have a similar distribution of j_p(≲6R_e) as the ePN.S ETGs (see Fig. 8). Our main results are:

(i) The kinematic diversity of the ETG stellar halos already explored in Pulsoni et al. (2018) is also visible, and quantified, through their large variety of λ(a) and j_p(≤a) profiles (Figs. 3 and 4). FRs can have rapid rotation at large radii or more or less steeply decreasing λ(a). Elliptical FRs have reduced halo rotational support compared to S0s, which often contain extended disks or rapidly rotating stellar halos. SRs generally have increasing λ with radius but the rotational support at large radii is still modest. The contribution to j_p from off-major-axis rotation increases with radius and stellar mass (Fig. 5), with the largest values (≳60%) being in the halos of the SRs.

(ii) The large radial extent of the PN data show that ETGs do not contain large mass fractions of high j_* in their stellar halos. For the FRs, aperture j_p(≤a) profiles are approximately converged within the radial range of ePN.S data (typically 6R_e), but this is not the case for SRs (Fig. 4). We therefore used the simulated ETGs for estimating the total j_* for the ePN.S galaxy sample from their measured j_p(≲6R_e). In the TNG100 FRs, j_p(≲6R_e) has likewise essentially converged; however, for the SRs j_p increases by up to 40% beyond 6 Re (Fig. 9). For the projection factor C_i = j_t/j_p, relating the (galaxy-integrated) projected j_p to the three-dimensional sAM j_t, we also used the TNG100 ETGs. We find that C_i can be inferred from the observed projected ellipticity and rotational support with little scatter except for nearly face-on FRs (Fig. 11). The reconstructed j_t values are our estimates of the total stellar sAM j_* of the ePN.S ETGs.

(iii) We find that j_p and j_* increase with M_* and confirm the well-known dependence of j_* on morphology (Figs. 7 and 12): at fixed M_*, ellipticals have lower sAM than S0s, which have lower sAM than spirals. However, we do not detect a significant difference between FRs and SRs: SRs lie at the high-M_*, high-j_* end of the increasing trend with mass.

(iv) The total mass-weighted j_p (and j_*) for ETGs are 0.07 dex lower than light-weighted sAM, due to color gradients. Correcting for IMF gradients using results from ETG stellar population studies in the literature further reduces values by 0.2 dex relative to the light-weighted case, due to the additional stellar mass in the central regions. This mostly reduces the normalization of the j_* − M_* relation.

(v) A power-law fit to the distribution of ellipticals and S0s in the mass-weighted j_* − M_* diagram gives a slope of 0.55 ± 0.17 for the S0s and 0.76 ± 0.23 for the ellipticals. Remarkably, similar to the dark matter halos, the most massive ETGs have the largest sAM. However, the normalization is systematically lower: the ePN.S ellipticals have a factor two lower sAM than S0s and nine times lower than spirals. The gap between ellipticals and spirals increases up to a factor of 13 in the IMF-corrected case. The scatter increases from spirals and S0s to ellipticals, where it reaches 0.3 dex in the mass-weighted and 0.37dex in the IMF-corrected case. This likely reflects the varied formation histories of ETGs as well as the scatter in the halo spin and in the stellar-mass-halo-mass relation.

(vi) In the theoretical framework in which stars inherit a fraction f_j of the primordial sAM of their host halo, j_h, we estimate that, at stellar masses in the range 10¹⁰ − 10^10.5 M_⊙, S0s retain f_j ∼ 25% of j_h, while ellipticals only retain f_j < 10%. For higher masses, this fraction strongly decreases with the stellar mass, by ∼1.5 orders of magnitudes at M_∗ ∼ 10¹² M_⊙ (Fig. 14). This strongly contrasts with the behavior of spirals which have f_j ≳ 0.8 weakly dependent on M_∗ > 10¹⁰ M_⊙ (Posti et al. 2019; Di Teodoro et al. 2023).

Our results clearly show that the stellar sAM of ETGs increase with stellar mass roughly as $M_{*}^{2/3}$ but they do not converge to those of the spirals in their stellar halos, confirming previous findings. However, the ratio of j_* between ellipticals and spirals has increased with our analysis of the ePN.S sample from a factor of ∼4.5 to a factor of nine. It is likely that, during their evolution, ETGs have lost a significant fraction of their sAM to the dark matter component through dynamical friction. In addition, they may never have accreted some of the high j baryons which could be retained in the baryonic component that did not condense into stars, such as the hot circumgalactic medium. Also, in the most massive systems, the high-j_* stellar component associated with the dark matter halo is the satellite galaxies that have not yet merged with the central ETG. Therefore, while the decreasing trend of f_j with M_* is consistent with the increased gas-poor merger rate and AGN activity at higher masses, the actual magnitude of f_j at M_∗ > 10¹¹ M_⊙ is uncertain, due to the uncertainties in the star formation efficiency and in how representative the galaxy j_* is of the sAM of the all the stars in the host halo.

Acknowledgments

We thank the anonymous referee for their comments which improved the clarity of our manuscript. We thank C. Spiniello and T. Parikh for helpful discussions on IMF variations in ETGs. We are also grateful to R. Ragusa and M. Spavone for kindly providing us their photometric data in tabular form. C.P. is grateful to F. Hofmann for his support.

References

Appendix A: Kinematic profiles for the new MUSE velocity fields

Fig. A.1. Kinematics profiles for the new MUSE velocity fields. The profiles from MUSE are shown with red symbols. Blue symbols show major axis slit data available in the literature (Kormendy & Illingworth 1982; Saglia et al. 2000); green circles show the PN kinematic profiles.

Fig. A.2. As in Fig. A.1, but for NGC1399, and we also show the velocity dispersion profile from Vaughan et al. (2018) based on different MUSE observations with open symbols.

In this section, we show the kinematic profiles extracted from the new MUSE velocity fields. The full two-dimensional kinematic data will be made available in the future (Ennis et al. in preparation). For each galaxy, we show the major axis velocity profiles

$$\begin{array}{c}\hfill V(a,\varphi =P{A}_{\mathrm{kin}})=V_{\mathrm{rot}}(a)+V_{\mathrm{c}3}(a),\end{array}$$(A.1)

that is Eq. (1) for ϕ = PA_kin, the azimuthally averaged velocity dispersion profiles σ(a) in elliptical annuli, and the misalignment between the kinematic position angle profile PA_kin(a) and the mean photometric position angle ⟨PA_phot⟩. We compare these kinematic profiles with slit data available in the literature and complement them with PN data at large radii. The different data sets agree well within the uncertainties.

NGC4594 contains a thin, dusty disk seen close to edge-on. As already discussed in Pulsoni et al. (2018), the disk dominates the major axis kinematics but not that of the PNe, which mainly follow the kinematics of the luminous bulge. This is demonstrated by the fact that the PN kinematics agrees well with the kinematic extracted from a slit aligned with the photometric major axis but offset by 30 arcsec, and therefore probing the bulge kinematics. The velocity dispersion profiles are azimuthally averaged, therefore MUSE data, PN, and offset slit are consistent with each other and systematically higher than the major axis velocity dispersion probing the disk.

Appendix B: The distribution of sAM in simulated TNG100 galaxies

In this section, we compare the sAM radial distribution of simulated TNG100 galaxies with observations. Previous works using IllustrisTNG galaxies found that their j_t − M_* relation as well as its dependence on morphology agrees well with the observations (Di Teodoro et al. 2023; Rodriguez-Gomez et al. 2022).

However, Pulsoni et al. (2020) found hints for a different radial distribution of angular momentum in the TNG100 galaxies compared to observations. The simulated galaxies in fact display shallower V_rot/σ profiles with radius, which tend to peak at larger radii (at a median 3R_e compared to the median 1.3R_e in ePN.S), but comparable mass-size relation. This difference in the V_rot/σ profile shapes could not be attributed to differences the sample selection or to resolution effects, as the effect was present also in galaxies from the higher resolution run of IllustrisTNG, TNG50.

Figure B.1 compares the TNG100 ETGs with the Atlas3D ETGs, a survey that targets a volume and magnitude limited sample of 260 morphologically selected ETGs (Cappellari et al. 2011). The Atlas3D ETGs have mostly red colors, nevertheless we apply the same color selection used for selecting the red sequence galaxies TNG100 to single out the few bluer objects (g − r colors from the NASA-Sloan Atlas⁸ or, when unavailable, B− colors from the Hyperleda⁹ catalog converted to g − r using the relations in Pulsoni et al. 2020). The Atlas3D stellar masses are derived from the 2MASS K-band magnitudes as described in Sect. 2.4 for the ePN.S galaxies. Their projected angular momentum j_p(≤ 1R_e) is derived by applying Eq. (6) to the Atlas3D velocity fields and fluxes, integrated within elliptical apertures of semi-major axis a = 1R_e (ellipticities and R_e from Cappellari et al. 2011 and Emsellem et al. 2011).

Fig. B.1. Projected sAM jp(≤ 1Re) in the TNG100 ETGs compared to Atlas3D galaxies. Colored symbols show the sample of ETGs, gray small symbols show all the other galaxies (namely the bluer and lower-mass systems). The ETGs are color-coded according to the angular momentum parameter λe integrated within 1Re. Arrows on the Atlas3D galaxies show values integrated out to a < 1Re, hence the show jp is likely a lower-limit estimate. Dashed black lines show the jt − M* relations for different bulge fractions β as derived by Fall & Romanowsky (2018). The different jp − M* relation between observed and simulated galaxies shows that the TNG100 ETGs have their angular momentum distributed at larger radii.

The j_p(≤ 1R_e)−M_* diagram for the Atlas3D ETGs is shown in the top left panel of Fig. B.1. This is compared with the j_p(≤ 1R_e)−M_* relation for the TNG100 ETGs observed along the z-axis of the simulation box: we show results for both j_p(≤ 1R_e) weighted with fluxes (top right panel) and weighted with stellar masses (bottom left panel). We also include with small gray symbols all the other TNG100 galaxies at z = 0 that are not part of the selected ETG sample. The location of the TNG100 ETGs appears shifted toward lower j_p(≤ 1R_e) by more than 0.5 dex compared to Atlas3D in both light- and mass-weighted cases. Even the blue star-forming systems in TNG (in gray) barely reach the j_p(≤ 1R_e) values of the Atlas3D ETGs. This figure demonstrates that the difference between simulated and observed ETGs does not originate from the sample selection: galaxies with j_p(≤1R_e) as high as in Atlas3D are just not produced by the simulation.

The bottom right panel shows that we need to integrate j_p out to ∼3R_e to move the TNG100 ETGs roughly into the same location as the Atlas3D galaxies. This result shows that even though the simulated TNG100 ETGs have total sAM converging to the observed values (Di Teodoro et al. 2023; Rodriguez-Gomez et al. 2022) and aperture values at 6R_e consistent with the ePN.S galaxies (Fig. 8), their j_p profiles increase more slowly with radius compared to observations in the central regions. In Pulsoni et al. 2020, we suggested that the shallower j_p(≤a) profiles in the TNG galaxies could be due to a too efficient conversion of the gas into stars that does not allow the gas to collapse to small enough radii before forming stars. The investigation of origin of this discrepancy is beyond the scope of this paper.

Appendix C: Comparison of stellar masses and specific angular momenta with previous work

Figure C.1 compares different stellar mass determinations, assuming a universal Chabrier IMF. The stellar masses used in this work are based on magnitudes in the K-band and mass-to-light ratio M_*/L_K from stellar population modeling. Fall & Romanowsky (2013) use a M_*/L_K ratio based on color, also calibrated on stellar population models. The galaxy-by-galaxy agreement is excellent at intermediate masses, however at high and low masses the values of Fall & Romanowsky (2013) are slightly off as their M_*/L_K - M_K relation is slightly less steep. Figure C.1 also shows the comparison with the stellar masses from Forbes et al. (2017) using IR observations and M_*/L_3.6μm based on age. In this case, the comparison with our determinations does not highlight any systematic effect. The mean variation between M_* values from the three methods for the same galaxies is 0.05 dex.

Fig. C.1. Comparison between the stellar masses adopted in this work, those determined by Fall & Romanowsky (2013), and the stellar masses derived by Forbes et al. (2017) based on Spitzer 3.6 μm photometry. The measurements shown assume a constant Chabrier IMF.

Figure C.2 compares the total sAM derived in this work (the light weighted j_t, light) with the values from Romanowsky & Fall (2012) for the same galaxies. The previous j_t measurements are systematically larger by a median 0.3 dex (a factor 2). This discrepancy can be at least partially explained by their assumption of cylindrical velocity field in galaxies with rotation concentrated along the major axis, overestimating j_t.

Fig. C.2. Galaxy-by-galaxy comparison of jt determinations of this work and Romanowsky & Fall (2012) for the subset of ePN.S galaxies in common. Galaxies without error-bars on the y-axis are those for which Romanowsky & Fall (2012) used the approximation jt ∝ ReVrot(2Re), see their Eq. 6. The discrepancies between measures can be partially explained by the assumption of cylindrical velocity fields in Romanowsky & Fall (2012), which lead to an overestimate of jt. Symbols as in Fig. 12.

To elucidate this argument, we consider the two E-FRs NGC0821 and NGC3377, for which Romanowsky & Fall (2012) used the same PN data set to trace the halo kinematics. We checked that the major-axis V_rot profiles extracted here from the two-dimensional stellar and PN kinematic data are in fact consistent with the major-axis profiles published in their Figs. 28 and 30. Yet the total sAM j_t measured by Romanowsky & Fall (2012) for these galaxies are more than a factor of two larger than the values obtained in the current analysis.

The reconstruction of j_t from the observable sAM involves additionally (i) an extrapolation from the radial coverage of the data to the whole extent of the galaxy, and (ii) a deprojection to the true three-dimensional j_t. In the case of these two relatively low-mass E-FRs, (i) is not critical. The correction from the measured to the galaxy-integrated j_p estimated from the TNG simulations is of order 10 − 20% (see Fig. 9), and the j_p profiles are essentially converged within the radial coverage of the PNe (see Fig. 4), as is also shown by the profiles in Figs. 28 and 30 of Romanowsky & Fall (2012). This means that the dominant difference in j_t between the two studies must come from the different definitions of j_p and the correspondingly different deprojection factors C_i.

The direct comparison of the j_p values between the two studies is nontrivial because of their different definitions. In this work, j_p is the projected sAM defined in Eq. 5. C_i is defined to be the ratio between the two-dimensional and the three-dimensional sAM (Eq. 12) and is adjusted to each galaxy based on its flattening and rotational support (Sect. 7.2). In Romanowsky & Fall (2012), $j_p^{RF}$ is a quantity derived from the major-axis kinematics which, if the galaxy has cylindrical rotation, is readily connectable to j_t through its own "deprojection factor" $C_i^{RF}$. They use an inclination averaged value which accounting for inclination bias in the morphological classification of galaxies, and galaxies classified as ellipticals are assigned a uniform $C_i^{RF}$ corresponding to a mean inclination of 0.72 radians, that is, about 40 degrees.

Using only the major-axis rotation of NGC0821 and NGC3377, and deriving j_p via Eq. 5 assuming no dependence of V_rot on the distance from the major axis y_n, as for an edge-on cylindrical rotator¹⁰, would increase the obtained j_p by a factor 1.3 and 1.9, respectively. Adopting a deprojection factor C_i for an average inclination of 40 degrees, gives a larger correction than the flattening of these galaxies would imply (a factor 3.5 rather than 2.4, see also Fig. 10). This returns j_t values of 281 kpc km/s for NGC0821 and of 378 kpc km/s for NGC3377, closer to the measurements of Romanowsky & Fall (2012) which are ${346}_{-165}^{+264}$ kpc km/s for NGC0821 and 330 ± 50 kpc km/s for NGC3377. This experiment shows that the cylindrical rotation assumption naturally leads to an overestimate of j_t and can, at least in part, explain the differences seen in Fig. C.2. Additionally the inclination averaged $C_i^{RF}$ may amplify this for galaxies with inclined disk components.

Returning to Fig. C.2, a notable exception is the merger NGC5128, for which we measure a higher j_t compared to Romanowsky & Fall (2012). This galaxy has a large contribution to j_t coming from minor axis rotation that is not accounted for in their major-axis-based measurements.

Appendix D: Specific angular momentum measured from velocity fields versus particles

Figures D.1 and D.2 compare measurements of j_p and Λ in simulated TNG100 ETGs from projected mean velocity fields and from discrete particles. The results are shown for the 1327 selected ETGs projected along the z-axis of the simulation box. Each data-point in the figures is a local measurement within a galaxy.

Fig. D.1. sAM jp derived from smoothed velocity fields and compared with jp from discrete particles in simulated TNG100 ETGs at each radius.

Fig. D.2. Parameter Λ derived from smoothed velocity fields and compared with Λ from discrete particles in simulated TNG100 ETGs at each radius.

The good agreement between the two methods allows us to derive j_p and Λ directly from particles without deriving the mean velocity fields, and quickly extend the calculation to 1327 × 100 random line-of-sight projections. It also implies that we can consistently derive Λ from the ePN.S mean velocity fields.

Figure D.3 shows the Λ(≤ a) profile of an example galaxy, NGC3115. We show the comparison between light-weighted and mass-weighted Λ(≤ a), derived using the blue and IR photometry respectively, and the comparison with the traditional λ parameter defined by Emsellem et al. (2007).

Fig. D.3. Angular momentum parameter profiles in NGC3115 quantifying its rotational support. Different lines show the aperture λ(≤ a) profile, the light-weighted Λlight(≤ a) profile, and the mass-weighted Λmass(≤ a) profile.

All Tables

All Figures

Fig. 1. Kinematic data of NGC1023. Top: fitted kinematic parameters in elliptical bins as a function of major axis distance a for NGC1023 (see Eq. (1)). Results on the Atlas3D velocity field are shown in red, and those on the ePN.S smoothed velocity field are in green. Blue symbols show rotational velocities and kinematic position angles from kinemetry on SLUGGS data (Foster et al. 2016). These agree within errors with the Vrot and PAkin values measured from Atlas3D. ePN.S values also agree within errors except for the first PN velocity point at 70 arcsec, where the PN systematically underestimate the mean velocity (Pulsoni et al. 2018, Sect. 7.2). In the radial range marked with vertical lines, we therefore use the Vrot(a) and PAkin(a) profiles from SLUGSS, and Vs3(a) and Vc3(a) from the linear interpolation of the values fitted on the Atlas3D and ePN.S velocity fields. Bottom: the two-dimensional rotation velocity field of NGC1023 reconstructed from Atlas3D, SLUGGS, and ePN.S., using Eq. (1). The velocity field is rotated so that the photometric major axis (PAphot = 83.3 deg) is horizontal. The final mean velocity field is given by the Atlas3D data in the center (within the black square), the reconstructed two-dimensional field from SLUGGS at intermediate radii (within the black ellipse), and the ePN.S data at large radii. The PN positions are highlighted by open black circles.

In the text

	Fig. 2. Photometric data of NGC1023. Top: photometric profiles for NGC1023 extrapolated to large radii. Bottom: the galaxy image of NGC1023 reconstructed on a grid of coordinates from the photometric profiles above.
In the text

	Fig. 3. Local λ profiles of FRs (top panel) and SRs (bottom panel). S0 galaxies are shown with dashed lines, elliptical FRs with solid lines. The two mergers NGC1316 and NGC5128 are highlighted in gray.
In the text

Fig. 4. Light-weighted aperture sAM profiles for fast and slow rotators (left panels): ellipticals and S0s are distinguished by solid and dashed lines, respectively. The two mergers, NGC1316 and NGC5128, are shown with dotted gray lines. The right panels show the median jp(≤R)/jp(≤6Re) (top) and the median cumulative angular momentum Jp(≤R)/Jp(≤6Re) profiles (bottom) for galaxies divided into elliptical FRs, S0s, and SRs. The two recent major mergers NGC1316 and NGC5128 are highlighted with gray lines.

In the text

	Fig. 5. Median contribution to the local j from misaligned rotation as a function of major axis distance for different groups of galaxies as in the legend. The median profiles have been smoothed to highlight the radial trends. The profiles for the two mergers NGC1316 and NGC5128 are shown separately.
In the text

Fig. 6. Light-weighted and mass-weighted jp(≤a) profiles in four example galaxies, two FRs (top) and two SRs (bottom). For each galaxy, we show the profiles weighted by the blue-band fluxes, by the IR fluxes, by the mass-to-light ratio profile given by the colors, and by the IR fluxes corrected for IMF-driven gradients in the mass-to-light ratio. The error bar shows the error on jp derived from the errors on the colors and the dispersion in the mass-to-light ratio versus color relations.

In the text

Fig. 7. Projected sAM as a function of stellar mass for the ePN.S galaxies. Lenticulars, fast, and slow rotators ellipticals are shown with different colors and symbols as in the legend. We show with progressively decreasing sizes jp, light vs. M*, jp, mass vs. M*, and jp, mass + IMF vs. M*, var IMF. Arrows connect different measures for the same galaxies. Dashed lines show the power-law fit to jp, light vs. M*; solid lines show the fit to jp, mass vs. M*; dotted lines show the fit to jp, mass + IMF vs. M*, var IMF. Purple lines show the fits to the elliptical galaxies (E-FRs+SRs), while green lines show the fits to the lenticulars.

In the text

Fig. 8. Comparison between the projected sAM of each ePN.S galaxy jp,ePN.S and the median of the distribution of jp(≤6Re) of the ePN.S analogs among the TNG100 ETGs, divided by the sum in quadrature of 1σ-scatter of distribution of the TNG analogs and the error on the ePN.S jp,ePN.S. The three histograms correspond to the three jp determinations for the ePN.S galaxies: the light-weighted jp, light, the mass-weighted jp, mass, and the jp, mass + IMF corrected for IMF gradients. The jp values for the TNG100 galaxies are mass-weighted. Vertical dashed lines show the median of the three distributions.

In the text

Fig. 9. Median difference between the total, galaxy-integrated jp and jp ≤ 6Re, divided by jp, as a function of stellar mass in the TNG100 ETGs. Galaxies are divided into FRs (blue lines) and SRs (red lines). Solid thick lines are for jp integrated out to 15Re; the shaded areas report the quartiles of the distribution. The dashed lines are for jp integrated out to the outermost bound particle, while thin solid lines are for apertures of 10 and 12.5Re.

In the text

	Fig. 10. Distribution of Ci for different inclinations measured on 100 random line-of-sight projections of the TNG100 ETGs, as labeled in the legend. The legend also lists the mean, the median, and the standard deviation of each distribution. The gray histogram show the distribution of Ci for all inclinations.
In the text

Fig. 11. Distribution of Ci in the Λ − ε plane, where Λ and ε are measured at the location of maximum rotational support for the entire sample of simulated ETGs. Galaxies are divided in low-mass FRs (left), high-mass FRs (center), and SRs (right). The first row of panels shows the median μ(Ci) in each bin. The second row shows the root-mean-square scatter around the mean, σ(Ci). The third row shows the ratio between scatter and median. The location of the ePN.S galaxies in the Λlight − ε plane is also shown with open circles. For each ePN.S ETG, the median Ci value of its “analogs”, selected to have similar Λ, ε, and M*, is used to “de-project” jp (see Sect. 8).

In the text

Fig. 12. jt − M* plane for the ePN.S galaxies (left) and comparison with previous works (right). Left panel: total sAM of the ePN.S galaxies as a function of their stellar mass. For clarity we only show the mass-weighted jt, mass values data-points. Solid lines show the power-law fit to jt, mass vs. M*; dashed lines show the fit to the mass-weighted jt, light vs. M*; dotted lines show the fit to the values corrected for IMF gradients jt, mass + IMF − M*, var IMF using the mean mass excess profiles from Bernardi et al. (2023) as discussed in Sect. 4.2.2. Green lines refer to S0 galaxies, purple lines refer to the ellipticals. Right panel: dashed green and purple lines show the power-law fits to the ePN.S lenticulars and elliptical galaxies, respectively, as in the left panel. Dotted black lines show the jt − M* relations for different bulge fractions β as derived by Fall & Romanowsky (2018). Dashed and solid black lines show the fits of Romanowsky & Fall (2012) and Fall & Romanowsky (2013) to their sample of S0s and ellipticals. In both panels, the magenta line shows the results of Di Teodoro et al. (2023) on their samples of spiral galaxies.

In the text

	Fig. 13. Comparison with simulations. The distribution of ePN.S galaxies on the jt − M* plane (mass-weighted measurements with Chabrier IMF) is compared to the fit on their simulated analogs from TNG100 and to results from other works investigating ETGs in cosmological-hydrodynamical simulations (see text).
In the text

Fig. 14. Ratio of the stellar to halo sAM as a function of stellar mass fj(M*) for ETGs assuming the SMHM relation of Moster et al. (2013). Ellipticals and S0s are shown with thick purple and green lines, respectively. Colored bands show the scatter around the median given by the Monte-Carlo simulations described in the text. The thinner purple lines show the median fj(M*) profiles using fM from Dutton et al. (2010) for ETGs, Behroozi et al. (2013, 2019) for quiescent galaxies as labeled. For comparison, we also show the results of Posti et al. (2019) for spiral galaxies with a magenta line.

In the text

	Fig. A.1. Kinematics profiles for the new MUSE velocity fields. The profiles from MUSE are shown with red symbols. Blue symbols show major axis slit data available in the literature (Kormendy & Illingworth 1982; Saglia et al. 2000); green circles show the PN kinematic profiles.
In the text

	Fig. A.2. As in Fig. A.1, but for NGC1399, and we also show the velocity dispersion profile from Vaughan et al. (2018) based on different MUSE observations with open symbols.
In the text

Fig. B.1. Projected sAM jp(≤ 1Re) in the TNG100 ETGs compared to Atlas3D galaxies. Colored symbols show the sample of ETGs, gray small symbols show all the other galaxies (namely the bluer and lower-mass systems). The ETGs are color-coded according to the angular momentum parameter λe integrated within 1Re. Arrows on the Atlas3D galaxies show values integrated out to a < 1Re, hence the show jp is likely a lower-limit estimate. Dashed black lines show the jt − M* relations for different bulge fractions β as derived by Fall & Romanowsky (2018). The different jp − M* relation between observed and simulated galaxies shows that the TNG100 ETGs have their angular momentum distributed at larger radii.

In the text

	Fig. C.1. Comparison between the stellar masses adopted in this work, those determined by Fall & Romanowsky (2013), and the stellar masses derived by Forbes et al. (2017) based on Spitzer 3.6 μm photometry. The measurements shown assume a constant Chabrier IMF.
In the text

Fig. C.2. Galaxy-by-galaxy comparison of jt determinations of this work and Romanowsky & Fall (2012) for the subset of ePN.S galaxies in common. Galaxies without error-bars on the y-axis are those for which Romanowsky & Fall (2012) used the approximation jt ∝ ReVrot(2Re), see their Eq. 6. The discrepancies between measures can be partially explained by the assumption of cylindrical velocity fields in Romanowsky & Fall (2012), which lead to an overestimate of jt. Symbols as in Fig. 12.

In the text

	Fig. D.1. sAM jp derived from smoothed velocity fields and compared with jp from discrete particles in simulated TNG100 ETGs at each radius.
In the text

	Fig. D.2. Parameter Λ derived from smoothed velocity fields and compared with Λ from discrete particles in simulated TNG100 ETGs at each radius.
In the text

	Fig. D.3. Angular momentum parameter profiles in NGC3115 quantifying its rotational support. Different lines show the aperture λ(≤ a) profile, the light-weighted Λlight(≤ a) profile, and the mass-weighted Λmass(≤ a) profile.
In the text