MaNGA DynPop – IV. Stacked total density profile of galaxy groups and clusters from combining dynamical models of integral-field stellar kinematics and galaxy-galaxy lensing

We present the measurement of total and stellar/dark matter decomposed mass density profile around a sample of galaxy groups and clusters with dynamical masses derived from integral-field stellar kinematics from the MaNGA survey in Paper I and weak lensing derived from the DECaLS imaging survey. Combining the two data sets enables accurate measurement of the radial density distribution from several kpc to Mpc scales. Intriguingly, we find that the excess surface density derived from stellar kinematics in the inner region cannot be explained by simply adding an NFW dark matter halo extrapolated from lensing measurement at a larger scale to a stellar mass component derived from the NASA-Sloan Atlas (NSA) catalogue. We find that a good fit to both data sets requires a stellar mass normalization about 3 times higher than that derived from the NSA catalogue, which would require an unrealistically too-heavy initial mass function for stellar mass estimation. If we keep the stellar mass normalization to that of the NSA catalogue but allow a varying inner dark matter density profile, we obtain an asymptotic slope of 𝛾 gnfw = 1 . 82 + 0 . 15 − 0 . 25 and 𝛾 gnfw = 1 . 48 + 0 . 20 − 0 . 41 for the group bin and the cluster bin respectively, significantly steeper than the NFW case. We also compare the total mass inner density slopes with those from TNG300 and find that the values from the simulation are lower than the observation by about 2 𝜎 level.


