The central dark matter fraction of massive early-type galaxies

Dark matter (DM) is predicted to be the dominant mass component in galaxies. In the central region of Early-type galaxies it is expected to account for a large amount of the total mass, although the stellar mass should still represent the majority of the mass budget, depending on the stellar Initial Mass Function (IMF). We discuss latest results on the DM fraction and mean DM density for local galaxies and explore their evolution with redshifts in the last 8 Gyr of the cosmic history. We compare these results with expectations from the $\Lambda$CDM model, and discuss the the role of the IMF and galaxy model, through the central total mass density slope. We finally present future perspectives offered by next generation instruments/surveys (Rubin/LSST, Euclid, CSST, WEAVE, 4MOST, DESI), that will provide the unique chance to measure the DM evolution with time for an unprecedented number of galaxies and constrain their evolutionary scenario.


INTRODUCTION
Early-type galaxies (ETGs) are the most massive galaxy systems in the Universe and represent the final stage of galaxy evolution. As such, they carry the fossil record of the assembly of stellar and dark matter through time, and, because they are bright and massive, can be studied in details from low to high redshifts. They are the result of a complex alchemy of physical processes that can be grouped in two major categories: galaxy merging (Trujillo et al. 2006(Trujillo et al. , 2007Cenarro and Trujillo 2009;Naab et al. 2009;Hopkins et al. 2010) and feedback mechanisms (e.g. Dekel and Burkert 2014;Dekel and Birnboim 2006). The variety of evolutionary mechanisms they have experienced make them a crucial test bench of the galaxy formation theories and the ultimate probe of the hierarchical evolutionary scenario. In particular, they can be used to trace the assembly of both the luminous and the dark components across time and check if the observations are consistent with predictions from simulations. For example, the total stellar-to-dark mass ratio of ETGs depends strongly on the galaxy mass, and observational and theoretical studies have suggested this to be connected to the overall star formation efficiency (Benson et al. 2000;Marinoni and Hudson 2002;Napolitano et al. 2005;Mandelbaum et al. 2006;van den Bosch et al. 2007;Conroy and Wechsler 2009;Moster et al. 2010;Alabi et al. 2016).
The ETG centers should experience, on small scales, the effect of the mechanisms acting at the large scales and show similar correlations for the central masses as the ones involving their total masses. In the last two decades, increasing evidences have been collected, showing that the central DM fraction (typically within one effective radius, R e hereafter) is higher in larger and more massive galaxies (e.g. Hyde and Bernardi 2009b;Tortora et al. 2009;Ruszkowski and Springel 2009;Auger et al. 2010a;Napolitano et al. 2010; Thomas et al. 2011;Tortora et al. 2012).
This positive correlation with mass seems almost insensitive to the adopted galaxy mass profile or initial mass function, IMF (e.g., Cardone et al. 2009;Cardone and Tortora 2010;Cardone et al. 2011), but it can become uncertain in case a non-ΛCDM scenario, with mass following the (non-homologous) light distribution, is adopted (e.g., Trujillo et al. 2004;Tortora et al. 2009Tortora et al. , 2012. In general, a general consensus about the existence of such a trend is not yet reached as there are sparse evidences of an anti-correlation of the dark matter fraction with mass (e.g., Grillo et al. 2009;Grillo 2010;Grillo and Gobat 2010).
Among the largest source of uncertainty to quantify the stellar and DM mass budget in the central galactic regions, the IMF remains the most difficult to account for. Direct constraints from gravity sensitive spectral lines (see Spiniello et al. 2012; La Barbera et al. 2013) are strongly required, albeit quite expensive in terms of observational requirements. In absence of these constraints, the simplest approach is to assume a "universal" IMF. Unfortunately, this can introduce unreasonable fluctuations of the stellar mass involved in the models, i.e. by a factor of 2 or more (e.g. assuming a Chabrier 2001 or a Salpeter 1955 IMF or even super-Salpeter IMF, e.g. Tortora et al. 2009). This can strongly affect the conclusions on the central DM fraction in these extreme cases. On the other hand, "non-universality" of the IMF, i.e. the systematic variation of the IMF with mass (or velocity dispersion), from a bottom-lighter (i.e., 'lower-mass') IMF for low mass ETGs to a bottom-heavier (i.e., 'higher-mass') IMF in massive galaxies, has accumulated convincing evidences in the last decade (e.g., Treu et al. 2010;Conroy and van Dokkum 2012;Cappellari et al. 2012, Spiniello et al. 2012Goudfrooij and Kruijssen 2013;La Barbera et al. 2013;Tortora et al. 2013;Martín-Navarro et al. 2015;Li et al. 2017;Domínguez Sánchez et al. 2019), despite, also in this case, contradicting evidences pointing to a bottom-light IMFs for massive systems have been reported (Smith et al. 2015). This non universal IMF can, in principle, incorporate most (if not all) the "apparent" DM fraction trend with mass (e.g., Thomas et al. 2011;Tortora et al. 2013).
A still almost unexplored area, to test the correlation of the dark-matter fraction with mass, is represented by high-redshift (i.e. z > 0.3) galaxies. For these systems, the difficulty of obtaining extended kinematics over large samples of galaxies is a barrier for large statistical studies. Here, one possibility is the use of strong gravitational lensing , Tortora et al. 2010bSonnenfeld et al. 2013). However, at the moment the number of available systems with all the necessary data to perform accurate mass models is still restricted to small samples and small redshift windows, i.e. z ∼ < 0.7 (e.g. Shajib et al. 2020). Another option is to use simple dynamical analysis based on integrated aperture kinematics, as the one provided by large sky surveys, like SDSS/BOSS (e.g. Thomas et al. 2013), GAMA (e.g. Liske et al. 2015) and LAMOST (Napolitano et al. 2020). It has been proven that in these cases Jeans equation analysis can produce results consistent with more sophisticated dynamical modeling , see also Section 2).
The first systematic studies of the evolution of the central DM fraction with redshift has been performed by Beifiori et al. (2014) and Tortora et al. (2014b), which provided evidences that high−z ETGs are less DM dominated than their local counterparts. Only recently, we have updated our previous analysis to include the effect of the IMF in the dynamical analysis and discussed in more details the results in the context of the hierarchical scenario (Tortora et al. 2018b).
The paucity of high−z studies remains a strong limitation in the understanding of the DM content in galaxies. In fact, by studying correlations at higher-redshift we could possibly separate the effect of the DM assembly and the IMF non-universality more clearly, when one or either effects start to emerge or the freedom on the choice of some parameters (e.g. age, metallicity of stars, concentration of the DM haloes, etc.) is lower.
In this paper we present some updated results about the DM content of the central regions of massive ETGs, compare them with other observations and simulations, and discuss some perspectives in the context of the upcoming sky surveys. Despite in the last decade there has been a significant production of studies addressing these specific measurements in ETGs (Padmanabhan et al. 2004;Cappellari et al. 2006Cappellari et al. , 2013bHyde and Bernardi 2009b;Tortora et al. 2009Tortora et al. , 2010aTortora et al. , 2012Tortora et al. , 2014bThomas et al. 2009Thomas et al. , 2011Ruszkowski and Springel 2009;Auger et al. 2010a;Cardone et al. 2009;Cardone and Tortora 2010;Napolitano et al. 2010;Beifiori et al. 2014;Shu et al. 2015;Nigoche-Netro et al. 2016, 2019Xu et al. 2017;Lovell et al. 2018), only few reviews on dark matter in galaxies have touched this specific subject (e.g. Courteau et al. 2014;Salucci 2019), with little emphasis on the insight the dark matter fractions can provide about the ETG structure and evolution.
The paper outline is as follows. In Section 2 we present the procedure we have used to determine the DM content and datasamples. We correlate the DM fraction with galaxy parameters in Section 3, deriving hints on the DM and total mass density slopes in Section 4. The evolution of DM fraction with redshift and the interpretation within the galaxy evolution scenario is provided in Section 5. Alternatives to the ΛCDM are briefly touched on in Section 6. A summary of the results and future prospects are presented in Section 7.
1. The projected luminosity profile I(r) is parameterized by a Sersic (1968) model. 2. The total cumulative (deprojected) dynamical mass profile M (r) is assumed to follow different distributions: a) a constant-M/L profile M (r) = Υ 0 L(r), where L(r) is the cumulative luminosity profile, e.g. from the curve of growth of the surface brightness profile I(r), and Υ 0 is the constant stellar mass-to-light ratio, or b) a singular isothermal sphere (SIS), where M (r) ∝ σ 2 SIS r. In both cases the parameterized profiles have a single free parameter (Υ 0 or σ SIS , respectively). The use of SIS model, in particular, is motivated by evidences from strong gravitational lensing, showing that the total mass density of massive ETGs is close to a power law with slope equal to −2 (e.g. Koopmans et al. 2006. Gavazzi et al. 2007). More complex two-component models can also be adopted, e.g. using DM profiles from N-body simulations, as the NFW double power-law (Navarro et al. 1996), its adiabatically contracted version (e.g. Gnedin et al. 2004), or cored DM models (e.g. Burkert 1995). Exploring different kinds of (more general) models (e.g. Tortora et al. 2007;Cardone et al. 2009Cardone et al. , 2011 can finally help probing how the total mass distribution slope can change in terms of mass and galaxy parameters (e.g. Tortora et al. 2014a).
(3) is projected along the line-of-sight to obtain the projected velocity dispersion: where is the projected density profile.
This is a provisional file, not the final typeset article 5. σ los is integrated within a fixed aperture R Ap to obtain the aperture velocity dispersion, σ Ap using the Equation: where L(R) = R 0 2πsI(s) ds is the luminosity within the projected radius R, and s is the generic coordinate along the los. To avoid lengthy calculations, we have adopted the compact formulae for σ Ap calculated in Equation B7 of Mamon and Lokas (2005) (note that the correct version of Equation B7 is reported in Mamon and Lokas 2006). 6. The modelled σ Ap derived above is fitted to the observed velocity dispersion, σ, by varying the mass model free parameters until the best-fit is achieved. In this paper, we will adopt one-parameter mass models. As matter of fact, even in the case of multi-parametric models, we will reduce these to a single parameter, by using realistic correlations among the model parameters. In particular, as one-parameter models we will adopt the SIS or constant-M/L profiles (e.g. Tortora et al. 2009Tortora et al. , 2012. As a multi-parameter profile we will use NFW+Sérsic for DM+stars, where we fix the NFW dark mass parameters by using independent information from the literature, e.g. the concentration and virial mass by adopting realistic c vir -M vir and also force the the virial mass to obey an observed virial-to-stellar mass relation, M vir -M , where the stellar mass, M , can be a free parameter (the impact of these assumptions have been extensively discussed in Tortora et al. 2013Tortora et al. , 2014a. In the latter case, this allows us to have freedom on the stellar mass-to-light ratio and evaluate the IMF normalization. In this paper we will adopt as reference results the ones derived in Tortora et al. (2013), where the DM halo is set using the c vir -M vir correlation from N-body simulations  and M vir -M correlations from abundance matching results in Moster et al. (2010). This reduction of the parameter space is obliged, in our case, by having a single observable for each galaxy. However, the robustness of such approaches have been demonstrated in more complex analysis using spatially resolved kinematics (e.g. Cappellari et al. 2013b) and combining dynamics with strong and weak gravitational lensing (e.g. Auger et al. 2010b).
7. The resulting best-fit mass profile then provides the total spherical M dyn (r).
The dynamical procedure described above does not take into account the contribution of the black hole (BH) mass, orbital anisotropy and rotation. We have estimated that these latter, in most of the cases, do not impact the M dyn (r) estimates for more than few per cents (see more in Tortora et al. 2009Tortora et al. , 2012Tortora et al. , 2018b.

Central dark matter
The main prediction of the virial formula (Eq. 1) and Jeans analysis (Section 2.1) is M dyn (r), which, under the reasonable assumption of no cold gas and dust in the galaxy centers, is made up only from the stellar mass, M * and the dark mass, M DM , at a given radius r. Hence, if an accurate and independent estimate of M * (r) is available, one can quantify the dark mass just by M DM (r) = M dyn (r) − M * (r) in the galaxy regions probed by the observed σ. Thus, the DM content can be characterized by computing the following quantities. First, the de-projected DM fraction, f DM (r) = 1 − M (r)/M dyn (r) or, equivalently, the total-to-stellar mass ratio M dyn (r)/M (r), which avoids the negative values arising in the DM fractions when M (r) > M dyn (r). Second, the deprojected average DM density, ρ DM = (M dyn (r) − M (r))/( 4 3 πr 3 ), to probe the local average density of DM, indipendently of the slope of the mass density profile from the center to r. Hereafter, we will also refer to M dyn /M (r = R e ) ( ρ DM (r = R e ) ), i.e. the DM fraction (density) computed within a de-projected radius equal to the projected R e , as the "central" DM fraction (density).

Datasamples
In order to apply the method illustrated in Section 2.1 and derive the quantities in Section 2.2, one need datasets for which accurate surface photometry (e.g. high quality imaging in one or more optical and/or near-infrared, NIR, bands) and accurate stellar masses (e.g from optical + NIR multi-band stellar population synthesis models or from spectroscopy) are available. In case one wants to study the evolution with redshift of the DM content, these datasets should provide uniform measurements from z = 0 to some higher redshits. For our analysis we rely on two datasamples with these properties.
Our reference sample is made of massive and red galaxies collected in Tortora et al. (2018b), including stellar masses and structural parameters derived from Kilo Degree Survey (KiDS) imaging (de Jong et al. 2015(de Jong et al. , 2017. This sample encompasses a quite large range of redshifts and it is suitable to study the dark and luminous properties of galaxies from z ∼ 0.1 to z ∼ 0.7. The dataset consists of ∼ 9700 selected ETGs with optical and NIR photometry information. The sample is complete in stellar mass, obtained by the SED fitting of the KiDS optical bands and assuming a Chabrier (2001), at log M /M > 11.2 and redshifts z < 0.7. Structural parameters (e.g., effective radius R e and Sérsic index n) are obtained by a PSF-convolved Sérsic fit of the KiDS galaxy images in g-, rand i-bands (using 2DPHOT, La Barbera et al. 2008). In our calculations, these are rest-framed, as described in Tortora et al. (2018b). Spectroscopic redshifts and central velocity dispersions are taken from SDSS-DR7 (z < 0.2 Abazajian et al. 2009) and BOSS (z > 0.2, Thomas et al. 2013). We refer the reader to Tortora et al. (2018b) for more details about data selection and analysis.
We also use an external properly selected sample of local galaxies with analogous photometry and spectroscopy information, as control sample at z = 0. This is made of ∼ 4300 giant ETGs drawn from the SPIDER project (La Barbera et al. 2010;Tortora et al. 2012). It includes stellar masses derived from SED fitting their optical and NIR photometry (Swindle et al. 2011) using a Chabrier (2001 IMF. This dataset also includes galaxy structural parameters, derived from g through K wavebands, and the SDSS central-aperture velocity dispersions, SPIDER ETGs are defined as luminous bulge-dominated systems, featuring passive spectra in the central SDSS fibre aperture (La Barbera et al. 2010).
There are other datasets that have been used to perform similar studies. In the following, we will compare our results with this independent literature, when available.

CORRELATION WITH STRUCTURAL PARAMETERS AND MASS PROBES
To characterize the central DM content in massive ETGs, we start looking for correlations among M dyn /M and galaxy parameters.

Results for KiDS and SPIDER datasamples
We will consider first the results found using the KiDS dataset and compare these with results obtained on the SPIDER dataset at z ∼ 0. In Figure 1 we show M dyn /M as a function of effective  Figure 1. The total-to-stellar mass ratio M dyn /M within effective radius, R e , assuming a Chabrier IMF, as a function of R e (left panel), total stellar mass M (middle panel) and velocity dispersion within effective radius σ e (right panel). Lines correspond to medians, and shaded regions or error bars to 25-75th percentile. Red squares with error bars are for the whole KiDS sample. Dark blue line and light blue region are for SPIDER galaxies with M > 10 11.2 M . Red and blue dashed lines are medians for KiDS and SPIDER galaxies adopting a Salpeter IMF. Green squares and error bars are for SLACS lenses from Auger et al. (2010a). The cyan lines in the middle panel are from IllustrisTNG simulated galaxies in Lovell et al. (2018): solid line is from the reference IllustrisTNG, while dashed one is for galaxies which, simulated within the full-physics IllustrisTNG simulation, are placed in their corresponding DM haloes simulated within the DM-only simulation. In the right panel, black open square and error bars are for the results obtained applying a Schwarzschild's orbit superposition technique to ETGs in Thomas et al. (2011). radius, stellar mass and velocity dispersion, We fix the IMF to the Chabrier one, which implies that the measured M dyn /M trends correspond to a variation in the DM content. Red symbols represent the estimates of the KiDS sample with redshift z < 0.7 and masses log M /M > 11.2. Dashed blue lines with light blue shaded regions represent the SPIDER sample, assuming g-band structural parameters. Error bars and the shaded regions are the 25-75th per cent quantiles. All the correlations are significant at more than 99 per cent. The best-fitted slopes of the correlations are also reported. Figure 1 shows a tight and positive correlation between (the logartithm of) M dyn /M and R e with a slope α = 0.72. This can be interpreted as a physical aperture effect, where a larger R e subtends a larger portion of a galaxy DM halo (e.g. Napolitano et al. 2010Napolitano et al. , 2011. In practice, larger R e are found in galaxies with larger stellar mass; however, being these massive ETGs characterized by a steep halo-to-stellar mass relation (e.g. Moster et al. 2010), the halo DM mass is increasing too. The net effect is the observed positive correlation between M dyn /M and R e . A similar correlation can also be found for the Sérsic index (see Tortora et al. 2018b), due to the positive correlation existing between the n−index and effective radius (e.g. Tortora et al. 2012). However, if we plot M dyn /M as a function of R e and bin the galaxies in terms of the Sérsic index, then we see that there is a negligible dependence on n (this result is not shown in the figure for brevity). Because smaller effective radii correspond to higher stellar densities, this correlation with R e also translates to a sharp anti-correlation between DM content and central average stellar density, which has been reported for the first time in Tortora et al. (2012) and subsequently confirmed in Tortora et al. (2018b). Galaxies with the smallest R e (∼ 5 kpc) have the smallest DM fraction (∼ 35 per cent), while the largest galaxies (R e ∼ 35 kpc) present the largest DM content (∼ 85 per cent).  Figure 2. The total-to-stellar mass ratio M dyn /M within rest-frame effective radius, R e , as a function of rest-frame effective radius R e , assuming a Chabrier IMF (left panel) and Salpeter and variable IMF (right panel). Red and blue symbols are for KiDS and SPIDER, respectively. Following Figure 1, red squares with errors bars and solid lines with shaded regions adopt a fixed IMF (Chabrier and Salpeter in the two panels). Short-dashed-dotted is for IMF variable from Tortora et al. (2013) and Tortora et al. (2014a) with log M /M > 11.2, while long-dashed-dotted with log M /M > 10.5. The gray region is from Illustris TNG (Wang et al. 2020), extracted from Figure 15 in Shajib et al. (2020); note that in that figure they only show results for R e ∼ < 10 kpc. Orange triangle is from the lensing analysis in Shajib et al. (2020), in a good agreement with Auger et al. (2010a) results, plotted in green.
Moving to the correlations of the M dyn /M with "mass" parameters, we find a shallow correlation with M , with an average M dyn /M ∼ 3 (i.e. 67 per cent of DM) and slope of 0.11. We also find that M dyn /M correlates with σ e (M dyn /M ∝ σ e 0.89 ). It can be shown that the strong correlations with R e and those with σ e translates to a strong positive correlation with M dyn (Tortora et al. 2018b).
The impact of the IMF is also shown: M dyn /M for a Salpeter IMF are plotted as dashed lines for both the KiDS and SPIDER samples, obtained by multiplying the stellar mass by a factor of 1.8 ). The change of IMF move downward the M dyn /M of a factor ∼ 0.25 dex. On average, the DM fractions are positive in these very massive galaxies, with f DM close to 0 only in the smallest galaxies.
These trends are consistent with the ones found at z ∼ 0 using SPIDER galaxies (blue lines with shaded regions; Tortora et al. 2012), also confirming previous results reported in literature at z ∼ 0 (e.g., Cappellari et al. 2006;Hyde and Bernardi 2009a;Tortora et al. 2009Tortora et al. , 2012Napolitano et al. 2010;Nigoche-Netro et al. 2016;Lovell et al. 2018), or at intermediate redshift (Tortora et al. 2010b;Auger et al. 2010a;Tortora et al. 2014b). The evolution with redshifts will be discussed in the next Sections.

Comparison with the literature
We compare our results with M dyn /M estimates from gravitational lensing of the SLACS sample with log M /M > 11.2 and average redshift of z ∼ 0.2 (Auger et al. , 2010a. For a fair comparison, we have homogenized the lensing results, and derived medians and 25-75th percentiles. These are shown in Figure 1 with green symbols. A good agreement is clearly seen for the M dyn /M -R e correlation, while we notice that at fixed M and σ e SLACS M dyn /M are smaller of ∼ 0.1 − 0.2 dex than the median KIDS relation. However, we need to remark that, at fixed M , the SLACS sizes are smaller than the ones of the KiDS sample by ∼ 0.15 dex, while velocity dispersion are ∼ 0.03 dex higher, which implies than that SLACS M dyn /M are smaller of ∼ 0.1 dex within their R e . The smaller sizes of SLACS galaxies are also clear from the M dyn /M -R e correlation, where SLACS galaxies have sizes concentrated towards smaller values, with respect to the range of sizes of SPIDER and KiDS datasamples (see also Tortora et al. 2018b).
In the middle panel of Figure 1 we also compare our M dyn /M with the results from the state-ofthe-art IllustrisTNG simulations from Lovell et al. (2018), which adopt a Chabrier IMF (cyan lines). The agreement with both SPIDER and KiDS galaxies is quite good, considered the uncertainties in observations and simulations.
In the right panel of the same figure we also show that the M dyn /M -σ e correlation for KiDS galaxies, and in particular for the SPIDER sample, is fairly consistent with the DM fractions from Schwarzschild's orbit superposition models in axisymmetric potentials in Thomas et al. (2011) applied to a sample of 16 COMA ETGs.
In Figure 2 we expand the comparison of the M dyn /M -R e with further theoretical and observational results, and investigate the impact of a non-Universal IMF. In particular, in the left panel we confirm one of the results presented in Figure 1 thanks to the comparison with Lovell et al. (2018): our DM fractions are fairly well consistent with the expectations from Wang et al. (2020), who also use IllustrisTNG. IllustrisTNG sample also includes galaxies with log M /M < 11.2 and have smaller sizes, for this reason we find an overlap with the smallest galaxies in our sample. This result provide a further confirmation that the prescriptions adopted by IllustrisTNG are able to realistically provide a quite good agreement with observations. Moreover, the comparison with state-of-the-art observations is presented in the right-hand panel of Figure 2, where we also compare with the recent results from Shajib et al. (2020), who determine the mass density slope and DM fraction of a sample of SLACS lenses, using their strong lensing data, velocity dispersions and weak lensing constraints. These results are plotted as orange symbols. Shajib et al. (2020) show that most of the totality of SLACS galaxies are well fitted by a NFW profile with a Salpeter IMF. For this reason, we converted our M dyn /M to Salpeter-based versions in this panel. As a consistency check we also overplot the results from Auger et al. (2010a), using the same IMF and find a striking consistency of the two analyses.

The intruder: the impact of Initial Mass Function
In Figure 1 we have shown that, on average the DM fractions are positive, independently of the IMF adopted. However the f DM estimates present single cases where galaxies have negative f DM < 0, implying unphysical M dyn (R e ) < M (R e ). In the KiDS dataset, only the ∼ 6 per cent of the galaxies show negative DM fractions, if a Chabrier IMF is adopted. On the other hand, using a Salpeter IMF, produce, by definition, smaller DM fraction, and for ∼ 23 per cent even negative (unphysical) values. The fraction of such negative f DM is larger at higher redshifts: this means that if the Salpeter IMF is an inappropriate choice, this is even worse at higher redshift. This is a well known critical effect that has been discussed in previous works (see e.g. Tortora et al. 2009;Napolitano et al. 2010;Tortora et al. 2012). Despite a significant fraction of these negative f DM is, in principle, compatible with observational scatter in M and M dyn (see Napolitano et al. 2010), in order to avoid unphysical results, one is left with no complete freedom on the assumption of the IMF to adopt. In particular, a higher stellar M/L normalization is unphysical for those systems which tend to have smaller f DM (e.g. the ones with smaller sizes and dynamical masses, larger stellar densities, higher redshift, etc.).
These indications cope with a large number of observational studies that in the last decade have suggested that IMF in ETGs is different from the one estimated by star counts in our Milky Way, where the standard forms for the IMF were identified (e.g. Kroupa 2001;Chabrier 2003). Now, this hypothesis of IMF universality is questioned by different lines of observational evidence, using completely independent data, as spectral features, galaxy dynamics and gravitational lensing (e.g., Treu et al. 2010;Conroy and van Dokkum 2012;Cappellari et al. 2012Cappellari et al. , 2013aSpiniello et al. 2012;Goudfrooij and Kruijssen 2013;La Barbera et al. 2013;Tortora et al. 2013Tortora et al. , 2014aTortora et al. ,c, 2018aMcDermid et al. 2014;Martín-Navarro et al. 2015;Corsini et al. 2017;Li et al. 2017). Among the others, using SPIDER data, we have found a larger stellar mass in the most massive galaxies (high velocity dispersion) than that provided by a Milky-Way Chabrier IMF, that can be translated to a bottom-heavy IMF (Salpeter-like) or a larger dwarf-to-giant star ratio (Tortora et al. 2013). On the contrary, in ETGs with a low-velocity dispersion, the IMF resembles the one that is found in the Milky-Way. We have found that this very tight correlation with velocity dispersion is safe independently of the galaxy model adopted, and also in alternative gravity scenario, in which DM is not included, as MOND (Tortora et al. 2014c) and Emergent Gravity (Tortora et al. 2018a). Despite these positive strong trends with σ and stellar mass (when this last takes into account of the dynamically inferred IMF), milder trends are found with R e , n and M when a fixed IMF is adopted for this latter (Tortora et al. 2014a).
If, a strongly varying IMF as a function of velocity dispersion is translated into a central DM fraction, we find a fairly universal DM fraction, consistent with f DM ∼ 0.2 for a NFW profile DM halo 1 with no further baryonic effects (e.g. adiabatic contraction, AC, Gnedin et al. 2004), or ∼ 0.4 if one accounts the physics of baryon collapse, e.g. via AC (Tortora et al. 2013). The positive correlation of DM fraction with mass has been typically invoked as the driver for the "tilt" of the ETG fundamental plane ), but it turns out that this latter can also be explained, also if partially, by a realistic DM halo model and a non-universal IMF.
Here, for the first time, in the right-panel of Figure 2, we can show that, instead, the positive trend between M dyn /M and R e still survives if a non-universal IMF is adopted. And this is expected considering that, unlike the trend with σ, IMF varies mildly with R e . Short-dashed-dotted is for IMF variable from Tortora et al. (2013) and Tortora et al. (2014a) with log M /M > 11.2, while long-dashed-dotted with log M /M > 10.5. The M dyn /M values found in this way are slightly larger than those made with SPIDER and KiDS, and the results from lensing; however, at least a part of these discrepancies can be also ascribed by the fact that K-band is used, and thus using g-band the curves should move towards right.

Tortora & Napolitano
Dark matter in early-type galaxies  tricky because obtaining inferences on the total DM halo from the measurements in the center is impracticable. The only possible way is to assume some standard recipe for the DM halo and try to match the predictions from these recipes with the central DM inferences (see e.g. Tortora et al. 2009Tortora et al. , 2012Tortora et al. , 2014bNapolitano et al. 2010).
There are different mechanisms that can affect the central DM content and decouple it from the overall DM. In fact, the central f DM can also reflect the local conditions and thus the environment density at the time of initial halo collapse. Moreover, it is well known that baryons interact with DM, changing its distribution (as in the case of AC, Gnedin et al. 2004). From a practical point of view, the DM fraction is somewhat more strongly dependent on the particular values of R e for the stars rather than the DM properties directly. These effects imply strong degeneracies among the parameters in the galaxy models (see e.g. Napolitano et al. 2010), and providing a quantitative comparison to cosmological theory is a necessary but difficult task.
One direct way to investigate the DM halos and study their properties is to define its average density within some small radius, i.e., ρ DM . Following Napolitano et al. (2010) and Tortora et al. (2010b), we focus on correlations of the average DM density with galaxy size and stellar mass for SPIDER and KiDS datasamples. Figure 3 shows that ρ DM strongly anti-correlates with R e . Considering the aperture effect and assuming DM halo homogeneity, this implies that we are measuring a mean DM density profile with radius. In fact, as firstly discussed in Napolitano et al. (2010), if we assume a power-law for the DM profile, i.e. ρ DM (r) ∝ r α , with α negative, then, with a little of algebra, for α > −3 one finds M DM ∝ r α+3 and finally ρ DM (r = R e ) ∝ R e α . Hence, the trend of ρ DM (R e ) with R e , calculated on a sample of galaxies, provides the average slope of the DM profile in the central regions of that galaxy sample. For an NFW halo, near R e is α ∼ −1.1, −1.3. AC makes the DM cuspier, with α ∼ −1.6, −1.9. For KiDS and SPIDER, the slope of this correlation is ∼ −1.8, −1.9, consistently with what originally shown by Napolitano et al. (2010), using a completely different datasample. Therefore, this steep slope can be indicative of cuspier-than-NFW halos, perhaps as induced by AC. But this result can be driven by the choice of a Chabrier IMF. Hence, we cannot exclude that a standard cuspy NFW model and a larger amount of stars due to a Salpeter IMF (e.g., Shajib et al. 2020). We have verified that using an alternative constant-M/L profile yields similar results still fully consistent with a cuspy halo (see also a detailed discussion in Napolitano et al. 2010).
A similar negative correlation is also obtained when the average densities are plotted as a function of stellar mass. The most massive galaxies are characterized by systematically lower values of the DM density; in fact, at increasingly higher masses, size get bigger and bigger, lowering the DM average density. In both the panels, the more distant galaxies in the KiDS sample have, at fixed R e or M larger average DM densities, also because of the systematically smaller sizes in these galaxies. The results for SPIDER galaxies are almost unchanged if the lower-mass cut of log M /M = 11.2 dex is removed (dashed lines).
In the bottom panels of Figure 3 we focus on the local SPIDER datasample and investigate the impact of different assumptions for galaxy model. Differently from the previous plots, we use here the K-band structural parameters to calculate the DM fractions. Consistently with the results in Tortora et al. (2009) and Tortora et al. (2012), we find that if we model the total mass distribution with a constant-M/L (long-dashed line), then the DM content is decreased with respect to an isothermal profile (short-dashed line). We notice that the trend with radius is steeper than the SIS profile when the constant-M/L is adopted, pointing to a slope of ∼ −2.3.
We can now do a step forward, by comparing these correlations with the average DM densities obtained fitting a NFW. We recall that, in this case, we leave the IMF free to change, and assume a realistic c vir -M vir and M vir -M correlations, to constrain the shape of the total mass profile at ∼ < R e . In this case we find that in order to match the DM distribution of a NFW profile the total mass profiles of the most massive galaxies are well reproduced by an isothermal profile, while lower mass galaxies are better fitted by a constant-M/L, this implying a steeper mass density in the central regions. This result implies a DM non-homology in the central galaxy regions, i.e. the total mass density slope in ETGs is not universal and is a function of the stellar mass, a result that we had suggested in Tortora et al. (2009), and demonstrated in more details in Tortora et al. (2014a). These results have been confirmed by independent observations and simulations (e.g., Remus et al. 2013, Dutton and Treu 2014, Wang et al. 2020. The steepest slopes are found at M ∼ 3 × 10 10 M , while isothermal profiles are common in the most massive galaxies. This trend with mass may be related to a varying role of dissipation and galaxy mergers with galaxy mass (see also Tortora et al. 2019).

EVOLUTION WITH REDSHIFT
One firm evidence of the hierarchical scenario in a ΛCDM cosmology is the size and mass growth of galaxies with time, after the Big Bang (Daddi et al. 2005;Trujillo et al. 2006Trujillo et al. , 2007Trujillo et al. , 2011Saglia et al. 2010;Tortora et al. 2014b, 2018b, Roy et al. 2018. These observations all seem to exclude a simple monolithic-like scenario, according to which the bulk of the stars is formed in a single dissipative event followed by a passive evolution, while they seem to support a scenario where galaxy evolution is mainly driven by merging, even though of different kind. While the understanding of the stellar component of galaxies has been supported by a multiplicity of observational probes, basically measuring the properties of the light distribution in galaxies (e.g., La Barbera et al. 2010, Tortora et al. 2016, 2018cRoy et al. 2018 and reference therein) and their stellar populations (e.g., Kauffmann et al. 2003;Gallazzi et al. 2005;Renzini 2006;Sánchez et al. 2012;Maraston et al. 2013), evidences about the growth of the DM with time come mainly from theory (e.g., Conroy and Wechsler 2009;Moster et al. 2010;Lovell et al. 2018) but there are yet little analyses trying to systematically study the evolution of size and mass of DM halo with redshift.

Evolution of the dark matter fractions at fixed mass
In the previous paragraphs we have seen that the M dyn /M , as well as the f DM , are tightly correlated to observables like the velocity dispersion and the (stellar or dark) masses. When studying the evolution of the galaxy DM fractions with redshift, one needs to estimate how much of this evolution is driven by the correlated quantities and how much by the actual mass assembly of galaxies. In Figure 4 we show the dependence on redshift of R e , σ e , M dyn /M and ρ DM , at fixed stellar mass, using the KiDS sample as reference 2 . With respect to Tortora et al. (2018b), we show here the results for all the galaxies with log M /M > 11.2. First, we clearly see the effect of the size growth of galaxies, as the galaxy effective radii are smaller at earlier epochs and increase toward z = 0 (Daddi et al. 2005;Trujillo et al. 2006Trujillo et al. , 2007Buitrago et al. 2008;van der Wel et al. 2008 see also Roy et al. 2018, for further details about size evolution in KiDS galaxies). To quantify this size-redshift relation, we can use a standard formula R e = R e,0 (1 + z) α to fit the data. In the case of no progenitor bias correction (i.e., assuming that the same kind of galaxies have evolved from one redshift bin to another, red filled squares with error bars), we find a slope, α = −2.4. If the progenitor bias is taken into account (i.e. galaxies observed today as passive systems might be active in earlier redshift bin, open squares with dashed red lines), the slope becomes shallower, i.e. α = −1.8. This is because active/disk systems tend to have larger sizes than passive/spheroids systems at high−z, hence increasing the overall sizes at high−z and diluting the size-redshift relation. Overall, these measured slopes are steeper than other literature results for spheroid-like systems with M > 10 11 M (solid black line in the top panels in Figure 4; Trujillo et al. 2007;Conselice 2014). However, at low−z we find a good agreement with the R e s from the local sample from SPIDER.
In the panel b) of Figure 4, we plot the effective velocity dispersion, σ e , as a function of the redshift. Here the evolution with redshift looks shallower than the one shown by the R e , with higher−z galaxies having slightly larger velocity dispersions than the lower−z ones (Cenarro and Trujillo 2009;Posti et al. 2014). We can quantify the σ e −z evolution using the relation σ e = σ e,0 (1 + z) α . The estimated slopes are α = 0.48 (without progenitor bias) and α = 0.73 (with progenitor bias). In general, we find a good agreement with local (La Barbera et al. 2010;Tortora et al. 2012), intermediate-z (Beifiori et al. 2014) and higher-z (Saglia et al. 2010;Tortora et al. 2014b) measures.
The total-to-stellar mass ratio (assuming a Chabrier IMF) is plotted in the panel c) of Figure 4. In this case we see that galaxies show a dominating DM content in their R e s at lower redshift (i.e. 75-80% of DM at z ∼ 0.2), while the f DM are smaller at higher-z (40-50% at z ∼ 0.6). In this 2 We refer the reader to Tortora et al. (2018b) for the analysis of systematics and for a full description of data selection and analysis. Open black square with error bar is median and 25-75th percentiles for SPIDER galaxies, while open black square is the median for the SPIDER sample when the progenitor bias is taken into account. The black solid line in the top panel is taken from the average R e /R e,0 -z trends (with R e,0 = R e (z ∼ 0)) for spheroid-like galaxies in Trujillo et al. (2007) normalized to R e,0 = 15 kpc. Shaded gray region and green line are our predictions using the merging model of Hopkins et al. (2009) and the "puffing-up" scenario from Fan et al. (2008), respectively. The cyan line in the bottom-left panel is for galaxies with a stellar mass of ∼ 10 11 M from IllustrisTNG simulations (Lovell et al. 2018). See explanations in the main text.
case we can use the relation M dyn /M = (M dyn /M ) 0 (1 + z) α to quantify the dependence of this quantity on the redshift and find α = −2.4 (without progenitor bias) and α = −1.3 (with progenitor bias).
Going to the theoretical interpretation of these observed trends, we can compare them with the predictions from two different scenarios, typically adopted to explain the galaxy size evolution. First, the merging scenario (MS, hereafter), which predicts that size growth is driven by the accretion of matter as sizes of the merger remnants are larger than those of their progenitors. The merging model of Hopkins et al. (2009) also predicts that the velocity dispersion decreases as a consequence of the size growth as the relation σ (z) ∝ (1 + γ) −1/2 γ + R e (0)/R e (z) holds, where the parameter γ sets the DM contribution to the potential relative to that of the baryonic mass. This γ parameter is expected to vary between 1 and 2 (which are the best fitted values for M ∼ 10 11 and ∼ 10 12 M , respectively). Second, the "puffing-up" scenario (PS, hereafter) from Fan et al. (2008), which predicts that galaxies grow by the effect of quasar feedback, which removes cold gas from the central regions, quenching the star formation and increasing the size of the galaxy. This model predicts also that velocity dispersion increases with increasing redshift as σ (z) ∝ R e −1/2 .
In order to check these predictions against the observed trends we need to derive the expected redshift dependences of velocity dispersions and M dyn /M on the redshift, starting from the R e -z relation from KiDS median values in panels a). We have used the best-fit R e −z relations discussed above and inserted into the σ (z) equation to derive the expected velocity dispersion in the two schemes. To obtain the M dyn /M (z) relation we need to translate the predicted velocity dispersions into a M dyn /M . This is done by solving the spherical Jeans equation (see Eq. 2), which includes 1) the 3D luminous density profile, and 2) the total potential, as a function of redshift. For the light distribution we have assumed a Sérsic profile with n = 4 for simplicity (i.e. a pure de Vaucouleurs) and the effective radius given by our interpolated R e (z) relation as defined above. For the total potential we have used the SIS model. By imposing that the velocity dispersion inferred by the Jeans equation (averaged within R e ) equals the σ(z) in the two scenarios, MS and PS, we can obtain M dyn and M dyn /M as a function of redshift. We remark here that in this calculation the relevant information we are interested on is the trend with redshift and not the normalization, which can be adjusted by hand, since in the σ (z) formulae the normalization factor is unspecified.
The predicted trends for σ e and M dyn /M are plotted in the panels b) and c) of Figure 4. Here, the PS predicts a too strong evolution (with a variation of ∼ 100 km s −1 in the redshift window analyzed), which unfits the KiDS results for both the σ e and M dyn /M . On the other hand, the milder evolution predicted from MS matches closer the observations (Cenarro and Trujillo 2009). However, even though the agreement with the velocity dispersion seem very good, the PS model predict a too shallow M dyn /M -z trend with respet to the observed one.
To complement this analysis, we also compare the M dyn /M -z trend with IllustrisTNG expectations from Lovell et al. (2018). Although the agreement with IllustrisTNG is quite good in terms of normalization, as confirmed by the analysis of the trends with M and R e (Figures 1 and 2), the simulations predict a change with redshift which is milder than what we find with KiDS data.
For the first time, in this paper we also present the evolution with z of the average DM density. This is shown in panel d) of Figure 4, where for simplicity we limit to show the results for KiDS and SPIDER datasamples, without comparing with any toy-model. The average DM density within R e reaches its largest values at z ∼ 0.6 (< log ρ DM (R e )/(M /kpc 3 ) >∼ 7.75 ) and decreases systematically dowm to < log ρ DM (R e )/(M /kpc 3 ) >∼ 7 at the lowest redshifts. This trend is explained by the fact that sizes are larger at smaller redshifts, and the impact on the denominator is stronger than that of the numerator where DM mass is larger at smaller z.

