# Revisiting the Bulge–Halo Conspiracy. I. Dependence on Galaxy Properties and Halo Mass

Francesco Shankar, Alessandro Sonnenfeld, Gary A. Mamon, Kyu-Hyun Chae, Raphael Gavazzi, Tommaso Treu, Benedikt Diemer, Carlo Nipoti, Stewart Buchan, Mariangela Bernardi, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Francesco Shankar, Alessandro Sonnenfeld, Gary A. Mamon, Kyu-Hyun Chae, Raphael Gavazzi, et al.. Revisiting the Bulge–Halo Conspiracy. I. Dependence on Galaxy Properties and Halo Mass. The Astrophysical Journal, 2017, 840 (1), pp.34. 10.3847/1538-4357/aa66ce. hal-01554680

HAL Id: hal-01554680

https://hal.science/hal-01554680v1

Submitted on 10 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Revisiting the Bulge–Halo Conspiracy. I. Dependence on Galaxy Properties and Halo Mass

Francesco Shankar<fnref n="1" />, Alessandro Sonnenfeld<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, Gary A. Mamon<fnref n="5" />, Kyu-Hyun Chae<fnref n="6" />, Raphael Gavazzi<fnref n="5" />, Tommaso Treu<fnref n="2" />,<fnref n="3" />,
Benedikt Diemer<fnref n="7" />, Carlo Nipoti<fnref n="8" />, Stewart Buchan<fnref n="1" />, Mariangela Bernardi<fnref n="9" />, Ravi Sheth<fnref n="9" />, and Marc Huertas-Company<fnref n="10" />,<fnref n="11" />

*Received 2016 September 9; revised 2017 February 28; accepted 2017 March 7; published 2017 May 2*

## Abstract

We carry out a systematic investigation of the total mass density profile of massive (log $M_{\star}/M_{\odot} \gtrsim 11.3$) early-type galaxies and its dependence on galactic properties and host halo mass with the aid of a variety of lensing/dynamical data and large mock galaxy catalogs. The latter are produced via semi-empirical models that, by design, are based on just a few basic input assumptions. Galaxies with measured stellar masses, effective radii, and Sérsic indices, are assigned, via abundance matching relations, host dark matter halos characterized by a typical ΛCDM profile. Our main results are as follows. (1) In line with observational evidence, our semi-empirical models naturally predict that the total, mass-weighted density slope at the effective radius $\gamma'$ is not universal, steepening for more compact and/or massive galaxies, but flattening with increasing host halo mass. (2) Models characterized by a Salpeter or variable initial mass function (IMF) and uncontracted dark matter profiles are in good agreement with the data, while a Chabrier IMF and/or adiabatic contractions/expansions of the dark matter halos are highly disfavored. (3) Currently available data on the mass density profiles of very massive galaxies (log $M_{\star}/M_{\odot} \gtrsim 12$), with $M_{\text{halo}} \gtrsim 3 \times 10^{14} M_{\odot}$, favor instead models with a stellar profile flatter than a Sérsic one in the very inner regions ($r \lesssim 3–5$ kpc), and a cored NFW or Einasto dark matter profile with median halo concentration a factor of $\sim 2$ or $\lesssim 1.3$, respectively, higher than those typically predicted by N-body numerical simulations.

**Key words:** cosmology: theory – galaxies: evolution – galaxies: statistics

## 1. Introduction

Early-type galaxies constitute a family of objects of remarkable regularity, captured by a series of tight scaling relations such as the fundamental plane (e.g., Djorgovski & Davis 1987; Dressler et al. 1987), the color–magnitude relation (e.g., Bower et al. 1992; Mei et al. 2012) and the relation between the mass of the central black hole and the global properties of the host galaxy (e.g., Ferrarese & Merritt 2000; Gebhardt et al. 2000; Shankar 2009; Shankar et al. 2016). This regularity is reflected in their mass structure: the total mass density profile of massive<fnref n="12" /> galaxies is well approximated, around the scale of the half-light radius, by a power-law $\rho(r) \propto r^{-\gamma}$ with slope close to isothermal ($\gamma \approx 2$) and little scatter across the population (e.g., Koopmans et al. 2006, 2009; Gavazzi et al. 2007; Barnabè et al. 2011). This compelling observation is named the bulge–halo conspiracy: stars and dark matter have separate density profiles that are significantly different from an isothermal profile, yet they conspire to produce a total density profile close to isothermal (e.g., Treu & Koopmans 2004).

Recent studies based on strong lensing and stellar dynamics have shown how the density profile of massive early-type galaxies is not exactly universal: the density slope $\gamma'$ correlates

with projected stellar mass density, being steeper for more compact objects, and anti-correlates with redshift, being shallower for higher redshift objects with respect to local systems of the same mass and size (e.g., Auger et al. 2010a; Ruff et al. 2011; Bolton et al. 2012; Sonnenfeld et al. 2013; Dye et al. 2014; Tortora et al. 2014a, 2014b). The slope of the density profile appears to also correlate with halo mass: Brightest Cluster Galaxies (BCGs) have on average shallower density profiles than field galaxies of similar stellar mass (e.g., Newman et al. 2013a, 2013b). The residual scatter left once these scaling relations are taken into account is as small as 6% on the mass density slope (Sonnenfeld et al. 2013), highlighting once more the high degree of self-similarity across the population of early-type galaxies.

The origin of such a regularity is not fully understood. If massive, central early-type galaxies are predominantly grown by mergers (Naab et al. 2009; Shankar et al. 2010a, 2010b, 2015; Bai et al. 2014; Bellstedt et al. 2016; Buchan & Shankar 2016), stochastic processes may tend to alter pre-existing scaling relations, or at least increase the scatter in the population (e.g., Nipoti et al. 2009a, 2009b; Shankar et al. 2013), unless additional mechanisms are present such as gas dissipation (e.g., Robertson et al. 2006; Oser et al. 2010; Remus et al. 2013, 2017; Sonnenfeld et al. 2014). More recent numerical and semi-analytic studies suggest that even collisionless mergers could create tight correlations close to those

<footnote n="1">¹Department of Physics and Astronomy, University of Southampton, Highfield, SO17 1BJ, UK; F.Shankar@soton.ac.uk</footnote>

<footnote n="2">²Physics Department, University of California, Santa Barbara, CA 93106, USA</footnote>

<footnote n="3">³Physics and Astronomy Department University of California Los Angeles CA 90095-1547, USA</footnote>

<footnote n="4">⁴Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="5">⁵Institut d'Astrophysique de Paris (UMR 7095: CNRS & UPMC, Sorbonne Universités), 98bis Bd Arago, F-75014 Paris, France</footnote>

<footnote n="6">⁶Department of Physics and Astronomy, Sejong University, 209 Neungdong-ro Gwangjin-Gu, Seoul 05006, Korea</footnote>

<footnote n="7">⁷Institute for Theory and Computation, Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge MA 02138, USA</footnote>

<footnote n="8">⁸Department of Physics and Astronomy, Bologna University, via Berti-Pichat 6/2, I-40127 Bologna, Italy</footnote>

<footnote n="9">⁹Department of Physics and Astronomy, University of Pennsylvania, 209 South 33rd St, Philadelphia, PA 19104, USA</footnote>

<footnote n="10">¹⁰GEPI, Observatoire de Paris, CNRS, Université Paris Diderot, 61, Avenue de l'Observatoire F-75014, Paris France</footnote>

<footnote n="11">¹¹Université Paris Denis Diderot, F-75205 Paris Cedex 13, France</footnote>

<footnote n="12">¹²This is not necessarily the case if one considers a broader mass range (e.g., Shu et al. 2015).</footnote>
observed (e.g., Taranu et al. 2013, 2015; Shankar et al. 2014b). Remus et al. (2013) specifically investigated the evolution of the total mass density slope in a set of galaxies in a cosmological simulation, and were able to naturally produce slopes close to isothermal consistent with observations. Despite recent progress (e.g., Dubois et al. 2013; Remus et al. 2013, and references therein), however, the still unclear sub-grid physics and large computational power required to run high-resolution cosmological simulations do not allow for a proper quantitative comparison between detailed models and observations.

A powerful, complementary approach to cosmological simulations and semi-analytic models is to use semi-empirical models (e.g., Dutton & Treu 2014). The latter, by design, are an extremely effective way of making use of the known properties of galaxies together with minimal theoretical inputs to make testable predictions on a set of observables, and set unique, independent constraints on galaxy evolution processes.

In this work, we use the observed correlations between photometric properties of massive galaxies (i.e., effective radius, stellar mass), together with basic assumptions on their host dark matter halos, as inputs to create large mock galaxy catalogs. We then test the validity of the (few) input assumptions and parameters, by comparison with data at different environments spanning from the field to the cluster halo mass scale. In the following papers of this series, we will specifically include velocity dispersion in our mocks as an additional probe for models, and study the evolution of the mass density profile with redshift.

The structure of this work is the following. In Section 2, we describe the method used to construct mock samples of galaxies and relative measurements of the density slope and dark matter fraction. In Section 3, we present the predictions for the dependence of the total mass density profile on galactic properties, and environment in Section 4. In Section 5, we discuss our results and conclude in Section 6.

In the following, we will adopt the reference cosmology with parameters $\Omega_m = 0.30$, $\Omega_b = 0.045$, $h = 0.70$, $\Omega_\Lambda = 0.70$, $n = 1$, and $\sigma_8 = 0.8$, to match those usually adopted in the observational studies on the stellar mass function and lensing considered in this work. We will also by default adopt a Salpeter (Salpeter 1955) as our reference Initial Mass Function (IMF), though we will carefully discuss the implications on our model predictions when switching to a Chabrier (Chabrier 2003) or even variable (Cappellari et al. 2012, 2015) IMF.

## 2. Method

### 2.1. The Semi-empirical Model

The aim of this work is to compare semi-empirical models based on $\Lambda$CDM halo relations coupled to observed galaxy scaling relations, to mass modeling analysis principally from strong lensing. Our semi-empirical approach makes use of mostly observational quantities with only a few, basic theoretical inputs. In this respect, it is extremely powerful because it does not rely on any physical assumption for evolving galaxies, e.g., via mergers and/or in situ star formation.

Our procedure to build mock galaxies relies on the following steps.

1. We take large catalogs of central, early-type galaxies with measured stellar masses, projected effective radii, and Sérsic indices (Sérsic 1963).

2. To each galaxy, we assign a host dark matter halo via abundance matching relations.

3. We associate to each central galaxy a 3D de Vaucouleurs (de Vaucouleurs 1948) or Sérsic stellar profile (Sérsic 1963), depending on which data set we compare to. Each profile depends on its specific projected effective radius and Sérsic index (equal to four for a de Vaucouleurs profile), as detailed in Appendix B (e.g., Lima Neto et al. 1999).

4. A Navarro et al. (1997, NFW hereafter) or Einasto (Einasto 1965) profile characterizes instead each host dark matter halo.

5. For each galaxy, we then build the scale-dependent, total mass density profile.

Although we will adopt the above as our reference model, we will also extensively discuss the impact on our predictions induced by relaxing one or more of the previous assumptions. We will consider changes in both the dark matter and stellar profiles, in the IMF, in the stellar mass-halo mass mapping, and also in dark matter concentration.

In detail, for our mocks, we extract galaxies from the Meert et al. (2015) sample derived from the Sloan Digital Sky Survey (SDSS) DR7 spectroscopic sample (Abazajian et al. 2009) in the redshift range of $0.05 < z < 0.2$, and with a probability $P(E) > 0.85$ of being elliptical galaxies based on the Bayesian automated morphological classifier by Huertas-Company et al. (2011). The vast majority of the early-type galaxies in this sample are labeled as “central” galaxies according to the Yang et al. (2007) host halo catalog. Stellar masses are obtained from integrated light profiles multiplied by the color-dependent mass-to-light ratios of Bell et al. (2003) renormalized to a Chabrier IMF. In the following, depending on the data set to compare with, we will consider three different light profiles, a de Vaucouleurs+Exponential (the “cmodel,” see Bernardi et al. 2010), a Sérsic+Exponential, and pure Sérsic profiles (see, e.g., Bernardi et al. 2013). Making use of directly measured physical quantities to build the mock has the advantage that no assumptions have to be made about, e.g., intrinsic correlations and scatters (e.g., Gaussian) among galaxy sizes, stellar masses, and Sérsic indices. It is important to stress that SDSS is characterized by an average seeing of ~1.5 arcsec, which would correspond to projected scales of $R \sim 3$ kpc at the average redshift of $z \sim 0.1$. This observational limit is small enough not to bias the measurement of the stellar profile on scales of $R \gtrsim 0.2 R_\odot$ for galaxies $M_{\star} \gtrsim 4 \times 10^{11} M_\odot$, though it mostly prevents the possible detection of a stellar “core” or flattening below this scale. In the following sections, we will discuss if extrapolations of a power-law Sérsic profile to very small scales is supported by strong lensing and kinematic data.

On the dark matter side, our reference will be the NFW profile with mass density

$$ \rho(r) = \frac{\rho_0}{(r/r_s)(1 + r/r_s)^2}, \quad (1) $$

where $r_s$ is the scale radius, and $\rho_0$ is a reference density. Where relevant, we will also consider other analytic forms for the dark matter profile, namely a “cored” NFW profile (Newman et al. 2013a)

$$ \rho(r) = \frac{b\rho_0}{(1 + br/r_s)(1 + r/r_s)^2}, \quad (2) $$

with a reference value of $r_{\text{core}} = r_s/b = 14$ kpc as the one adopted by Newman et al. (2013a), and an Einasto
profile, expressed as (see, e.g., Mamon & Łokas 2005a, their Appendix A)

$$ \rho(r) = \rho_{-2} \exp(2\mu) \exp\left[-2\mu\left(\frac{r}{r_{-2}}\right)^{1/\mu}\right]. \quad (3) $$

In Equation (3), $\rho_{-2}$ is the local mass density at $r_{-2}$, the radius at which the logarithmic slope of the density profile is equal to -2, which is equal to $r_s$ for an NFW profile. We will further discuss below the values of $\mu$ adopted in this work.

Following several previous attempts in the literature (e.g., Dutton & Treu 2014; Oguri et al. 2014), we also consider NFW profiles modified by the inclusion of adiabatic contraction and/or expansion. As in Dutton et al. (2007), assuming spherical collapse without shell crossing, we express the actual relation between the final and initial radius of the mass distribution as $r_f = \Gamma^\nu r_i$. The parameter $\Gamma = r_f/r_i$ is the contraction factor (Blumenthal et al. 1986), which can be numerically calculated via the conservation equation

