Illustris early-type galaxies

The inner structure of early-type galaxies in the Illustris simulation

Dandan Xu, Volker Springel, Dominique Sluse, Peter Schneider, E-mail: Dandan.Xu@h-its.org    Alessandro Sonnenfeld, Dylan Nelson, Mark Vogelsberger, Lars Hernquist
Heidelberg Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany
Zentrum für Astronomie der Universität Heidelberg, Astronomisches Recheninstitut, Mönchhofstr. 12-14, 69120 Heidelberg, Germany
STAR Institute, Quartier Agora - Allée du six Août, 19c B-4000 Liège, Belgium
Argelander-Institut fr Astronomie, Universitt Bonn, Auf dem Hgel 71, 53121 Bonn, Germany
Kavli Institute for the Physics and Mathematics of the Universe of Tokyo, 5-1-5 Kashiwanoha Kashiwa, 277-8583 Japan
Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, Postfach 1317, D-85741 Garching, Germany
Department of Physics, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, United States
Harvard Astronomy Department, 60 Garden Street MS 46, Cambridge, MA 02138, United States
Abstract

Early-type galaxies provide unique tests for the predictions of the cold dark matter cosmology and the baryonic physics assumptions entering models for galaxy formation. In this work, we use the Illustris simulation to study correlations of three main properties of early-type galaxies, namely, the stellar orbital anisotropies, the central dark matter fractions and the central radial density slopes, as well as their redshift evolution since . We find that lower-mass galaxies or galaxies at higher redshift tend to be bluer in rest-frame colour, have higher central gas fractions, and feature more tangentially anisotropic orbits and steeper central density slopes than their higher-mass or lower-redshift counterparts, respectively. The projected central dark matter fraction within the effective radius shows a very mild mass dependence but positively correlates with galaxy effective radii due to the aperture effect. The central density slopes obtained by combining strong lensing measurements with single aperture kinematics are found to differ from the true density slopes. We identify systematic biases in this measurement to be due to two common modelling assumptions, isotropic stellar orbital distributions and power-law density profiles. We also compare the properties of early-type galaxies in Illustris to those from existing galaxy and strong lensing surveys, we find in general broad agreement but also some tension, which poses a potential challenge to the stellar formation and feedback models adopted by the simulation.

keywords:
gravitational lensing: strong - galaxies: haloes - galaxies: structure - cosmology: theory - dark matter.
pagerange: The inner structure of early-type galaxies in the Illustris simulationReferencespubyear: 2016

1 Introduction

Cosmic structures grow in a hierarchical fashion whereby small halos merge to form larger ones. This prevailing picture of structure assembly is predicted by the cold dark matter (CDM) cosmological theory. Early-type galaxies are in some sense the end products of the corresponding galaxy merging and accretion processes (e.g., White & Rees 1978; Davis et al. 1985; Kauffmann et al. 1993; Cole et al. 1994; Kauffmann 1996), and thus provide an interesting testing ground of the CDM cosmological model.

Due to their association with mergers, early-type galaxies tend to live in high-density environments and have quite old stellar populations. Though traditionally thought to be structureless and “red and dead”, mounting evidence has shown over the past decades that there is considerable richness and complexity in the origin and evolution of their mass-size relations, star formation activities, and central density profiles, etc. Utilizing for example their fundamental plane relations and/or their gravitational lensing effects, early-type galaxies are also widely used as probes of the high-redshift Universe, making them a powerful and important tool in modern astrophysics and cosmology.

Among the various properties of early-type galaxies, their central dark matter fraction and central density profiles are particularly closely tied to their formation and evolution histories. The CDM model predicts universal NFW (Navarro et al. 1997) profiles for the dark matter distribution in halos over a wide range of mass scales. However, on galaxy-scales, baryonic matter strongly dominates the central regions of dark matter halos. As a result of dissipation, baryons follow more centrally concentrated density distributions, which in turn changes the central mass fraction of dark matter (e.g., Treu & Koopmans 2004; Koopmans et al. 2006; Napolitano et al. 2010; Barnabè et al. 2011; Ruff et al. 2011; Cappellari et al. 2013; Sonnenfeld et al. 2015; Oguri et al. 2014) and also modifies the inner dark matter slopes, making them steeper than the NFW prediction (e.g., Sonnenfeld et al. 2012; Grillo 2012; Johansson et al. 2012; Remus et al. 2013; Cappellari et al. 2013; Oguri et al. 2014).

Regarding the total density profiles in the central regions, the most intriguing fact is that the sum of dark and baryonic matter approximately follows an isothermal profile, i.e. , even though neither the dark nor the baryonic matter exhibit an isothermal distribution individually. Evidence for such a profile comes from stellar kinematical studies (e.g., Binney & Tremaine 2008; Cappellari et al. 2015), strong and weak lensing observations (e.g., Rusin et al. 2003; Koopmans et al. 2006, 2009; Gavazzi et al. 2007; Barnabè et al. 2009, 2011; Auger et al. 2010b; Ruff et al. 2011; Bolton et al. 2012; Sonnenfeld et al. 2013), as well as X-ray studies of early-type galaxies (e.g., Humphrey et al. 2006; Humphrey & Buote 2010). Theoretically, the formation of such a total density distribution is speculated to occur through a two-phase process, where active (central) star formation and adiabatic contraction in an early stage is followed by dissipationless mergers and accretion later on. The former steepens the central density slopes while the latter in general makes them shallower. The observed central density slopes and their evolutionary trend, therefore, put stringent constraints on both the CDM structure formation model and models for baryonic physics processes.

Galaxy-scale strong gravitational lensing is among the major tools to probe galaxies out to high redshifts. It robustly measures the projected mass within the Einstein radius, typically at within a few kpc from the centre of a lensing galaxy. Traditionally, the strong lensing technique has also been combined with stellar kinematics, which provides the mass measured within some different aperture radius (Romanowsky & Kochanek 1999). The combination thus allows measurements of the radial density slopes and dark matter fractions. The method has been put into good use for many existing strong lensing surveys, e.g., the Lenses Structure and Dynamics Survey (LSD; e.g., Treu & Koopmans 2004); the Sloan Lens ACS Survey (SLACS; Koopmans et al. 2006, 2009; Bolton et al. 2008a; Auger et al. 2010b), the BOSS Emission-Line Lens Survey (BELLS; Brownstein et al. 2012; Bolton et al. 2012), and the Strong Lensing Legacy Survey (SL2S; Ruff et al. 2011; Gavazzi et al. 2012; Sonnenfeld et al. 2013, 2015). To date, hundreds of strong-lensing early-type galaxies have been well studied out to redshift . The average central density slopes have been found to be approximately consistent with isothermal radial profiles with a small intrinsic scatter. This isothermal behaviour seems to have evolved very little in the past 7 Gyrs (i.e. since ).

To date various theoretical approaches, including semi-analytical models and N-body simulations have been exploited in order to address open issues regarding the mass-size relations, the central dark matter fractions, as well as the total density slopes of early-type galaxies (e.g., Nipoti et al. 2009a, b; Johansson et al. 2012; Remus et al. 2013; Dubois et al. 2013; Sonnenfeld et al. 2014). These studies have made important progress in finding the missing links in the framework of the formation and evolution theories and in understanding potential systematic biases of the observational techniques. However, these modelling techniques often lacked a self-consistent treatment of baryonic physics in a cosmological context, which limited their predictive power.