The evolution of Size-and Dark matter-mass relations: constraints on merger scenario
In the previous section we have found evidence that simple recipes for the size growth and the velocity dispersion evolution in a merging scenario reproduce better the observed evolution of the dark-to-stellar mass ratio, at least for massive ETGs. In this section we want to check further the merging scenario and see whether we can gain more insight on the mechanisms driving the size and mass growth in massive ETGs. In particular, we want to focus on two scaling relations discussed in previous sections, the R e -M and M dyn /M -M , and figure how galaxies evolve in this parameter space in response of the joint evolution of size and mass predicted by different kinds of merging events.
Indeed, dissipationless major mergers from simulations of elliptical galaxies have predicted that the DM fraction within a certain physical radius decreases mildly after the merger (Boylan-Kolchin et al. 2005). But they have also shown that the DM fraction within the final R e is greater than the DM fraction within the initial R e , because the total mass within R e , M tot (R e ), changes, after the merger, more than M (R e ). The problem has been investigated in detail with hydrodynamic simulations by Hilz et al. (2013) who find that the equal-mass mergers produce a smaller size increase of multiple minor mergers. In particular, the variation of R e with respect to the initial radius, R e /R 0 , in terms of the variation of M with respect to the initial stellar mass, M /M 0 , is found to be R e /R 0 ∝ (M /M 0 ) 0.91 for the equal-mass merger and ∝ (M /M 0 ) 2.4 for the minor mergers.
We can test all these predictions against the R e -M and M dyn /M -M relations (assuming a fixed Chabrier IMF) for KiDS galaxies at different redshift, as shown in Figure 5. Blue, green and red lines are for KiDS galaxies in three redshift bins 0.1 < z ≤ 0.3, 0.3 < z ≤ 0.5 and 0.5 < z ≤ 0.7.
Here we see that lower redshifts galaxies are larger in size and contain more DM in their cores at all values of M , confirming the trends in Figure 4. The trend is weaker if we consider the progenitor bias, which affects mostly the lowest redshift bin (dashed lines in Figure 4). We will check if mass, size and total-to-stellar mass evolution in KiDS galaxies can be explained, consistently, through major or minor mergers.
To reproduce these trends we need to account for the dark mass inside the effective radius, as we have to reproduce the M dyn /M , where M dyn = M + M DM and the evolution of the R e /R 0 vs. M /M 0 is incorporated in the relations discussed above. To do that, we have constructed some simplified toy-models assuming that M DM ∝ M vir r η around R e , with η ∼ 2 for a standard NFW and η ∼ 1.2 for a contracted NFW, hereafter AC+NFW (according with Boylan-Kolchin et al. 2005). We start assuming that the systems participating to merging (i.e. the progenitors) all have the same M vir /M , which is reasonable for most of the stellar mass range covered by our sample, i.e. δM vir ≈ δM . However, we cannot exclude for the minor merging case that the virial mass changes at a different rate of the stellar mass, since the main galaxy is merging with another galaxy with a different M vir /M .
Having all ingredients set, we have considered the evolution tracks related to the two different merging types: for the major merging we take the average evolution of R e in terms of M evolution for the equal-mass merging, i.e. R e /R 0 ∝ (M /M 0 ) 0.91 ; for minor merging we take R e /R 0 ∝ (M /M 0 ) 2.4 .
In Figure 5 we show the evolution tracks for a progenitor galaxy with log M /M = 11.4. Dots and squares are for NFW and AC+NFW profiles, respectively. In the left panels, the major merging tracks are shown as black lines with dots/squares indicating the events corresponding to masses δM × M ,0 , with the mass increments δM = 1, 2, 4, ..., assuming that δM vir = δM . The first dot/square at log M /M = 11.4 corresponds to the initial progenitor galaxy. The second dot/square is the result of one major merger, which doubles the initial mass of the progenitor, while the third dot/square corresponds to a second major merger with mass 4 times the mass of the initial progenitor and 2 times the mass of the remnant of the first merging event. The minor merging tracks are shown in the right panels by black lines. Dots/squares indicate remnants with masses M ,0 + δM × M ,0 where δM = 0, 0.2, 0.4, 0.6, ..., and we use two different increment laws for M vir . In the top-right From these different tracks we can see that major mergers can be excluded, since the predicted evolution in size of a galaxy in the highest redshift bin (with z ∼ 0.6) and with stellar mass 10 11.4 M is parallel to the size-mass relation in this same redshift bin (see left panels in Figure 5).
On the contrary, the same galaxy can evolve to z ∼ 0.2 experiencing few (5 or 6) minor mergers, which accrete ∼ 100 per cent of the initial stellar mass M ,0 . The example galaxy is not evolving on the top of the z ∼ 0.2 R e -M and M dyn /M -M relations if we consider that DM and star accrete at the same rate, for both NFW and AC+NFW (dots and filled squares in Figure 5). The M dyn /M -M evolution is too steep, thus after ∼ 2 minor mergers the example galaxy would end up on the z ∼ 0.2 observed M dyn /M -M , this number is not consistent with what is found analyzing the R e -M evolution, which would require 4 − 6 mergers to transform the z ∼ 0.6 galaxy in a typical galaxy observed at z ∼ 0.2. The number of minor mergers needed to transform the z ∼ 0.6 galaxy in a bigger and more DM-dominated galaxy is found if the DM mass is accreting with a lower rate and only for the AC+NFW mass profile. This is possible if the main progenitor is merging with lower mass galaxies with smaller M dyn /M . This result is consistent with the indications provided in terms of DM profile by the steep DM density slope found and discussed in Section 4.
This result suggests that our massive galaxies at z ∼ 0.6 have to merge with a population of less massive galaxies with lower total-to-stellar mass ratios. Abundance matching studies predict that galaxies with M < 10 11.4 M have, on average, smaller M dyn /M (see e.g. Moster et al. 2010).