$$ M_i(r_i) r_i = M_f(r_f) r_f, \quad (4) $$

with $M_i$ the initial (baryonic plus dark matter) mass distribution, and

$$ M_f(r_f) = M_{\text{star}}(r_f) + (1 - f_{\text{gal}}) M_i(r_i) \quad (5) $$

the final one, with $f_{\text{gal}}$ the ratio between the galaxy stellar (possibly plus gas) mass and the host halo mass (see, e.g., Barausse 2012). The parameter $\nu$ is set to zero for no profile modification, while it assumes positive or negative values for contraction or expansion, respectively.<fnref n="13" /> It has already been shown that extreme contraction ($\nu \gtrsim 1$) and/or expansion ($\nu \lesssim -0.5$ tend to be disfavored by present data (Dutton & Treu 2014). We will therefore limit our analysis to non-maximal effects, such as the one by Gnedin et al. (2004) for contraction ($\nu = 0.8$), and a slightly milder expansion with $\nu = -0.3$.

Once a total mass density profile for a galaxy has been constructed, one direct quantity usually adopted in the literature to compare with lensing data has been $\gamma$, i.e., the local logarithmic slope of the spherically averaged total density profile $\rho(r) \propto r^{-\gamma(r)}$. In more recent years, the quantity $\gamma'$ has been more often adopted. The latter is defined as the mass-weighted slope of the total density profile within a radius $r$, and it is computed as (Dutton & Treu 2014)

$$ \begin{aligned} \gamma'(r) &= -\frac{1}{M(<r)} \int_0^r \frac{d \log \rho}{d \log x} 4\pi x^2 \rho(x) dx \\ &= -\frac{1}{M(<r)} \left( \rho(r) 4\pi r^3 - 3 \int_0^r \rho(x) 4\pi x^2 dx \right) \\ &= 3 - \left. \frac{d \log M}{d \log x} \right|_{x=r} . \end{aligned} \quad (6) $$

The mass-weighted slope $\gamma'(r)$, specifically computed within the effective radius $R_e$, has been shown to well approximate the slope measured in joint strong lensing and stellar kinematics studies (Sonnenfeld et al. 2013). For power-law density profiles, $\gamma' = \gamma$. In the following, unless otherwise noted, we

will always refer to the mass-weighted slope $\gamma'$ computed at $R_e$. Sonnenfeld et al. (2013, see their Figure 6) and Dutton & Treu (2014) have shown that the mass-weighted $\gamma'$ computed as in Equation (6) is a good proxy of $\gamma'$ derived from strong lensing and dynamical measurements (see Sonnenfeld et al. 2013 for a full discussion).

## 2.2. The Stellar Mass–Halo Mass and Effective Radius–Stellar Mass Relations

In this work, we focus on central massive galaxies only, and thus accordingly adopt the mean central stellar mass–halo mass relation to map galaxies into halos. The complexities behind establishing a secure mapping for satellites have been under intense study for several years (e.g., Neistein et al. 2011; Rodríguez-Puebla et al. 2012; Watson & Conroy 2013), and are not vital to the present study focused on massive galaxies with $\log M_{\text{star}}/M_{\odot} \gtrsim 11.5$, for which the fraction of satellites is very low, at the percent level (e.g., Guo et al. 2014).

The empirical stellar mass–halo mass correlations adopted as inputs of our empirical models are derived from abundance matching between the stellar mass and (sub)halo mass functions

$$ \Phi(>M_{\text{star}}, z) = \Phi_c(>M_{\text{halo}}, z) + \Phi_s(>M_{\text{halo}}, z) \quad (7) $$

with $M_{\text{halo}} = M_{200c}$ the halo masses defined within $r_{200c}$, such that the average density within $r_{200c}$ is 200 times the critical density of the universe at redshift $z$. The $\Phi_c(>M_{\text{halo}}, z)$ term refers to the host halo mass function of Tinker et al. (2008), with the subhalo term $\Phi_s(>M_{\text{halo}}, z)$ with mass at infall $M_{\text{halo}}$ taken from Behroozi et al. (2013).

For the concentration $c_{200c} = r_{200c}/r_s$, we use the median $c_{200c}-M_{\text{halo}}$ relation as derived by Diemer & Kravtsov (2015) for our cosmology, and include a log-normal scatter of 0.16 dex. As we discuss in Appendix A, the Diemer & Kravtsov (2015) fit is, at least at $z = 0$, in very good agreement with the Bullock et al. (2001) model, as revised by Macciò et al. (2008) and Dutton & Macciò (2014). We find that the exact mass-dependent normalization and scatter of the concentration–mass relation (e.g., Dutton & Treu 2014), have a relatively small impact on most of our conclusions, except possibly at cluster scales, as discussed below. Note that we assume the scatter in halo concentration to be uncorrelated with stellar mass or galaxy size (see discussions in, e.g., Katz et al. 2017; Papastergis & Shankar 2016). We have checked that at fixed galaxy stellar mass, size, and host halo mass a systematic variation of 0.16 dex in halo concentration induces a relatively small variation of ~4% in the implied $\gamma'$ computed at the effective radius.

The stellar mass function on the left side of Equation (7) depends on the assumed light profile. We adopt the results by Bernardi et al. (2013), who provide stellar mass functions for all SDSS galaxies characterized by cmodel, Sérsic+Exponential, and pure Sérsic magnitudes (all based on a Chabrier IMF). As recently discussed by a number of groups, adopting Sérsic profiles naturally results in more integrated light than de Vaucouleurs profiles, and thus, at fixed mass-to-light ratio, larger abundances of massive galaxies (see, e.g., Bernardi et al. 2013, 2016a, 2016b, 2017; D'Souza et al. 2015; Thanjavur et al. 2016). This in turn yields steeper stellar-to-halo mappings, i.e., more massive mean stellar masses at fixed halo mass, in better agreement with direct

<footnote n="13">¹³ Dutton & Treu (2014, and references therein) discuss that $\nu$ may also be correlated with the shift $\Delta \log M_{\text{star}}$ from a Chabrier IMF. This effect is rather small, and we ignore it in the present study, though we discuss the effects of changing IMF, host halo masses, etc.</footnote>
measurements of massive BCG (see, e.g., Kravtsov et al. 2014; Shankar et al. 2014a). For each light profile we explore the impact of varying the input IMF, from a Chabrier to a Salpeter and variable IMF. For the former, we simply add 0.25 dex to each stellar mass (Bernardi et al. 2010). For the latter, we take the velocity dispersion-dependent dynamical/ strong lensing mass-to-light ratios from Cappellari (2016, their Figure 19), expressed as a difference with respect to the stellar (Salpeter) mass-to-light ratio

$$ \log(M/L)_{\text{dyn}} - \log(M/L)_{\text{Salp}} = -0.0576 + 0.364 \times \log\left(\frac{\sigma_e}{200 \text{ km s}^{-1}}\right) \quad (8) $$

with an intrinsic scatter of 0.11 dex.

The resulting mapping between stellar mass and host halo mass<fnref n="14" /> is well reproduced by the following two-power-law relation

$$ M_{\text{star}} = M_{\text{star}}^{0} \left( \frac{M_{200c}}{M_{200c}^{0}} \right)^{\alpha} \left[ 1 + \left( \frac{M_{200c}}{M_{200c}^{0}} \right)^{\beta} \right]^{-1}, \quad (9) $$

with $M_{200c} = M_{\text{halo}}$, valid in the approximate range of $10^{10} \lesssim M_{\text{star}}/M_{\odot} \lesssim 10^{12}$ (Chabrier IMF).

We always assume an intrinsic Gaussian scatter of 0.15 dex in $\log M_{\text{star}}$ at fixed halo mass, which is a valid approximation especially at high stellar masses and at low redshifts, in agreement with a number of diverse observational and theoretical studies (e.g., Guo et al. 2014; Shankar et al. 2014a; Gu et al. 2016; Tinker et al. 2016, and references therein). To include scatter, we follow the methodology outlined in Shankar et al. (2014a). We first fit Equation (9) to the results of abundance matching without scatter, and then vary the parameter $\beta$ until the implied mock galaxy–halo samples (inclusive of scatter) reproduce well the input stellar mass function within the observational uncertainties. The final values of the parameters of Equation (9) are given in Table 1 for a Chabrier, a Salpeter, and a variable IMF, and for each IMF we consider three different light profiles: a de Vaucouleurs, a Sérsic+Exponential, and a Sérsic profile. For each combination of light profile and IMF, the implied mean stellar mass–halo mass relation is calculated by inserting in Equation (7) the relevant light profile-dependant stellar mass function (Bernardi et al. 2013) corrected for an appropriate mass-to-light ratio.

In turn, when assigning a host halo mass to a galaxy in our SDSS sample with measured stellar mass, we make use of the mean halo mass at fixed stellar mass relation, i.e., the “inverse” of Equation (9). The latter mapping is computed by first creating a large halo mock catalog extracted from the halo mass function, assigning galaxies directly applying Equation (9), with a scatter of 0.15 dex and parameters given in Table 1, and then self-consistently computing the implied mean and scatter in halo mass as a function of stellar mass. As a consistency check, we make sure that our inverse relations and related scatters properly reproduce the halo mass function from random galaxy catalogs extracted from the input stellar mass functions.

<table><caption>**Table 1**
Stellar Mass–Halo Mass Relations (Equation (9))</caption><thead><tr><th>IMF</th><th>log M<sup>0</sup><sub>star</sub></th><th>log M<sup>0</sup><sub>200c</sub></th><th>α</th><th>β</th></tr></thead><tbody><tr><td>Chabrier-SerExp</td><td>10.68</td><td>11.80</td><td>2.13</td><td>1.68</td></tr><tr><td>Chabrier-Sersic</td><td>10.71</td><td>11.82</td><td>2.13</td><td>1.66</td></tr><tr><td>Chabrier-deVac</td><td>10.63</td><td>11.80</td><td>2.17</td><td>1.77</td></tr><tr><td>Salpeter-SerExp</td><td>10.98</td><td>11.84</td><td>2.14</td><td>1.72</td></tr><tr><td>Salpeter-Sersic</td><td>11.03</td><td>11.87</td><td>2.09</td><td>1.66</td></tr><tr><td>Salpeter-deVac</td><td>10.93</td><td>11.84</td><td>2.15</td><td>1.77</td></tr><tr><td>varIMF-SerExp</td><td>10.93</td><td>11.83</td><td>2.03</td><td>1.60</td></tr><tr><td>varIMF-Sersic</td><td>10.96</td><td>11.85</td><td>2.05</td><td>1.60</td></tr><tr><td>varIMF-deVac</td><td>10.83</td><td>11.81</td><td>2.04</td><td>1.64</td></tr></tbody></table>

**Note.** Parameters of the relations at $z = 0.1$ between the central stellar mass and host halo mass (200 times the critical density) from top to bottom, for a Chabrier, Salpeter, and variable IMF, and for three different light profiles, a Sérsic-Exponential, a Sérsic, and a de Vaucouleurs profile. A constant intrinsic scatter of 0.15 dex in stellar mass at fixed halo mass is assumed in all models. For each model, the most appropriate stellar mass function is adopted in the abundance matching routine (Equation (7)).

In Figure 1, we show the mean stellar mass–halo mass (top), effective radius–stellar mass (middle), and velocity dispersion–stellar mass (bottom) relations assuming a Salpeter (left) or variable IMF (right). As labeled, each panel includes three models corresponding to three different light profiles, de Vaucouleurs (cmodel, black, solid lines with their $1\sigma$ dispersions marked by gray areas), Sérsic–Exponential (blue, dotted–dashed lines), and Sérsic (red, long-dashed lines).

It is clear from the top panels that all of our models, and in particular the ones based on Sérsic profiles, predict mean stellar mass–halo mass relations with a high-mass end slope of $\delta = \alpha - \beta \sim 0.45$ (see Table 1), steeper than the Moster et al. (2013; purple, long-dashed line) or Behroozi et al. (2013) relations, characterized by slopes of $\delta \sim 0.30–0.35$, which would imply up to a factor of $\sim 3$ lower stellar masses at fixed halo mass (see, e.g., Kravtsov et al. 2014; Shankar et al. 2014a; Buchan & Shankar 2016). All of our models are in broad agreement with the Sonnenfeld et al. (2013) data (blue squares), with the Sérsic-based stellar mass–halo mass relations particularly well aligned with the Newman et al. (2013b, 2015) data (green triangles and red circles, respectively), which are also characterized by Sérsic-type light profiles. All in all, the variable IMF models are characterized by scaling relations that are very similar to the ones extracted from pure Salpeter models. At best, they might predict stellar masses slightly lower by $\sim 0.1$ dex at low halo masses. Such a high degree of similarity is expected in the range of high stellar masses considered in this work.

Our models’ results are significantly steeper than the Moster et al. (2013) relation mainly because ours are based on the Bernardi et al. (2013) stellar mass functions, which predict up to two orders of magnitude more massive galaxies than the Li & White (2009) stellar mass function, adopted as the reference one in the Moster et al. (2013) analysis. The Li & White (2009) stellar mass function adopted Petrosian magnitudes (Petrosian 1976), which have long been known to underestimate the light in the most massive galaxies (e.g., Bernardi et al. 2010). Moreover, the Petrosian magnitudes used by Li & White (2009) were based on the SDSS pipeline, which suffers from bad sky estimates, particularly in crowded fields where massive galaxies are typically found (see, e.g., Bernardi et al. 2010, 2014, and

<footnote n="14">¹⁴ Note that Equation (9) is the mean relation for the entire galaxy population. Dutton et al. (2010), however, showed that the stellar mass–halo mass relation of early-type galaxies is practically indistinguishable from the one characterizing the total population for $M_{\text{star}} \gtrsim 2 \times 10^{11} M_{\odot}$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Main galaxy scaling relations adopted in the semi-empirical modeling assuming a Salpeter (left panels) or variable (right panels) IMF. Top: stellar mass vs. halo mass (defined as 200 times the critical density) relation derived from abundance matching between the stellar mass and host halo mass functions assuming an intrinsic scatter of 0.15 dex in log M<sub>star</sub> at fixed halo mass. The black solid, blue dotted-dashed, and red long-dashed lines refer, respectively, to the relations inferred from the Bernardi et al. (2013) cmodel, Sérsic-exponential, and Sérsic stellar mass functions as described in the text (see also Shankar et al. 2014a), while the purple long-dashed line is the relation derived by Moster et al. (2013). Middle: two-dimensional, circularized effective radius R<sub>e</sub> vs. stellar mass as derived by Bernardi et al. (2014) from cmodel, Sérsic-exponential, and Sérsic light profiles, as labeled. Bottom: average line-of-sight velocity dispersion within a circular aperture of radius equal to half of the effective radius σ<sub>e/2</sub> as a function of stellar mass for the cmodel, Sérsic-exponential, and Sérsic light profiles. In all the left panels, we compare, where possible, with the Salpeter-based data by Sonnenfeld et al. (2013), Newman et al. (2013a, 2015; blue squares, green triangles, and red circles, respectively). Gray shaded regions in all panels mark the 1σ dispersions around the cmodel-de Vaucouleurs curves.</figcaption></figure>

references therein). More recent work by D'Souza et al. (2015)
and Bernardi et al. (2016a) confirms that the Li & White (2009)
photometry was indeed biased faint. Finally, Li & White (2009)
adopted stellar mass-to-light ratios based on templates from
Blanton & Roweis (2007), which are inappropriate for massive

galaxies (e.g., Bernardi et al. 2013, 2016a, 2016b, 2017). Bernardi et al. (2017) have also carefully demonstrated that the difference with the SDSS pipeline photometry does not depend on whether or not one includes intracluster light, but rather it is truly missing light from the body of the galaxy.
The middle panels of Figure 1 report instead the two-dimensional effective radius $R_e$ versus stellar mass for the Salpeter (left) and variable IMF (right), and for our three reference stellar light profiles, as labeled. It is clear that Sérsic-based profiles provide significantly more extended light profiles at fixed stellar mass, and thus proportionally larger effective radii up to a factor of ≳2–3 higher than what is inferred from pure de Vaucouleurs fits. Within uncertainties, our SDSS-cmodel data compare reasonably well with the effective radii and stellar masses of the galaxies from the Sonnenfeld et al. (2013),<fnref n="15" /> which are based on de Vaucouleurs profiles. Even closer is the match between our Sérsic-based size–mass relations with the data by Newman et al. (2013b, 2015), whose half-light radii and total luminosities are extracted from pseudo-isothermal and Sérsic profiles, respectively. Note that relatively tight correlations such as those between stellar mass and halo mass, and size versus stellar mass, should naturally induce additional dependencies between galaxies and their dark matter halos. We find, for example, that as a consequence of the relation between galaxy size and stellar mass, and the assumption of a monotonic stellar mass–halo mass relation, our models produce a roughly linear correlation between effective radius and virial radius, in broad agreement with the one independently measured by Kravtsov (2013).

For completeness, the bottom panels of Figure 1 plot the line-of-sight velocity dispersion within half of the effective radius<fnref n="16" /> as a function of stellar mass, in the same format as for the other panels, for Salpeter (left) and variable (right) IMF and for the three reference light profiles. In this case, we find the velocity dispersions in the Sonnenfeld et al. (2013) data to be slightly higher by ~0.05 dex, at fixed stellar mass or, alternatively, stellar masses to be lower by ~0.1 dex at fixed velocity dispersion. Irrespective of this, it is indeed quite remarkable that, despite possible differences and biases in the selections, and in the surface brightness profile fitting and/or background subtraction algorithms, the data by Sonnenfeld et al. (2013) and Newman et al. (2013b, 2015) show with respect to our SDSS catalogs offsets of at most ≲0.1 dex in effective radius, stellar mass or velocity dispersion.

It could be argued from the apparent (systematic) discrepancy in the (mean) velocity dispersion–stellar mass relations (lower, bottom panel of Figure 1), that the reference sample of lensed galaxies is biased high in velocity dispersion at fixed stellar mass with respect to SDSS galaxies. We disfavor a substantial bias in the SLACS galaxies as the galaxies do not appear proportionally more compact at fixed stellar mass (middle, left panel in Figure 1). Indeed, this (relatively small) discrepancy could be mostly ascribed to residual differences between the cmodel and de Vaucouleurs light profiles, and/or between mass-to-light ratios (see, e.g., Bernardi et al. 2013, 2017, and references therein). Auger et al. (2010a) claimed a closer match between SLACS galaxies and SDSS; however, their comparison was based on the velocity dispersion–stellar mass relation calibrated by Hyde & Bernardi (2009), who adopted mass-to-light ratios systematically slightly lower than the ones adopted by Bernardi et al. (2013). Treu et al. (2006)

also did not find any difference between SLACS and SDSS galaxies at fixed velocity dispersion. More relevantly to what follows, we verified that all of our main results are robust against the level of mild systematics seen in the left panels of Figure 1. For instance, we checked that reducing stellar masses in our SDSS mocks by ~0.1 dex does not affect any of our conclusions. Even when selecting only the SDSS galaxies above the mean cmodel velocity dispersion–stellar mass relation to better match the SLACS sample, also yields similar results. The latter behavior is not unexpected given that, on average, we verified that massive galaxies in all our models share a very weak dependence of $\gamma'$ on velocity dispersion (see also Xu et al. 2016; Poci et al. 2017).

In the following, we will adopt the scaling relations discussed in this section as inputs for our mock galaxy–halo catalogs. For self-consistency reasons, we will specifically adopt the cmodel and single-Sérsic scaling relations when comparing with the Sonnenfeld et al. (2013) and Newman et al. (2013b, 2015) data, respectively. Unless otherwise stated, from now on all our stellar masses $M_{\text{star}}$ will refer to Salpeter stellar masses, to make contact with the observations considered in this work, though the input IMF may not necessarily be a Salpeter IMF. In practice, when computing $\gamma'$ in some models, we will assume a different input IMF, Chabrier or variable IMF, but when comparing to the data we will convert our mock stellar masses into the values an observer would measure assuming a Salpeter IMF. This, we remind the reader, corresponds to a positive offset of 0.25 dex (e.g., Bernardi et al. 2010) in stellar mass in case the input IMF is a Chabrier IMF, or a velocity-dependent offset as given in Equation (8) for a variable IMF. All other photometric properties, namely effective radii and Sérsic indices, are instead clearly independent of the choice of input IMF but mainly rely on the choice of assumed light profile.

### 3. Dependence of $\gamma'$ on Galactic Properties

In this section, we compare the predictions of our models to observational data at $z < 0.3$. For simplicity, we start by analyzing the average predictions of the semi-empirical models in Section 3.1, and in Section 3.2 we proceed by analyzing further dependencies on galaxy properties. In this work, we will not discuss the dependence of $\gamma'$ on velocity dispersion, which shares with $\gamma'$ a strong correlation in the observational errors (e.g., Sonnenfeld et al. 2013). We will instead mainly focus on dependencies of $\gamma'$ on independently measured quantities such as effective radius and stellar mass. We will get back to velocity dispersions when modeling the velocity dispersion–stellar mass relation in a separate work. In Section 4, we will discuss environmental trends as quantified by the host halo mass. We expect the latter to be an important variable because massive galaxies at a given stellar mass could live in quite diverse halos.

#### 3.1. Average Trend

The top panels of Figure 2 show the predicted 3D mass density profile for galaxies with stellar masses (Salpeter) in the range of $11.3 < \log M_{\text{star}}/M_{\odot} < 11.7$ (left panels) and $11.7 < \log M_{\text{star}}/M_{\odot} < 12.2$ (right panels). The bottom panels of Figure 2 show the ratio of the total mass density profiles with respect to a pure power law, a profile of the type $\rho(r) \propto r^{-2.2}$. The latter choice is motivated by the findings of Cappellari

<footnote n="15">¹⁵ The lenses in Sonnenfeld et al. (2013) come from the SL2S and SLACS surveys (Auger et al. 2010a; Gavazzi et al. 2014).</footnote>

<footnote n="16">¹⁶ Our SDSS velocity dispersions are defined within an aperture of $R_e/8$ and have been corrected to $R_e/2$ following the mean aperture correction given by Cappellari et al. (2006), based on well-resolved two-dimensional spectroscopic data from SAURON.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Top: predicted 3D density profiles for the combined galaxy–dark matter halo system, for subsamples of galaxies from the Monte Carlo catalogs with stellar masses in the range of $11.3 < \log M_{\text{star}}/M_{\odot} < 11.7$ (left) and $11.7 < \log M_{\text{star}}/M_{\odot} < 12.2$ (right). The solid black line is the total density profile predicted for a de Vaucouleurs stellar density profile, with the blue dotted–dashed, and red long-dashed lines marking the relative contributions of the stellar and dark matter components, respectively. The gray bands mark instead the $2\sigma$ dispersion in mass density at a fixed scale. For comparison, also shown are the mean profiles for Sérsic and Hernquist profiles (purple, dotted–dashed, and green, dashed lines, respectively). The pink very thin stripe in the left panel marks the best-fit empirical fit by Cappellari et al. (2015) $\rho \propto r^{-2.2}$, strictly valid for galaxies with $10.2 < \log M_{\text{star}}/M_{\odot} < 11.7$. Bottom: ratios of the (total) mass density profiles and the Cappellari et al. (2015) $\rho \propto r^{-2.2}$ fits. For lower mass galaxies, the total mass profiles are remarkably consistent with a profile $\rho(r) \propto r^{-\gamma(r)}$ with $\gamma(r) \sim 2.2$, especially around the de Vaucouleurs and Sérsic effective radii (vertical, dotted, and dashed lines, respectively).</figcaption></figure>

et al. (2015, see also Serra et al. 2016) who, confirming and extending previous results, infer from two-dimensional stellar kinematics a total mass density profile of $\rho(r) \propto r^{-\gamma}$ with $\langle\gamma\rangle = 2.19 \pm 0.03$, valid in the range of $0.1R_e$ to $4R_e$ and $10.2 \lesssim \log M_{\text{star}}/M_{\odot} \lesssim 11.7$ (their results are reported with pink stripes in the left panels of Figure 2).

It is, first of all, evident from Figure 2 that the inner regions of massive galaxies are dominated by the stellar component (blue dotted–dashed lines), with the dark matter (red long-dashed lines) gradually taking over at larger scales $r \gtrsim R_e$, as already pointed out over a decade ago (e.g., Borriello et al. 2003; Mamon & Łokas 2005a, 2005b). Nevertheless, the total mean mass density profile (black solid lines) is remarkably well approximated by $\gamma \approx 2.2$ (e.g., Koopmans et al. 2009; Cappellari et al. 2015), especially in the region around the effective radius, with a slight tendency toward a somewhat warped profile (Mamon & Łokas 2005b; Chae et al. 2014). This prediction of the semi-empirical model is in very good agreement with direct strong lensing and dynamical measurements (e.g., Cappellari et al. 2015). It is important to note that the Sérsic effective radii in massive galaxies tends to become progressively larger than the one based on de Vaucouleurs profiles, thus naturally probing regions with flatter density profiles (right panels).

Figure 2 also includes two Hernquist profiles, with their core radii linearly correlated to the 3D half-mass radii (Hernquist 1990). The latter are, in turn, derived from the 2D circularized effective radii (vertical black dotted and orange dashed lines) adopting Equation (13). The green dashed and purple dashed lines show, respectively, the resulting Hernquist profiles for the same stellar masses but different core radii derived from a de Vaucouleurs and Sérsic 2D effective radius and index (equal to four for the de Vaucouleurs profile). It is clear that, irrespective of the input effective radius, the Hernquist profiles are disfavored with respect to a de Vaucouleurs or, even better, a pure Sérsic profile (purple dotted–dashed lines), tending to drop too quickly in the inner few kiloparsecs with respect to the data by Cappellari (2015, left bottom panel).

## 3.2. Dependence on Stellar Mass and Effective Radius

While probing the total mass density profile of massive galaxies has been a subject of intense study by a number of teams (e.g., Koopmans et al. 2009; Barnabè et al. 2011; Ruff et al. 2011; Bolton et al. 2012), only recently the galaxy lensing samples have become large enough to enable the study of its dependence on galactic properties (e.g., Auger et al. 2010a; Cappellari et al. 2013; Chae et al. 2014; Shu et al. 2015; Serra et al. 2016).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Posterior probability distribution function for the model parameters of Equation (10), normalized at $z = 0.1$. The red, yellow, and pink regions mark the $1\sigma$, $2\sigma$, and $3\sigma$ contour levels for each pair of the parameter space extracted from the Sonnenfeld et al. (2013) data. The data are compared to the predictions from the semi-empirical models based with a Salpeter and Chabrier Initial Mass Functions (black and green symbols, respectively). The circles, triangles, and squares refer, respectively, to models with standard profiles, with adiabatic expansion, and with adiabatic contraction, as labeled.</figcaption></figure>

Using strong lensing, Sonnenfeld et al. (2013) have empirically parameterized the dependence of the mass-weighted $\gamma'$ on galactic properties as (their Table 4, second column)

$$ \langle \gamma' \rangle = \gamma_0 + \alpha (z - 0.3) + \beta (\log M_{\text{star}} - 11.5) + \xi \log(R_e/5), \quad (10) $$

with $\gamma_0 = 2.08_{-0.02}^{+0.02}$, $\alpha = -0.31_{-0.10}^{+0.09}$, $\beta = 0.40_{-0.15}^{+0.16}$, and $\xi = -0.76_{-0.15}^{+0.15}$, with an average dispersion around the median of $\sigma_{\gamma'} = 0.12_{-0.02}^{+0.02}$.

The data tend to favor a positive increase of $\gamma'$ with increasing stellar mass, at fixed effective radius. This can be broadly understood in terms of the progressively less significant impact of the (steeper) stellar component when moving to lower mass galaxies, inducing a parallel flattening of the overall profile and thus of $\gamma'$. When fixing the stellar mass and varying the effective radius, $\gamma'$ is instead expected to steepen with decreasing size, as the profile gets progressively dominated by the steeper stellar component. Similarly, $\gamma'$ is expected to steepen even faster with mean (stellar) surface density at fixed stellar mass, since $\Sigma = M_{\text{star}}/(2\pi R_e^2)$.

In essence, the dependence of $\gamma'$ on stellar mass and/or effective radius is connected to probing different portions of the total mass density profile, from the inner steeper profile of the stellar component to the outer flatter profile of the dark matter component, as well as varying the relative contributions of the stellar and dark matter components (see Figure 2). Although the trends of $\gamma'$ summarized in Equation (10) are broadly expected, their precise measurements provide invaluable constraints to galaxy evolutionary processes. To this purpose, in this section, we will compare detailed predictions of our

semi-empirical models to the Sonnenfeld et al. (2013) data and specifically to Equation (10).

Figure 3 reports the probability distribution function (PDF) measured by Sonnenfeld et al. (2013) evaluated at the average redshift of our mock SDSS sample, $z = 0.1$. The red, yellow, and pink regions bracket, respectively, the $1\sigma$, $2\sigma$, and $3\sigma$ contour levels for each pair of the four parameters $\xi$, $\beta$, $\gamma_0$, and $\sigma_{\gamma'}$. In order to make the closest comparison to the data, we treat the outputs of each rendition of our semi-empirical model as in the observations. We apply the same fitting procedure used in Sonnenfeld et al. (2013) to infer the distribution of $\gamma'$ across the population of galaxies, which we model as a Gaussian distribution with a mean given by Equation (10) and with dispersion $\sigma_{\gamma'}$. The fit produces a posterior PDF for the parameters $\gamma_0$, $\beta$, $\xi$, and $\sigma_{\gamma'}$. For each model realization, we plot the peak of the posterior PDF in Figure 3. Due to the large sample size of our mock data sets, statistical uncertainties on these parameters are very small compared to the uncertainty on the same parameters in the Sonnenfeld et al. (2013) observations, and are therefore omitted in the figure. The circles, triangles, and squares in Figure 3 refer, respectively, to models with standard profiles, with adiabatic expansion, and with adiabatic contraction, as labeled, computed as discussed in Section 2.1.

The first apparent feature characterizing the Salpeter IMF model is that a standard, uncontracted NFW dark matter profile (black circles) is consistent with the data at the $1-1.5\sigma$ level in all the parameters' spaces identified by Sonnenfeld et al. (2013). Adiabatic contraction will instead tend to decrease $\beta$ and increase $\xi$, significantly worsening the match to the data (black squares). The endpoint of adiabatic contraction is in
fact to attract dark matter toward more central regions (see Figure 9). By increasing the relative contribution of dark matter at all scales weakens the dependence of $\gamma'$ on stellar mass, thus effectively decreasing the value of $\beta$ in Equation (10). At the same time, adiabatic contraction will also weaken the dependence of $\gamma'$ on effective radius, rendering $\xi$ less negative, as the profile tends to steepen at larger scales increasing $\gamma'$. Conversely, adiabatic expansion (black triangles) tends to progressively lower the dark matter contribution in the inner regions allowing for a more dominant role of the steeper stellar component, thus increasing the zero point $\gamma_0$, increasing $\beta$, and decreasing $\xi$. Figure 3 shows that even milder adiabatic contractions or expansions, as the ones adopted here, tend to be ruled out at the $\gtrsim 3\sigma$ level in at least two or more pairs of parameters. Adiabatic expansion, in particular, appears to be excluded at high significance (black triangles).

The second important point to mention is that a Chabrier IMF (green symbols in Figure 3) is also disfavored by the data, irrespective of any adiabatic contraction/expansion included in the models. The mean $\gamma_0$, in particular, appears to be $>3\sigma$ away the values constrained by the data, irrespective of the exact dark matter profile assumed (green symbols in the middle panels). The model with a Chabrier IMF, in fact, predicts lower stellar masses at fixed halo mass, which in turn reduces the contribution of the steeper stellar component and induces a flatter total mass density profile (see Figure 2), especially around the effective radius, in tension with the strong lensing data by Sonnenfeld et al. (2013, middle panels in Figure 3).

In Appendix C, we also show a series of other systematic variations in the inputs of our reference semi-empirical model (Section 2.1). First off, the PDFs for a Salpeter and variable IMF are quite similar and thus they are both consistent with the Sonnenfeld et al. (2013) strong lensing data (Figure 13). This is partly expected given that the input scaling relations appear to be quite similar for both choices of IMF (Figure 1). In future work, we will explore additional observational proxies such as the velocity dispersion-stellar mass relation to hopefully break these types of degeneracies. Assuming a Hernquist stellar profile (Figure 14) also has a relatively mild effect on the model predictions, as expected given that it mainly flattens out the mass density profile within the regions inside the effective radius (Figure 12), leaving the total mass density profile almost unaltered around $R_e$. Nevertheless, the model characterized by a Hernquist profile in the stellar component tends to be only marginally consistent with the data at the $\sim 2-3\sigma$ level. Moreover, a shallower stellar profile in the inner regions of intermediate-mass early-type galaxies tends to be in tension with strong lensing and 2D dynamical measurements, which point to steep total mass density profiles (e.g., Auger et al. 2010a; Cappellari et al. 2015), as we showed in Figure 2. Introducing a core in the dark matter (Figure 15) severely worsens the match with the data in all parameters, mainly inducing steeper mass density profiles thus larger $\gamma'$ values (middle panels). In this case, adiabatic contraction may help to bring model predictions closer to the data, though discrepancies remain at the $\sim 3\sigma$ level, especially in terms of the very weak dependence on stellar mass ($\beta \sim 0$) and also significantly weaker in size (green squares). As expected, also models characterized by combinations of Hernquist plus cored dark matter profiles tend to be significantly less aligned with the data, being excluded at the $\sim 3\sigma$ level, irrespective of any adiabatic contraction or

expansion (Figure 16). Switching to the Moster et al. (2013) stellar mass-halo mass relation (Figure 17) has a relatively small impact to model predictions, as the Sonnenfeld et al. (2013) data are dominated by intermediate-mass galaxies, and not BCGs, where the difference becomes most prominent (Figure 1). Nevertheless, the adiabatic contraction or expansion tend to still be excluded at the $\sim 3\sigma$ level.

### 4. Dependence of $\gamma'$ on Environment

So far, we have been discussing trends of $\gamma'$ with galaxy properties, ignoring the possible effects of the larger-scale environment, as labeled by the host dark matter halo mass. In order to add this important dimension to our analysis, we need to extend our data sets adopted as comparison. As discussed by Treu et al. (2009), the SLACS survey, which is the foundation of the Sonnenfeld et al. (2013) sample, is dominated by massive galaxies not in cluster-scale halos<fnref n="17" /> (only $\sim$17% of the original sample resides at the center of massive groups and clusters). A stacked weak lensing analysis of a subset of those galaxies also indicated a mean host halo mass of $M_{\text{halo}} \sim 1.6 \times 10^{13} M_{\odot}$ (Gavazzi et al. 2007), with an intrinsic scatter of a factor of $\sim 2$ (Newman et al. 2015). The Sonnenfeld et al. (2013) sample also tends to run out of galaxies above $M_{\text{star}} \sim 9 \times 10^{11} M_{\odot}$. We thus complemented our reference observational sample with 10 galaxies from Newman et al. (2015), specifically selected to be very massive galaxies with means of $M_{\text{star}} \sim 10^{12} M_{\odot}$ (Salpeter IMF) at the center of massive groups with an average mass of $M_{\text{halo}} \sim 10^{14} M_{\odot}$, and seven very massive galaxies by Newman et al. (2013b), all at the center of clusters and with stellar masses of $M_{\text{star}} \gtrsim 10^{12} M_{\odot}$ within clusters of mass $4 \times 10^{14} \lesssim M_{\text{halo}}/M_{\odot} \lesssim 2 \times 10^{15}$, with total mass density profiles obtained via a combination of weak and strong lensing, resolved stellar kinematics within the BCG and X-ray kinematics (Table 9 in Newman et al. 2015). All halo masses are always self-consistently defined as $M_{200c}$ (see Section 2).

#### 4.1. Dependence on Host Halo Mass

To start off with, we compare in Figure 4 the dependence of $\gamma'$ versus (Salpeter) stellar mass for the whole Sonnenfeld et al. (2013) lensing sample (gray band), with the “group” galaxy sample by Newman et al. (2015, green triangles), and with the cluster galaxy sample by Newman et al. (2013a, red circles). The Sonnenfeld et al. (2013) data indicate a very weak dependence of $\gamma'$ on stellar mass. In particular, at masses around $M_{\text{star}} \sim 10^{12} M_{\odot}$, the three samples predict very different total mass density slopes, ranging from values of $\gamma' \sim 2.2$ to values as low as $\gamma' \sim 1.1-1.2$. To better understand the origin of such huge discrepancies in these different data sets, we compare with the predicted mean $\gamma'$ for several different models, a standard de Vaucouleurs profile for a Chabrier, Salpeter, and variable IMF (blue dotted-dashed, red long-dashed, and black solid lines, respectively), and a Sérsic profile with a Salpeter and variable IMF (purple dotted and green dashed lines, respectively), as labeled. The first point to note is that all models predict a more or less rapid drop of $\gamma'$ with increasing stellar mass, irrespective of the input stellar profile or IMF (left panel). However, if in our

<footnote n="17">¹⁷ We point out that when comparing with Sonnenfeld et al. (2013) in Figure 3 or in the Appendices, we do not apply any halo mass cut in our mock galaxies, though we verified that retaining only galaxies with, say, $\log M_{\text{halo}}/M_{\odot} < 13.5$, does not alter our main results.</footnote>
<figure><img src="image_7.png" /></figure>

**Figure 4.** Predicted $\gamma'$ dependence on (Salpeter) stellar mass for several different models, as labeled, a standard de Vaucouleurs profile for a Chabrier, Salpeter, and variable IMF (blue dotted-dashed, red long-dashed, and black, solid lines, respectively), and a Sérsic profile with a Salpeter and variable IMF (purple, dotted and green, dashed lines, respectively). All models assume an NFW profile for the dark matter component. The left and right panels show, respectively, predictions for the same models including all host dark matter halos and restricted to $\log M_{\text{halo}}/M_{\odot} < 13.5$. The data, the same in both panels, are extracted from the Sonnenfeld et al. (2013) sample (gray band), from the “group” galaxy sample by Newman et al. (2015, green triangles), and from the cluster galaxy sample by Newman et al. (2013a, red circles). Galaxies of similar stellar mass can show quite different total mass density profiles most probably induced by differences in structural properties and in the environments they live in.

<figure><img src="image_8.png" /></figure>

**Figure 5.** Predicted $\gamma'$ dependence on effective radius for galaxies with $\log M_{\text{star}}/M_{\odot} > 11.9$ for three different models, as labeled: a standard Sérsic+NFW (black, solid lines), a Hernquist+NFW (blue, dotted-dashed lines), and a Sérsic+NFW adopting the stellar mass–halo mass relation by Moster et al. (2013, purple, dotted lines). All models assume a Salpeter IMF. The left, middle, and right panels show predictions for the same models assuming dark matter halo profiles uncontracted, with Adiabatic Expansion, and Adiabatic Contraction, respectively. Filled, blue squares, green triangles, and red circles are data from the $z < 0.3$ Sonnenfeld et al. (2013) sample, from the “group” galaxy sample by Newman et al. (2015), and from the cluster galaxy sample by Newman et al. (2013b), all with $\log M_{\text{star}}/M_{\odot} > 11.9$ and Salpeter IMF. Galaxies of similar stellar mass and effective radius can still show quite different total mass density profiles most probably induced by the different environments they live in.

mocks we select only “field” galaxies (right panel), i.e., we disregard galaxies at the center of massive groups and clusters above, say, $\log M_{\text{halo}} = 13.5$ (which is around the upper limit probed by Sonnenfeld et al. 2013), we tend to find steeper profiles in high stellar mass galaxies induced by the lower contribution of the host dark matter halo characterized by flatter profiles. This simple cut in host halo mass is sufficient to reconcile de Vaucouleurs model predictions with the Sonnenfeld et al. (2013) data (gray band). It can be seen that both the de Vaucouleurs models with Salpeter and variable IMF match observations, with the Chabrier IMF model predicting generally lower $\gamma'$, in line with what is discussed in reference to Figures 3 and 13.

However, even after excluding massive host halos, the Sérsic-based models tend to fall to lower $\gamma'$ values at large stellar masses. This key difference with respect to the predictions of de Vaucouleurs-based models could be mostly ascribed to the

progressively larger effective radii in more massive galaxies with Sérsic profiles (middle panels of Figure 1). Having larger radii, the Sérsic-based models probe larger and flatter portions of the mass density profile, in which the dark matter starts dominating, as emphasized in the right panel of Figure 2. We conclude that the Newman et al. galaxies show significantly flatter total mass density profiles with respect to Sonnenfeld et al. mainly because they are central galaxies in groups and clusters, thus they are physically more extended and better described by Sérsic profiles (Bernardi et al. 2011a, 2011b, 2014), and also characterized by significantly more massive host dark matter halos.

To further probe the dependence on structural properties, in Figure 5, we plot $\gamma'$ against the effective radius for the most massive galaxies in the Sonnenfeld et al. (2013) catalog, and the group and cluster galaxy samples of Newman et al. (2013b, 2015), shown by blue squares, green triangles, and red circles, respectively. We include three semi-empirical models in this
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Predicted dependence of $\gamma'$ on host halo mass for five different models, a standard Sérsic+NFW (black, solid lines), a Sérsic+cored NFW (red, long-dashed lines), Hernquist+NFW (blue, dotted-dashed lines), a Hernquist+cored NFW (orange, dashed lines), and a Sérsic+NFW adopting the stellar mass-halo mass relation by Moster et al. (2013, purple, dotted lines). All models assume a Salpeter IMF. The left, middle, and right panels show predictions for the same models assuming dark matter halo profiles uncontracted, with Adiabatic Expansion and Adiabatic Contraction, respectively. As in Figure 5, filled blue squares, green triangles, and red circles in the top panels are data from Sonnenfeld et al. (2013), Newman et al. (2015), and Newman et al. (2013b). The bottom panels show the predicted trends for a narrow bin in stellar mass and effective radius, as labeled. The predicted total mass density profile has a relatively strong dependence on halo mass with $\gamma' \propto M_{\text{halo}}^{0.2-0.3}$, mostly because a higher host halo mass induces an overall flatter mass density profile and thus lower $\gamma'$.</figcaption></figure>

figure, characterized by standard stellar Sérsic and dark matter NFW profiles (black solid lines), stellar Hernquist and dark matter NFW profiles (blue dotted-dashed lines), and a stellar Sérsic and dark matter NFW based on the stellar mass-halo mass relation by Moster et al. (2013, purple dotted lines). For consistency with the data, in Figure 5 (as well as in Figure 6), we show results for the subsample of mock galaxies with $log M_{\text{star}} > 11.9$, which is comparable to the lower limit in stellar mass probed in the Newman et al. samples, and still low enough to retain a statistically relevant number of galaxies in our mocks. For the same reason, we also restrict the Sonnenfeld et al. (2013) data to only galaxies with $log M_{\text{star}} > 11.9$. Despite the variety of mass profiles and stellar-dark matter mappings the predicted $\gamma'$ distributions show quite similar behaviors, differing only in normalization by just a modest amount of $\lesssim 10\%$. More interestingly, the predicted mean $\gamma'$ substantially decreases with increasing size, as the (flatter) dark matter component becomes more dominant at larger scales (see discussion in Section 3.2). Our model predictions, at least the ones without adiabatic contraction, are in fact, in broad agreement with the data, and the predicted trend is by itself sufficient to explain a large portion of the difference between the $\gamma'$ values in the Sonnenfeld et al. (2013) sample and in the group galaxies of Newman et al. (2015).

However, the models appear to still be highly inconsistent with the $\gamma'$ measurements by Newman et al. (2013b) for BCGs (red circles). Their median $\gamma'$ lays significantly below our estimates, and even below the data by Newman et al. (2015) for galaxies with similar stellar masses and effective radii but residing in group-scale halos. The mean offset in stellar mass

between the Newman et al. (2015) group sample and the Newman et al. (2013b) cluster sample is, in fact, only ~0.17 dex, which is too modest to explain the strong observed systematic difference in $\gamma'$ purely in terms of stellar mass. This finding is consistent with the hypothesis that halo mass may play a non-negligible role in shaping the overall mass density profiles of very massive galaxies at the center of clusters (see also Kochanek & White 2001).

To explore this issue further, we present in Figure 6 the predicted dependence of $\gamma'$ on host halo mass from our reference models, compared to data from Sonnenfeld et al. (2013) for field/group galaxies (blue squares), Newman et al. (2015) for group central galaxies (green triangles), and Newman et al. (2013b) for BCGs (red circles). It is first of all apparent, as also pointed out by Newman et al. (2015), that $\gamma'$ depends on environment/halo mass. While some of the data appear more mixed in the $\gamma'–R_e$ plane (Figure 5), they tend to separate well when $\gamma'$ is expressed as a function of host halo mass. Predictions for a series of models (with Salpeter IMF) are included in Figure 6: standard stellar Sérsic and dark matter NFW profiles (black, solid lines), stellar Sérsic and dark matter cored NFW profiles (red, long-dashed lines), stellar Hernquist and dark matter NFW profiles (blue, dotted-dashed lines), stellar Hernquist and cored NFW profiles (orange, dashed lines), and stellar Sérsic and dark matter NFW profiles, adopting the stellar mass-halo mass relation by Moster et al. (2013, purple, dotted lines).

All models, especially the ones with no adiabatic contraction, predict a drop of $\gamma'$ with increasing halo mass in line with the data. Our reference model with uncontracted NFW dark matter profiles,
in particular, predicts (left panel) a dependence on halo mass of the type

$$ \gamma' \propto M_{\text{halo}}^p \quad (11) $$

with $p = -0.30 \pm 0.01$, in good agreement with Newman et al. (2015) who find $p = -0.33 \pm 0.07$ when combining different samples of lensed galaxies in the field, in groups, and in clusters (each single data set, by itself, does not show any dependence in halo mass). The value of $p$ we find is strictly valid for $\log M_{\text{halo}}/M_{\odot} \gtrsim 13$, with a possible gradual flattening below this mass scale. To make a fairer comparison with the data, we again consider in Figure 6 only galaxies from our semi-empirical models with $\log M_{\text{star}}/M_{\odot} > 11.9$. However, we have verified that the dependence of $\gamma'$ on halo mass is quite independent of the exact cuts in stellar mass or effective radius. As such, it could be considered to be an addition to Equation (10). We note, however, that observational calibrations of the slope $p$ require uniform and precise measurements of strong and weak lensing galaxies over a large range of host halo masses. Non-negligible errors in halo mass could, in fact, wash out some of the trends with halo mass (e.g., Shankar et al. 2014b). The dependence of $\gamma'$ on halo mass is a genuine effect and only partly induced by the increase of effective radius with stellar mass. This is clearly seen in the bottom panels of Figure 6 in which we plot the same models as in the upper panels but restricted in the narrow bins of stellar mass and effective radius, $12.05 < \log M_{\text{star}}/M_{\odot} < 12.15$ and $1.2 < \log R_e/\text{kpc} < 1.3$. In this case, the trend of $\gamma'$ with halo mass gets somewhat flatter but is still highly significant $p = -0.21 \pm 0.02$.

While Equation (11) has a relatively weak dependence on galaxy properties, it has a significant dependence on the choice of input NFW dark matter halo profile. The left, middle, and right panels of Figure 6 show predictions for the same models assuming dark matter halo profiles uncontracted, with Adiabatic Expansion and Adiabatic Contraction, respectively. It is clear that irrespective of the chosen stellar and/or dark matter mass profiles, and/or stellar mass-halo mass mapping, the predicted values of $\gamma'$ tend to decrease at broadly the same pace with increasing halo mass, with only a small offset in normalization of $\lesssim 15\%$. This behavior is mainly induced by structural non-homology, with the relative density distribution of dark matter becoming gradually more dominant in cluster-sized halos. Despite decreasing with increasing halo mass, the predicted values of $\gamma'$ still tend to be steeper than those actually measured in halos above $M_{\text{halo}} \gtrsim 3 \times 10^{14} M_{\odot}$ by Newman et al. (2013b). Including adiabatic expansion in the models yields a stronger decrease of $\gamma'$ with effective radius, mirroring the flattening of the mass profile at $r \gtrsim R_e$ (see Section 3.2), in apparent better agreement with the data. Still, we do not properly reproduce the total mass density profiles derived by Newman et al. (2013b) for dominant galaxies in rich clusters.

Figure 7 reports a closer comparison between model predictions and the total mass density profile measured by Newman et al. (2013b). In all top panels of Figures 7 and 8, the black long-dashed line and gray area mark, respectively, the Newman et al. (2013b) mean total mass density profile and $1\sigma$ dispersion, while the light green area roughly marks their $1\sigma$ dispersion in the BCG stellar profile. The bottom panels of

Figures 7 and 8 plot instead the ratios between model predictions and data for each model. In all panels, the predicted stellar and dark matter contributions to the total mass density profiles are reported by blue triple dotted-dashed, and red long-dashed lines, respectively. The left panel in Figure 7 shows predictions for the reference stellar Sérsic and dark matter NFW profiles. There are two main issues. The reference model is much steeper than the data at small scales, severely overpredicting them at $r \lesssim 4$ kpc by up to a factor of $\sim 3-4$. At intermediate scales of $10 \lesssim r \lesssim 300$ kpc, the models instead significantly underpredict the data by a systematic factor of $\sim 2$.

Switching to a Hernquist profile in the stellar component naturally produces profiles at small scales much better aligned with the data (middle panel of Figure 7). The resulting stellar profile is indeed in good agreement with the stellar profile inferred by Newman et al. (2013a, light green areas), predicting a density of $\rho_{\text{star}} \sim 3 \times 10^8 M_{\text{star}}/\text{kpc}^3$ at $r \sim 2$ kpc, and $\rho_{\text{star}} \sim 1.5 \times 10^7 M_{\text{star}}/\text{kpc}^3$ at $r \sim 10$ kpc. A Hernquist profile is a simple, but clearly not a unique, solution. Adding, for example, an appropriate “core” radius (e.g., Trujillo et al. 2004; Newman et al. 2013b) to the stellar profile could yield similar improvements. A flatter core in the inner regions of very massive galaxies, especially at the center of clusters, has indeed been reported several times in the literature (e.g., Kormendy & Bender 2009; Del Popolo 2012; Postman et al. 2012; Huang et al. 2013; Oldham & Auger 2016). On the other hand, simply assigning a lower Sérsic index of, say, $n \lesssim 3$ to the most massive galaxies, though flattening the inner density profile (see Appendix B), would also decrease the predicted stellar mass density profiles in more external regions around a few to 10 kiloparsecs, worsening the match to the galaxy stellar profiles measured by Newman et al. (2013b). Moreover, detailed light profile fits to massive galaxies have revealed their Sérsic indices to be as large as $n \gtrsim 5$ (e.g., Huertas-Company et al. 2013; Bernardi et al. 2014).

Despite the possible improvements at small scales, our reference models still significantly underpredict the mass density observationally inferred by Newman et al. (2013b) at intermediate scales beyond $\gtrsim 10$ kpc, where the mass profile is expected to become fully dominated by the dark matter component. We find that additionally allowing for a boost of a factor of two in the median halo concentration helps to resolve the mismatch with the data, as shown in the right panel of Figure 7. Increasing the concentration in fact decreases the core radius of $r_s = r_{200c}/c_{200c}$ at fixed halo mass and virial radius $r_{200c}$, thus naturally increasing the mass density in the inner regions, with the end result of flattening the total mass density profile, especially around the effective radius (vertical, dotted lines in Figure 7). However, we note that an increased concentration also tends to predict a larger mass density at $r \lesssim 10$ kpc. To avoid the latter effect, when adopting a higher concentration, we must also allow for a core in the dark matter profile, for which we adopt Equation (2). We checked that also adiabatic expansion could reduce the dark matter in the inner regions at comparable levels to those of a core.

Is there any observational justification for increasing concentrations over the predictions of N-body simulations? An increase in dark matter concentration at cluster scales is indeed possible (e.g., Comerford & Natarajan 2007; Auger et al. 2013, and references therein), though, at present, good evidence for real clusters to be on average more concentrated is still not convincing (e.g., Ettori
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Total mass density profile predicted from three models, the reference Sérsic+NFW (left), Hernquist+NFW (middle), and a Hernquist+cored NFW with $r_{\text{core}} = 14$ kpc and a mean concentration higher by a factor of two (right). The top panels compare the model outputs with the data by Newman et al. (2013b). In each panel, the black, long-dashed line along with the gray area mark, respectively, the mean profile and its $1\sigma$ dispersion, while the light green area roughly brackets the $1\sigma$ dispersion in the average stellar profiles of the Newman et al. (2013b) central galaxies. The bottom panels show the ratios between model predictions and data. The stellar and dark matter contributions to the total mass density profiles are reported by blue, triple dotted-dashed, and red, long-dashed lines, respectively. An inner flatter stellar component and higher dark matter concentrations tend to yield better fits to the data. The vertical dotted lines in all panels mark the position of the mean effective radius in the mock galaxies.</figcaption></figure>

et al. 2010; Dutton & Macciò 2014; Meneghetti et al. 2014; Merten et al. 2015; Shan et al. 2015; Umetsu et al. 2016, and references therein). Newman et al. (2013b), and more recently Amodeo et al. (2016) at $z \gtrsim 0.4$, empirically constrained the median concentration at $M_{\text{halo}} \sim 10^{15} M_{\odot}$ from their samples to be $c \sim 5–6$, in broad agreement with our inferred values of $c \sim 2 \times 3 \div 4 \sim 6 \div 8$, which are consistent with the numerical results by Prada et al. (2012, but see Meneghetti & Rasia 2013). Evidence for larger concentrations for red, massive galaxies consistent with what is presented in Figure 7 also come from the kinematics of satellites, though mostly inferred from isolated galaxies (Wojtak & Mamon 2013). Giocoli et al. (2014, see also Giocoli et al. 2016) discuss how selection effects may bias the measured concentration–mass relation at cluster scales, yielding up to 30% higher concentrations than what is predicted by N-body simulations (see also Auger et al. 2013; Sereno et al. 2015; Lieu et al. 2017). Several studies have also highlighted the possibility of biases arising from halo triaxiality (e.g., Foëx et al. 2014). An alignment of the major axis and the line of sight may cause an apparent higher concentration, that may be evident in lensed groups and clusters. Although the Newman et al. (2013b) cluster sample is relatively small, the selection is well defined, and in fact Newman et al. (2013b) find that only in one case there is substantial elongation along the line of sight (A383), which could artificially boost the measured concentration.

As numerical simulations and direct observations do not tend to unanimously support dark matter concentrations in cluster-sized halos much higher than what predicted by numerical simulations, we have explored in Figure 8 the effects of switching from an NFW to an Einasto profile, which is

characterized by an additional parameter $\mu$ (Equation (3)), and it has been shown to be an accurate fit to massive dark matter halos by a number of groups (e.g., Navarro et al. 2004, 2010; Stadel et al. 2009; Ludlow et al. 2014; Del Popolo 2014; Dutton & Macciò 2014; Klypin et al. 2016). The left panel of Figure 8 shows the results of a standard Sérsic+unperturbed dark matter Einasto profile with $\mu = 4$ (or $\alpha = 1/\mu = 0.25$), which is the average value extracted from numerical N-body simulations of very massive clusters (e.g., Navarro et al. 2004; Gao et al. 2008; Dutton & Macciò 2014). This basic model tends to present broadly the same shortcomings as the standard NFW-based model presented in the left panel of Figure 7. However, the match is substantially improved at all scales when simply allowing for a Hernquist stellar profile and a median concentration just 30% higher than the reference numerical value (middle panel). A higher value of $\mu = 7$ (right panel) further increases the dark matter density at intermediate scales without the need to invoke any increase in median concentration, but the match to the data in this case is poorer.

In summary, the present available lensing data for galaxies at the center of cluster-scale halos (Newman et al. 2013b), tend to favor models with profiles in their inner stellar component at $r \lesssim 3–5$ kpc, flatter than what is predicted by extrapolations of their outer Sérsic profiles. The same data also tend to align better with models characterized by an NFW profile with a core and a median concentration in the dark matter component twice the one predicted by dark-matter-only N-body simulations, or by an Einasto profile with $\mu \sim 4$ and a dark matter concentration just $\lesssim 30\%$ higher than the simulation results. The latter results are in nice agreement with Auger et al. (2013)
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Same format as Figure 7 showing models characterized by a Sérsic+Einasto with $\mu = 4$ (left), a Hernquist+Einasto with $\mu = 4$ and a mean concentration higher by a factor of 1.3 (middle), and a Hernquist+Einasto with $\mu = 7$ (right). The match to the data in this case does not require a prominent core nor a much higher median dark matter concentration.</figcaption></figure>

who also found 26 group- and cluster-scale strong gravitational lenses to be only ~26% more overconcentrated than similar-mass halos from dark matter simulations. We stress that the relatively small comparison observational sample prevent us to generalize our findings to all galaxy clusters. The idea that Einasto models may be more accurate representations of the mass density profiles of especially the most massive halos is in line with a number of analytic and numerical studies (e.g., Zhao et al. 2003; Gao et al. 2008; Lapi & Cavaliere 2011; Cen 2014; Nipoti 2015; Angulo et al. 2016; Klypin et al. 2016; Ludlow & Angulo 2017). The common theme in these works is the connection between the final shape of the system, and thus of its parameter $\mu$, with its mass accretion history, in a way that $\alpha = 1/\mu$ and “curvature” tend to increase with fast growth, as it might be the case, on average, for the inner regions of very massive halos.

## 4.2. Dark Matter Fractions

Up to this point, we have mainly discussed the outputs of semi-empirical models with respect to $\gamma'$, but clearly additional observables can and must be used to place tighter constraints on galaxy evolution. In particular, scaling relations with velocity dispersions (e.g., Dutton et al. 2013), evolution of the fundamental plane (e.g., Treu et al. 2001; Shankar & Bernardi 2009; Shankar et al. 2013; Oldham et al. 2017), evolution of the virial relations (e.g., Peralta de Arriba et al. 2014), building up of age/metallicity gradients (e.g., Montes et al. 2014), are all relevant issues to test models.

While we aim to explore at least some of these additional issues in a separate work, we here begin by comparing in Figure 9 our predicted 3D ratio of dark to total (dark plus stellar) mass as a function of scale (radius normalized to the effective radius) for all

galaxies in the range of $11.3 < \log M_{\text{star}}/M_{\odot} < 11.5$ (left panel) and $6 < R_{\text{e}}/\text{kpc} < 7$ and $\log M_{\text{star}}/M_{\odot} > 12$ (right panel), compared to the lensing-based dark matter fractions derived from the Sonnenfeld et al. (2015), thick blue lines, left panel) and Newman et al. (2013a, thick solid and dotted red lines, right panel) galaxy samples. Individual dark matter fractions in the Newman et al. (2013b) cluster central galaxies are calculated by integrating the density profiles of their Figure 3. We then take the mean of the dark matter fractions at each radius to produce the solid and dotted thick red lines reported in the right panel of Figure 9. The black solid, orange long-dashed, and purple dotted-dashed lines show predictions for the uncontracted NFW profile, with adiabatic contraction, and adiabatic expansion, respectively. Both panels show predictions for our standard stellar Sérsic and dark matter NFW profiles. It is apparent that at cluster scales (right panel) the dark matter fractions rise relatively more rapidly than for lower mass galaxies in the lower mass halos (left panel), irrespective of the details in the stellar and/or dark matter profiles. This trend is mostly induced by the fact that more massive galaxies, characterized by progressively higher effective radii, naturally will include more dark matter mass in their central regions.

Our semi-empirical models, within the uncertainties, align with the observationally based dark matter fractions, and we verified this broadly continues to hold true even with flat inner stellar profiles, and/or with cores in the dark matter component. In this respect, dark matter fractions appear to be less constraining for models, though adiabatic contraction or expansion continue to be disfavored by the data irrespective of the exact inputs in the semi-empirical models, in line with Oguri et al. (2014). Our predicted dark matter fractions are in full agreement with the $f_{\text{dm}}(R_{\text{e}}) \sim 20\%–30\%$ values inferred by Cappellari et al. (2013, see also Wojtak & Mamon 2013) for $M_{\text{star}} \sim 2–3 \times 10^{11} M_{\odot}$, or on larger scales with Alabi et al.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Ratio of dark matter mass and total (dark+stellar) mass as a function of scale normalized to the effective radius for galaxies in the range of 11.3 < log M<sub>star</sub>/M<sub>☉</sub> < 11.5 (left) and 6 < R<sub>e</sub>/kpc < 7 and log M<sub>star</sub>/M<sub>☉</sub> > 12 (right). The black solid, orange long-dashed, and purple dotted-dashed lines show predictions for the uncontracted NFW profile, with Adiabatic Contraction and Adiabatic Expansion, respectively. Both panels show predictions for standard Sérsic +NFW profiles. Overall, the data tend to disfavor any adiabatic contraction or expansion.</figcaption></figure>

(2017), who infer in the SLUGGS survey $f_{dm} (<5R_e) \gtrsim 60\%$ from globular cluster kinematics data. At larger masses, our models predict $f_{dm}(R_e) \sim 60\%–70\%$ for $M_{star} \gtrsim 10^{12} M_{☉}$, consistent with the extrapolation of the Cappellari et al. (2013) analytic fit, $f_{dm}(R_e) \sim 1.3 + 0.24(\log M_{star} - 10.6)^2$, to the most massive galaxies in our mocks.

For completeness, to visualize the full dependence of the mass density profiles on the relative contributions between stellar and dark matter components, in Figure 10, we plot $\gamma'$ as a function of dark matter fractions $f_{dm}(<R_e)$ within one effective radius for galaxies with $M_{star} \gtrsim 3 \times 10^{11} M_{☉}$ as predicted by our reference semi-empirical model with uncontracted NFW (solid and dotted lines), with adiabatic contraction (orange, dashed line), and with adiabatic expansion (purple, dotted-dashed line). All models predict a steady drop of $\gamma'$ with increasing dark matter fraction, as expected given that, as discussed with respect to Figure 9, larger $f_{dm}(<R_e)$ tend, on average, to be linked to larger effective radii with more prominent dark matter contributions, and thus flatter $\gamma'$, also in line with Figure 5. However, models with adiabatic expansion tend to predict a substantially faster drop with $f_{dm}(<R_e)$ than the other two models. Our uncontracted NFW predictions nicely match those from the dynamical modeling on ATLAS3D and SPIDER galaxies by Tortora et al. (2014a, empty squares with error bars) who also assume an underlying NFW profile. On the other hand, while model predictions agree with the dark matter fractions extracted from Sonnenfeld et al. (2015, blue squares), they are highly inconsistent with the dark matter fractions derived from Newman et al. (2013a, red circle) for galaxies in clusters. As already emphasized for Figure 5, galaxies characterized by similar stellar mass, effective radius, and even dark matter fractions tend to significantly differ in their total mass density slopes, an additional sign that host halo mass and its (outer) structure are systematically different.

# 5. Discussion

## 5.1. Constraints from the Dependence of $\gamma'$ on Size and Stellar Mass

We have carried out a detailed comparison of a set of semi-empirical models with a diverse set of observational data. The

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Dependence of $\gamma'$ on dark matter fraction within one effective radius for galaxies with $M_{star} \gtrsim 3 \times 10^{11} M_{☉}$ as predicted by our reference Sérsic +NFW semi-empirical model with uncontracted NFW (solid and dotted lines), with adiabatic contraction (orange dashed line), and with adiabatic expansion (purple dotted-dashed line). Empty squares are the results of the dynamical modeling on ATLAS3D and SPIDER by Tortora et al. (2014a), blue squares from Sonnenfeld et al. (2015), and red circle from Newman et al. (2013b) for galaxies in clusters. It is clear that the latter subsample of galaxies is significantly off with respect to galaxies of similar dark matter fraction, indicating that galaxies in rich clusters are systematically different.</figcaption></figure>

models are based on assigning galaxies, with measured stellar mass, effective radius, and Sérsic index, to dark matter halos via abundance matching techniques. By design, our approach is based on a minimal number of parameters and underlying assumptions, mainly related to the underlying dark matter mass distribution, which is relatively well constrained by high-resolution N-body simulations (Prada et al. 2012; Dutton & Macciò 2014; Diemer & Kravtsov 2015), and the IMF, which we usually set equal to a Salpeter or even variable IMF, having shown in Figure 3 that a Chabrier IMF is ruled out by the data in massive galaxies, in line with a number of independent studies based on surface gravity-sensitive absorption lines in stellar population synthesis models (e.g., Conroy & van
Dokkum 2012; La Barbera et al. 2013; Spiniello et al. 2014), and/or stellar kinematics also coupled to strong lensing measurements (e.g., Auger et al. 2010a; Treu et al. 2010; Spiniello et al. 2011; Thomas et al. 2011; Barnabè et al. 2013; Cappellari et al. 2012, 2013; Posacki et al. 2015).

Our semi-empirical models thus represent ideal tools to
probe the consistency among independent data sets, to reliably
test the predictions of dark matter numerical simulations, to pin
down the role of adiabatic contraction/expansion, and
environment in shaping the total mass density profiles of
massive galaxies. The semi-empirical models successfully
reproduce the total mass density profile of massive galaxies
with $M_{\text{star}} \gtrsim 3 \times 10^{11} M_{\odot}$ (Salpeter IMF), along with its
dependence on stellar mass and effective radius, at least in
non-cluster environments. We find that even moderate
modifications to the halo profile in terms of adiabatic
contraction/expansion are disfavored at high significance. In
the next steps (F. Shankar et al. 2017, in preparation), we will
mainly focus on the evolution of $\gamma'$ as a function of redshift, as
well as the additional constraints that can be gained by a
comprehensive Jeans-based calculation of velocity dispersions
(e.g., Chae et al. 2014).

Our results on the dependence of $\gamma'$ on galaxy properties are
in line with several previous works. We discuss some of these
here below.

Xu et al. (2016) have recently compared the properties of
simulated early-type galaxies in the Illustris hydrodynamic
cosmological simulation with the same SLACS data as in this
paper. Their simulations broadly align with observations,
suggesting a flattening of the total mass density profile with
increasing effective radius, though generally shallower than
what the data suggest.

Poci et al. (2017) have also recently confirmed from ATLAS3D that the slopes of the total mass density profiles are approximately isothermal, and depend on various galactic properties, most significantly on mass surface density, while the correlation with velocity dispersion becomes negligible above $\log \sigma/\text{km s}^{-1} \gtrsim 2.1$.

Tortora et al. (2014a) have explored the density slope
dependence in early-type galaxies over quite a broad range of
stellar mass finding clear signs of non-universality. They
suggested that the relation with effective radius is the most
robust to changes in the underlying modeling. They find, in
fact, that the largest galaxies have a density slope close to
isothermal, while the smallest ones are characterized by a
steeper profile, in line with our findings.

Dutton & Treu (2014) have also addressed the question of
dependencies of $\gamma'$ on galactic properties, though adopting only
one stellar mass–halo mass mapping. In agreement with our
work, they disfavor models with significant adiabatic contrac-
tion or expansion, and also predict an overall non-universal $\gamma'$
dependence on galactic properties.

Auger et al. (2010b) compared data on 53 massive early-type
lensed galaxies with a variety of semi-empirical models with
varying baryon efficiency, adiabatic contraction, and stellar
IMF. They highly disfavored light IMFs, such as a Chabrier of
Kroupa, as well as strong adiabatic contraction, in nice
agreement with the results presented in this work. Their
models also suggest that only mass-dependent IMFs could be
consistent with no adiabatic contraction, a conclusion not
necessarily supported by our findings (see, e.g., Figure 13).

Jiang & Kochanek (2007) also carried out semi-empirical
galaxy–halo modeling, treating the stellar mass-to-light ratio as a
free parameter. They analyzed surveys of gravitational lenses
with velocity dispersion measurements assuming combinations
of Hernquist and NFW profiles. Their total mass density slopes
are peaked around the isothermal values with non-negligible
variations around the mean. They also found that models without
adiabatic contraction favor a stellar fraction of around 3%, in
agreement with weak lensing estimates (e.g., Mandelbaum
et al. 2006) and the most recent abundance matching results by
Shankar et al. (2014a) considered in this work (Figure 1).

Dutton et al. (2011) modeled the scaling relations of early
and late-type galaxies, finding that massive galaxies character-
ized by a Salpeter IMF are not consistent with halo contraction.
Napolitano et al. (2010) also emphasized that that the lower
central dark matter fractions inferred in the oldest early-type
galaxies are inconsistent with adiabatic contraction.

By analyzing a suite of simulated spheroids formed out of
binary mergers, Remus et al. (2013) concluded that the exact
value of the density slope depends on the amount of gas
involved in the merger, with more dissipative gas-rich mergers
producing steeper profiles. This would naturally imply some
non-homology, as more dissipative mergers would produce
more compact galaxies at fixed stellar mass with steeper slopes,
in line with what observed. Dissipative processes may, how-
ever, not be the only cause for the dependence of $\gamma'$ on galactic
properties. Indeed, as discussed in Section 3.2, in basic ΛCDM
models with uncontracted halos, the trend of $\gamma'$ with effective
radius is simply controlled by the rate of progressive change
from the steeper stellar component to the flatter dark
matter one.

Remus et al. (2017, see also Dubois et al. 2013; Sonnenfeld
et al. 2014) also recently emphasized that the decreasing
importance of dissipative processes toward lower redshifts and
for more massive systems is mostly responsible for shaping the
strong flattening of $\gamma'$ with increasing dark matter fractions
$f_{\text{dm}}(<R_e)$ within one effective radius as measured by Tortora
et al. (2014a). As detailed in Figure 10, galaxies in clusters do,
however, show significantly smaller $\gamma'$ distributions at fixed
$f_{\text{dm}}(<R_e)$, supporting a scenario in which also the larger-scale
environment plays a significant role in defining the total mass
density profiles of early-type, massive galaxies.

## 5.2. Constraints from the Dependence of $\gamma'$ on Halo Mass

We showed that the mass-weighted total mass density slope
at the effective radius $\gamma'$ is a strong decreasing function of the
host halo mass, ranging from values close to, or higher than,
the canonical isothermal value of $\gamma' \sim 2$ in the field/low-
mass groups, to values closer to $\gamma' \sim 1.5$ at the cluster scale.
This behavior is mainly induced by the increasing contrib-
ution of the dark matter component around the effective
radius, thus inducing $\gamma'$ to progressively flatten out (see
Figure 2).

A closer comparison to the data interestingly reveals that for
very massive galaxies with $\log M_{\text{star}}/M_{\odot} \gtrsim 11.9$ in cluster-
scale environments (see Figure 7), the stellar profile at small
scales $r \lesssim 3–5$ kpc should be significantly flatter than what
was predicted by an extrapolation of a Sérsic profile. This is in
line with several independent lines of evidence (Kormendy &
Bender 2009; Postman et al. 2012), as mentioned in
Section 4.1.
In general, several physical processes are thought to play some role in shaping the density profiles of massive galaxies, especially in overdense regions. Laporte & White (2015), for example, have presented a series of collisionless N-body resimulations to follow the growth of central cluster galaxies and their dark matter halos since redshift $z \lesssim 2$. In line with what is discussed in this work, Laporte & White (2015) report a steepening at small scales in the stellar component not mirrored in the data. They propose that the action of repeated super-massive black hole mergers could create enough large cores of a few kiloparsecs to alleviate the tension between models and observations (see also, e.g., Merritt 2006).

However, other in situ mechanisms may be able to induce a permanent flattening of the central matter distributions, such as the action of efficient feedback (e.g., Peirani et al. 2008; Martizzi et al. 2013) from central active galactic nuclei (AGNs). On the other hand, an "in situ" process such as AGN feedback, a phenomenon more and more supported by direct and indirect observations (e.g., Shankar et al. 2016, 2017; Barausse et al. 2017, and references therein), should produce a stellar core in all massive galaxies, irrespective of their environment.

Laporte & White (2015) also discuss that multiple dissipationless mergers can effectively flatten out the inner distributions, with the dark matter component producing a shallower central cusp, and gravitational heating redistributing material in the outer regions (see also, e.g., El-Zant et al. 2004, 2016; Nipoti et al. 2004; Tonini et al. 2006; Lapi & Cavaliere 2009). This type of evolution is in line with our proposed modifications presented in Figure 7 of a cored NFW profile to faithfully align model predictions to observations in the range of $3 < R/kpc < 20$.

Finally, we note that other non-baryonic processes mainly linked to the dark sector, such as self-interacting dark matter, may contribute to the diversity in total mass density profiles inferred in this work (e.g., Kaplinghat et al. 2016; Di Cintio et al. 2017).

Chae et al. (2014) fitted two-component models to local SDSS nearly spherical galaxies. Adopting a velocity dispersion-dependent variable IMF, they found a mean slope for the total mass density profiles of $\langle\gamma_e\rangle = 2.15 \pm 0.04$, computed in the radial range of $0.1R_e < r < R_e$, in full agreement with our findings. Chae et al. also went forward in testing the dependence on galactic properties including effective radius and halo mass. Their least-square fit total mass density slope between $0.1R_e$ and $R_e$ also decreases with effective radius, though more weakly than the one reported in Figure 5. Chae et al. (2014) also infer a dependence of the total mass density slope on host halo mass weaker than the one measured here, though they also show that a much stronger dependence is naturally obtained when the slope of the mass density profile is fitted on scales larger than their effective radius.

## 6. Conclusions

In this work, we have set up a semi-empirical approach to create large catalogs of mock galaxies and host dark matter halos at $z = 0.1$. To this purpose, we have selected a large sample of galaxies from SDSS with a probability $P(E) > 0.85$ of being ellipticals, and with measured stellar masses, effective radii, and Sérsic indices. We then assigned host dark matter halos to each galaxy in the sample making use of

the most recent stellar mass-halo mass relations (Shankar et al. 2014a), as well as exploring the impact of different mappings. We finally computed for each galaxy the total mass density profile and compared to a variety of data sets from strong and weak lensing as cataloged and analyzed by Sonnenfeld et al. (2013), Newman et al. (2015), and Newman et al. (2013b), and also, where appropriate, with dynamical measurements from Cappellari et al. (2015). Our main results, for galaxies with masses above $M_{star} \gtrsim 2-3 \times 10^{11} M_{\odot}$ (Salpeter IMF), can be summarized as follows.

1. In line with observational evidence, the semi-empirical models naturally predict a non-universal mass weighted total density profile at the effective radius $\gamma'$, independently of the exact input assumptions. Specifically, the models confirm the observational evidence that $\gamma'$ decreases for increasing galaxy size at fixed stellar mass, and for lower mass galaxies at fixed effective radius.

2. The strong lensing data disfavor at $\gtrsim 3\sigma$ a Chabrier IMF and are consistent with a Salpeter or variable IMF.

3. Data also disfavor at $\gtrsim 2-3\sigma$ departures from a Sérsic stellar profile and an uncontracted NFW dark matter profile in the inner regions.

4. Our modeling predicts the total density profile to be roughly isothermal or steeper ($\gamma' \gtrsim 2$) in low-mass halos of $M_{halo} \lesssim 10^{13} M_{\odot}$, but shallower ($\gamma' \sim 1.5$) in cluster environments, with $M_{halo} \gtrsim 3 \times 10^{14} M_{\odot}$. This is mainly a manifestation of structural non-homology, with the relative density distribution of stars and dark matter varying systematically from isolated galaxies to central galaxies in clusters. Despite this substantial decrease with increasing halo mass, the predicted $\gamma'$ are still significantly steeper than the total mass density slopes currently measured in BCGs, which suggest values closer to $\gamma' \sim 1.1$.

5. Adequately reproducing the full total mass density profiles of massive galaxies at the center of massive halos with $M_{halo} \gtrsim 3 \times 10^{14} M_{\odot}$ requires a stellar mass profile flatter than Sérsic within $r \lesssim 3-5$ kpc. An improved match at intermediate scale $10 \lesssim r \lesssim 300$ kpc can instead be obtained by assuming: (1) either a cored NFW dark matter profile, coupled to a median halo concentration a factor of $\sim 2$ higher than the one predicted by N-body numerical simulations; (2) or an Einasto profile with $\mu \sim 4$ ($\alpha \sim 0.25$), coupled to a median halo concentration only a factor of $\lesssim 1.3$ higher than numerical predictions.

F.S. acknowledges valuable discussions with A. Newman, C. Kochanek, A. Dutton, and S. Dye. R.G. acknowledges support from the Centre National des Etudes Spatiales (CNES). T.T. acknowledges support from the Packard Foundation through a Packard Research Fellowship. K.H.C. was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2016R1D1A1B03935804). We thank two anonymous referees for several useful comments and suggestions that greatly improved the quality and presentation of our results.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Concentration vs. halo mass for our cosmology computed from the Diemer & Kravtsov (2015) model (solid line) along with its $1\sigma$ uncertainty (gray region), also plotted at $z=1$ (dotted line). Also shown for comparison the results of the Bullock et al. (2001) concentration model as revised by Macciò et al. (2008) and Dutton & Treu (2014; red and purple long-dashed lines) for the same cosmology, using $K_{200c} = 4$ (see the text for details).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Schematic sketch of the change in stellar mass profile with decreasing Sérsic index $n$, as labeled, for a galaxy at fixed stellar mass and effective radius. The profile tends to flatten out with decreasing $n$ inducing a proportional flattening in $\gamma'$. For comparison, the predicted Hernquist profile is also shown.</figcaption></figure>

## Appendix A
### The Concentration–Mass Relation

The reference concentration-halo mass relation we always use in this paper has been derived by Diemer & Kravtsov (2015). Their results are shown in Figure 11, for our chosen cosmology (Section 1), at $z=0$ (solid line), with the gray area marking the $1\sigma$ uncertainty region of 0.16 dex, and at $z=1$ (dotted line).

We compare the results by Diemer & Kravtsov (2015) with the analytic model put forward by Bullock et al. (2001) and further revised by Macciò et al. (2008). The latter first look for the redshift of collapse $z_c$ of a given halo as the redshift at which the characteristic mass is equal to a fraction $F=0.01$ of the halo mass at the observation redshift $z$. It then computes the concentration for halos defined as 200 times the critical density as

$$c_{200c} = K_{200c} \left[ \frac{H(z_c)}{H(z)} \right]^{2/3} \quad (12)$$

where $H^2(z) = H_0^2[\Omega_\Lambda + \Omega_m(1+z)^3]$, and we set $K_{200c} = 4$ in line with the numerical results by Dutton & Macciò (2014).

The Bullock et al. (2001) and Macciò et al. (2008) analytic model (red and purple long-dashed lines) is in very good agreement with the Diemer & Kravtsov (2015) results, but mostly at low redshifts, while it departs from it at increasing redshifts (dotted line; see also Klypin et al. 2016).

## Appendix B
### The Sérsic 3D Profile

Given a stellar mass, a projected effective radius $R_{2D}$, and a Sérsic index $n$ assigned to an SDSS galaxy from photometric image analysis, we work out the full 3D Sérsic profile following, in the specific, the analytic fits by Lima Neto et al. (1999), though other groups have provided useful approximations (e.g., Prugniel & Simien 1997; Mamon & Łokas 2005a). We first convert the input, empirical 2D radius to a 3D one

using their approximation

$$r = R_{2D}(1.356 - 0.0293m + 0.0023m^2), \quad (13)$$

with $m=1/n$. The 3D mass density profile at any radius $r$ reads then as

$$f(x) = f_0 x^{-p} \exp(-x^m), \quad (14)$$

where

$$p = 1.0 - 0.6097m + 0.05463m^2 \quad (15)$$

and $x=r/a$, $\ln(a) = \ln(R_{2D}) - k$, and $k = [0.6950 - \ln(m)]/m - 0.1789$. The normalization $f_0$ is then fixed by combining Equations (9) and (18) in Lima Neto et al. (1999) as

$$f_0 = \frac{M_{\text{star}}}{4\pi a^3 \Gamma[(3-p)/m]} \quad (16)$$

The cumulative mass $M(<r)$ within a given radius $r$ can be analytically expressed as

$$M(<r) = 0.5M_{\text{star}} \frac{\gamma[(3-p)/m, x^m]}{\gamma[(3-p)/m, (r/a)^m]} \quad (17)$$

with $\gamma$ here indicating the (lower) incomplete gamma function.

An example of the Sérsic profile for different values of the Sérsic index $n$ at fixed stellar mass and effective radius is given in Figure 12, where we also compare it to the Hernquist profile.

## Appendix C
### Alternative Models

We here present a series of measured PDFs as extracted from different semi-empirical models and compared to the data by Sonnenfeld et al. (2013). Each figure has the same format as Figure 3, with the reference Salpeter IMF model, with stellar de Vaucouleurs and dark matter uncontracted NFW profiles, compared to a model with a variable IMF
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Same format as Figure 3 with the Salpeter IMF model compared to the variable IMF one. Both models perform.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Same format as Figure 3 with the standard de Vaucouleurs model compared to a Hernquist stellar profile model. The latter model is less favored by the data, in line with what is inferred from kinematic data (Cappellari et al. 2015, see Figure 2).</figcaption></figure>

(Figure 13), a Hernquist stellar profile (Figure 14), a core in the dark matter profile as in Newman et al. (2013b) of size $r_{core} = 14$ kpc (Figure 15), a combination of Hernquist plus cored ($r_{core} = 14$ kpc) dark matter profiles (Figure 16), or

assuming a reference model characterized by the Moster et al. (2013) stellar mass–halo mass relation (Figure 17). Most (but not all) of these systematic variations to the reference model tend to be disfavored by the data at different significance.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Same format as Figure 3 with the standard NFW model compared to a cored NFW model as in Newman et al. (2013b). The latter model is disfavored, at least at the stellar and halo mass scales probed by Sonnenfeld et al. (2013).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Same format as Figure 3 with the de Vaucouleurs+NFW standard model compared to a stellar Hernquist+cored NFW model. The latter model is disfavored in intermediate-mass galaxies with respect to a standard model with steeper stellar and dark matter profiles.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Same format as Figure 3 with the standard model compared to a model characterized by a flatter stellar mass-halo mass relation (Moster et al. 2013). The latter model tends to be disfavored with respect to a model characterized by a steeper stellar mass-halo mass relation (Figure 1).</figcaption></figure>

References

Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543

Alabi, A. B., Forbes, D. A., Romanowsky, A. J., et al. 2017, arXiv:1701.05904

Amodeo, S., Ettori, S., Capasso, R., & Sereno, M. 2016, A&A, 590, A126

Angulo, R. E., Hahn, O., Ludlow, A., & Bonoli, S. 2016, arXiv:1604.03131

Auger, M. W., Budzynski, J. M., Belokurov, V., Koposov, S. E., & McCarthy, I. G. 2013, MNRAS, 436, 503

Auger, M. W., Treu, T., Bolton, A. S., et al. 2010a, ApJ, 724, 511

Auger, M. W., Treu, T., Gavazzi, R., et al. 2010b, ApJL, 721, L163

Bai, L., Yee, H. K. C., Yan, R., et al. 2014, ApJ, 789, 134

Barausse, E. 2012, MNRAS, 423, 2533

Barausse, E., Shankar, F., Bernardi, M., Dubois, Y., & Sheth, R. K. 2017, arXiv:1702.01762

Barnabè, M., Czoske, O., Koopmans, L. V. E., Treu, T., & Bolton, A. S. 2011, MNRAS, 415, 2215

Barnabè, M., Spiniello, C., Koopmans, L. V. E., et al. 2013, MNRAS, 436, 253

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 149, 289

Bellstedt, S., Lidman, C., Muzzin, A., et al. 2016, MNRAS, 460, 2862

Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, MNRAS, 436, 697

Bernardi, M., Meert, A., Sheth, R. K., et al. 2016a, arXiv:1604.01036

Bernardi, M., Meert, A., Sheth, R. K., et al. 2016b, MNRAS, 455, 4122

Bernardi, M., Meert, A., Sheth, R. K., et al. 2017, MNRAS, tmp https://doi.org/10.1093/mnras/stx176

Bernardi, M., Meert, A., Vikram, V., et al. 2014, MNRAS, 443, 874

Bernardi, M., Roche, N., Shankar, F., & Sheth, R. K. 2011a, MNRAS, 412, 684

Bernardi, M., Roche, N., Shankar, F., & Sheth, R. K. 2011b, MNRAS, 412, L6

Bernardi, M., Shankar, F., Hyde, J. B., et al. 2010, MNRAS, 404, 2087

Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734

Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, ApJ, 301, 27

Bolton, A. S., Brownstein, J. R., Kochanek, C. S., et al. 2012, ApJ, 757, 82

Borriello, A., Salucci, P., & Danese, L. 2003, MNRAS, 341, 1109

Bower, R. G., Lucey, J. R., & Ellis, R. S. 1992, MNRAS, 254, 601

Buchan, S., & Shankar, F. 2016, MNRAS, 462, 2001

Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559

Cappellari, M. 2016, ARA&A, 54, 597

Cappellari, M., Bacon, R., Bureau, M., et al. 2006, MNRAS, 366, 1126

Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2012, Natur, 484, 485

Cappellari, M., Romanowsky, A. J., Brodie, J. P., et al. 2015, ApJL, 804, L21

Cappellari, M., Scott, N., Alatalo, K., et al. 2013, MNRAS, 432, 1709

Cen, R. 2014, ApJL, 790, L24

Chabrier, G. 2003, PASP, 115, 763

Chae, K.-H., Bernardi, M., & Kravtsov, A. V. 2014, MNRAS, 437, 3670

Comerford, J. M., & Natarajan, P. 2007, MNRAS, 379, 190

Conroy, C., & van Dokkum, P. G. 2012, ApJ, 760, 71

Del Popolo, A. 2012, MNRAS, 424, 38

Del Popolo, A. 2014, Journal of Cosmology and Astro-Particle Physics, 7, 019 de Vaucouleurs, G. 1948, AnAp, 11, 247

D'Souza, R., Vegetti, S., & Kauffmann, G. 2015, MNRAS, 454, 4027

Di Cintio, A., Tremmel, M., Governato, F., et al. 2017, arXiv:1701.04410

Diemer, B., & Kravtsov, A. V. 2015, ApJ, 799, 108

Djorgovski, S., & Davis, M. 1987, ApJ, 313, 59

Dressler, A., Lynden-Bell, D., Burstein, D., et al. 1987, ApJ, 313, 42

Dubois, Y., Gavazzi, R., Peirani, S., & Silk, J. 2013, MNRAS, 433, 3297

Dutton, A. A., Conroy, C., van den Bosch, F. C., et al. 2011, MNRAS, 416, 322

Dutton, A. A., Conroy, C., van den Bosch, F. C., Prada, F., & More, S. 2010, MNRAS, 407, 2

Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359

Dutton, A. A., Macciò, A. V., Mendel, J. T., & Simard, L. 2013, MNRAS, 432, 2496

Dutton, A. A., & Treu, T. 2014, MNRAS, 438, 3594

Dutton, A. A., van den Bosch, F. C., Dekel, A., & Courteau, S. 2007, ApJ, 654, 27

Dye, S., Negrello, M., Hopwood, R., et al. 2014, MNRAS, 440, 2013

Einasto, J. 1965, TrAlm, 5, 87

El-Zant, A. A., Freundlich, J., & Combes, F. 2016, MNRAS, 461, 1745

El-Zant, A. A., Hoffman, Y., Primack, J., Combes, F., & Shlosman, I. 2004, ApJL, 607, L75

Ettori, S., Gastaldello, F., Leccardi, A., et al. 2010, A&A, 524, A68

Ferrarese, L., & Merritt, D. 2000, ApJL, 539, L9

Foëx, G., Motta, V., Jullo, E., Limousin, M., & Verdugo, T. 2014, A&A, 572, A19

Gao, L., Navarro, J. F., Cole, S., et al. 2008, MNRAS, 387, 536

Gavazzi, R., Marshall, P. J., Treu, T., & Sonnenfeld, A. 2014, ApJ, 785, 144

Gavazzi, R., Treu, T., Rhodes, J. D., et al. 2007, ApJ, 667, 176

Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJL, 539, L13

Giocoli, C., Bonamigo, M., Limousin, M., et al. 2016, MNRAS, 462, 167

Giocoli, C., Meneghetti, M., Metcalf, R. B., Ettori, S., & Moscardini, L. 2014, MNRAS, 440, 1899

Gnedin, O. Y., Kravtsov, A. V., Klypin, A. A., & Nagai, D. 2004, ApJ, 616, 16
Gu, M., Conroy, C., & Behroozi, P. 2016, arXiv:1602.01099

Guo, H., Zheng, Z., Zehavi, I., et al. 2014, MNRAS, 441, 2398

Hernquist, L. 1990, ApJ, 356, 359

Huang, S., Ho, L. C., Peng, C. Y., Li, Z.-Y., & Barth, A. J. 2013, ApJ, 766, 47

Huertas-Company, M., Aguerri, J. A. L., Bernardi, M., Mei, S., & Sánchez Almeida, J. 2011, A&A, 525, A157

Huertas-Company, M., Shankar, F., Mei, S., et al. 2013, ApJ, 779, 29

Hyde, J. B., & Bernardi, M. 2009, MNRAS, 394, 1978

Jiang, G., & Kochanek, C. S. 2007, ApJ, 671, 1568

Kaplinghat, M., Tulin, S., & Yu, H.-B. 2016, PhRvL, 116, 041302

Katz, H., Lelli, F., McGaugh, S. S., et al. 2017, MNRAS, 466, 1648

Klypin, A. V., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340

Kochanek, C. S., & White, M. 2001, ApJ, 559, 531

Koopmans, L. V. E., Bolton, A., Treu, T., et al. 2009, ApJL, 703, L51

Koopmans, L. V. E., Treu, T., Bolton, A. S., Burles, S., & Moustakas, L. A. 2006, ApJ, 649, 599

Kormendy, J., & Bender, R. 2009, ApJL, 691, L142

Kravtsov, A., Vikhlinin, A., & Meshscheryakov, A. 2014, arXiv:1401.7329

Kravtsov, A. V. 2013, ApJL, 764, L31

La Barbera, F., Ferreras, I., Vazdekis, A., et al. 2013, MNRAS, 433, 3017

Lapi, A., & Cavaliere, A. 2009, ApJL, 695, L125

Lapi, A., & Cavaliere, A. 2011, AdAst, 2011, 903429

Laporte, C. F. P., & White, S. D. M. 2015, MNRAS, 451, 1177

Li, C., & White, S. D. M. 2009, MNRAS, 398, 2177

Lieu, M., Farr, W. M., Betancourt, M., et al. 2017, arXiv:1701.00478

Lima Neto, G. B., Gerbal, D., & Márquez, I. 1999, MNRAS, 309, 481

Ludlow, A. D., & Angulo, R. E. 2017, MNRAS, 465, L84

Ludlow, A. D., Navarro, J. F., Angulo, R. E., et al. 2014, MNRAS, 441, 378

Macciò, A. V., Dutton, A. A., & van den Bosch, F. C. 2008, MNRAS, 391, 1940

Mamon, G. A., & Łokas, E. L. 2005a, MNRAS, 362, 95

Mamon, G. A., & Łokas, E. L. 2005b, MNRAS, 363, 705

Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., & Brinkmann, J. 2006, MNRAS, 368, 715

Martizzi, D., Teyssier, R., & Moore, B. 2013, MNRAS, 432, 1947

Meert, A., Vikram, V., & Bernardi, M. 2015, MNRAS, 446, 3943

Mei, S., Stanford, S. A., Holden, B. P., et al. 2012, ApJ, 754, 141

Meneghetti, M., & Rasia, E. 2013, arXiv:1303.6158

Meneghetti, M., Rasia, E., Vega, J., et al. 2014, ApJ, 797, 34

Merritt, D. 2006, ApJ, 648, 976

Merten, J., Meneghetti, M., Postman, M., et al. 2015, ApJ, 806, 4

Montes, M., Trujillo, I., Prieto, M. A., & Acosta-Pulido, J. A. 2014, MNRAS, 439, 990

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Naab, T., Johansson, P. H., & Ostriker, J. P. 2009, ApJL, 699, L178

Napolitano, N. R., Romanowsky, A. J., & Tortora, C. 2010, MNRAS, 405, 2351

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039

Navarro, J. F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 402, 21

Neistein, E., Li, C., Khochfar, S., et al. 2011, MNRAS, 416, 1486

Newman, A. B., Ellis, R. S., & Treu, T. 2015, ApJ, 814, 26

Newman, A. B., Treu, T., Ellis, R. S., et al. 2013a, ApJ, 765, 24

Newman, A. B., Treu, T., Ellis, R. S., & Sand, D. J. 2013b, ApJ, 765, 25

Nipoti, C. 2015, ApJL, 805, L16

Nipoti, C., Treu, T., Auger, M. W., & Bolton, A. S. 2009a, ApJL, 706, L86

Nipoti, C., Treu, T., & Bolton, A. S. 2009b, ApJ, 703, 1531

Nipoti, C., Treu, T., Ciotti, L., & Stiavelli, M. 2004, MNRAS, 355, 1119

Oguri, M., Rusu, C. E., & Falco, E. E. 2014, MNRAS, 439, 2494

Oldham, L., Auger, M., Fassnacht, C., et al. 2017, MNRAS, 465, 3185

Oldham, L. J., & Auger, M. W. 2016, MNRAS, 457, 421

Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312

Papastergis, E., & Shankar, F. 2016, A&A, 591, A58

Peirani, S., Kay, S., & Silk, J. 2008, A&A, 479, 123

Peralta de Arriba, L., Balcells, M., Falcón-Barroso, J., & Trujillo, I. 2014, MNRAS, 440, 1634

Petrosian, V. 1976, ApJL, 209, L1

Poci, A., Cappellari, M., & McDermid, R. M. 2017, MNRAS, 467, 1397

Posacki, S., Cappellari, M., Treu, T., Pellegrini, S., & Ciotti, L. 2015, MNRAS, 446, 493

Postman, M., Lauer, T. R., Donahue, M., et al. 2012, ApJ, 756, 159

Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, MNRAS, 423, 3018

Prugniel, P., & Simien, F. 1997, A&A, 321, 111

Remus, R.-S., Burkert, A., Dolag, K., et al. 2013, ApJ, 766, 71

Remus, R.-S., Dolag, K., Naab, T., et al. 2017, MNRAS, 464, 3742

Robertson, B., Bullock, J. S., Cox, T. J., et al. 2006, ApJ, 645, 986

Rodríguez-Puebla, A., Drory, N., & Avila-Reese, V. 2012, ApJ, 756, 2

Ruff, A. J., Gavazzi, R., Marshall, P. J., et al. 2011, ApJ, 727, 96

Salpeter, E. E. 1955, ApJ, 121, 161

Sereno, M., Giocoli, C., Ettori, S., & Moscardini, L. 2015, MNRAS, 449, 2024

Serra, P., Oosterloo, T., Cappellari, M., den Heijer, M., & Józsa, G. I. G. 2016, MNRAS, 460, 1382

Sérsic, J. L. 1963, BAAA, 6, 41

Shan, Y., McDonald, M., & Courteau, S. 2015, ApJ, 800, 122

Shankar, F. 2009, NewAR, 53, 57

Shankar, F., & Bernardi, M. 2009, MNRAS, 396, L76

Shankar, F., Bernardi, M., Sheth, R. K., et al. 2016, MNRAS, 460, 3119

Shankar, F., Bernardi, M., & Sheth, R. K. 2017, MNRAS, 466, 4029

Shankar, F., Buchan, S., Rettura, A., et al. 2015, ApJ, 802, 73

Shankar, F., Guo, H., Bouillot, V., et al. 2014a, ApJL, 797, L27

Shankar, F., Marulli, F., Bernardi, M., et al. 2010a, MNRAS, 405, 948

Shankar, F., Marulli, F., Bernardi, M., et al. 2010b, MNRAS, 403, 117

Shankar, F., Marulli, F., Bernardi, M., et al. 2013, MNRAS, 428, 109

Shankar, F., Mei, S., Huertas-Company, M., et al. 2014b, MNRAS, 439, 3189

Shu, Y., Bolton, A. S., Brownstein, J. R., et al. 2015, ApJ, 803, 71

Sonnenfeld, A., Nipoti, C., & Treu, T. 2014, ApJ, 786, 89

Sonnenfeld, A., Treu, T., Gavazzi, R., et al. 2013, ApJ, 777, 98

Sonnenfeld, A., Treu, T., Marshall, P. J., et al. 2015, ApJ, 800, 94

Spiniello, C., Koopmans, L. V. E., Trager, S. C., Czoske, O., & Treu, T. 2011, arXiv:1103.4773

Spiniello, C., Trager, S., Koopmans, L. V. E., & Conroy, C. 2014, MNRAS, 438, 1483

Stadel, J., Potter, D., Moore, B., et al. 2009, MNRAS, 398, L21

Taranu, D., Dubinski, J., & Yee, H. K. C. 2015, ApJ, 803, 78

Taranu, D. S., Dubinski, J., & Yee, H. K. C. 2013, ApJ, 778, 61

Thanjavur, K., Simard, L., Bluck, A. F. L., & Mendel, T. 2016, MNRAS, 459, 44

Thomas, J., Saglia, R. P., Bender, R., et al. 2011, MNRAS, 415, 545

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Tinker, J. L., Brownstein, J. R., Guo, H., et al. 2016, arXiv:1607.04678

Tonini, C., Lapi, A., & Salucci, P. 2006, ApJ, 649, 591

Tortora, C., La Barbera, F., Napolitano, N. R., et al. 2014a, MNRAS, 445, 115

Tortora, C., Napolitano, N. R., Saglia, R. P., et al. 2014b, MNRAS, 445, 162

Treu, T., Auger, M. W., Koopmans, L. V. E., et al. 2010, ApJ, 709, 1195

Treu, T., Gavazzi, R., Gorecki, A., et al. 2009, ApJ, 690, 670

Treu, T., Koopmans, L. V., Bolton, A. S., Burles, S., & Moustakas, L. A. 2006, ApJ, 640, 662

Treu, T., & Koopmans, L. V. E. 2004, ApJ, 611, 739

Treu, T., Stiavelli, M., Bertin, G., Casertano, S., & Møller, P. 2001, MNRAS, 326, 237

Trujillo, I., Burkert, A., & Bell, E. F. 2004, ApJL, 600, L39

Umetsu, K., Zitrin, A., Gruen, D., et al. 2016, ApJ, 821, 116

Watson, D. F., & Conroy, C. 2013, ApJ, 772, 139

Wojtak, R., & Mamon, G. A. 2013, MNRAS, 428, 2407

Xu, D., Springel, V., Sluse, D., et al. 2016, arXiv:1610.07605

Yang, X., Mo, H. J., van den Bosch, F. C., et al. 2007, ApJ, 671, 153

Zhao, D. H., Mo, H. J., Jing, Y. P., & Börner, G. 2003, MNRAS, 339, 12