MaNGA DynPop – VI. Matter density slopes from dynamical models of 6000 galaxies versus cosmological simulations: the interplay between baryonic and dark matter

We try to understand the trends in the mass density slopes as a function of galaxy properties. We use the results from the best Jeans Anisotropic Modelling (JAM) of the integral-field stellar kinematics for near 6000 galaxies from the MaNGA DynPop project, with stellar masses 10 9 M ⊙ ≲ 𝑀 ∗ ≲ 10 12 M ⊙ , including both early-type and late-type galaxies. We use the mass-weighted density slopes for the stellar 𝛾 ∗ , dark 𝛾 DM and total 𝛾 T mass from the MaNGA DynPop project. As previously reported, 𝛾 T approaches a constant value of 𝛾 T ≈ 2 . 2 for high 𝜎 e galaxies, and flattens for lg ( 𝜎 e / km s − 1 ) ≲ 2 . 3 galaxies, reaching 𝛾 T ≈ 1 . 5 for lg ( 𝜎 e / km s − 1 ) ≈ 1 . 8. We find that total and stellar slopes track each other tightly, with 𝛾 T ≈ 𝛾 ∗ − 0 . 174 over the full 𝜎 e range. This confirms the dominance of stellar matter within 𝑅 e . We also show that there is no perfect conspiracy between baryonic and dark matter, as 𝛾 ∗ and 𝛾 DM do not vary inversely within the 𝜎 e range. We find that the central galaxies from TNG50 and TNG100 simulations do not reproduce the observed galaxy mass distribution, which we attribute to the overestimated dark matter fraction, possibly due to a constant IMF and excessive adiabatic contraction effects in the simulations. Finally, we present the stacked dark matter density profiles and show that they are slightly steeper than the pure dark matter simulation prediction of 𝛾 DM ≈ 1, suggesting moderate adiabatic contraction in the central region of galaxies. Our work demonstrates the power of stellar dynamics modelling for probing the interaction between stellar and dark matter and testing galaxy formation theories.