A COMMENT ON ALTERNATIVES TO COLD DARK MATTER
In this paper we have presented a dynamical analysis that allows us to derive constraints on the total mass in the central regions of ETGs from simple aperture velocity dispersion measurements (e.g. from fiber or multi-slit spectroscopy). As starting working hypothesis we have assumed a strict Cold Dark Matter (CDM) paradigm (e.g. galaxies are hosted in cuspy NFW haloes with or without baryonic effects like adiabatic contraction, see Sect. 4). In this scenario, our results point to a significant DM fraction within R e , which strongly depends on the assumptions made on the IMF. However, both the lack of direct evidence of the DM particle constituent and the presence of mismatches between predictions from the CDM and observations (e.g. the "missing satellite" problem, Simon and Geha 2007, or  Neither scenario, though, has found full observational or theoretical support. There are simulations adopting different DM flavours, or more complex cosmologies based on dark energy models or other alternative gravity theories. However, to our knowledge, these are not providing predictions for central DM fractions in massive ellipticals we can test against the data . From the modeling point of view, there have been a few attempts to incorporate in dynamical analyses of ETGs, some of these alternative theories. E.g., we have demonstrated that the outcomes of particular f (R) theories 3 , predicting an effective total potential including a Yukawa modification to the standard Newtonian potential, fit quite well the extended kinematic data of ellipticals . Moreover, in Tortora et al. (2014c) and Tortora et al. (2018a), we have shown that the central velocity dispersion of massive ETGs, can be interpreted within such alternative scenarios to DM, if a non-universal IMF is considered. Interestingly, within these alternative theories, we find taht a bottom-heavy (bottom-light) IMF is needed in the ETGs with the largest (smallest) σ, similarly to what is found within a standard ΛCDM paradigm.

CONCLUSIONS AND FUTURE PERSPECTIVES
In this paper we have summarised a simple but accurate dynamical method based on Jeans equations (see Section 2.1), to estimate the DM content within 1 effective radius in Early-Type Galaxies, making use of high quality (i.e. good seeing and high resolution) surface photometry and aperture internal kinematics (i.e. velocity dispersion). The method has been proved to reproduce results consistent with more accurate techniques (e.g. Schwarzschild's orbital superposition, Thomas et al. 2011;Tortora et al. 2018b), if limited to the central effective radii, measuring DM fractions across cosmic time and constraining physical processes driving this evolution.
We have focused, in particular, on the study of the M dyn /M and f DM for massive ETGs (log M /M > 11.2, Chabrier IMF) as a function of redshift, discussing the impact of galaxy models and IMF. We base our conclusions on the results obtained with two reference samples of ETGs: one made of z ∼ 0 local galaxies from the SPIDER project ) and one assembled cross-matching KiDS photometry with SDSS and BOSS spectroscopy (Tortora et al. 2018b). We list below the main results.
• The DM fraction strongly correlates with different galaxy parameters, e.g. the effective radius, the total mass within the effective radius, the Sérsic index, the mean stellar density, while it has a milder correlation with stellar mass and velocity dispersion, if the IMF is taken constant (see Section 3). In particular we have found that more massive and larger galaxies have a larger amount of DM. However, looking at the M dyn /M vs. M relation, a significant part of the scatter of the M dyn /M vs. M relation comes from the variation of the M dyn /M with redshift at fixed mass (see below), hence containing crucial information on the galaxy assembly processes. These DM fractions are quite consistent with independent literature, as combined dynamical+lensing analysis (Auger et al. 2010a;Shajib et al. 2020) or simulations (e.g. Lovell et al. 2018;Wang et al. 2020).
• Different lines of evidences suggest that IMF is not universal and change with galaxy parameters. M dyn /M analysis presented in this paper offers a simple way to reach this conclusion since higher stellar M/L normalizations, as those produced by a Salpter IMF, produce unphysical negative f DM for those systems which tend to have smaller sizes and masses (see also Tortora et al. 2013 andTortora et al. 2014b). If from one side the large variation of the IMF with velocity dispersion make DM fraction constant with this parameter, the milder correlation with sizes, produce a strong variation of DM fraction with size.
• The central mean DM density show indications of cuspy DM haloes in both local and higherredshift galaxies. Moreover, we have confirmed the evidences of a DM non-homology first discussed in Tortora et al. (2009), showing that observations prefer nearly isothermal total mass profiles at large masses and sizes and constant-M/L profiles at lower masses. This is consistent with a non-universality of the total mass density slope with mass. In fact, lower-mass galaxies, down to M ∼ 3 × 10 10 M , present steeper profiles, and this can be explained by a varying role of dissipation and merging in terms of mass (Dutton and Treu 2014;Tortora et al. 2014aTortora et al. , 2019.
• DM fractions of massive ETGs decrease with redshift: at a given mass galaxies at higher redshift tend to have smaller M dyn /M than galaxies at lower redshift. On the contrary, because of the smaller R e at larger redshift, the average DM density is larger at higher redshift. We have demostrated that this is an effect of minor merging of progenitor galaxies with lower-mass galaxies characterized by a smaller total-to-stellar mass ratio. Indeed, minor merging, unlike major merging, produce the correct growth of the galaxies moving from the highest redshift bin to the lower-one and still supply the necessary stellar mass for the stellar mass growth, but possibly a faster growth of the DM, which, combined with larger sizes, contributes to a faster growth of the overall M dyn /M .
The results discussed are very promising in view of the plethora of new observations that will increase data for massive ETGs up to z = 1 and above. The simple Jeans modelling described here has the great advantage to be applicable to large galaxy samples with elementary measurements, that will be soon under the reach of all sky surveys either from ground (photometry: e.g. Vera Rubin observatory/LSST, Ivezić et al. 2019;spectroscopy: StePS@WEAVE, Costantin et al. 2019, 4MOST, de Jong 2011, DESI, DESI Collaboration et al. 2016 or from Space (photometry: Euclid, Laureijs et al. 2011;CSST, Zhan 2018). The disadvantage is that the methods cannot constrain the stellar orbital anisotropy (which is assumed in the analysis to be reasonably constantly equal to zero in the very galaxy centers, Tortora et al. 2009) and is limited to the inference of the total mass; thus, the DM halo properties are model dependent (e.g. Tortora et al. 2014a). This makes the method equivalent to "gravitational lensing" only analyses (see e.g. Auger et al. 2010a). Hence, the two methods are often used in combination to break mutual degeneracies and/or derive constraints on the DM density slope (e.g., Shajib et al. 2020). However, if one is interested to collect information on the DM content of very large samples, eventually up to millions of galaxies in the era of next generation spectroscopic surveys, we have demonstrated in this paper that the Jeans method has a great potential in providing DM estimates that can be used to test models for dark and stellar mass assembly (see e.g. Section 5.2), with lower computation times with respect to more complex approaches. Furthermore, it can be easily adapted to test alternative theories of gravity or different DM flavor predictions.
However, on smaller samples, there are two different ways to obtain more precise constraints on the DM distribution. One consists in the combination of the Jeans analysis with strong gravitational lensing and stellar population. This will be possible with new samples of strong lenses that will be discovered in future wide-field observations (e.g. from Euclid and Rubin Observatory). In particular, Euclid is expected to find ∼ 170, 000 potential galaxy-scale gravitational lenses within the 15, 000 sq. deg. of the survey (; Petrillo et al. 2019). For a large part of these systems a measure of spectroscopic redshifts and velocity dispersions will be available, providing us and exceptional dataset for the determination of DM fractions. A second one is offered by repeated observations on individual galaxies observations from multiple instruments/surveys (see e.g. BOSS and LAMOST: Napolitano et al. 2020). In particular, the possibility of relying on different fiber/slit apertures will provide multiple constraints on the central mass profile of these galaxies, that can be eventually further improved if combined with lensing and/or high-spatial resolution from adaptive optics (e.g. MAVIS@VLT, McDermid et al. 2020).
Finally, the brightest possible future to push forward this kind of analysis will consist of deep spectroscopic surveys measuring velocity dispersions for lower-mass and higher-z ETGs, complementing what is only possible now at very high masses (M > ∼ 10 11 M ) and z ∼ < 0.7 (e.g. Tortora et al. 2018b). Collecting new observations for ETGs with M ∼ 3 × 10 10 M will allow us to map, as a function of cosmic time, the dynamical properties and the DM content of those galaxies which are characterized by the largest baryonic content and star formation efficiency. This "bimodality" mass scale, corresponding to a virial mass of ∼ 10 12 M , emerges from different kinds of observations and typically separate passive from star-forming systems, separating two mass regimes where galaxies are driven by different kinds of physical processes (e.g. Napolitano et al. 2005;Dekel and Birnboim 2006;Dekel et al. 2019;Tortora et al. 2019 and references therein). Going further below this characteristic mass scale, these observations, combined with current and future cosmological simulations, will allow us to describe the DM assembly of galaxies across a wide range of masses and redshifts within a single coherent galaxy formation framework.