In this regard, the latest generation of cosmological hydrodynamical simulations of galaxy formation represents a significant step forward and enables a much closer comparison between theoretical predictions and observations (e.g., Wellons et al. 2015, 2016; Remus et al. 2016). Among these new simulations is the Illustris Project (Vogelsberger et al. 2014a, b; Genel et al. 2014; Sijacki et al. 2015; Nelson et al. 2015), which provides an ideal tool for such purposes. Run with the accurate moving-mesh hydro solver AREPO (Springel 2010), the Illustris simulation took into account a wide range of baryonic processes, resolved the formation of 40 000 galaxies of different morphology types, and managed to reproduce many fundamental properties of observed galaxies.

In this paper, we report a variety of properties of early-type galaxies in the highest resolution simulation of the Illustris project111The highest resolution Illustris run covers a cosmological volume of (106.5 Mpc) and has a dark matter mass resolution of and an initial baryonic mass resolution of , resolving gravitational dynamics down to a physical scale of pc.. In particular, we investigate the dependencies and the redshift evolution since of (1) the stellar orbital anisotropies, (2) the central dark matter mass fractions, and (3) the central radial density slopes over the past Gyrs. The main aim of this work is to unveil correlations between these galaxy properties and to link them to underlying physical processes. Also, we are interested in identifying differences between simulation predictions and observations, through which one can establish systematic biases of observational techniques and the interpretations of the measurements.

The paper is organized as follows. In Sect. 2, we describe in detail how the light distributions of the simulated galaxies were determined and how the observed galaxy properties were measured. In Sect. 3, we report general properties of the selected galaxies, including galaxy and total matter morphologies (Sect. 3.1), the mass-size-velocity dispersion relations (Sect. 3.2), and fundamental plane properties (Sect. 3.3). We then present the dependencies and the redshift evolution of the stellar orbital anisotropies in Sect. 4, those of the central dark matter fraction in Sect. 5 and those of the central density slopes in Sect. 6. Finally, a discussion and our conclusions are given in Sect. 7. We note that all the galaxy properties reported in this paper have been made publicly available from the Illustris website. In the Appendix, we give a brief summary of the content of this online catalogue.

In this work, we adopted the same cosmology as used in the Illustris simulation, i.e., a matter density of = 0.27, a cosmological constant of = 0.73, a Hubble constant and a linear fluctuation amplitude . These values are consistent with the Wilkinson Microwave Anisotropy Probe (WMAP)-9 measurements (Hinshaw et al. 2013).

2 Illuminating galaxies and measuring galaxy properties

The galaxies simulated in Illustris are identified as gravitationally bound structures of gas cells, dark-matter particles, and stellar particles using the subfind algorithm (Springel et al. 2001; Dolag et al. 2009). In this section, we describe how we calculated in a post-processing procedure observational properties for the simulated galaxies.

It is noteworthy to recall that when a smaller halo is accreted onto a bigger structure, its dark and baryonic matter at the outskirt can be tidally stripped while sinking to the centre of the host. As a result, a (galaxy) halo in a group or cluster environment can be composed of a tightly-bound central region and a loosely-bound outskirt that extends to large radii. Observationally, the measured “galaxy” properties are only accounting for the former. However, the latter component, which in the context of galaxy clusters is also known as “intracluster light”, can make up for of the total stellar mass and luminosity (e.g., Lin & Mohr 2004; Zibetti et al. 2005; Puchwein et al. 2010). This stellar mass is still bound gravitationally to the galaxy and hence normally included in the raw measurement of SUBFIND.

Often, a simple radial cut-off radius has been used in numerical simulations to deal with this problem and to calculate properties that are associated with the tightly-bound galaxy component. For example, in Puchwein et al. (2010, see their eq. 1), a halo mass-dependent radial cut was applied to the central brightest cluster galaxies (BCGs) more massive than . Schaye et al. (2015) applied a three-dimensional (3-D) radial cut of 30 physical kpc to all galaxies from the EAGLE project. Such a choice was found to significantly affect the measured properties of the tightly-bound galaxy component that has a stellar mass more than . Such a radial cut was also able to reproduce the observed galaxy stellar mass function that was derived using the frequently applied Petrosian apertures (see in Schaye et al. (2015) for detailed discussion). In this work, we used a similar strategy. In calculations of projection-independent intrinsic galaxy properties, such as the stellar mass, a 3-D radial cut of 30 kpc is adopted. For the luminosity-based and thus projection-dependent properties, we used a 2-D radial cut of 30 kpc from the centre of light for a given galaxy projection. Any resolution element still bound to the system but located or projected outside this radius is excluded from the corresponding calculations.

2.1 From stellar particles to galaxy light

We calculated the emission properties of individual simulated galaxies from the luminosities of their constituent stellar particles. Each stellar particle of is treated as a coeval single stellar population that has the Chabrier (2003) initial mass function (IMF). Given the star formation time and metallicity of each stellar particle, a “raw” luminosity in a given bandpass can thus be derived using the Bruzual & Charlot (2003) stellar population synthesis (SPS) model galaxev.

This raw galaxy light can be processed by dust through absorption and scattering at shorter wavelengths and re-emission at longer wavelengths. We implemented such a dust attenuation process through a simple semi-analytical approach as follows. Assuming that a galaxy can be approximated as a uniformly mixed slab of stars, gas and dust, the amount of extinction of each stellar particle is given by

(1)

where and are the “observed” (dust-attenuated) and the raw (dust-free) luminosities, respectively. is the optical depth, which depends on the wavelength and is caused by both dust absorption and scattering along the line of sight.

In order to predict the amount of dust extinction for each stellar particle, a mesh that covers the (projected) central region of a simulated galaxy, from to times the half-stellar mass radius in either dimension, was used to tabulate the distribution of . The -mesh can be derived from the neutral hydrogen (HI) distribution of the gas cells using a semi-analytical prescription as follows.

A redshift- and metallicity-dependent optical depth due to the dust absorption process is traditionally modelled as (e.g., Guiderdoni & Rocca-Volmerange 1987; Devriendt et al. 1999; Devriendt & Guiderdoni 2000)

(2)

where is the solar-neighbourhood extinction curve (Cardelli et al. 1989), is the gas metallicity (of the galaxy) and is the measured value for the Sun. is the average neutral hydrogen column density. A power-law index for Å  and for Å  was found by Guiderdoni & Rocca-Volmerange (1987) for the metallicity dependence. In addition, Kitzbichler & White (2007) found that can reproduce measurements of Lyman-break galaxies at . We also adopted this value for .

To also account for dust scattering, we used an approximate solution from Calzetti et al. (1994), in which the total effective optical depth is given by

(3)

where is the albedo, defined as the ratio between the scattering and the extinction coefficients, and and are the weighing factors for the isotropic and the forward-only scattering, respectively.

For each cell of the -mesh, a mean neutral hydrogen column density was first calculated by scattering the (fractional) cold hydrogen masses of each gas cell onto the mesh using a SPH smoothing technique. was then converted into via Eqs. (2) and (3), assuming being the effective wavelength of a given bandpass. For stellar particles that are projected within the mesh coverage, the exact (at the position of a given stellar particle) was then interpolated from the -mesh. For those outside, was assumed. The total light distribution of a simulated galaxy in its viewing direction in a given bandpass was then computed from , the “observed” (dust-attenuated) luminosities of the constituent stellar particles.