INTRODUCTION
The prevailing Cold Dark Matter cosmological model proposes that structures in the Universe form hierarchically, where smaller dark matter halos form first and grow into larger ones through accretion and merger processes.If there were no baryonic matter, the evolution of the density profile of pure dark matter can be accurately traced using N-body numerical simulations.Various studies have shown that in a universe of cold dark matter, the density of dark matter at the inner regions of halos changes with radius according to a  −1 law (Navarro et al. 1996(Navarro et al. , 1997;;Springel et al. 2008;Gao et al. 2012;Wang et al. 2020b).
However, in the real Universe, baryonic matter condenses at the centre of dark matter haloes, forming galaxies, and the galaxies coevolve with the dark matter halo.Thus, the total material density profile at the center of a dark halo is the result of a combination of both ★ E-mail: ranl@bao.ac.cn baryonic and dark matter.Galaxies of different types possess varying stellar mass density profiles.For instance, early-type galaxies (ETGs) have stellar density slope with  * ∼ 2.3, assuming a mass-followslight case, while the stellar density slope (slope of the luminosity profile) of late-type galaxies (LTGs) is noticeably flatter, with some galaxies even approaching  * ∼ 1.2 (Li et al. 2019).In general, it is clear that if the halo has density slope of dark matter  DM and star  * , the total density slope  T can vary between  * and  DM when the dark matter fraction  DM varies between 0 and 1 in the same radial range.As a result,  T of a galaxy can vary significantly depending on the ratio of dark matter to baryonic matter at the center of the dark halo.
Apart from this ratio, the interplay between baryonic matter and dark matter during the galaxy formation process can also alter the overall material density profile in the center of the dark halo.Firstly, the collapse and accumulation of baryonic matter in the center of the dark halo may induce an 'adiabatic contraction', making the distribution of dark matter more concentrated (Blumenthal et al. 1986;Gnedin et al. 2004;Gustafsson et al. 2006).Secondly, the formation and subsequent evolution of stars can potentially impact the structure of the dark matter halo.For example, the feedback effects of stars and supermassive black holes can swiftly expel gas from the central region, thereby changing the gravitational potential well of the dark halo, resulting in a flatter central density profile of the dark halo (e.g.Read & Gilmore 2005;Pontzen & Governato 2012;Benítez-Llambay et al. 2018;Bose et al. 2019).
Deviations from 'cold' dark matter properties can also cause differences in the density profiles of dark matter halos' inner regions.In self-interacting dark matter model, scattering among the dark matter particles can generate a core-like profile in the centre of dark matter halo (Kaplinghat et al. 2016;Tulin & Yu 2018;Robertson et al. 2019;Bondarenko et al. 2021;Andrade et al. 2022;Eckert et al. 2022).
Observationally, density profiles can be measured using dynamics and gravitational lensing.For spiral galaxies, the rotation curves of the gas disk can provide an excellent probe of galaxy mass distribution (e.g.Bosma 1978;Faber & Gallagher 1979;Rubin et al. 1980).Recently, Tortora et al. (2019) analyzed the rotation curves of 175 spiral galaxies from the SPARC project (Lelli et al. 2016), calculating the mass-weighted density slopes for these galaxies.They found that the density profiles of these galaxies become steeper as galaxy mass increases.At the low-mass end (less than 10 8 M ⊙ ), the total matter density profile slope of the galaxies is nearly 1, while at the high-mass end (greater than 10 10.5 M ⊙ ), the density profile slope of the galaxies is close to 2, indicating an isothermal density profile.
Galaxy-galaxy strong gravitational lensing is typically applied to ETGs of high velocity dispersion.The Sloan Lens Advanced Camera for Surveys (SLACS; Bolton et al. 2006) searched for lensing candidates within the SDSS survey, and through subsequent observations with the Hubble Space Telescope, high-resolution images were obtained for lens modeling.For galaxies with stellar masses greater than 10 11 M ⊙ , they discovered that the average inner density profile slope of total mass ⟨ T ⟩ = 2.078 (Koopmans et al. 2009;Auger et al. 2010), which is very close to predictions made by the isothermal model.This result has been confirmed by subsequent similar strong gravitational lensing analyses (Sonnenfeld et al. 2013;Li et al. 2018b;Etherington et al. 2023) and the near-isothermal feature of ETGs is sometimes termed as the 'bulge-halo conspiracy' (see Section 4.2 for details).
Compared to strong gravitational lensing and gas dynamics, stellar kinematics can be more widely applied to different types of galaxies over a broad range of galaxy masses.Cappellari et al. (2015) was the first to systematically apply the Jeans Anisotropic Modelling (JAM) method to analyze the dynamic data of 14 large-mass ETGs, reporting a nearly 'universal' average total density profile slope of ⟨ T ⟩ = 2.19, with small scatter.Using the JAM method, Poci et al. (2017) analyzed 260 ETGs from the ATLAS 3D survey (Cappellari et al. 2013a).They found for galaxies with velocity dispersion higher than a threshold lg( e [km s −1 ]) ∼ 2.1, ⟨ T ⟩ = 2.193, but for low velocity dispersion galaxies,  T shows a decreasing trend with decreasing velocity dispersion.
Thanks to the SDSS-IV MaNGA survey (Bundy et al. 2015), a new generation of Integral Field Unit (IFU) survey, a broad range of galaxy types over a wide mass range have obtained stellar kinematic data usable for dynamic modeling.Li et al. (2019) analyzed MaNGA galaxies in SDSS Data Release 14 (SDSS DR14; Abolfathi et al. 2018), which covers a mass range from 10 9 M ⊙ to 10 12 M ⊙ , including both early-and late-type galaxies.Using the mass distribution models obtained with the JAM method from Li et al. (2018a), they computed the total density slope of these galaxies.For ETGs with  e > 100 km s −1 , they obtained similar results to previous studies, with the average mass-weighted total density slope ⟨ T ⟩ = 2.24.Moreover, they showed not only for ETGs, but for all types of galaxies as a whole, how the total density slope decreases with decreasing velocity dispersion below a velocity dispersion threshold.Additionally, they showed that central galaxies in clusters have a flatter  T than satellite galaxies.
In this study, we utilize the mass distribution models derived from the MaNGA DynPop project (Zhu et al. 2023(Zhu et al. , 2024;;Lu et al. 2023b,a;Wang et al. 2024), which analyze the full data release of MaNGA survey that contains the IFU data for over 10,000 nearby galaxies.This dataset cover a wide range of galaxy types and stellar masses, from 10 9 to 10 12 M ⊙ (Wake et al. 2017).The MaNGA DynPop project takes advantage of this IFU data for galaxy dynamical modeling (Zhu et al. 2023, Paper I, hereafter).In Lu et al. (2023b, Paper II, hereafter), we obtain the stellar population properties of MaNGA galaxies under the Salpeter (Salpeter 1955) initial mass function (IMF) assumption, including the stellar mass-to-light ratio used here to calculate the stellar mass  * .In Zhu et al. (2024, Paper III, hereafter), we calculate the total density slope for 6000 MaNGA galaxies.Our results once again confirm that that for ETGs, the average mass-weighted total density slope is slightly steeper (approximately 2.2) than the isothermal case.This finding is consistent with the previously observed trend in the variation of galaxy total density slopes with velocity dispersion, as reported by Li et al. (2019).In Wang et al. (2024, Paper IV, hereafter), we select central galaxies of clusters and groups that have weak gravitational lensing measurements, and jointly constrain their density profile slopes using both dynamics and weak gravitational lensing.We show that the inner total density profiles of central galaxies in clusters and groups are close to isothermal profile, and find that such density profiles require the stellar mass of the galaxies to be significantly higher than the stellar mass derived from stellar population synthesis models assuming a Chabrier (Chabrier 2003) IMF.
In this paper, we will use the galaxy mass distribution models obtained from JAM to more comprehensively analyze the evolution of galaxy density slope with galaxy properties, the connection between the slopes of total and stellar density, along with the dark matter fraction, and the comparison of observational results with hydrodynamic numerical simulations.Additionally, we will also analyze the stacked density profile of galaxies.
The structure of this paper is as follows.In Section 2 , we briefly introduce the MaNGA project and the simulation data we use for comparison.Section 3 provides an brief overview of the dynamical modeling methods used in MaNGA DynPop, as well as the method we use to compute the mass-weighted density slope in this paper.Our results are presented in Section 4. In Section 5, based on the results of this study, we discuss the bulge-halo conspiracy scenario, the implications of our comparison with simulations, and the impact of the quality flag in our results.The final section, Section 6, summarizes our conclusions.

MaNGA galaxies
This paper investigates density slopes of galaxies in the final data release of the MaNGA survey (SDSS DR17;Abdurro'uf et al. 2022), which includes unprecedentedly more than 10000 nearby galaxies.For further information of the MaNGA survey, the readers are referred to papers listed below: an overview of MaNGA (Bundy et al. 2015), SDSS-IV technical summary (Blanton et al. 2017), MaNGA instrumentation (Drory et al. 2015), sample design (Wake et al. 2017), observing strategy (Law et al. 2015), spectrophotometric calibration (Smee et al. 2013;Yan et al. 2016a), survey execution and initial data quality (Yan et al. 2016b) and an overview of SDSS telescope is included in Gunn et al. (2006).
We adopt the mass-weighted density slopes (see Section 3.2 for details) derived in Paper I, which construct accurate mass models using the Jeans Anisotropic Modelling (JAM; Cappellari 2008Cappellari , 2020) ) method and the MaNGA stellar kinematics extracted from the Data Analysis Pipeline (DAP; Belfiore et al. 2019;Westfall et al. 2019).A brief introduction to the dynamical modelling is presented in Section 3.1.We select galaxies with an acceptable visual modelling quality (i.e.Qual ⩾ 1) to ensure that the density slopes are reliable.For these approximately 6000 galaxies, we require the difference in the mass-weighted total density slope between JAM cyl and JAM sph modelling (see Section 3.1) is smaller than three times of 0.079, which is the observed root-mean-square (rms) scatter of this dynamical property among different model assumptions for galaxies with Qual = 1 (see Fig. 13 and Table 3 of Paper I), to ensure the reliability of our conclusion.In result, our final MaNGA sample includes 5688 galaxies.For more information of modelling qualities and the dynamical modeling, readers are referred to Paper I.

IllustrisTNG galaxies
We compare the measured mass-weighted density slopes with theoretical prediction derived from The Next Generation Illustris Simulations (IllustrisTNG, TNG hereafter), which is a suite of state-ofthe-art magneto-hydrodynamic cosmological galaxy formation simulations carried out in large cosmological volumes with the movingmesh code arepo (Springel 2010).The TNG collaboration has publicly released the general properties of galaxies, which are available for use (Nelson et al. 2019a).In this paper, we use the highest resolution version of two simulations from the suite: the TNG 50 (Nelson et al. 2019b;Pillepich et al. 2019), which is a full-physics version with a box size of 51.7 Mpc and has a mass resolutions for both baryonic and dark matter component of  baryon = 8.5 × 10 5 M ⊙ and  DM = 4.5 × 10 5 M ⊙ , respectively; the TNG100 (Pillepich et al. 2018;Marinacci et al. 2018;Naiman et al. 2018;Nelson et al. 2018;Springel et al. 2018), which is a full-physics version with a cubic box of 110.7 Mpc side length and has a mass resolution for baryonic and dark matter component of  baryon = 1.4 × 10 6 M ⊙ and  DM = 7.5 × 10 6 M ⊙ , respectively.The gravitational softening length for dark matter and stellar particles is  softening = 0.288 kpc for TNG50 and  softening = 0.738 kpc for TNG100.Galaxies residing within their host dark matter halos are identified using the subfind algorithm (Springel et al. 2001;Dolag et al. 2009).To match the MaNGA observation, we select central TNG galaxies from snapshot 99 (corresponding to redshift  = 0) with stellar mass  * ⩾ 10 9 M ⊙ .Combined with the selection criterion based on galaxy size (see Section 3.2), our final TNG sample ended up with 1733 TNG50 galaxies and 6107 TNG100 galaxies.

Dynamical modelling
In Paper I, we use the axisymmetric Jeans Anisotropic Modelling (JAM; Cappellari 2008Cappellari , 2020) ) method to derive the dynamical quantities for the whole MaNGA sample and provide the values of each quantity for eight models.Among the eight models, we adopt two assumptions of velocity ellipsoids (cylindrically-aligned model JAM cyl and spherically-aligned model JAM sph ) and four mass models (for each velocity ellipsoid assumption) that have different assumptions on the dark matter distribution.The surface brightness of each galaxy is obtained from Multi-Gaussian Expansion (MGE; Emsellem et al. 1994;Cappellari 2002) fitting to the SDSS −band images, and then is deprojected to obtain the luminosity distribution of the kinematic tracers.The total mass distribution (or equivalently the total gravitational potential) is described as the luminosity density multiplied by a spatially constant stellar mass-to-light ratio and the dark matter mass distribution.For a given gravitational potential, the modelled velocity second moments derived from the axisymmetric Jeans equations are compared to the observed one to determine the best-fitting parameters.The JAM method has been tested on mock galaxies generated using cosmological hydrodynamical simulation, which has shown that the density slopes of galaxies can be recovered robustly (Li et al. 2016).
We mainly use the JAM cyl + generalized Navarro-Frenk-White (gNFW) dark model in this work, which is the most flexible mass model provided in Paper I and is supposed to provide the most accurate measurements of density slopes.The gNFW profile (Wyithe et al. 2001) is described as where   is the characteristic radius,   is the characteristic density, and  is the inner density slope.For  = 1, this function reduces to the NFW profile.The density slopes of other mass models can be used to access the systematic uncertainties, which are listed in Table 3 of Paper I. Specifically, for our sample with Qual ⩾ 1, the rms scatter of the mass-weighted total density slopes between different models is smaller than 0.079, indicating the robustness of density slopes measurements.We adopt three times of this value in Section 2.1 to filter galaxies as a guarantee of reliability in JAM modeling.In the following sections, unless explicitly stated otherwise, all the dynamical properties of MaNGA galaxies utilized are from the catalogue of Paper I.

Mass-weighted density slope
With galaxy density profile approximated by power-law  ∝  − , one can follow Dutton & Treu (2014) to calculate the mass-weighted density slope within the effective radius  e as For MaNGA galaxies, the projected circularized half-light radius  e ('Re_arcsec_MGE' in the catalogue of Paper I) is derived from MGE fitting of MaNGA galaxy's r-band image.For TNG galaxies, we did not directly use the 3D half-mass radius provided in the TNG subhalo catalog but assumed that the projection direction was along the x-axis and based on the surface mass density of the stellar particles in the y-z plane, we calculated the projected radius of a circle enclosing half of the galaxy stellar particles as  e .Throughout this paper, we use  x to denote the mass-weighted density slope, where x represents the material components we are considering.We use 'T', '*' and 'DM' to label total, stellar and dark matter, respectively.The mass-weighted density slopes of MaNGA galaxies come from the catalogue of Paper I. For TNG galaxies, we divide galaxy's radius from  softening to 100 kpc into 50 equal intervals in logarithmic space to obtain the density and enclosed mass profile and interpolate them as function of radius.When considering the total mass, we combine the masses of star, dark matter, and gas particles while when focusing on the mass of dark matter, we sum up the masses of dark matter and gas particles.In cosmological N-body simulations, the gravitational softening length  softening is introduced into the Newtonian gravity formula to prevent close encounters between particles.In order to mitigate the impact of unrealistic gravitational softening on the density profiles of simulated galaxies (see details in Zhang et al. 2019), it is necessary to exclude the innermost region of the simulated galaxies.Similar to Wang et al. (2020a), we set 0.3  e as the lower integral limit (in Wang et al. 2020a, they choose 0.4  e to include more ETG samples) and only choose simulated galaxies with 0.3  e >  softening .We should aslo change Eq. 2 into and substitute 0.3  e into  min .For TNG50, 92 per cent (1733) central galaxies satisfy this criterion.But for TNG100, due to the larger  softening , only 50 per cent (6107) central galaxies satisfy this criterion.However, we have verified (although we do not show) that the results presented in the remainder of this paper remain unchanged, regardless of the inclusion or exclusion of these TNG100 galaxies with smaller size.Furthermore, if we change the filtering criteria to 0.3  e > 2.5 softening , it will not affect our conclusion.

Total density slope
In Fig. 1, we present the scaling relations between the mass-weighted total density slopes  T and the velocity dispersion within effective radius  e or the stellar mass  * . e ('Sigma_Re' in the catalogue of Paper I) is defined as the square root of luminosity-weighted second velocity moments within the elliptical half-light isophote ('Rmaj_arcsec_MGE' in the catalogue of Paper I) and is calculated as where   ,   and   are the flux, stellar velocity, and stellar velocity dispersion in the -th IFU spaxel.For each galaxy,  * of Salpeter IMF is calculated as where ( * /) SPS (<  e ) (comes from 'ML_int_Re' in the catalogue of Paper II) is the r-band stellar mass-to-light ratio obtained from the stacked spectrum again within the elliptical half-light isophote, while  (comes from 'Lum_tot_MGE' in the catalogue of Paper I) represents the r-band total luminosity derived from MGE models.
To compare with the simulations which assume Chabrier IMF, we then subtract the  * by 0.25 dex (Bernardi et al. 2010).As demonstrated in Figure 8 of Paper III,  T increases rapidly with velocity dispersion from 1.6 to 2.2 for galaxies with lg( e [km s −1 ]) ≲ 2.28 (≈ 190 km s −1 ).Above this velocity dispersion, the median  T remain relatively constant at 2.2, slightly steeper than the slope of isothermal profile (that is, 2).This trend can be well described by equation ( 13) of Paper III, which reads: with { 0 ,   , , , } = {2.18,189, 11.13, −0.02, 0.34} the bestfitting parameters. T increases with  * as well.A mild transition can be found at  * ≈ 10 11  ⊙ , above which the median value of  T is constantly 2.2.
Here we make a comparison with the results of Li et al. (2019).While the overall trend of the  T - e relation is similar to the findings in Li et al. (2019), the transition of the slope decrease occurs at a lower value of  e in Li et al. (2019), at around 100 km s −1 , and the trend difference is less apparent in the plot because logarithmic coordinates are used.We infer this discrepancy is mainly due to the difference of measurements in  e for different data releases.With the updated MaNGA DAP (Law et al. 2021), it can be found that the velocity dispersion in SDSS DR17 had been underestimated in previous MaNGA data releases, especially for low-mass galaxies.Therefore, the same galaxy at the final data release will have a higher velocity dispersion and in turn have a steeper total density profile predicted by dynamical modelling.As a result, the transition of the T - e relation in ETGs can be found in Table 1 of Paper III, and it shares the same form as the one for the full sample.The parameter values for this equation are {  0 ,   , , , } = {2.24,150, 397.85, −0.03, 0.11}.The turnover occurs earlier at  e ≈ 150 km s −1 .The gradient of the decreasing slope with decreasing velocity dispersion becomes gentler, reflecting the fact that ETGs generally exhibit steeper profiles.
We also plot the mass-weighted stellar density slope  * derived under the gNFW model in Fig. 1.Interestingly, the trend of the mass-weighted density slope of the stellar component is very similar to that of the total mass.We fit the  * −  e relation using a function of the same form as for  T with {  0 ,   , , , } = {2.33,190, 9.98, 0.03, 0.25} the best-fitting parameters.The mean difference between the fitted values of  * and  T is consistently 0.174 over the full  e range with the standard deviation of the differences is roughly 0.016, indicating that the total and stellar slopes closely follow a similar trend.As shown in Section 3.4 of Paper III and Fig. 4, the JAM modelling generally predicts galaxies with low dark matter fraction within  e ,  DM (<  e ).Especially for ETGs, 90 per cent of them have  DM (<  e ) < 23 per cent.Therefore, the similarity of  * and T may be a direct result of the stellar component dominating the total mass budget within the effective radius.
From Fig. 1 we can see that gNFW and NFW model assumption derive almost identical median scaling relations on  T (except for small  e and  * end), reflecting the reliability of JAM modelling.In the following figures, unless otherwise specified, we all display the results of MaNGA galaxies under JAM cyl + gNFW model.
In Fig. 2, we present  T as a function of halo mass  200 (left) and  e (right) for central galaxies.We adopt the central-satellite classification results and  200 from an updated version 1 of the SDSS DR4 group catalog (Yang et al. 2007) (SDSSGC, hereafter).Following SDSSGC, there are 3831 central galaxies (the galaxy with the largest stellar mass in one group) in our MaNGA samples, of which 2922  2019) is plotted by pink dashed line.The median and the 1- scatter are showed by red diamonds with error bars.The black contours show the kernel density estimate for the sample distribution.The slope value of the isothermal profile is depicted with a black horizontal dash-dot line.
have halo mass that are considered reliable and are plotted in left panel of Fig. 2. The median  T exhibits a subtle increase from 1.9 to about 2.2 among central galaxies with  200 ≲ 10 13 M ⊙ , and remains consistently fixed at 2.17 when  200 ≳ 10 13.5 M ⊙ .The behavior of  T - e for central galaxies is the same as that for full sample, except this relationship has a slightly more gentle turn for central galaxies.Li et al. (2019) found that, for different values of  e , the mean  T of satellite galaxies is overall about 0.1 higher than that of central galaxies, suggesting that this may be due to differences in the galactic formation background.Here, we further confirm this phenomenon with SDSS DR17.
We also compare the results of TNG50 and TNG100 in Fig. 2. The  e of TNG galaxies is calculated as the velocity dispersion of stellar particles along the X-axis within the 2D half mass radius described in 3.2.The line widths of the median value at high  e bins reflect the lack of massive galaxies in TNG due to the limited volume of the simulation box.It can be seen that the  T of TNG50 galaxies with  200 > 10 13 M ⊙ and  e > 150 km/s are generally in good agreement with MaNGA, while results from TNG100 show obvious differences with them: for TNG100 galaxies with  200 > 10 13 M ⊙ ( e > 150 km/s), their median  T decreases gradually from 2.03 (2.22) to 1.76 (1.76).However, this decreasing trend is consistent with Remus et al. (2017), where ETGs selected from the Magneticum Pathfinder simulations (Dolag et al. 2016) were studied.At low mass end (at  e ∼ 100 km/), both TNG50 and TNG100 overestimate the total density slope, by ∼0.4 for TNG50 and ∼0.37 for TNG100.To understand why do the TNG100 and TNG50 differ at high mass end, we further plot the  e −  200 ,  e − e and  e −  * relations in Fig. 3.The range of the horizontal axis for the comparison is constrained to be approximately the same as that in Fig. 2. The lg  e of TNG100 is obviously higher by at least 0.2 dex than TNG50, which results the observed lower value of  T for TNG100 in high-mass end in Fig. 2. In the related study Paper IV, we measures the total density slope for a subset of central galaxies in groups and clusters by combining their stellar kinematics with weak gravitational lensing.For galaxies in the group bin with the best-fitting result at lg  200 [M ⊙ ] of 13.2,  T of their mean density profile is 2.15 +0.04 −0.05 , which is in excellent agreement with the findings presented in this paper.However, for galaxies in the cluster bin with the best-fitting result at lg  200 [M ⊙ ] of 13.92,  T of their mean density profile is 1.95 +0.08 −0.09 , which is lower than the result in this paper by about 0.2 for the same mass range.The source of this difference primarily stems from two aspects.First, the data points in Fig. 2 represent the stacked  T , which is obtained by first averaging the density profiles of galaxies within a specific mass bin and then calculating the density slope as Eq. 2. Through our tests we found that for pure JAM modelling results of MaNGA, the stacked  T in the lg  200 [M ⊙ ] range of 13.5 − 14.8 is lower than the median value of galaxy density slopes by 0.1.In Sec 4.3, we will further observe differences between calculating  T from averaged density profiles and obtaining median (or mean)  T values from the galaxies in the same bin.Second, the residual portion of the discrepancy can be attributed to differences in sample selection.Paper IV include galaxies with Qual⩾0, whereas the results presented in this paper is based on galaxies with Qual⩾ 1.We have exclusively selected galaxies with Qual ⩾ 1 to ensure the reliability of our density slope measurements.Paper IV employed mass measurements within a fixed radius (approximately 10-20 kpc) derived from Paper I, in conjunction with gravitational lensing and this calculation of total mass for MaNGA galaxies is relatively reliable even for galaxies with Qual = 0.The discrepancy in the measurement of  T between this paper and Paper IV yields an intriguing implication: the total density profiles of galaxies are contingent upon their dynamical state.Galaxies with Qual ⩾ 0 exhibit relatively complex dynamical states and flatter density profiles.Furthermore, in Sec.5.3 we present the sample distribution characteristics of Qual on  T values to illustrate the systematic differences in the matter distribution among galaxies of different JAM modelling quality.If we calculate the stacked  T for the same group-and cluster-central galaxies, including those Qual = 0 sample, but solely based on the JAM modelling results, we also predict the consistent  T , as Figure 9 of Paper IV shows.We also compare our observation with recent measurements of total density slope.Following the images captured by the Hubble Space Telescope, Etherington et al. (2023) employed the pixel-based strong lensing modelling technique, as detailed in Nightingale et al. (2021), to derive lensing-only measurements of the logarithmic total density slope  for 42 early-type galaxies (ETGs) within the SLACS sample, which are originally presented in Bolton et al. (2008).By jointly combining the weak and strong lensing constraints from 22 SLACS galaxies, Gavazzi et al. (2007) measured the virial-to-stellar mass ratio  119 / * = 54 +28 −21 . * for 36 galaxies, assuming a Chabrier IMF, is available in Auger et al. (2010), and we use these values to estimate the halo mass,  200 , for the SLACS sample.We employ a mass-concentration relation from Dutton & Macciò (2014) to convert  119 to  200 .The  values obtained through pure lensing measurements exhibit a slight shallower trend when compared to the median result from JAM, but these measurements are generally consistent with the JAM results.
For the high-mass end, we compare our result with ten massive galaxies from Newman et al. (2015) at the center of massive groups with an average mass of  200 ∼ 10 14 M ⊙ , and seven very massive galaxies from Newman et al. (2013a,b) at the center of clusters have halo mass  200 = 0.4 − 2 × 10 15 M ⊙ , with total mass density profiles obtained via a combination of weak and strong lensing, resolved stellar kinematics and X-ray kinematics.The mass-weighted total density slopes  T of these massive galaxies reach ∼ 1.7 for galaxy groups and ∼ 1.2 for clusters.Comparing these results to our study, we find that the  T values for the ten group-scale central galaxies are notably lower than the mean value of MaNGA.However, they are in good agreement with the results from TNG100.Nevertheless, through semi-empirical approach, Shankar et al. (2017) found that some of their model predictions can match the variation of  T from SLACS galaxies to group-scale central galaxies in Newman et al. (2015), but the density slopes of cluster-scale galaxies in Newman et al. (2013a) are still too low to be explained.They also demonstrated that the dependence of  T on halo mass is a genuine effect, only partially influenced by the increase in effective radius with stellar mass, reflecting the structural non-homology, where the relative density distribution of stars and dark matter vary systematically from isolated galaxies to central galaxies in clusters.The dependence of the total density slope on the environment cannot be ignored.Further investigations are needed to understand the origin of the discrepancy and to improve the accuracy of the total density slope measurements for galaxy groups and clusters.