INTRODUCTION
In the modern lambda cold dark matter (ΛCDM) cosmogony, the structure forms hierarchically.Small dark matter haloes form early and then larger dark matter halos form through halo-halo merging and continuous dark matter accretion (Frenk & White 2012).N-body cosmological simulations show that in such an ΛCDM universe, the dark matter haloes obtain a self-similar universal structure of mass distribution (e.g.Navarro et al. 1996b;Navarro et al. 1997;Springel et al. 2008;Gao et al. 2011) that the spherically averaged dark matter mass profile follows () ∝  −1 at the inner part and () ∝  −3 at ★ E-mail: chunxiang_wang@sina.cn † E-mail: ranl@bao.ac.cn the outer part.In a real universe, galaxies form at the cente of dark matter haloes and co-evolve with dark matter haloes.The baryonic process during the galaxy formation and evolution can modify the mass distribution of the dark matter halo, especially at the inner part, through baryonic condensation and contraction (Blumenthal et al. 1986;Gnedin et al. 2004a;Gustafsson et al. 2006;Duffy et al. 2010;Schaller et al. 2015), the expulsion of gas during the feedback process (e.g.Navarro et al. 1996a;Read & Gilmore 2005;Pontzen & Governato 2012).The impact of these processes is complex and the net effect is not clear yet.Accurate measurement of the total mass density profile, as well as the decomposed stellar and dark matter density profiles, can provide a unique tool to probe the galaxy formation process.
The density profiles of real galaxies were measured mainly using dynamical modeling of the gas or stellar kinematics (see review by Courteau et al. 2014), or using gravitational lensing (see review by Treu 2010).
The presence of gas in a nearly circular motion in the equatorial plane of spiral galaxies made these the first targets for studies of rotation curves using either ionized gas at optical wavelengths (e.g.Rubin et al. 1980) or neutral Hi in the radio (e.g.Bosma 1978).These studies found that the rotation curves of spiral galaxies tend to be flat at large radius (beyond about 4 projected half-light radius  e ) and provided one of the first convincing evidence for the presence of dark matter inside galaxies (Faber & Gallagher 1979).The general flatness of rotation curves of spiral galaxies at large radius was later confirmed by numerous papers (e.g.Martinsson et al. 2013).
If one assumes the total density of galaxies to be wellapproximated by power laws of the form  T () ∝  − T , there is a simple relation between the logarithmic slope  vel of the circular velocity and the one  T of the density (Binney & Tremaine 2008, eq. 2.61) (1) This implies that the observed flat rotation curves ( vel ≈ 0) of spiral galaxies suggest their total densities of spiral galaxies to be close to isothermal ( T ≈ 2) at large radius.Unlike spiral galaxies, early-type galaxies (ETGs) generally do not contain extended Hi discs.For ETGs the total density slopes were measured via strong gravitational lensing or stellar dynamics.The requirement for galaxies to act as a lens implies that samples of lenses galaxies tend to be ETGs with large masses and effective velocity dispersion  e .For these lens ETGs the mean power-law slope out to a median radius of  e /2 was found to be ⟨ tot ⟩ = 2.078, or nearly isothermal, in the project Sloan Lens Advanced Camera for Surveys (SLACS) with a sample of 73 strong galaxy lenses (Koopmans et al. 2009;Auger et al. 2010b).A similar value was reported in other strong lensing projects at higher redshift  ∼ 0.5 (Li et al. 2018b).
Unlike strong gravitational lensing, stellar dynamical modeling can be applied to statistically-representative and larger samples of galaxies, without the need for the galaxies to be lenses.Moreover, one can sample different radius, not just close to the lens.However, stellar kinematics is challenging to measure out to large radius due to the low surface brightness.Various, early studies examined the density profiles of individual ETGs (e.g.Weĳmans et al. 2008Weĳmans et al. , 2009;;Forestell & Gebhardt 2010;Morganti et al. 2013;Napolitano et al. 2014).Cappellari et al. (2015) conducted the first systematic investigation of the density profile of ETGs extending to large radius.They used the Jean Anisotropic Models (JAM, Cappellari 2008) to analyze extended stellar kinematics for 14 massive ETGs from the SAGES Legacy Unifyinig Globulars and GalaxieS (SLUGGS) survey (Brodie et al. 2014) out to a median radius of 4 e and found they are well described by power laws with the logarithmic slope of the total density profile slightly steeper than isothermal ⟨ tot ⟩ = 2.19 with a scatter of just   = 0.11.This sample was later extended to 21 ETGs by Bellstedt et al. (2018) who found a very similar value ⟨ tot ⟩ = 2.24 with a scatter of   = 0.05.This mean value and the small scatter were confirmed for a sample of 16 massive ETGs with Hi gas discs out to a median radius of 6 e by Serra et al. (2016), who found a remarkably close mean slope ⟨ tot ⟩ = 2.18.Subsequently, Poci et al. (2017) used JAM to model a much larger sample of 260 ETGs from the ATLAS 3D survey (Cappellari et al. 2011), but only out to a median radius of 1 e .They found that above a stellar dispersion lg( e /km s −1 ) ≈ 2.1, these galaxies have a mean total density slope of ⟨ tot ⟩ = 2.193, with an observed scatter of   = 0.17, in ex-cellent agreement with, the previous values, and consistently slightly larger than the lensing value.However, at lower  e a trend was discovered with the slope decreasing with  tot .Moreover, the decrease of  tot was found to correlate better with  e than with stellar mass (Cappellari 2016, fig. 22c).
The Sloan Digital Sky Survey (SDSS) Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey (SDSS-MaNGA, Bundy et al. 2015) substantially expanded the size of galaxy samples with integral-field stellar kinematics and included both ETGs and spirals.Using these data Li et al. (2019) derived the mass-weighted total density slope out to a median radius of 1.5 e with JAM for about 2K nearby galaxies, and again accurately reproduced the mean total inner density slope of ⟨ tot ⟩ = 2.24 for ETGs with  e ≳ 100kms −1 , below which the density slope decreases with  e as previously noted.However, their study was the first to model with a consistent method a large sample of both spiral galaxies and ETGs.They found a much larger variation and much clearer trend in the total slopes, with  tot decreasing, namely becoming more shallow, for spiral galaxies vs ETGs.
We revisited the  tot trends in Zhu et al. (2023a), using JAM modeling of a sample of about 6K galaxies with the best data quality, extracted from the final data release of the MaNGA survey.We again confirmed the nearly constant mean slope ⟨ tot ⟩ = 2.20 above lg( e /km s −1 ) ≈ 2.2 but also found a clear variation of  tot at fixed  e , with the slopes decreasing for younger ages.
Studies on galaxy groups and galaxy clusters(e.g.Newman et al. 2015Newman et al. , 2013) ) that combine stellar kinematics and weak gravitational lensing have shown that the total density slope may also be shallower than 2, reaches ∼ 1.7 for galaxy groups and ∼ 1.2 for galaxy clusters.
Unlike the total density distribution, the decomposed dark matter and baryonic density profiles usually cannot be measured reliably.The decomposition of the two density components is challenging when observational measurements are available only within about 1 e , due to the model degeneracy between the total stellar mass and the dark matter inner density slope.It is important to have observational data from multiple scales, which helps to obtain constraints on the mass of the dark matter halo, and shrinking the freedom of dark matter density profile (e.g.He et al. 2020;Newman et al. 2015).Sonnenfeld et al. (2012) derived dark matter density slope  = 1.7 ± 0.2 for SDSSJ0946+1006, the 'Jackpot' lens, by combining the stellar kinematics and the double Einstein ring data of the system, implying a scenario of dark matter contraction at the halo center.A value  = 1.602 ± 0.079 syst , also consistent with halo contraction, was obtained for the Milky Way using APOGEE and Gaia data out to 5 e in what represents the most accurate determination for any galaxy, due to the availability of full six-dimensional stellar phase space information (Nitschai et al. 2021).On the other hand, Yang et al. (2020) derived decomposed mass models by combining stellar kinematics at inner one effective radius (1 e ) and Hi kinematics within 5 e , and obtained a dark matter inner density slope  dm = 0.6 +0.3 −0.2 for NGC2974, a bright nearby ETG, which is much shallower than the standard Navarro-Frenk-White (NFW) value (1).Newman et al. (2015) presented the dark matter inner density slope for a sample of galaxy groups and clusters, they found that the galaxy groups have a mean  dm ∼ 1, consistent with NFW, but the clusters have a mean  dm ∼ 0.5, significantly lower than the NFW prediction.Recently, Sartoris et al. (2020), however, derive a  dm = 0.99 ± 0.04 for Abell S1063.Overall, the measurement cases of decomposed density profiles are not enough and the results are controversial, thus independent observations are valuable to clarify the situation.
In this work, we performed measurement of mass density profiles for a sample of galaxy groups and clusters by combining stellar kine-matics from the integral field unit (IFU) data of the SDSS-MaNGA survey and the weak gravitational lensing measurement with Dark Energy Camera Legacy Survey (DECaLS; Dey et al. 2019) shear catalogue.MaNGA data provides us the mass distribution information within the effective radius of the group/cluster central galaxy, and the stack galaxy-galaxy lensing measurement can constrain the mean density profile of the dark matter halo from ∼ 100 kpc to a few Mpc.Combining the two data sets, we derive the total mass density profile and decomposed stellar/dark matter profile for the selected galaxy groups and clusters.
The structure of this paper is organized as follows.In Section 2, we describe the data set.In Section 3, we show the methods used in this work, including the dynamical model of the galaxies in Section 3.1 and the gravitational lensing measurement in Section 3.2.In Section 4, we present our results.Discussions and conclusions are shown in Section 5. Throughout the paper, we adopt a flat ΛCDM cosmological model from the Planck 2015 results (Planck Collaboration et al. (2016), Ω m = 0.3075, H 0 = 67.74kms −1 Mpc −1 ).