2.2 Measurement of galaxy centres, ellipticities and orientation angles

Each simulated galaxy (together with its dark matter halo) has a “centre” that was calculated as the position of the particle with the minimum gravitational potential found using subfind. The observed galaxy centre is light-based and projection-dependent, we therefore defined a two-dimensional (2-D) galaxy centre (, ) as:

(4)

where and are the - and -coordinates of the -th stellar particle in the plane of a given galaxy projection. is the “observed” (dust-attenuated) luminosity, and is the total luminosity within a given aperture.

The orientation and ellipticity (or axis ratio) of a galaxy were measured through luminosity-weighted second moments, which are defined as:

(5)

The axis ratio of a galaxy projection is then given by:

(6)

The orientation angle is given by:

(7)

2.3 Luminosity and effective radius measurement

In order to calculate the galaxy luminosity and effective radius within non-circular apertures, we first assumed that the 2-D surface brightness distribution in a given viewing projection follows a series of elliptical isophotes that are well described by and measured using Eqs. (6) and (7) within three-times the half stellar-mass radius from the centre of the galaxy. A “summed” galaxy luminosity is calculated by directly summing up for all constituent stellar particles that are projected within 30 kpc.

A “direct” effective radius is determined as the geometric mean of the semi-major and semi-minor radii of the elliptical isophote which encloses half of . We used and as approximate estimates of the intrinsic luminosity and size of a simulated galaxy. They are, however, different from those used in observations. In order to make a fair comparison, we followed one of the observational conventions to derive a “model” luminosity and a “model” effective radius , i.e., by fitting a Sersic profile (Sérsic 1963) to the radial distribution of the elliptical isophotes that are assumed to closely trace the 2-D surface brightness distribution of a given galaxy projection.

The surface brightness of a Sersic profile at radius is given by:

(8)

where is the effective radius of the Sersic profile, and is the Sersic index. The factor can be determined by satisfying . We adopted the values of from Ciotti & Bertin (1999) for and those of MacArthur et al. (2003) for . The luminosity of a Sersic distribution enclosed within a radius of is given by:

(9)

where is the incomplete gamma function, . Note that in the formulae above, the radius is defined as the geometric mean of the semi-major and semi-minor radii of the elliptical isophotes.

For each galaxy projection, the binned radial profile between 0.05 and 3.0 was fitted by the Sersic model using a minimum fitting approach. In general, the radial surface brightness distributions of the simulated galaxies well follow Sersic profiles. The model effective radius of a given galaxy projection was set to be the best-fitting Sersic effective radius. The model luminosity was then given by the best-fitting Sersic luminosity within .

We measured the light distribution of each simulated galaxy in a variety of optical filter bandpasses in order to compare with latest observations. Note that galaxy luminosities (magnitudes) and effective radii are band-dependent (also see La Barbera & de Carvalho 2009). Hereafter, the effective radii specifically refer to the Sersic model effective radius that was derived in the rest-frame Johnson -band. For galaxies at measured in this band, the luminosity ratios have a (stellar-mass weighted) mean of with a standard deviation of 0.16, and the effective radii ratios have a (stellar-mass weighted) mean of with a 0.33 scatter. These ratios vary little across the investigated redshift range.

In Fig. 1 we present the dust-corrected luminosity functions (LF) of the Illustris galaxies at intermediate redshifts (), measured in the Johnson- band. Plotted in the same panel are the Schechter function fits to the observed LFs from DEEP2 and COMBO-17 surveys (Faber et al. 2007). The red and blue galaxy samples were selected from the simulation according to eq. 1 in Faber et al. (2007). The Illustris galaxy LFs have been found to roughly match observational results in various bands and within a wide range of redshifts (also see Vogelsberger et al. 2014b; Hilbert et al. 2016). We note, in particular, that the effect of dust attenuation is significant especially for blue galaxy sample, which, without taking dust into account, would contribute to at least half of the total galaxy count at any given luminosity up to . In comparison, the observed LFs have this transition at two magnitudes higher, where the blue and red galaxy samples have equal contributions to the total LF and below which (in the brighter end), red galaxies would dominate the total LF. As can be seen, applying dust attenuation to the simulated galaxies has pushed this transition magnitude much closer to the observation. At the bright end, there are larger measurement uncertainties, and the derived magnitudes can differ strongly depending on the assumed light profiles (e.g., Bernardi et al. 2013).

Figure 1: Galaxy luminosity functions at intermediate redshifts (), measured in Johnson- band. Simulation data are given by square symbols, where black, red and blue colours represent the total, “red”, and “blue” galaxy samples, respectively. The latter two are rescaled by a factor of 0.1 for graphical clarity. The red and blue galaxy samples were selected from the simulation according to eq. 1 in Faber et al. (2007). The dotted and the dashed lines give single power-law Schechter function fits to the observed galaxy luminosity functions from the COMBO-17 and DEEP2 surveys, respectively.

2.4 Galaxy type classification

The method that we used for galaxy classification is similar to the practice of the Sloan Digital Sky Survey (SDSS). Using the rest-frame SDSS , and filters simultaneously, we fitted both de Vaucouleurs profiles (de Vaucouleurs 1948) and exponential profiles to the radial surface brightness distributions of the elliptical isophotes (between 0.05 to 3.0). If the former provides a better fit, then the galaxy is classified as an early-type (elliptical) galaxy, whereas if the latter fits better, then it is considered a late-type (disk) galaxy.

Fig. 2 shows histograms of the best-fitting Sersic indices of the simulated early- and late-type galaxies at redshift . As expected, the former have larger () Sersic indices, while the latter have . As an illustration, we select two typical galaxies with different type classifications and show their synthesized light distributions in the top panel of Fig. 3, where the left and right sub-panels display an elliptical and a disk galaxy at , respectively. The image was made by combining the surface brightness distributions in the rest-frame SDSS , and filter bandpasses, within a region from the light centres of the corresponding galaxy projections.

Figure 2: Histograms of the best-fitting Sersic indices of galaxies whose radial surface brightness distributions are better fitted by de Vaucouleurs profiles (red) and of those whose radial distributions are better fitted by exponential profiles (blue).
Figure 3: The top row presents synthesized images of an elliptical (left panel) and a disk (right panel) galaxy from the simulation at redshift . The image was made by combining the surface brightness distributions in the rest-frame SDSS , and filter bandpasses, within a region from the galaxy centres. The middle and bottom panels show radial surface brightness distributions of the early- and late-type galaxies, respectively. In these panels, black dots are the binned data measured in the SDSS band, blue lines give the best-fitting Sersic profiles (Sersic indices of 2.35 and 1.27, respectively), and red lines show the best-fitting de Vaucouleurs and exponential profiles of the two galaxies, respectively.

The lower panels of Fig. 3 show the radial surface brightness distributions of the early- (middle) and late-type (bottom) galaxies above. In these panels, black dots are the binned data measured in the SDSS band, blue lines give the best-fitting Sersic profiles (Sersic indices of 2.35 and 1.27, respectively), while red lines show the best-fitting de Vaucouleurs and exponential profiles for the two galaxies, respectively.