Dark matter and stellar components
To inspect the reason of a relatively shallow  T in galaxies of Il-lustrisTNG simulation, we investigate the decomposed dark matter and stellar components respectively.The central dark matter fraction  DM (<  e ) is defined as the mass ratio of dark matter over the total mass enclosed within the 2D effective radius  e for both MaNGA and TNG galaxies.In Paper III (Figure 11) we already have looked into the  DM (<  e )- * relation.Here we show again the relation in Fig. 4 to compare it with the results from TNG.For MaNGA galaxies, according to the suggestion in Table 2 of Paper I, we add the filtering condition |  DM,cyl −  DM,sph | < 0.1, to improve the accuracy of measuring the dark matter fraction.Although the trend of the  DM (<  e )- * relation is roughly same, the difference between MaNGA and TNG is obvious: simulated galaxies tend to have much more dark matter within their central region and this difference is greater at the large  * end, especially for TNG100.
In Fig. 5 we compare the  * - * relation between MaNGA and TNG galaxies.It is evident that the numerical simulations predict much higher stellar density slopes compared to the observations.Therefore, the shallower overall density profiles of the simulated galaxies, as depicted in Fig. 2, can be attributed entirely to the disparity in the dark matter fraction between the simulations and observations, rather than differences in the stellar density distribution.Fig. 4 and Fig. 5 also demonstrate that although TNG50 appears to better predict the overall density profiles at the high-mass end, this is merely an incidental effect of the slightly flatter stellar density profiles and relatively lower dark matter fraction combined.It is not a robust indication of TNG50 outperforming TNG100 in modeling the total density profiles.
Previous lensing and dynamical observations find that massive elliptical galaxies have a mass density profile close to isothermal ( ∝  −2 ,  T = 2) within an effective radius (Koopmans et al. 2009;Auger et al. 2010;Tortora et al. 2014;Poci et al. 2017;Li et al. 2018bLi et al. , 2019) ) or well beyond the effective radius (Gavazzi et al. 2007).This phenomenon is commonly referred to as the 'bulge-halo conspiracy' (e.g., Dutton & Treu 2014).Given that neither baryons nor dark matter exhibit an isothermal density profile on their own, it suggests that there is an interaction or 'conspiracy' between baryons and dark matter that collectively results in the formation of an isothermal density profile.
We intend to use our results from MaNGA to elucidate this issue.Specifically, we examine the characteristics of the distribution of stars and dark matter in galaxies with respect to  T and  e in Fig. 6 and Fig. 7, with particular focus on galaxies that exhibit a higher velocity dispersion while maintaining a total density slope of ∼2, to investigate how their stellar and dark matter distributions vary.This will allow us to explore whether there is a matching mechanism between luminous and dark component that contributes to the 'bulge-halo conspiracy'.We utilize  * and  DM (<  e ) to characterize the distribution of stellar and dark matter, respectively.For galaxies within a specified  e value, a smaller effective radius corresponds a steeper density profile within  e , leading to a higher .To mitigate the impact of the diversity in  e among galaxies and ensure a fair comparison, we employ the median  e - e relation to derive an interpolation function,  e ( e ).Then we recompute  and label them as ' within  e ( e )'.For brevity, we will continue to refer to them as  T and  * in the following text.
To capture the overall trend of the data distribution, the colored scattering points have been smoothed using the loess technique, a Locally Weighted Regression method developed by Cleveland (1979), and implemented in the loess Python procedure 2 created by Cappellari et al. (2013a).For full sample, the dispersion of  T is not random but instead depends on the stellar density profile and the proportion of dark matter.Given the stellar mass, a higher density slope of the stellar component corresponds to a higher total density slope of the galaxy.Similarly, a higher proportion of dark matter leads to a flatter total density profile.There is no sign that stellar and dark matter components balance each other to yield a  T ∼ 2, as predicted by the bulge-halo conspiracy theory.But ETGs exhibit a markedly different scenario: across the entire  e range, the median value of  T ∼ 2 remains constant at 2.  slope at the same  e and consider it a fairer comparison for  T .This appears to corroborate the conspiracy theory regarding massive elliptical galaxies, if we disregard that  T is around 2.2 (rather than 2) and that their  * and  DM (<  e ) still exhibit a systematic layered structure.