OBSERVATIONAL DATA
In this project, we measure the lensing signal around galaxy groups and galaxy clusters in the overlapping region of the MaNGA and DECaLS survey, the latter of which provides the source catalogue of weak lensing measurements.We describe the data sets in this section.

Lens galaxies
We select our lens sample by matching the ETGs in MaNGA with the central galaxies of groups/clusters in the SDSS DR7 group catalogue (SDSSGC; Yang et al. 2005;Yang et al. 2007).MaNGA (Bundy et al. 2015) is a multi-object IFU spectroscopy survey that makes use of the 2.5-meter Sloan Foundation Telescope and the Baryon Oscillation Spectroscopic Survey (BOSS) spectrograph (Smee et al. 2013).
The BOSS spectrograph provides continuous coverage between 3600 Å and 10300 Å with a spectral resolution R∼2000.The full sample has been selected at low redshift (0.01 <  < 0.15) to follow a flat distribution of stellar mass across the range of 10 9 M ⊙ − 10 12 M ⊙ .
We draw our lens sample from the final data release of MaNGA, which contains 10010 unique galaxies in the SDSS Data Release 17 (DR17, Abdurro'uf et al. 2022).The MaNGA data cubes are obtained by the spectrophotometric calibration (Yan et al. 2016) and the Data Reduction Pipeline (DRP; Law et al. 2016).For each data cube, the spaxels are Voronoi-binned (Cappellari & Copin 2003) to reach a target S/N ∼10.The MaNGA Data Analysis Pipeline (DAP; Westfall et al. 2019), which uses ppxf (Cappellari 2017) and a subset of MILES library (Sánchez-Blázquez et al. 2006), MILES-HC, extracts the stellar kinematics for each binned spectrum.The stellar kinematics from DAP products publicly available since SDSS Data Release 17 1 provide the spatial distribution of the projected stellar velocity and stellar velocity dispersion for each galaxy.Zhu et al. (2023b, hereafter Paper I) used the MaNGA stellar kinematics to construct accurate dynamical models of these galaxies and we use the result from the models in this paper.
The following criteria are used in the lens selection process, • Θ M,S < 1.0 arcsec.
1 https://www.sdss.org/dr17/manga/manga-data Here, the Θ M,S is the angular separation between the MaNGA galaxy and the central galaxies in SDSSGC, where the central galaxies are defined as the galaxy with the largest stellar mass, and Δ is the difference between the corresponding redshifts from the two catalogues.Paper I assigns a quality grade (Qual = −1, 0, 1, 2, 3) to each MaNGA galaxy.The higher the grade is, the better the model fitting quality is.We discard those galaxies with Qual = −1 whose kinematic/dynamic properties can not be trusted.For galaxies with Qual ⩾ 0, we further require the lg  T (<  e ) gnfw,sph and lg  T (<  e ) gnfw,cyl to be consistent within 1 ( 0.13 √ 2 dex, derived from lts_linefit2 software by fitting galaxies with Qual=0 without clipping), where the  T (<  e ) gnfw,sph and  T (<  e ) gnfw,cyl are the total mass within a sphere of radius  e measured by dynamical models assuming either a spherically-aligned "JAM sph "+ "gNFW" model or cylindrically-aligned velocity ellipsoid "JAM cyl "+ "gNFW" model, respectively.We refer readers to Paper I for model details.
In this work, we only use groups/clusters whose central galaxy is ETG defined in Domínguez Sánchez et al. (2022).We removed the late-type galaxies to suppress the mis-center effect (Gao & White 2006;Leauthaud et al. 2010;Shan et al. 2017;Wang et al. 2018).We further require the lens sample to overlap with the survey footprint of DECaLS, and we remove lenses with redshift  < 0.02 to ensure an efficient measurement of galaxy-galaxy lensing.We divide these central galaxies into 3 bins according to their assigned halo mass  200m in the SDSSGC, where  200m is defined as the total mass enclosed in  200m within which the mean density is 200 times of the mean matter density of the universe at the redshift of the halo.The mass range and the number of lenses in different bins are listed in Table 1.
In Fig. 1, we show the stellar mass and effective radius of our lens galaxies and the full sample of MaNGA ETGs, where the values are derived from the NASA-Sloan Atlas (NSA) catalogue3 (Blanton & Roweis 2007;Blanton et al. 2011).The  e is the circularized effective radius which is calculated from the multi-Gaussian expansion (MGE; Emsellem et al. 1994;Cappellari 2002) formalism of the galaxy band luminosity distribution. e is finally scaled by a factor of 1.35 to match the values determined from 2MASS (Skrutskie et al. 2006) plus RC3 (de Vaucouleurs et al. 1991).In Fig. 1, the three subsamples with increasing halo mass are represented by blue, orange, and green dots.The full sample of MaNGA ETGs is depicted by gray dots.Contours show the full sample of MaNGA ETGs (solid lines) and the combined there sub-samples used in this study (dashed lines) at 30%, 60%, and 90% probability levels.One can see that there is only a slight difference between the contours of these two samples, indicating that the selected sources represent well the population of ETGs in the MaNGA full dataset.
In this work, we focus on the median and the high mass bins, which we will refer to as the group bin and the cluster bin respectively.The observational lensing error of the lowest mass bin is too large to derive meaningful results, as shown in the Appendix.Distributions of the number of member galaxies from SDSSGC in the group bin and the cluster bin are shown in Fig. 2.

Source galaxies
In this work, we use the Dark Energy Camera Legacy Survey4 (DE-CaLS, see Dey et al. 2019) shear catalogue which has been utilized in Table 1.Posterior constraints of free parameters of different models.The first three columns, show the halo mass range, the number of lenses, and the data set and model.The following four columns show the fitting parameters, halo mass, concentration, stellar mass normalization, dark matter fraction within the mean value of  e , and reduced  2 /.The best-fit value represents the peak of the one-dimensional marginalized posterior distribution of one parameter.We determine a horizontal line that intersects the marginalized distribution of the parameter at two points, such that the probability between these two intersections sums up to 68%.These two intersections correspond to the upper and lower limits of the 1 interval.The photo- of each source galaxy in DECaLS DR8 shear catalogue is taken from Zou et al. (2019), where the redshift of a target galaxy is derived with its k-nearest-neighbor in the SED space whose spectroscopic redshift is known.The photo-z is derived using 5 photometric bands: three optical bands, , , and  from DECaLS DR8, and two infrared bands, W1, W2, from Wide-Field Infrared Survey Explorer.By comparing with a spectroscopic sample of 2.2 million galaxies, Zou et al. (2019) shows that the final photo-z catalogue has a redshift bias of Δ norm = 2.4×10 −4 , the accuracy of  Δ norm = 0.017, and outlier rate of about 5.1%.

Dynamical Model
In this paper, we make use of dynamical models lens galaxies derived in Paper I, where the JAM (Cappellari 2008(Cappellari , 2020) ) is performed for the final SDSS-MaNGA data release.The JAM method has been applied to the mock stellar kinematic data from cosmological simulations and is demonstrated to be robust in recovering the total mass profile (Lablanche et al. 2012;Li et al. 2016).Using the JAM method, Paper I predicts second velocity moments maps with an assumed parametric mass distribution model and fits them to the observed one extracted from MaNGA.
To investigate the systematics introduced by different theoretical assumptions and different parametric forms, Paper I uses eight mass models to fit the observed stellar kinematics.They adopt two extreme assumptions on the orientation of velocity ellipsoid, i.e.JAM cyl (cylindrically-aligned) and JAM sph (spherically-aligned), and for each orientation case, they adopt four parametric mass distributions, including one mass-follow-light mass model and three stellar-darkmatter two-component models.
In this work, we do not fit stellar kinematics and weak lensing signal simultaneously, mainly because the dynamical systems cannot be stacked linearly as the galaxy-galaxy lensing signal.Instead, we use the average value of  dyn (<  0 ), the JAM derived 3D total mass enclosed in a radius of  0 , for different mass ranges, as the observational constraints.We note that an alternative approach would consist of constraining the outer density to follow the average values for the halo mass when fitting the dynamical models of individual galaxies with JAM.
The primary sample in SDSS-MaNGA has IFU data within 1.5  , within which the dynamical total mass estimation are most reliable, thus we set  0 to be the 1.5 ⟨  ⟩ of the central galaxies. e is the effective radius.For the group and the cluster bins,  0 are 13.57and 20.10 kpc respectively.
For each halo mass bin, we calculate the  dyn (<  0 ) using the "JAM cyl "+ "gNFW" model of Paper I which is the most flexible model in Paper I. We estimate  ,dyn , the uncertainties of  dyn (<  0 ) as follows.First, we use the bootstrap method to estimate the statistical uncertainties  stat using the total mass estimated by "JAM cyl "+ "gNFW" model.Then, we estimate the systematic error  sys of  dyn (<  0 ) induced by imperfect model assumptions

Galaxy-Galaxy Lensing
We measure the stacked galaxy-galaxy lensing signal in 10 logarithmic radial bins from 0.1 to 10 Mpc, where the excess surface density, ΔΣ() is derived as where Σ(< ) is the mean density within the radius , the Σ() is the azimuthally averaged surface density at radius  (e.g.Miralda-Escude 1991; Wilson et al. 2001;Leauthaud et al. 2010),  ls t is the tangential shear, and where the critical surface density Σ crit can be written as where  s ,  l , and  ls are respectively the angular diameter distance between the observer and the source, the observer and the lens, and the source and lens, and  is the constant of light speed in the vacuum.The weight factor  n is introduced to account for intrinsic scatter in ellipticity and shape measurement error of each source galaxy (Miller et al. 2007(Miller et al. , 2013)), defined as  n = 1/( 2  +  2 e ), where   = 0.27 is the intrinsic ellipticity dispersion derived from the whole galaxy catalogue (Giblin et al. 2021). e is the error of the ellipticity measurement defined in Hoekstra et al. (2002).To suppress the dilution effect from the photo- uncertainties of the source galaxies, we remove the lens-source pairs with  s −  l < 0.1.
We apply the correction of multiplicative bias to the measured excess surface density as where where  is the multiplicative bias as described in Sec.2.2.The lensing signal is multiplied by boost factor () = ()/ rand (), which is the ratio of the number density of sources relative to the number around random points, in order to account for dilution by sources that are physically associated with lenses, and therefore not lensed (Mandelbaum et al. 2005;Mandelbaum et al. 2006).Throughout this work, we use the Super W Of Theta (SWOT) code 6 (Coupon, J. et al. 2012) to calculate the excess surface density.

Joint constraint
The  2 of lensing can be written as We model the stacked excess surface density ΔΣ() as where the first term represents the contribution of the central baryonic component, the second term represents the contribution of host projected dark matter haloes of the galaxies, and the third term represents the contribution of neighboring halos, namely the 2-halo terms.
In many previous lensing analyses, the contribution of the baryonic component is often modeled as a point mass, the value of which is set to the sum of the stellar mass of the galaxies.However, the approximation is not accurate within around 2  .Following Paper I, we assume the mass distribution of the baryonic component follows the r-band light distribution, the form of which is derived by fitting the Multi-Gaussian Expansion (MGE; Emsellem et al. 1994;Cappellari 2002) to the r-band image of the galaxies in SDSS DR17.The excess surface density of the stellar component can be written as where where Σ * ,nsa is the stellar mass distribution which follows the r-band MGE brightness distribution derived from Paper I, with normalization fixed to the stellar mass  * ,nsa , which adopts a Chaberier IMF. nsa is a normalization parameter which describes the probable mismatch between Σ * ,nsa and the ground truth.In our "ggl only model" and "ggl+dyn" model (defined below), we fix  nsa = 1.
In our fiducial model, we use the NFW profile (Navarro et al. 1996b) to model the dark halo component: where  s is the scale radius where the local logarithmic slope d ln  d ln  = −2, which can be derived from dark matter halo radius through the concentration parameter  200 =  200 / s .
The excess surface density ΔΣ NFW () can be calculated by integrating the three-dimensional density profile along the line of sight, which we assume aligned with the  axis, as follow: For the NFW model, we use the analytical expression for ΔΣ presented in Oaxaca Wright & Brainerd (1999) to perform the model fitting.The 2-halo term ΔΣ 2h can be written as 6 http://jeancoupon.com/swot where the tangential shear profile due to the neighboring halos (Oguri & Hamana 2011) is: where  2 is the second-order Bessel function,  m () is the mass density at ,  A () is the angular diameter distance,  m () is the linear matter power spectrum,  h () is the halo bias derived by Tinker et al. (2010).
As mentioned in subsection 3.1, we do not stack kinematic data as we do for the lensing data, instead, we incorporate the dynamical data of MaNGA galaxies by adding a term of  2 dyn as Therefore, the total  2 of the joint analysis can be written as The averaged stellar mass within  0 is given by where  * ,JAM (<  0 ) = ⟨/⟩ * ,JAM  (<  0 ) is the stellar mass of JAM model within  0 .⟨/⟩ * ,JAM is the stellar mass-to-light ratio (/) of JAM model and (<  0 ) is the luminosity within  0 derived from Paper I. Different data sets and mass models are described as follows.
• ggl only: only the weak gravitational lensing is used in fitting.The stellar mass normalization  nsa = 1, and the model has two free parameters,  200 and  200 .
• ggl+dyn: both dynamical data and weak gravitational lensing signal are used in fitting.The stellar mass normalization  nsa = 1.
• ggl+dyn+free ml: both dynamical data and weak gravitational lensing signal are used in fitting.The mass model has three free parameters,  200 ,  200 , and the stellar mass normalization  nsa .
• ggl+dyn+gnfw: similar to ggl+dyn, but the gNFW model describe is used instead of NFW model for the host halo component.The stellar component is fixed as  nsa = 1.
We adopt the Markov chain Monte Carlo (MCMC) technique to perform the modeling fitting and calculate the posterior distribution.We use 24 chains of 300,000 steps with the MCMC ensemble sampler Emcee7 (Foreman-Mackey et al. 2013).
Studying galaxy properties in sub-sample bins with stacked galaxy-galaxy gravitational lensing method may introduce biases (Sonnenfeld & Leauthaud 2018).We performed a test using halo masses that are identical to the sub-samples used in this work and assumed they follow an NFW profile distribution and a massconcentration relation.We then stacked the signals and fit them using the NFW model.The results showed a bias between the best-fit model parameters ( 200 ,  200 ) and the true mean value.However, the true mean values were consistent with the fitting results within the 1 error range.Therefore, the bias introduced by the stacked method will not change our conclusion.

RESULTS
We present the gravitational lensing signal and the best-fit models in Figs. 3 and 4 for the two halo mass bins.We also show the decomposed components of the best-fit model, namely the dark matter halo (dashed), stellar component (dotted), and the two halo term (dashdot).Each figure contains four different panels, showing the fitting results of four different models.The posterior distribution of the model parameters can be found in Figs. 5, 6, and 7.The best-fit parameters with uncertainties are listed in Table 1 and Table 2.We derive mean halo mass of lg( 200 [M ⊙ ]) = 13.07 +0.13 −0.14 and 13.87 +0.09 −0.12 for the group and the cluster bins respectively (ggl+dyn model).The best-fit values of the halo mass are stable among different mass model assumptions, with differences within 1 error, which reflects the fact that the total mass of a halo is mainly constrained by the weak lensing data which measures the density distribution at a larger scale.
On the other hand, the total density profiles at the inner 100 kpc region depend strongly on whether stellar kinematics is included in the fitting.If we use weak lensing alone to constrain the mass model with an NFW profile and the stellar mass normalization fixed with NSA catalogue, we get  200 = 2.73 +3.80  −1.70 for 10 13 − 10 14  ⊙ bin, and  200 = 4.54 +2.39  −1.66 for > 10 14  ⊙ bin, while the best-fit model predicts an amplitude of ΔΣ significantly lower than the result from MaNGA stellar kinematics.If we use the same mass model to match both inner stellar kinematic data and the weak lensing data (ggl+dyn case), the values of  200 raise to 13.26 +3.02  −1.52 , and 7.01 +1.90 −1.49 for the two bins respectively, significantly higher than that predicted by Nbody cosmological simulations (Duffy et al. 2008) as well as other weak lensing measurements (see Fig. 8).
The discrepancy between the weak lensing and stellar kinematics can be alleviated by allowing the variation of stellar mass contribution, by relaxing the assumption of a universal IMF, or the inner density slope of dark matter halo.In the bottom panels of Figs. 3 and  4, we show the fitting results of these two new mass models.One can find that with additional free parameters, the mass model can fit both data set well.
In Fig. 6, we show the posterior distribution for the mass model with a free stellar mass normalization parameter.By combining the weak lensing and stellar kinematic data, we find best-fit values  nsa = 3.34 +0.52  −0.73 for the group bin and  nsa = 2.87 +1.12 −1.38 for the cluster bin.In both cases, an / ∼ 3 times higher than that given by the NSA catalogue is preferred.In the figure, we also mark the mean stellar normalization derived from the Paper I catalogue using the red vertical line, where the decomposed stellar mass distribution is derived from stellar kinematics alone.For the group bin, the best-fit value of stellar mass normalization of this work is slightly higher than that derived from Paper I, but consistent at ∼ 1 level.For the cluster bin, the two results agree well with each other.
In this project, we make the assumption that the stellar mass surface density is proportional to the MGE surface brightness model and normalize it to the stellar mass derived from the NSA catalogue.However, it should be noted that the total luminosity used to derive the stellar mass by the NSA catalogue and that from the MGE model may differ slightly, which can also contribute to the normalization factor  nsa value but this difference does not affect the results of "ggl+dyn+ml free" model, thus won't change the conclusion of this paper.
Allowing the variation of stellar mass normalization has a direct impact on the best-fit value of the concentration parameter for the group bin.When using a free  nsa , the concentration parameter decreases from  200 = 13.26+3.02  −1.52 to  200 = 2.32 +3.78 −1.52 which is close to the value from using galaxy-galaxy lensing alone.The bestfit value of concentration for the cluster mass bin also decreases, but not significantly.In Fig. 8, we compare concentration derived from this work with mass-concentration relations measured in weak lensing surveys and numerical simulations (Shan et al. 2017;Mandelbaum et al. 2008;Duffy et al. 2008;Klypin et al. 2016;Umetsu et al. 2014;Covone et al. 2014;Merten et al. 2015).Here, we correct the 2D concentration to 3D concentration with the relation  2D () =  3D () × 1.630 −0.018 provided by Giocoli et al. (2012), who found that halo triaxiality and substructures within the host halo virial radius can bias the observed 2D mass-concentration relation.The concentrations are all scaled to  = 0 by dividing (1+) −0.67 , assuming the redshift evolution from Klypin et al. (2016).One can find that the concentration derived with stellar mass normalization fixed to the NSA catalogue ("ggl+dyn" ) is significantly higher than the previous lensing observation and the prediction of numerical simulation, while the results derived with "ggl+dyn+free ml" agree with these previous measurements within 1.
Different models also predict different dark matter fraction for the central region.Table 1 lists the dark matter fraction within the average  e of the sample for different models as  dm (< ⟨ e ⟩).For both mass bins, the "ggl+dyn" model yields a high dark matter fraction, with  dm (< ⟨ e ⟩) = 0.66 +0.02 −0.03 for the group bin and  dm (< ⟨ e ⟩) = 0.68 +0.05 −0.05 for the cluster bin.In contrast, the "ggl+dyn+free ml" model predicts much lower values, with  dm (< ⟨ e ⟩) = 0.06 +0.14 −0.04 for the group bin and  dm (< ⟨ e ⟩) = 0.23 +0.24  −0.11 for the cluster bin, which agrees with the values derived from stellar kinematic data alone in Paper I.
In Figs. 3 and 4, we also investigate whether a steeper asymptotic inner density slope can also explain the observational data.Since the inner density slope degenerates strongly with stellar mass normalization parameter, we choose to set the latter to two fixed values during the fitting, the NSA value, where  nsa = 1; the value derived with JAM model from Paper I, where  nsa = 2.57 and 2.49 for the group and the cluster bins respectively.We present the best-fit asymptotic density slope,  gnfw in Table 2.If we fix the stellar mass normalization to that derived from the NSA catalogue, we obtain a best-fit value of density slope  gnfw = 1.82 +0.15  −0.25 and  gnfw = 1.48 +0.2 −0.41 for the group and the cluster bins respectively.The case of the NFW profile ( = 1) is at the 10th percentile and 23th percentile of the posterior distribution of  gnfw for the group and cluster bins.If we choose to fix the stellar mass normalization to the mean value derived from JAM instead, we obtain  gnfw = 1.57+0.16  −0.43 ( = 1 is the 21th percentile ), and  gnfw = 1.21 +0.28 −0.59 ( = 1 is 53th percentile).For the gNFW profile, the inner density slope degenerates with scale radius  s (Dutton & Treu 2014;He et al. 2020), so we also calculate mass-weighted density slope within  e to quantify the shape of the density profile, where where () is the mass density and the  ( e ) is the mass within the radius R e .The best-fit values of  are shown in Table 2, and in the left panels of Fig. 9, we show the posterior distribution of massweighted density slope of dark matter component.For the group bin,  dm (<  e ) disfavor the NFW model prediction if we set  nsa = 1, but the two agree within 1, if  nsa is set to the values from Paper I. For the cluster bin,  dm (<  e ) is not tightly constrained, thus the predictions from models with different  nsa setting broadly agree with each other.
In Fig. 9, we also show the posterior distribution for the massweighted total density slope.As expected, the best-fit total density slope does not change with the choice of  nsa .For the lower mass bin, the total density slope  tot (<  e ) is 2.12 +0.05 −0.09 , slightly steeper than the singular isothermal case, while the slope for the cluster mass bin is flatter (1.93 +0.06 −0.10 ), which may due to the increasing dark matter fraction in the cluster scale haloes.
In Fig. 9, we also present the mass-weighted density slopes obtained from the cosmological hydrodynamical simulation, specifically the IllustrisTNG project (referred to as TNG300), with red stars.We select central galaxies from the TNG300 simulation (snap-shot=99,  = 0) whose FoF halo mass  200 falls within the 1 range of the halo masses in our lens sub-samples.The effective radius of the central galaxies in the  band, projected along the -axis, is obtained from Genel et al. (2018).The enclosed mass  ( e ) and the mass density ( e ) in Equ.18 are directly calculated from the simulation data without any mass density model fitting.For both cluster and group mass bins, the mass-weighted average dark matter density slopes from TNG300 are consistent within the 1 confidence levels with those predicted by mass models with differently fixed  nsa .The mass-weighted total density slope is consistent with prediction of TNG300 at 2 confidence levels when  nsa = 1, while the prediction of TNG300 is 2 lower than that of our model fitting when  nsa is set to the values from Paper I.
In Fig. 9, the pink stars represent the stacked mass-weighted den-sity slopes of the solely stellar kinematic data with JAM model fitting in Paper I. Results obtained solely from stellar kinematic fitting are consistent within a 1 range with those when  nsa is set to the values from Paper I.

DISCUSSION AND CONCLUSIONS
In this paper, we presented the joint analysis of the stellar kinematic data and the galaxy-galaxy lensing data for a sample groups and clusters in the overlapping region between the MaNGA survey and the DECaLS DR8 imaging survey.Combining the two data sets allows us to derive the mean halo mass and measure the radial density profile from 10 kpc to Mpc for the two sample bins.Intriguingly, we find the excess surface density derived using stellar kinematics cannot be naturally explained by adding an NFW halo derived from galaxy-galaxy lensing alone to a stellar mass component whose normalization is fixed to the value derived by the NSA catalogue.
To match the high surface density derived from stellar kinematic data, the best-fit NFW profile requires concentration parameter  200 = 13.26+3.02  −1.52 (group bin) and  200 = 7.01 +1.9 −1.49(cluster bin), which is significantly higher than that predicted by the massconcentration relation derived by Duffy et al. (2008) using cosmological numerical simulation.
By setting the stellar mass normalization as a free parameter, we obtain a better fit to the observational data and draw the concentration parameter to the normal amplitude.We find that observational data  allows us to put constraints on the stellar mass normalization  nsa .
In both the group and the cluster bins, the fitting favors a stellar mass normalization ∼ 3 times higher than that given by the NSA catalogue.
Much of the difference between the observation required stellar mass normalization and that given by the NSA catalogue can be alleviated by using a more bottom-heavy initial mass function (IMF).The stellar mass of the NSA catalogue is derived by fitting stellar population templates to the broadband photometry data with a Chabrier initial mass function (Chabrier 2003).Migrating the IMF from a Chabrier IMF to the Salpeter IMF (Salpeter 1955) on average can cause a 0.25 dex higher results of total stellar mass.Many recent observations show that the IMF may indeed vary as a function of velocity dispersion for the ETGs (see a review of Smith 2020), and that the massive ETGs have a heavier IMF than that of Milky Way (Auger et al. 2010a;Conroy & van Dokkum 2012;Cappellari et al. 2012Cappellari et al. , 2013;;Spiniello et al. 2014;Li et al. 2017;Lu et al. 2023b).Gnedin et al. 2004b) when stellar mass is obtained with a bottom-light Kroupa IMF (Kroupa 2001).In order to fully explain the observation without AC, the stellar masses would need to increase by a factor of two, meaning that a bottom-heavy Salpeter IMF would be required.
In this work, we derive a  nsa = 3.34 +0.52 −0.73 and  nsa = 2.87 +1.12 −1.38 for the group and cluster bin respectively, which are slightly larger, but consistent with the value from pure stellar kinematic results of the MaNGA observation (vertical solid line in the histogram).
Given the condition of the data, in this paper, we did not consider the / gradient of galaxies and assumed that the / does not change with radius.However, if galaxies inherently have a stellar / gradient (for example, Sonnenfeld et al. (2018) using the Bayesian hierarchical modeling method identified a gradient in galaxy /), then our results on stellar mass and dark matter decomposition would be biased.If the stellar / of a galaxy decreases with the galaxy's radius, then in the outer regions of the galaxy, our model would overestimate the dark matter mass fraction, consequently underestimating the dark matter density slope.Nevertheless, our article primarily focuses on central galaxies in galaxy groups, which are generally massive elliptical galaxies.According to Li et al. (2018a) and Lu et al. (2023a), these galaxies usually have a flat / gradient.
We finally explore whether the observational data can also be explained by using a steeper dark matter profile.If we fix the stellar mass normalization to the NSA value, the data requires a steep inner density profile for a gNFW dark matter model, with  gnfw = 1.82 +0.15  −0.25 , and a mass-weighted density slope  dm (<  e ) = 1.83 +0.13  −0.22 for the group bin, higher than that predicted by the NFW model.The best-fit massweighted total density slope does not depend strongly on the choice of stellar component normalization, which is  tot (<  e ) = 2.12 +0.05 −0.09 and  tot (<  e ) = 1.93 +0.06 −0.10 for the two bins respectively.The total density slope is steeper than the values from TNG300 by about 2 level.
The next generation weak lensing survey, such as China Space Station Telescope (CSST) (Zhan 2011(Zhan , 2021) ) mission, and the Eu-clid (Laureĳs et al. 2011) mission will all provide a weak lensing source sample at least one magnitude larger than the current DECaLS sample.Combining the stellar kinematic data and these incoming weak lensing data may eventually break the degeneracy between stellar mass normalization and dark matter inner density slope, helping us understand the interplay between dark and light at the very center of dark matter haloes.