3 General properties of the galaxy catalogue

Apart from the above mentioned luminous properties, we derived a wide range of projection-dependent lensing and dynamical properties, such as the Einstein radii , dark matter fractions , stellar velocity dispersions , orbital anisotropies , and various mass density slope estimators . All these calculations were carried out for galaxies (regardless of galaxy types) with stellar masses (corresponding to more than 10 000 stellar particles) at redshifts in the range with an interval spacing . We defined artificial source redshifts at accordingly, in order to calculate the Einstein radii222The Einstein radius is found as the radius within which the mean surface density is equal to the lensing critical density , where , and are the angular diameter distances to the lens, to the source, and from the lens to the source, respectively.. These choices for are motived by the observed lens-source redshift distributions. All the calculated properties have been catalogued and are publicly available from the Illustris website (www.illustris-project.org), and a detailed description of the different catalogue fields can be found in the Appendix.

In the following part of the paper, we aim at presenting the statistical properties of the simulated early-type galaxies and comparing them to those resulted from the SDSS early-type galaxy survey (Hyde & Bernardi 2009a), and to those from recent strong lensing surveys, i.e., the SLACS and SL2S surveys, as well as those from the Cosmological Monitoring of Gravitational Lenses projects (COSMOGRAIL; Sluse et al. 2012). These surveys provide comprehensive observational samples predominantly composed of isolated early-type galaxies within a wide redshift range, similar to the one studied here.

We are aware that the lensing selection effect has always been a complication when comparing simulation samples to observations (or using observed samples to interpret physical properties of galaxies). Detailed work in this regard can be found in, e.g., Mandelbaum et al. (2009). For the present study, we do not apply any sophisticated selection criteria. We selected Illustris galaxies at each given redshift according to the following two simple criteria: (1) central and early-type galaxies, and (2) the stellar line-of-sight velocity dispersions (measured within 0.5) satisfies . These criteria are mainly observationally motivated by the strong lensing galaxy surveys. As will be seen in Sect. 3.2, the criteria above resulted in galaxy samples that roughly reproduced the observed mass-size-velocity dispersion relations.

It is worth noting that when calculating properties of a central galaxy, any self-gravitationally-bound substructures (satellite galaxies) identified by subfind were excluded. To increase the sample size, we treated galaxies that are viewed along the three principal directions of the simulation box as independent. The selection criteria resulted in independent galaxy projections at each of the simulation redshifts investigated. In this section, we present the measured shapes of the luminous and dark matter distributions (Sect. 3.1), mass-size-velocity dispersion relations (Sect. 3.2) and fundamental plane relations (Sect. 3.3) of the selected galaxy samples at various redshifts.

3.1 Shape of luminous and total matter

It has always been a critical question “how well the light follows the matter”. In central regions of early-type galaxies, the projected total matter distributions are in general rounder333The situation is the opposite at larger radii for which weak lensing technique can be applied. At those radii, dark matter dominates the total matter distribution, which appears to be flatter than the light profile (e.g., Hoekstra et al. 2004). than the luminous (stellar) distributions. This can be seen from the top panel of Fig. 4, which shows the ratio between the luminous axis ratio and the total axis ratio , as a function of the central stellar velocity dispersion : the solid blue and red line indicates the median of the distribution measured within and , respectively; the dashed lines show the 90% boundaries of the distributions. Note that the ellipticity ratio decreases as the aperture size increases from to . This is mainly attributed to the fact that as the aperture size increases the total matter distribution becomes rounder, while the ellipticity of the stellar distribution only varies mildly. The ellipticity ratio distribution does not seem to evolve strongly with redshift, at higher redshift the scatter increases marginally.

Observationally, the shape of the total matter distribution of a (lensing) galaxy can be inferred via a lens modelling technique, while that of the stellar distribution can be obtained from direct imaging. In this regard, previous studies based on singular isothermal ellipsoidal (SIE) lens models found that the median ellipticity ratios of the galaxy samples from the SLACS, SL2S and COSMOGRAIL (see Koopmans et al. 2006; Gavazzi et al. 2012; Sluse et al. 2012; Shu et al. 2015) are very close to 1.0, which lies above both the median and the 90% upper boundary of the simulation distribution. Regarding this disagreement, we note that, as will be seen in Sect. 5, the central dark matter fractions in the Illustris early-type galaxies were found to be systematically higher compared to the observational results. This may also explain the systematically lower ellipticity ratios between the luminous and the total matter distributions of the simulated galaxies. However, it is also noteworthy to point out that recent lensing studies that have either applied non-parametric lens modelling (e.g., Bruderer et al. 2016) or adopted novel techniques to extract galaxy morphology data from adaptive optics observations (e.g., Rusu et al. 2016) found that the total matter distributions are systematically rounder than the stellar distributions in central regions of early-type galaxies, in agreement with our simulation results.

The orientation of the projected total matter distribution follows in general that of the light distribution very well. The bottom panel of Fig. 4 shows the misalignment angle between the orientation angles of the light distribution and of the total matter distribution, as a function of the galaxy axis ratio that was measured within . Larger scatter of for rounder galaxies is due to less clean measurements of their orientation angles. A marked evolution exists for the selected early-type galaxy samples at different redshifts, e.g., for a sub-sample of galaxies that have , the standard deviation of measured within decreases from at to at , and finally to at . This evolution coincides with a mild increase of the (projected) central baryonic fraction towards lower redshifts (see Fig. 19). In the intervening period, stars and dark matter become more mixed and thus better aligned with cosmic time. These results are consistent with strong lensing observations (e.g., Sluse et al. 2012; Dye et al. 2014, Rusu et al. 2016). In particular, for galaxies that have stellar axis ratios falling within similar ranges as above, Gavazzi et al. (2012) measured an rms scatter of for the SL2S sample () and Koopmans et al. (2006) reported an rms deviation of for the SLACS galaxy sample ().

Figure 4: Top panel: the ratio between the axis ratio measured for the light distribution and that for the total matter distribution, as a function of the central stellar velocity dispersion , for the selected early-type galaxy sample at . Bottom panel: the difference between the orientation angle of the light distribution and that of the total matter distribution, as a function of the galaxy axis ratio measured within , for the same galaxy sample. In both panels, blue and red represent the distributions measured within and , respectively. The solid and the dashed lines indicate the median and the 90% boundaries of the distributions, respectively.

3.2 Galaxy mass, size, and velocity dispersion

The relations between galaxy size, velocity dispersion and stellar mass are among the most basic scaling relations. Fig. 5 shows the and relations measured for the selected early-type galaxies at . In either panel, the black squares mark the simulation data, and the solid and the dashed red lines indicate the best linear fit to the data and the 90% boundaries of the distribution, respectively. Note that the lower cut in the stellar velocity distribution (in the upper panel) is due to our sample selection criteria.