Stacked density profile
Next, we stack (take the average of) the density profiles of the central galaxies in different mass bins, and calculate the slopes of the stacked density profiles of stars, dark matter, and total matter to study how the stacked density slope varies with mass.Firstly, we divide the MaNGA galaxies into sub-sample bins, based on  * ,  e and  200 , respectively.Then we obtain the mean density profiles of stars, dark matter, and total matter for each mass bin, which are plotted in Fig. 8, Fig. 9 and Fig. 10.In these plots, each panel represents one mass bin, and the dashed gray vertical line indicates the mean  e of galaxies in that bin.Note that, Yang07 only measured the reasonable dark matter halo mass of galaxies with  * bigger than ∼ 10 10 M ⊙ and the density profiles of galaxies in this mass bin corresponding to the upper left panel of Fig. 8 are almost absent in Fig. 9. Once again, we observe that the stellar matter dominates within  e , especially Median MaNGA TNG50 TNG100 Figure 5.  * as a function of  * for central galaxies.For MaNGA galaxies, the black solid line shows their median and the green shadow region indicates the 1- scatter.For TNG50 and TNG100, the medians and 1- scatters are showed by scattered dots with error bars.
in galaxies with larger masses, where the dark matter density in the central region is only 1/10 of the stellar density.We can see that within the mean  e , the slope of the stacked dark matter density profile in each mass bin is slightly steeper than the prediction of NFW profile, namely 1, while the slope of the stacked total matter density profile, except for the lowest  * bin, is very close to 2 in all other bins.
We investigate how the stacked density slope changes with galaxy properties.In Fig. 11, we show how the stacked density slope of different components in MaNGA and TNG galaxies varies as a function of the mean mass (velocity dispersion) of galaxies within each stacked mass (velodicty dispersion) bin.We used bootstrap resampling to obtain the uncertainty of stacked density slopes within each bin.Specifically, we perform a re-sampling with replacement of the galaxies in each bin.For each re-sampling, we select galaxies from the bin with replacement and stack their density and enclosed mass profiles to calculate the stacked density slope.The number of galaxies sampled in each re-sampling is equal to the total number of galaxies contained in that bin.We repeat this process 100 times to obtain the standard deviation of the stacked density slopes.
For a fair comparison, we also resample galaxies from TNG50 and TNG100 to match the  * ,  200 and e space distribution of galaxies of MaNGA galaxies in each subsample bin.After obtaining the re-sampled TNG galaxies, we computed the stacked density slopes.However, it should be noted that due to the limited volume of the simulation box, the massive galaxy subsample in TNG50 and TNG100 may be less representative.For the mass range of  200 > 10 13 M ⊙ , where TNG50 (TNG100) only has 24 (127) galaxies, while MaNGA has 1081 galaxies.
In Fig. 11, the stacked density slopes of different components are denoted as stacked  DM , stacked  * , and stacked  T , respectively.For MaNGA galaxies, overall, stacked  T increases with increasing  e and  * and become flatter in the highest two  e and  * bins.The behavior of stacked  T as a function of  200 is relatively flat, which is due to the lower limit cutoff of our dark halo mass at 10 12 M ⊙ , not including central galaxies with smaller  * or  e .
Different with the situation of median value in Fig. 1 that  * is uniformly higher by ∼ 0.2 than  T across nearly all MaNGA samples, stacked  * is higher than stacked  T by 0.36, 0.28, 0.21 and 0.24 in the four  * bins (by 0.38, 0.23, 0.20 and 0.24 in the four  e bins), respectively.We find that this is because at the low- * ( e ) end, the value of stacked  * (slope of the mean stellar profile) is higher by ∼ 0.2 than the median of  * while at the high- * (high- e ) end these two values match quite well, and the difference of these two values decreases with the increasing  * or  e .This results in a less pronounced growth trend of stacked  * with  * or  e compared to stacked  T .We attempt to explain the reasons for the different behaviors of stacked  * in different  * or  e bins in Appendix A. The dark matter density slope is slightly steeper than  ∼ 1, ranging between 1.1-1.4,presenting a flat inverted-U shape with changes in  * and  e , peaking at ∼ 1.8×10 10  ⊙ ( e ≈ 150 km/s), and tending towards  ∼ 1 at both the higher and lower mass ends.
There are significant differences in the values of stacked  DM and stacked  * between MaNGA and TNG galaxies, with TNG galaxies having higher values than MaNGA galaxies.Although both the decomposed stellar and dark matter components in the simulations show obvious higher density slopes, the relationship between total density slopes and stellar mass in TNG50 and TNG100 appears to produce a similar trend and slightly higher magnitude to observations.This is because the dark matter fraction in galaxies in the TNG simulations is evidently higher than those observed in MaNGA.The relatively flat dark matter density distribution substantially counterbalances the steep stellar density distribution, coincidentally resulting in a similar outcome to the observations.
For the relationship between  e and total density slope, and for the relationship between  200 and total density slope, the results given by TNG50 and TNG100 are quite different, and they are significantly different from the observations.TNG100 predicts significantly lower total density slopes at the high-velocity dispersion end compared to  TNG50, with a difference reaching 0.2-0.4.In fact, we find that for the central galaxies in large mass halos in TNG100, they have flatter stacked stellar and dark matter density profiles, and also have a larger dark matter fraction (Fig. 4).This is why it gives a flatter density profile.
It is worth noting that the TNG simulations and MaNGA observations also have different  e - * relationships.At the same stellar mass, MaNGA galaxies have a larger velocity dispersion.This is also a reason why the results are complex in the comparisons between TNG simulations and MaNGA in different density slope scaling relationships.