DATA AVAILABILITY
The catalogues of dynamical properties (Zhu et al. 2023b) and stellar population properties (Lu et al. 2023a) are publicly available on the website of MaNGA DynPop (https://manga-dynpop.github.io).Other data underlying this article will be shared on a reasonable request to the authors.

Figure 1 .
Figure 1.The relation between the NSA stellar mass and the effective radius  e of the selected lens galaxies in different group masses bins.Colored dots show the sample with different masses.The full sample of MaNGA ETGs are plotted using grey dots.Contours show the full sample of MaNGA ETGs (solid lines)and the combined there sub-samples used in this study (dashed lines) at 30%, 60%, and 90% probability levels.

Figure 2 .
Figure2.Distribution of the number of member galaxies in SDSSGC.We show that of the group bin and cluster bin in the left and right panels, respectively.
as  sys = √︂  2   , where   is standard deviation of  dyn (<  0 ) among 8 different mass models of Paper I for the th galaxy.Finally, we have  2 ,dyn =  2 stat +  2 sys .

Figure 3 .
Figure 3.The figure shows the observational lensing data for group mass bin of 10 13 M ⊙ <  200m < 10 14 M ⊙ .The dark blue circles with error bars show the measured ΔΣ (), and the blue squares with error bars show the predicted ΔΣ () by the mass model derived from MaNGA dynamical modeling alone.This 2D dynamical lensing signal is not involved in the model fitting and we just show it here for illustration.We show the best fit for different mass models in different panels, where the labels of the mass models are marked.The solid lines, dotted lines, dashed lines, and dashdot lines show the total signal, and the contribution of the stellar component, the host dark matter halo, and the two-halo term, respectively.

Figure 7 .
Figure 7.The figure shows the posterior distribution of ggl+dyn+gnfw model with different stellar mass model (blue: NSA; orange:JAM).The left panel shows the results of the 10 13 M ⊙ <  200m < 10 14 M ⊙ bin, and the right panel shows that of the  200m > 10 14 M ⊙ bin.The contours show 68% and 95% confidential levels. 200 is in units of M ⊙ .

Figure 8 .
Figure8.The comparison between our best-fit parameter result and the mass-concentration relation inShan et al. (2017).All the concentrations are corrected from 2D to 3D and rescaled to  = 0, assuming the redshift evolution fromKlypin et al. (2016).The blue, red, and orange circles with error bars show the fitting result of 'ggl only', 'ggl+dyn', and 'ggl+dyn+free ml' models.InShan et al. (2017), they binned the lens samples into two redshift bins, low- (0.2 <  < 0.4) and high- (0.4 <  < 0.6) shown with black triangles left and black triangles right.The black solid and dashed line shows the best-fit mass-concentration relation for the low- and high- subsamples.The grey shaded area shows the 1 uncertainty of this relation of high- subsamples.The grey dashed and dotted curves are the simulation predictions by Duffy et al. (2008) and Klypin et al. (2016).The grey symbols denote the lensing-based measurements of concentration and mass by Mandelbaum et al. (2008) with SDSS (median redshift ( ∼ 0.22); Covone et al. (2014) with CFHTLenS ( ∼ 0.36); and Umetsu et al. (2014) ( ∼ 0.35) and Merten et al. (2015) ( ∼ 0.40) with the CLASH cluster sample.The data points of Merten et al. (2015) were binned.

Figure 9 .
Figure 9. 2D distribution of the inner density slope at the scale of  e and dark matter halo mass  200 (in units of  ⊙ ).The results of 10 13 M ⊙ <  200m < 10 14 M ⊙ and  200m > 10 14 M ⊙ subsample bins are shown in the upper and lower panels, respectively.The dark matter inner density slope and total mass inner density slope at the scale of  e are shown in the left and right panels.In all panels, the blue lines show the fitting results using the NSA stellar mass in the model fitting and the orange lines show results using the JAM stellar mass.The red star in each panel is the corresponding mean value estimated from the TNG300 simulation.The black stars in the left panels are mass-weighted inner density slope of dark matter when we assume that an NFW profile has the best-fit dark halo mass of our sub-samples and follows the mass-concentration relation fromDuffy et al. (2008).Effective radius of the NFW profile are set equal to the average effective radius of the observed sub-samples.Pink stars show the stacked mass-weighted inner density slope of JAM model from Paper I.
2020SKA0110100), the science research grants from the China Manned Space Project (Nos.CMS-CSST-2021-B01, CMS-CSST-2021-A01), CAS Project for Young Scientists in Basic Research (No. YSBR-062), and the support from K.C.Wong Education Foundation.HYS acknowledges the support from NSFC of China under grant 11973070, Key Research Program of Frontier Sciences, CAS, Grant No. ZDBS-LY-7013 and Program of Shanghai Academic/Technology Research Leader.We acknowledge the support from the science research grants from the China Manned Space Project with NO.CMS-CSST-2021-A01, CMS-CSST-2021-A04.WWX acknowledges support from the National Science Foundation of China (11721303, 11890693, 12203063) and the National Key R&D Program of China (2016YFA0400703).JY acknowledges the support from NSFC Grant No. 12203084, the China Postdoctoral Science Foundation Grant No. 2021T140451, and the Shanghai Post-doctoral Excellence Program Grant No. 2021419.

Table 2 .
Posterior constraints of ggl+dyn+gnfw models.The first three columns, show the halo mass range, the number of lenses, and the stellar mass model.The following columns show the fitting parameters, halo mass, inner density slope, mass weighted dark matter density slope, mass weighted total density slope, and dark matter fraction within the mean value of  e .