For comparison, the solid blue and cyans lines represent the best linearly fitted relations for the SLACS galaxies (Auger et al. 2010b), assuming either Salpeter IMFs (Salpeter 1955) or Chabrier IMFs (Chabrier 2003), respectively. Recent studies suggest that a Salpeter IMF is more compatible with observational inferences for massive early-type galaxies (e.g., Auger et al. 2010a; Treu et al. 2010; Grillo & Gobat 2010; Barnabè et al. 2011; Spiniello et al. 2011; Oguri et al. 2014; Sonnenfeld et al. 2015), while a Chabrier IMF is more suitable for those at the lower end of the spectrum (e.g., Shu et al. 2015). In addition, the linearly fitted relations for the observed early-type galaxies from the SDSS survey (Hyde & Bernardi 2009a) are also given, as indicated by the solid orange lines.

As can be seen, there exists some slight underestimation of the central velocity dispersions for galaxies at the lower stellar-mass end. However, we note that the level of the disagreement between the simulation and the observations is comparable to that the systematic divisions among the observational relations themselves, which can be attributed to the uncertainty in the IMF and/or different sample biasing etc. Overall, the simulation broadly reproduces the observed and relations.

Figure 5: The (top panel) and (bottom panel) relations. The black squares mark the simulation data of the selected early-type galaxies at . The lower cut in stellar velocity distribution (in the top panel) is due to our selection criteria. The solid and the dashed red lines indicate the best linear fit to the data and the 90% boundaries of the distribution, respectively. For comparison, the solid blue and cyans lines represent the best linearly fitted relations for the SLACS galaxies (Auger et al. 2010b), assuming either Salpeter IMFs or Chabrier IMFs, respectively. The solid orange lines indicate the linearly fitted relations for the observed early-type galaxies from the SDSS survey (Hyde & Bernardi 2009a).

3.3 The fundamental plane and the mass planes

Elliptical galaxies are observed to tightly follow the so-called Fundamental Plane (hereafter FP; Faber et al. 1987; Djorgovski & Davis 1987; Dressler et al. 1987) and its mass counterparts – the Stellar Mass Plane (hereafter MP; Hyde & Bernardi 2009b) and Mass Plane (hereafter MP; Bolton et al. 2007). We examined how well the simulated galaxies trace these observed planes. To this end we adopted the same definitions as Auger et al. (2010b) and fitted the plane relations with the following form:

(10)
(11)
(12)
(13)

In Eq. (11), is the Sersic luminosity measured in the rest-frame Johnson- band. A measurement uncertainty of was assumed for each of the five quantities of the simulated galaxies. The fitting of the coefficients , and was done using of the lts_planefit program described in Cappellari et al. (2013), which combines the Least Trimmed Squares robust technique of Rousseeuw & Van Driessen (2006) with a least-squares fitting algorithm which allows for errors in all variables as well as intrinsic scatter.

Fig. 6 shows the FP relation of the selected early-type galaxy sample at from the simulation. The blue circles filled with black dots indicate individual galaxies. The red solid and dashed lines show and scatter, respectively, around the best-fitting relation, which is given by the black solid line. Galaxies indicated by green symbols are outside of the fitted relation. For the FP relation, we found best-fitting coefficients equal to , and . For the MP, we obtained , and ; and for MP, , and .

The fitting coefficient represents the normalizations of the plane relations; while and depict the slopes of the plane relations. In particular, coefficients and depend on the sample selection, observational bands, and fitting methods (see Bernardi et al. 2003). Observationally, ranges from 0.99 to 1.52, and ranges from to for the FP; for the MP, varies from 1.77 to 1.86, and ranges from to (Bernardi et al. 2003; Bolton et al.; Bolton et al. 2008b; Auger et al. 2010b; Cappellari et al. 2013). Within the fitting uncertainties, our best-fitting values for , and are not in stark disagreement with observational results. A more detailed study of the FP relation of the Illustris early-type galaxies will be presented in a separate paper (Li et al. in preparation).

Figure 6: The fundamental plane relation of the selected early-type galaxy sample at from the simulation. is and is . A measurement uncertainty of was assumed for each “observed” quantity. The blue circles filled with black dots indicate individual simulated galaxies. The red solid and dashed lines give the and bounds, respectively, from the best-fitting relation given by the black solid line. Galaxies indicated by green symbols are outside of the fitted relation. The fitting and plotting software lts_planefit used here was provided by Cappellari et al. (2013).

4 Stellar orbital anisotropies

One of the major uncertainties in interpreting kinematical data is the stellar orbital anisotropy, which cannot be directly measured or constrained with high accuracy (e.g., Li et al. 2016) as it is degenerate with galaxy density slopes, unless the latter can be determined independently. When combined with single-aperture stellar kinematics, strong lensing studies often have to assume zero orbital anisotropies (but see e.g., Barnabè et al. 2009, 2011 for 2-D kinematics studies of the SLACS lenses). In this section, we thus investigate the dependencies of the anisotropy parameter and its redshift evolution since to check the validity of this assumption.

For spherically symmetric systems, the anisotropy parameter can be written as (Binney & Tremaine 2008)

(14)

where , and are the second velocity moments, and , and are the velocity dispersions measured in the azimuthal , polar and radial directions of a spherical coordinate system. By definition, , where is the mean velocity. The second equality sign in Eq. (14) is valid for stationary non-rotating systems, where , and vanish. Note that observationally, as the measurements are carried out for the light components, , and by construction. We followed the same convention, measuring through instead of . In this sense, is constructed to measure the anisotropy of the velocity dispersion. corresponds to the “isotropic” case, and () describes a radially (tangentially) anisotropic orbital distribution.

For each simulated galaxy, measurements of were made (for stellar particles) within 3-D radii of and from the centres of galaxies. Fig. 7 shows the dependence of on for the selected galaxies at . The blue and red curves present the distributions of and , respectively. The solid and dashed lines indicate the medians and the 90% percentiles, respectively. The stellar orbits of more massive galaxies tend to be more radially anisotropic than those of their lower-mass counterparts. Observationally, Koopmans et al. (2009) applied two independent techniques to measure the logarithmic density slopes for SLACS early-type galaxies444The majority of the SLACS galaxies are located at lower redshifts () and have peaks around .. The combination of the two measurements provided a (weak) constraint on the orbital anisotropies, , consistent with the values we measured for their counterparts in the Illustris simulation.

It is interesting to note that the different behaviour of the average in low- and high-mass galaxies is possibly related to the recent star formation histories. To demonstrate this, we consider in Fig. 8 the central (cold) gas fractions versus for the same galaxy sample at . The median and 90% boundaries of the distribution are given by the solid and dashed lines, respectively. The blue and red curves show the distributions of the cold (HI) and total gas fractions, respectively, measured within a 3-D radius of from the galaxy centres. We can see that more massive galaxies that tend to have higher radial anisotropies also contain less cold gas in their central regions, while the less massive galaxies with higher tangential velocity contributions have on average higher central gas fractions. The former were also seen to have Johnson colours redder than the latter. We also note that galaxies at higher redshifts are markedly bluer and contain higher fractions of central cold gas than their lower redshift counterparts.

These correlations provide a consistent picture. As the cold gas is channelled down to the centre, star-formation activity preferably happens on tangential orbits as a consequence of gas accretion and rotational support. When the system (passively) evolves, more radial anisotropies emerge. In this case, one would expect that the stellar orbits at higher redshifts are more tangentially dominated than their lower-redshift counterparts.

Indeed, this can be clearly seen in Fig. 9, which shows the redshift evolution of since . In order to quantify this redshift evolution, we randomly assigned a galaxy from a given snapshot at with a redshift of , where . We then fit versus using a linear regression approach, which resulted in with a linear correlation coefficient and with . We verify that changing the range of from 0.01 to 0.1 makes no difference in the linear regression results.