Bulge-halo conspiracy
In Fig. 1 we illustrate that the median  T of ETGs remains at 2.2 only with  e ≳ 150 km s −1 .This is the same value originally reported as 'universal' for ETGs by Cappellari et al. (2015).Below this velocity dispersion,  T becomes more dispersed and decreases as  e decreases.This finding corroborates the conclusions from the dynamical studies of ETGs in previous ATLAS 3D and early MaNGA data releases (Poci et al. 2017;Li et al. 2019).In Fig. 6 and Fig. 7, we observe that when we keep  e fixed to calculate the  T of ETGs for a specific  e value, the median  T remains consistently at 2.2, regardless of changes in  e .This result seems to support the bulgehalo conspiracy theory.However, this situation does not persist if we employ different methods of ETGs classification.In the preceding sections, we have utilized the stringent selection criteria proposed by Domínguez Sánchez et al. (2022) to identify our ETG sample, which predominantly consists of elliptical (E) and lenticular (S0) galaxies.The primary criteria for this morphological classification of ETGs are P LTG < 0.5 and T-Type<0, effectively distinguishing ETGs from LTGs.
Following Genel et al. (2018) and Lu et al. (2021), the luminositymorphology parameter  dev / tot and the specific star formation rate (sSFR) is individually adopted by us. dev / tot is defined as the ratio of the de Vaucouleurs luminosity to the total galaxy luminosity and characterizes the galaxy's fraction of the elliptical component, while sSFR is defined as the ratio of the star formation rate (in M ⊙ /Gyr) to the galaxy stellar mass.ETGs typically have the elliptical morphology ( dev / tot > 0.5) and quenched star formation (lg sSFR/Gyr −1 ⩽ −2), while late-type galaxies are dominated by disks ( dev / tot < 0.5) and still have ongoing star formation (lg sSFR/Gyr −1 ⩾ −1.5). dev / tot comes from the MaNGA Deep Learning morphological catalogue, while sSFR is computed using the reliable ('QCFLAG'=0) SFR estimation derived from H luminosity ('log_SFR_H') and the photometric stellar mass ('log_Mass') sourced from the pyPipe3D analysis catalogue, which is based on SDSS DR17 (Sánchez et al. 2022).
We present the results for ETGs once again in Fig. 12, this time comparing the differences brought by different classification methods.Under the criterion of luminosity-morphology parameter, sSFR and morphological classification method, there are about 2.2k, 3k and 2.2k ETGs, respectively.Clearly, apart from the differences in quantity, these three classification methods also exhibit distinct distributions in the corrected  T - e parameter space for ETGs.Under the luminosity-morphology parameter classification, ETGs display the most scattered density slopes, with a median  T value dropping to 1.7 at lg  e = 1.8, nearly same as the changing trend of the full sample.The total density slope values of ETGs obtained under the sSFR classification exhibit a more compact distribution, with the median  T decreasing from 2.2 to 1.9 as  e decreases, but this trend is still more prominent compared to that observed under the morphological classification.
As mentioned in Lu et al. (2021), the luminosity-morphology parameter classification is not accurate for S0 galaxies, and different classification methods have their respective limitations.However, here we aim to illustrate that the distribution of stars and dark matter in ETGs always exhibits clear layering with respect to the variation in T and the difference lies in the range delineated by different classification methods in the galaxy  T −  e parameter space.Disregarding the influence of classification methods and solely focusing on the distribution of stars and dark matter within galaxies, for high-velocity-dispersion MaNGA galaxies, while the median  T at a given velocity dispersion approaches  T ∼ 2, our results for the separated stellar density profile and dark matter density profile reveal that these two components do not 'conspire' to form an isothermal total density profile.To the contrary, as the stellar density deviates more from the average level, the fraction of dark matter in the galactic center decreases, and the total density profile deviates further from the isothermal (Fig. 6 and Fig. 7).The outcome of isothermal like density slope for early type galaxies appears to be a coincidental phenomenon in the mass growth process of galaxy-dark halo systems within a given velocity dispersion bin.