We mention in passing that the observed correlation between and could strongly depend on the details of the adopted galactic wind and AGN feedback models, which, as shown in Genel et al. (2015), efficiently affect the gas distribution and determine the stellar angular momentum and thus orbital anisotropies. Observational constraints on the distribution of are crucial in establishing the validity of various feedback models.

Figure 7: The anisotropy parameter as a function of the central stellar velocity dispersion , measured for the selected early-type galaxies at . The blue and red curves present the distributions of and , respectively. The solid, dotted and dashed lines indicate the medians, the 68% and 90% boundaries of the distributions, respectively.
Figure 8: The central (cold) gas fractions versus for the early-type galaxy sample at . The median and 90% boundaries of the distribution are given by the solid and the dashed lines, respectively. The blue and red curves show the distributions of the cold (HI) and of the total gas fractions within a 3-D radius of , respectively.
Figure 9: The redshift evolution of the anisotropy parameter measured for the early-type galaxy samples at from the simulation. The blue and red curves represent the distributions of and , respectively. The solid, dotted and dashed lines indicate the best linear fit to the data, the 68% and 90% boundaries of either distribution, respectively. Linear regression resulted in with a linear correlation coefficient and with .

5 Projected central dark matter fractions

The projected central dark matter fraction of observed galaxies has often been constrained through combined measurements of the stellar and the total masses. A galaxy’s stellar mass can be obtained using the SPS method applied to multi-band photometric data, provided that the IMF is independently constrained. Constraints on the total mass may come from stellar kinematics and/or strong lensing measurements. As already discussed in Sect. 4, this is often complicated due to the lack of knowledge about stellar orbital anisotropies. In addition, it also needs to assume parameterized density profiles either of the total matter distribution (e.g., Koopmans et al. 2006; Auger et al. 2010b) or of the individual dark and luminous components (e.g., Cappellari et al. 2013; Sonnenfeld et al. 2015; Li et al. 2016).

In this section, we first compare the projected central dark matter fraction of the early-type Illustris galaxies to the ones derived for the SLACS (Auger et al. 2010b) and SL2S (Sonnenfeld et al. 2015) galaxies. To this end, we use the quantity , which is defined as the projected dark matter fraction within a fixed aperture of . The choice of a fixed aperture for comparison purposes instead of, say, , is to eliminate systematic differences due to possible sample bias.

Fig. 10 shows as a function of the central velocity dispersion . The blue curves give the distribution of the directly measured for the selected early-type galaxies at , whereas the red curves indicate the distribution of measurements for the same galaxies but obtained by modelling the total matter density distributions with power-law profiles (see Sect. 6.1.3 for details). The solid and the dashed lines give the median and the 90% boundaries of the simulation distributions. The black squares indicate measurements for the SLACS and SL2S galaxy samples, where the error bars show 1  error of the data. Specifically, a dark halo component was modelled by a NFW profile, and a de-projected best-fitting de Vaucouleur distribution was adopted to model the stellar distribution. The sum of the two was then used to fit both strong lensing and kinematics data (Sonnenfeld et al. 2015).

As can be seen from the figure, the distribution of the measurements derived under the power-law profile assumption (in red) has larger scatter than the true distribution (in blue). But both distributions indicate that dark matter on average contributes with to the (projected) total matter distributions in centres of early-type Illustris galaxies. This fraction is higher than suggested by the observational results. Such tension poses a potential challenge to the stellar formation and feedback models adopted by the simulation.

Observations also suggest that the projected central dark matter fraction within the effective radius is mass-dependent: the more massive galaxies are, the larger is their dark matter fraction. A noticeable positive correlation was also found between and , albeit with large scatter (e.g., Tortora et al. 2009; Napolitano et al. 2010; Humphrey & Buote 2010; Graves & Faber 2010; Auger et al. 2010b; Shu et al. 2015). Using a significantly larger statistical sample of early-type galaxies from the simulation, we investigated such dependences.

Fig. 11 presents versus (top panel) and versus (bottom panel) for the early-type galaxy sample selected at . Blue and red represent the fractions measured within a radius of and , respectively. The solid and the dashed lines give the best linear fit to the data and the 90% boundaries of the simulation distribution. Linear regression resulted in with a linear correlation coefficient ; with and with ; with .

As can be seen, for the early-type galaxy sample, the dependence on the stellar velocity dispersions is much weaker than on galaxy sizes. The latter shows a tight and clear positive correlation between the two quantities. Similar dependences were also found by Remus et al. (2016), where early-type galaxies selected from the Magneticum Pathfinder Simulations (Dolag et al. 2015) were studied.

These dependences suggest that the projected central dark matter fraction in terms of [or ] has a very mild mass dependence for our early-type galaxy samples. The clear positive correlation between [or ] and may purely be an aperture effect: the dark matter fraction drops with decreasing radius as baryons dominate more and more towards the galactic centre (also see Grillo 2010; and Fig. 3 of Xu et al. 2016).

We also studied the redshift dependence of the projected central dark matter fraction for the selected early-type galaxy samples in different redshift bins between to . The result is shown in Fig. 12, where the solid, dotted and dashed lines indicate the best linear fit to the data, the 68% and 90% boundaries of the distribution, respectively. Using linear regression approach, we found that both and are equal to and both fits have linear correlation coefficients .

Fair comparisons with observations (or studies using simulations of the same kind) would require applying identical sample selection criteria. Dye et al. (2014) reported a similar increasing trend of with redshift for a galaxy sample from the Herschel Astrophysical Terahertz Large Area Survey (see also Sonnenfeld et al. 2015, Fig. 6). The result from the simulation is not in stark contrast with observations. We note that the early-type galaxy samples at different redshifts were selected according to their , which also evolve with redshift for individual galaxies. Therefore, the redshift trends found here hold for a statistical sample defined as such, but not necessarily for individual galaxies (also see Remus et al. 2016).

Figure 10: The projected central dark matter fraction as a function of the central velocity dispersion . The blue curves present the distribution of the directly measured for the selected early-type galaxies at . The red curves indicate the distribution of the measurements for the same galaxies but obtained by assuming the total matter density distributions of galaxies to be power laws. The solid and the dashed lines give the median and the 90% boundaries of the simulation distributions. The black squares show measurements for the SLACS and SL2S galaxy samples, where the error bars indicate 1  error of the data. For the observed galaxy sample, a dark halo component was modelled by a NFW profile, and a de-projected best-fitting de Vaucouleur distribution was adopted to model the stellar distribution. The sum of the two was then used to fit both strong lensing and kinematics data (see Sonnenfeld et al. 2015 for details).
Figure 11: The projected central dark matter fraction as a function of the central velocity dispersion (top panel) and of the effective radius (bottom panel) for the selected early-type galaxy sample at . Blue and red represent the fractions measured within a projected radius of and , respectively. The solid and the dashed lines give the best linear fit to the data and the 90% boundaries of the simulation distribution, respectively. Linear regression resulted in with a linear correlation coefficient ; with and with ; with .
Figure 12: The redshift evolution of the central dark matter fractions (blue) and (red) measured for the selected early-type galaxy samples at . The solid, dotted and dashed lines indicate the best linear fit to the data, the 68% and 90% boundaries of either distribution, respectively. For the entire sample, both and are equal to and both fits have linear correlation coefficient .

6 Central matter density profiles

Figure 13: Matter density distribution of an early-type galaxy at (the same galaxy as presented in Fig. 3). The black, red and blue points indicate the density distributions of the total, dark matter and baryonic matter, respectively.

As already shown in Sect. 5, baryons contribute a large fraction of the total matter distribution in centres of early-type galaxies. Fig. 13 shows the density distribution of an early-type galaxy at (the same galaxy as presented in Fig. 3). The black, red and blue symbols indicate the density distributions of the total, dark matter and baryonic matter, respectively. The profile of baryons is much steeper than that of dark matter in the central region. Interestingly, the projected galactocentric radii where typical strong lensing and stellar kinematics data are available coincide with the radii where the radial profiles of dark matter and baryons intercept (e.g., see Fig. 3 of Xu et al. 2016). In this radial range (normally , corresponding to the inner few kpc), the slope of the total density profile depends on both components. Its quantification can be problematic because the sum of the two components does not necessarily obey simple global power-law distributions, i.e., no single slope can describe the overall distribution.

For a simulated galaxy, one can find approximate slope estimators within given radial ranges. For example, an average slope of the density profile between two radii and can be expressed using a power-law interpolation (the superscript “AV” refers to “average”):

(15)

One can also define as the local logarithmic slope of the power-law profile that best fits the radial density distribution between and (the superscript “PL” refers to “power-law”). In particular, this definition has been adopted in many studies on simulated galaxies (e.g., Nipoti et al. 2009a; Johansson et al. 2012; Remus et al. 2013; Li et al. 2016) when comparing to observations.

Another definition is a mass-weighted density slope (the superscript “MW” refers to “mass-weighted”, see Dutton & Treu 2014):

(16)

where the local slope is given by

(17)

and the 3-D enclosed total mass is given by

(18)

Note that for a matter density distribution that follows a perfect power law, i.e., , at all radii .

Despite different definitions, all the above-mentioned slope estimators quantify some intrinsic matter density distribution to first-order. The interpretation of the resulting measurements under these definitions is straightforward and model-independent. However, they cannot be applied to observed galaxies, unlike the brightness distribution, which can be directly measured as long as the galaxy is spatially resolved. The (total) matter density distribution can only be determined indirectly by dynamical methods, for example gravitational lensing or stellar kinematics, and through fitting parameterized models based on certain assumptions.

For galaxies at lower redshifts where 2-D kinematical data (e.g., the integral-field spectroscopic data) are available, one can implement sophisticated dynamical methods (e.g., Barnabè et al. 2011; Cappellari et al. 2015). Assuming parameterized density profiles within radial ranges (e.g., from to a few ) for data fitting purposes, the method allows simultaneous fitting to the matter distribution as well as the stellar orbital anisotropies. In particular, Li et al. (2016) investigated the validity of such techniques using the Illustris simulation. They found that although the orbital anisotropies cannot be accurately recovered and degeneracies exist between the dark matter and stellar components, the total mass distributions and their density slopes within are well recovered with 10% accuracy.

For galaxies at higher redshifts where only single-aperture kinematical data are available, simple approaches that use multiple mass measurements at different radii to make predictions about matter density slopes can be adopted. For example, in the SLACS (Auger et al. 2010b) and SL2S (Sonnenfeld et al. 2013, 2015) surveys, the central density slopes were derived for the observed lensing galaxies assuming spherical symmetry, power-law profiles, and isotropic orbital distributions (). The derived slopes could therefore suffer from more systematic biases than those using 2-D kinematics methods.

In order to make fair comparisons, one should apply the observational estimators also to the simulated samples. In this work, we adopted a simple approach along these lines that combines strong lensing and single-aperture kinematics for simulated early-type galaxies. In Sect. 6.1, we present the derived slopes, compare them with observational results, and discuss two associated major systematic effects. In Sect. 6.2, different slope estimators (as presented above) of the total matter density distributions are compared. In Sect. 6.3, we present the inner density slopes of the dark matter and stellar distributions of the early-type Illustris galaxies, and in Sect. 6.4, their cosmic evolution is discussed.

6.1 Total density slopes from combining strong lensing and single-aperture kinematics

For making fair comparisons of the total power-law slopes between the simulation and, in particular, the SLACS and SL2S survey results, we adopted a similar practice as used in these studies. Here we first briefly describe the main features of the method (details can be found in e.g., Koopmans et al. 2006; Auger et al. 2010b; Sonnenfeld et al. 2013).

The total matter distribution of a galaxy is assumed to be spherically symmetric, with a radial profile that follows a power law, i.e., . The stellar distribution is obtained by de-projecting the Sersic profile that best fits the surface brightness distribution (see Sect. 2.3). This latter component is assumed to be a massless tracer sitting in the gravitational potential of the former. The stellar orbital anisotropy , defined by Eq. (14), is assumed to be zero. Two “measurements” are made: (1) the mass projected within the Einstein radius (the strong lensing constraint); (2) the line-of-sight stellar velocity dispersion measured within a circular aperture of radius , same as used for the SLACS data (the stellar kinematics constraint). For any given power-law slope , measurement (1) constrains the normalization of the total matter density profile, with which the radial distribution of the stellar velocity dispersion can be derived by solving the spherical Jeans equation. The slope (where the superscript “LD” stands for “strong lensing and dynamics” and the subscript 0 refers to the assumption of ) that results in the best fit to measurement (2) is then taken as the power-law slope of the total density profile. We searched for within [1.2, 2.8] with a step of 0.02. This leads to the differences between the best-fitting and the “observed” velocity dispersion in most cases smaller than 1 km/s, and in all cases no larger than 2 km/s, much smaller than the observational uncertainty (2%-10%).

6.1.1 Comparisons to observations

Fig. 14 presents the total matter density slope as a function of (top panel), (middle panel), and (bottom panel). The solid and dashed lines indicate the median and the 90% boundaries of the distributions for the selected early-type galaxies between and . The blue and black squares with their error bars indicate measurements for the SLACS (Auger et al. 2010b) and SL2S (Sonnenfeld et al. 2013) galaxies, respectively.

The simulation reproduces the general observational trends and scatter. The dependencies are noticeable: on average decreases with increasing but increases with increasing and . In particular, higher-mass (higher-) galaxies on average have larger with smaller scatter in comparison to their lower-mass (lower-) counterparts. We note that, however, such a mass dependence is in part due to the correlation between the measurements of the velocity dispersion and the density slope, as well as other systematic biases (Sect. 6.1.2 and 6.1.3). In Sect. 6.2, mass dependences of different slope estimators are presented and discussed.

The measurements of for the simulated sample are overall shallower than observations. We note that although the selected Illustris galaxy sample shares a similar range of with the observed sample, they have rather different probability distributions. The former increases in number towards lower while the latter almost has a peak around . For this reason, we present quantitative comparisons made at given (within a small range): for galaxies with , (rms) for the Illustris galaxy sample, while (rms) for the combined SLACS and SL2S galaxy sample. For galaxies with , (rms) for the Illustris sample, and (rms) for the observation sample.