Discrepancy with hydrodynamical simulation
The scaling relationships for the total density profiles slopes of MaNGA galaxies obtained in this study cannot be replicated in the magneto-hydrodynamic cosmological TNG simulation.As shown in Fig. 2, for galaxies at high  200 and  e end, the TNG100 simulation underestimates the total density slope by ∼ 0.4, while at low  e end, both TNG100 and TNG50 significantly overestimate the total density slope of the galaxies.The differences can be partially attributed to the different  200 −  e and  e −  e relations between observations and TNG simulations (particularly evident in TNG100), but cannot be solely attributed to this factor.To correct the difference caused by  e , we derive interpolation functions for  e ( 200 ) and  e ( e ) based on the  e −  200 and  e −  e relations of MaNGA galaxies and recalculate the density slope with the fixed  e for TNG50 and TNG100 galaxies as discribed in Sec.4.2.The results of TNG50 are basically unchanged because the variation of its  e is similar to that of MaNGA.For TNG100, the difference on T with MaNGA at the high-mass end reduce from larger than 0.4 to smaller than 0.2, while in the low-mass range the difference still exists and is significant.The simulated trends in the  T - * (although we do not show in the text) and the  T - 200 relationships also differ markedly from the observations.The discrepancy becomes even larger when we separate and study the stellar density and dark matter components.Both the stellar density slope and the dark matter density slope are higher in numerical simulations.But why do the total density slopes of ob-served MaNGA galaxies, particularly those of larger mass, exceed those in the TNG100 simulation, despite the fact that their decomposed components exhibit lower slopes?The key lies in the fact that the proportion of dark matter in the central regions in the numerical simulations is significantly higher than in the observations.(referto Fig. 2, Fig. 4 and Fig. 11).The issue regarding the overestimation of the dark matter fraction in hydrodynamical simulations has been highlighted in previous studies.For the original Illustris simulation, Xu et al. (2017) found dark matter on average contributes about 40 -50 per cent to the (projected) total matter distributions at the centres of ETGs, which is higher than suggested by the observational results.Lovell et al. (2018) found a clear tension between the high dark matter fraction in TNG100 central galaxies and the low dark matter fractions within the central regions of the ATLAS 3D ellip- γ * within R e (σ e ) Figure 12.Same as Fig. 6, but we employ three distinct ETGs classification methods, indicated from left to right as the luminosity-morphology parameter, the specific star formation rate, and the morphological classification, each with its respective criterion listed in the corresponding tical galaxies.the present we confirm this discrepancy by utilizing a sample spanning a broader Theoretically, AGN feedback could alter the slope a galaxy's density profile.For instance, Peirani et al. (2019) discovered that incorporating AGN feedback in hydrodynamical simulations is essential for achieving improved agreement with observational values and the variation trends of  T .They also found the derived slopes are slightly lower than in the observational values when AGN is included because the simulated galaxies tend to be too extended, especially the least massive ones.Wang et al. (2020a) demonstrated that the kinematic wind feedback from AGNs could flatten the density profiles of ETGs and weakening the AGN feedback in numerical simulations might bring the resulting total density profiles closer to the observations of massive galaxies.However, given that the stellar density profiles of the simulated galaxies are already steeper than the observations as we see in Fig. 5, weakening the AGN feedback in the simulations would only exacerbate this difference.Therefore, it is hard to imagine that merely reducing AGN feedback could fully resolve the discrepancies we observed in this study.Indeed, a more critical factor causing the discrepancy between simulations and observations might lie in the assumption of an invariant IMF in the star formation processes within the TNG simulations.Within TNG, regardless of galaxy properties, the IMF consistently adheres to a bottom-light form as described by the Chabrier formula.However, previous studies suggest that the IMF of star formation could vary depending on the characteristics of the galaxy.Specifically, in ETGs with high velocity dispersion, the IMF may skew closer to a bottom-heavy situation (e.g.van Dokkum & Conroy 2010;Cappellari et al. 2012).Consequently, the stellar mass of high-velocitydispersion galaxies in simulations might be underestimated, leading to a higher derived dark matter fraction in the central regions.While the TNG simulations are calibrated against the observed stellar mass function of galaxies, this calibration also assumes a Chabrier IMF, which does not rectify the issue of galaxy mass underestimation.In Paper IV of our series of studies, we combined gravitational lensing observations with dynamical observations to reveal that the stellar mass estimates for central galaxies in groups and clusters, based on the Chabrier IMF, may be underestimated by approximately a factor of three.Furthermore, in Lu et al. (2023a) (Paper V), we found that the mass-to-light ratio of galaxies indeed aligns more closely with the predictions of a bottom-heavy model at the high-velocity-dispersion end.
It should be emphasized that our intention is not to assert that alterations in the IMF will necessarily result in a perfect agreement between observations and simulations.Instead, we aim to underscore the fixed IMF as a potentially significant factor contributing to disparities.In simulations, the recovery of the stellar mass function is often used as a benchmark for assessing simulation quality.However, when comparing with observations, simulations typically assume a fixed Chabrier IMF to compute the observed stellar mass function.This assumption makes it challenging to identify issues arising from a fixed IMF in such comparisons.Nevertheless, there have been simulations exploring the impact of variable IMFs on numerical simulation outcomes.For instance, Barber et al. (2018) introduced a variable IMF based on the initial interstellar medium (ISM) pressure and conducted simulations with fluctuating nucleosynthetic yields, star formation laws, and stellar feedback, all aligned with their locally varying IMFs.They concluded that a bottom-heavy IMF can effectively increase the stellar mass by introducing an excess of (dim) dwarf stars.Moreover, their research showed that conducting simulations with a variable IMF and making corresponding adjustments not only reproduces the observed correlation between mass-to-light ratio excess and central stellar velocity dispersion as reported in Cappellari et al. (2013b), but also achieves excellent agreement with the observational diagnostics initially employed to calibrate the subgrid feedback physics in the reference EAGLE model.
For dark halos less than 10 10 M ⊙ , the dark matter fraction in simulations is relatively close to the observed one.However, the stellar density profiles and dark matter density profiles in simulations are still steeper than observations (refer to the first row of Fig. 11).Lovell et al. (2018) showed that in TNG, the density profile of dark matter at 5-10 kpc might be affected by the adiabatic contraction of stars (Blumenthal et al. 1986;Gnedin et al. 2004), making it steeper than results from pure dark matter simulations.In this work, we found that the stellar density profiles of TNG100 and TNG50 in the low-mass range are significantly steeper than in MaNGA galaxies, implying that the dark matter distribution in TNG simulations would be affected by a stronger contraction of baryonic matter.This may be one of the sources of discrepancies between observations and simulations.
In the past, dark halos with  * less than 10 9 M ⊙ have been noted to have very flat dark matter halo density profiles, even flatter than the  DM = 1 prediction from the cold dark matter model (Bullock & Boylan-Kolchin 2017).Some studies suggest that this might be related to feedback from star formation, which can rapidly expel gas from the inner dark halo, thereby reducing the gravitational potential in the center of the halo and leading to a more flattened halo density profile in the center (e.g.Read & Gilmore 2005;Pontzen & Governato 2012;Benítez-Llambay et al. 2018;Bose et al. 2019).Other studies propose that self-interactions of dark matter might flatten dark matter density profiles (Kaplinghat et al. 2016;Tulin & Yu 2018;Robertson et al. 2019;Bondarenko et al. 2021;Andrade et al. 2022;Eckert et al. 2022).The galaxy masses involved in our study are approximately an order of magnitude higher than in these studies, so whether baryonic feedback still plays a significant role in these larger mass galaxies requires more detailed studies from future numerical simulations.

Sample selection
In Fig. 2, we also show that the total density slope derived from Paper IV by which included galaxies with a quality score of Qual ⩾ 0, whereas the results presented in this paper were based on galaxies with Qual > 0. If we include Qual = 0 galaxies in our analysis, we also predict a lower stacked  T for this mass range, consistent with the findings of Paper IV.We excluded galaxies with Qual = 0 in our study because of the large scatters in the total mass density slopes predicted by different mass models, even though these galaxies have a converged total mass estimate from dynamical modeling.However, Paper IV only used the total mass within a certain radius from dynamical modeling, which should not be affected by the scatter of density slope from different dynamical modeling models.
Fig. 13 shows the distribution of our galaxy samples with different Qual in the  T - * ,  T - e and  T - 200 diagrams, excluding samples with Qual = −1 since they are all untrustworthy and Qual = 0 since their density slopes are not reliable.Compared to the full sample, galaxies for different fitting quality groups exhibit distinct distributions.In most mass ranges, galaxies with Qual = 3 exhibit  T values that are 0.1 lower than galaxies with Qual = 1, while galaxies with Qual = 2 typically have  T values nearly all above 2.Although we do not show Qual = 0 galaxies, they have relatively shallower density slope than that of Qual ⩾ 1 galaxies.The reasons for the systematic differences in the matter distribution among galaxies with different modeling qualities need further investigation.