The notably shallower for the simulated galaxy sample is in fact consistent with lower central velocity dispersion and higher central dark matter fraction when comparing to the observational sample. This could indicate potentially inaccurate modelling of the involved baryonic physics processes in galaxy centres. We must note, however, the ad-hoc choice of a unique source redshift that was assigned to each lens redshift resulted in averagely smaller distribution of the simulated galaxy sample compared to the observation. As a mild increase of with increasing (i.e., steeper slopes at larger radii) was observed, such a systematic difference in would also cause a lower average of for the simulation sample.

Figure 14: The total matter density slope (see Sect. 6.2) as a function of (top panel), (middle panel), and (bottom panel). The solid and the dashed lines indicate the median and the 90% boundaries of the distributions from the selected early-type galaxies at . The blue and black squares show the observational results for the SLACS (Auger et al. 2010b) and SL2S (Sonnenfeld et al. 2013) samples, respectively, where the error bars indicate 1  error of the data.

6.1.2 Biases from the isotropic orbital assumption

Combining strong lensing and single-aperture stellar kinematics to derive the total density slopes, is commonly assumed due to a lack of sufficient observational constraints. To see the effect of a non-zero anisotropy parameter , we also calculated the total-density power-law slopes under the same assumptions as before, but using the true measured for the simulated galaxies. Fig. 15 shows the histograms of for the selected early-type galaxy sample at . Red and blue lines indicate galaxies that have and (measured within ), respectively. The former distribution (for radially anisotropic galaxies) peaks at , while the latter (for tangentially anisotropic galaxies) peaks at . This demonstrates that the true slopes of radially anisotropic systems tend to be overestimated assuming , while those of tangential ones tend to be underestimated (see also Koopmans 2006; 2009).

It is worth noting that, as shown in Fig. 7 and Fig. 9, markedly depends on stellar velocity dispersion and also evolves with redshift. This leads to possible systematic biases of the observational measurements. The total density slopes of higher- (lower-) mass galaxies or galaxies at lower- (higher-) redshifts are likely to be overestimated (underestimated) by when assuming .

Figure 15: Histograms of for the selected early-type galaxy sample at . Red and blue lines indicate galaxies that have and (measured within ), respectively.

6.1.3 Robustness of the power-law assumption

As shown by Xu et al. (2016), true profiles of realistic galaxies deviate from power-law distributions. In order to validate the robustness of the applied power-law assumption despite this fact, we first defined a curvature parameter for the 3-D enclosed mass distribution (defined by Eq. (18)) between and :

(19)

This curvature parameter directly quantifies the closeness of , and thus , to a power-law distribution between and . If the deviation of from unity is large, then the power-law approximation is poor. We set to be , which is the most relevant radial range for approaches that combine strong lensing and stellar kinematics.

We calculated for all the selected early-type galaxies from the simulation. Its dependence on is presented in the left-most panel of Fig. 16. The crosses indicate the simulation results at , the solid and the dashed lines indicate the median and the 90% boundaries of this distribution. In fact, both and show decreasing trends towards larger systems. In this case, the power-law approximation becomes worse for galaxies with central velocity dispersion . As shown by Xu et al. (2016), this leads to significantly biased determinations of the Hubble constant when power-law mass models are used to describe the lensing galaxies as they artificially break the so-called “mass sheet degeneracy” (e.g., Falco et al. 1985; Schneider & Sluse 2013).

The breakdown of the power-law approximation can also lead to biased estimates of other derived quantities, e.g., the 2-D enclosed mass distribution and the projected dark matter fraction , if they were derived assuming a power-law total density profile with slope (for , the dark matter mass is obtained by subtracting the observationally-constrained stellar mass from the total mass projected within a given aperture). In order to identify potential biases in these quantities, we further calculate two ratios: (1) between that is derived under the power-law assumption and the true mass distribution ; and (2) between that is derived under the power-law assumption and the true fraction . The “true” quantities are directly measured for simulated galaxies. Note that by construction, . However, this is not necessarily the case at other radii, unless the true distribution is indeed a power law.

We measured the two ratios within different aperture radii for the selected early-type Illustris galaxies. In Fig. 16, the middle and right-most panels show the enclosed mass ratio versus , and the projected dark matter fraction ratio versus , respectively (for the same galaxy sample at ). We note that the ratios deviate from unity due to a combination of two effects: (1) poor approximations of the power-law models for lower-mass galaxies; and (2) biased (power-law) slope estimates due to single-aperture kinematics data and the stellar isotropy assumption. These combined effects could lead to significantly biased estimates of the enclosed mass and the projected dark matter fraction, especially when contraints are made for lower-mass galaxies. The deviations from unity as well as the scatter also increase with redshift up to .

Figure 16: The curvature parameter (left panel), the enclosed mass ratio (middle panel) and the projected dark matter fraction (right panel), as a function of the central stellar velocity dispersion for the early-type galaxy sample at . In all three panels, the solid and the dashed lines indicate the median and the 90% boundaries of the distributions. The blue, red and orange curves present measurements made within different aperture sizes as specified in the panels. Note that the median values of the ratios measured within are close to 1.0, this is because the simulation sample has on average , and by construction . In addition, a small fraction of galaxies have negative (as can be seen from the right panel), which originates in the power-law assumption. For these galaxies, the power-law model that best fits both the strong lensing and kinematical constraints predicts total masses smaller than the stellar masses within the studied aperture.

6.2 Different estimators of the total density slopes

To investigate the variations among different total-slope estimators, we calculated several sets of slopes for the simulated galaxies within different radial ranges, namely, the simple power-law dynamical slopes and , the mass-weighted slopes and , the power-law fitted slopes and , and finally the average slopes and .

slope estimator    mean    standard deviation 1.83 0.24 1.80 0.23 1.95 0.14 1.99 0.14 2.08 0.27 2.07 0.26 2.08 0.26 2.06 0.22 1.49 0.22 1.58 0.17 2.74 0.30 2.87 0.26
Table 1: Stellar mass-weighted mean and standard deviation of the density slopes for the selected early-type galaxy sample at

Table 1 presents the stellar mass-weighted mean and standard deviation of each above-mentioned slope estimator for the selected early-type galaxy sample at . In Fig. 17, we plot, in particular, four typical slopes measured around and within , for the same galaxy sample, as a function of . The red, orange, green and blue curves indicate the distributions of , , and , respectively. The solid and the dashed line styles show the median and the 68% range of each distribution, respectively. The thicker dashed and dotted lines present the best linear fits to the data.

We see that within a large fraction of the velocity range investigated here, the median slopes are all about isothermal. In particular, the median and scatter of the intrinsic density slope estimators and are consistent with the observational results from 2-D kinematical data (e.g., Barnabè et al. 2011; Cappellari et al. 2015).

Interestingly, systematic discrepancies are seen among different slope estimators. For most galaxies, since the local slope (defined by Eq. (17)) tends to decrease, i.e., the profile turning flatter, towards smaller (see also Dutton & Treu 2014 for more discussion), one finds as expected. Between the two simple power-law dynamical slopes, the fact that on average (within a large span of ) can be explained by the stellar orbital anisotropies, because the majority of the galaxy sample has radial anisotropies (as can be seen from Fig. 7). For galaxies with , a moderate disagreement exists between the simple dynamical slope and the intrins