SUMMARY AND CONCLUSION
In this paper, we have investigated the distribution of galaxy mass by analyzing stellar kinematics modeling data from the MaNGA-DynPop project.Our study centers on approximately 6000 nearby galaxies, for which we employ reliable total density slopes obtained from the modeling work of Paper I. Furthermore, we explore the decomposed density profiles for both stellar and dark matter components of these galaxies.Our analysis has led to the following key findings: (i) Paper III presents the  T - e relationships derived from the analysis of these 6000 galaxies.The updated scaling relation reveals that, for galaxies with high  e values, the average mass-weighted total density slope ( T ) remains relatively constant at 2.2.However, as  e decreases below a specific threshold, the galaxy density profile becomes flatter.For the full sample, this threshold is approximately 189 km/s, while for ETGs, it is approximately 150 km/s.Here, we have also observed a strong resemblance between the trends in the variation of  T and the variation of  * .This similarity appears to be attributed to the dominance of stellar matter within  e of MaNGA galaxies.
(ii) We have compared our findings with central galaxies from TNG50 and TNG100 simulations.Our analysis reveals that the current magneto-hydrodynamic cosmological simulations do not yet accurately replicate the observed galaxy's central mass distribution.In particular, both  DM (<  e ) and  * tend to be overestimated in the simulations.Furthermore, the simulations fall short in reproducing the relationship between  T and galaxy's properties such as  * ,  200 , and  e .We attribute these discrepancies primarily to the simulations' tendency to over-predict the dark matter fraction, which could be linked to questionable assumptions like a constant IMF, excessive adiabatic contraction effects and feedback implementations.
(iii) We have observed that within a specific  e range, an increase in the stellar density slope corresponds to a higher total density slope.In this context, we have not observed dark matter compensating for the steepness of the stellar density profile to restore the density profile to a isothermal one ( T = 2).We conclude that there is no perfect conspiracy between baryonic matter and dark matter.
(iv) We have presented the stacked (the average) galaxy density profiles and calculated the changes in the stacked slopes with  * ,  e , and  200 .We find that in each sub-sample, the stacked dark matter density profile is slightly steeper than the pure dark matter simulation prediction of  −1 , which may indicate moderate adiabatic contraction in the central region of galaxy.
Our study underscores the effectiveness of stellar dynamics modeling as a valuable tool for investigating the interaction between stellar and dark matter, thereby constraining galaxy formation theories.In the future, we plan to further our investigations by integrating stellar dynamics techniques with both strong and weak gravitational lensing data to understand galaxy formation and evolution.
T - e relation is smoother than what is observed in SDSS-IV DR14.Based on the MaNGA Deep Learning morphological catalogue (Domínguez Sánchez et al. 2022) and adopt the most restrictive selection strategy listed in sec.3.4.1 of Domínguez Sánchez et al. (2022), we select 2180 ETGs from our MaNGA samples.The best-fitting equation for the 2. The trend of change is even flatter than the fitting result of  T - e relations plotted in Fig. 1 because here we have corrected for the varying impact of  e on the 12

Figure 2 .Figure 3 .
Figure 2. T as a function of  200(Left)  and  e (right) for central galaxies.The median of MaNGA, TNG50, and TNG100 is represented by green, red, and blue lines, respectively.The line width corresponds to the standard error of the median.The shaded region in the same color represents the 1- scatter.The four most massive galaxy of TNG50 are shown by red solid triangles in left panels.The data points with error bars in left panel represent the observational results fromEtherington et al. (2023) for SLACS galaxies,Newman et al. (2015) for central galaxies of group-and cluster-scale and Paper IV for central MaNGA galaxies of group-and cluster-scale combining JAM with weak lensing.The results from the C-EAGLE simulation(He et al. 2020) are also included.The dashed lines in the right panel represent the best-fitting results fromPoci et al. (2017) andLi et al. (2019), while the result for ETGs is from Paper III.The slope value of the isothermal profile is shown with a black horizontal dash-dot line.

Figure 4 .
Figure 4.  DM (<  e ) as a function of  * for central galaxies.For MaNGA galaxies, the black solid line shows their median and the green shadow region indicates the 1- scatter.For TNG50 and TNG100, the medians and 1- scatters are showed by scattered dots with error bars.

Figure 6 .Figure 7 .MFigure 8 .
Figure 6. T within  e (  * ) as a function of  * , colored by loess-smoothed  * within  e (  * ) (frac = 0.05), for full sample (left) and for ETGs.In each panel, the black line shows the median with the line width represents the stand error of the median, and the region enclosed by dashed lines indicates the 1- scatter.The blue dashed line marks the best fit for ETGs in Paper III.The black contours in each panel show the kernel density estimate for the sample distribution.The slope value of the isothermal profile is shown with a grey horizontal line.

Figure 9 .
Figure 9. Same as Fig. 8, but the mass bin is split by  200 .

Figure 10 .
Figure 10.Same as Fig. 10, but the mass bin is split by  e .

Figure 11 .
Figure 11.Top: stacked  DM ,  * and  T (from left to right) as a function of the mean  * for each stellar mass bin.The results of MaNGA, TNG50 and TNG100 are plotted by filled circles with error bars indicating the 1- scatter.Middle: the same as the top panels, but plotted as functions of the mean  200 for each halo mass bin.Bottom: the same as the top panels, but plotted as functions of the mean  e for each velocity dispersion bin.

Figure 13 .
Figure 13.The Qual distribution of MaNGA galaxies in  T - * ,  T - 200 and  T - e relation.Three different colored scatter points represent Qual = 1, 2, 3 groups, respectively.Their median values are indicated by solid lines with the line width represents the stand error of the median.The error bars enclose the 1- scatters.The slope value of the isothermal profile is shown with a black horizontal dash-dot line.
- e relations for the full sample and ETGs in Paper III are represented by blue and orange dashed lines, respectively.The result ofLi et al. ( T as a function of  e (left) and  * (right).Grey scattered points represent  T derived under the JAMcyl + gNFW assumption.The solid green line represents the median of  T , with the green shaded region denotes the 1- scatter.The medians of NFW  T are also plotted.The fitting result of  T