title: Dark matter halos of massive elliptical galaxies at $z \sim 0.2$ are well described by the Navarro-Frenk-White profile
authors: Anowar J. Shajib Tommaso Treu Simon Birrer Alessandro Sonnenfeld
text: 
  
  firstpage–lastpage
  Large-Range Frequency Tuning of a Narrow-Linewidth Quantum Emitter
      Richard J. Warburton
      December 30, 2023
  ==================================================================
  
  
  
  
  
  
  We investigate the internal structure of elliptical galaxies at z∼ 0.2 from a joint lensing–dynamics analysis. We model Hubble Space Telescope images of a sample of 23 galaxy–galaxy lenses selected from the Sloan Lens ACS (SLACS) survey. Whereas the original SLACS analysis estimated the logarithmic slopes by combining the kinematics with the imaging data, we estimate the logarithmic slopes only from the imaging data. We find that the distribution of the lensing-only logarithmic slopes has a median 2.08±0.03 and intrinsic scatter 0.13 ± 0.02, consistent with the original SLACS analysis. We combine the lensing constraints with the stellar kinematics and weak lensing measurements, and constrain the amount of adiabatic contraction in the dark matter (DM) halos. We find that the DM halos are well described by a standard Navarro–Frenk–White halo with no contraction on average for both of a constant stellar mass-to-light ratio (M/L) model and a stellar M/L gradient model. For the M/L gradient model, we find that most galaxies are consistent with no M/L gradient. Comparison of our inferred stellar masses with those obtained from the stellar population synthesis method supports a heavy initial mass function (IMF) such as the Salpeter IMF. We discuss our results in the context of previous observations and simulations, and argue that our result is consistent with a scenario in which active galactic nucleus feedback counteracts the baryonic-cooling-driven contraction in the DM halos.
  
   
  
  
  
  
  
  gravitational lensing: strong – galaxies: elliptical and lenticular, cD
  
  
  
  
  
  
  
  
  § INTRODUCTION
   
  
  
  
  
  
  Measurements of structural properties of elliptical galaxies can test prediction of galaxy formation theories in the cold dark matter (CDM) paradigm <cit.>. In  this paradigm, small haloes merge to hierarchically form larger halos. CDM-only N-body simulations predict that the dark matter is universally distributed according to the Navarro–Frenk–White (NFW) profile with a `cuspy' central slope, i.e., the 3D density scales as ρ∝ r^-1 in the inner region <cit.>. Such slopes have been observed in galaxy clusters <cit.>. However, shallower central density slopes also have been observed in some galaxy clusters and in dwarf and low-surface-brightness galaxies <cit.>. One possible explanation of these shallower slopes is given by alternative dark matter models, e.g., self-interacting dark matter and warm dark matter <cit.>. Another possible explanation, instead, is related to the astrophysics of galaxy formation and evolution. In fact, as massive elliptical galaxies are believed to be the end-product of the hierarchical merging and accretion processes, their mass density profiles are a sensitive probe of the physics of galaxy formation and evolution.
  
  
  In the galaxy formation process, baryons play an important role that can affect the central density slope of the dark matter distribution. If gas cools and inflows slowly towards the centre, then the dark matter distribution can be adiabatically contracted <cit.>. In contrast, the dark matter distribution can expand in dissipationless mergers or due to gas outflows driven by supernova feedback, stellar feedback, or active galactic nucleus (AGN) feedback <cit.>. Observational avenues to study these processes and their importance in galaxy formation have been limited. The combination of lensing and dynamics has been one of the most informative probes of the dark and luminous matter distributions in the inner region of galaxies and clusters <cit.>. The elliptical galaxies in the Sloan Lens ACS (SLACS) survey were found to have no contraction on average <cit.>, although individual galaxies can have contracted (or expanded) dark matter distributions <cit.>. As many or all of the above-mentioned baryonic processes happen at various points of the galaxy formation process, the degree of contraction or expansion depends on the relative importance of these baryonic processes.
  
  
  The interplay between baryon and dark matter in galaxy formation is also highlighted by the so-called `bulge-halo conspiracy' <cit.>. This conspiracy refers to the nearly isothermal profiles of the total matter distribution with small scatter (∼0.1–0.2) observed within half of the half-light radii to 100 half-light radii of elliptical galaxies (e.g., from strong and weak lensing: ; from stellar dynamics: ). 
  As neither of the baryonic and dark matter distributions follows a power law, fine-tuning between these two distributions is required to produce the isothermal distribution for the total mass. To understand the origin of this conspiracy through simulation, we can use two parameters: the dark matter fraction f_ dm within the inner region and the distribution of logarithmic slope γ for the total mass profile. Near the half-light radius, dark matter has as a logarithmic slope of γ≤ 1.5 and the baryonic distribution has a logarithmic slope of γ∼ 2.3 for a de Vaucouleurs profile <cit.>. Thus, fine-tuning in f_ dm is required to achieve γ∼ 2 for the combination of baryonic and dark matter. However, cosmological hydrodynamic simulations have been unable to match both the observed f_ dm and slope distribution. The γ distribution can be reproduced in simulations by having no or weak feedback, but this leads to an overestimated galaxy formation efficiency and underestimated f_ dm <cit.>. Conversely, reproducing the observed f_ dm requires strong feedback, but then the predicted γ distribution is too shallow. <cit.> similarly find that simulation with AGN feedback can predict the observed f_ dm in 0.4–8 × 10^13 M_⊙ halos at z=0, but underestimate the γ distribution. Without the AGN feedback, overestimated galaxy formation efficiency leads to underestimated f_ dm and overestimated γ distribution. More recently, <cit.> studied simulated elliptical galaxies from the Illustris hydrodynamic simulation that incorporates a number of baryonic processes <cit.>. These authors find higher f_ dm and lower average γ in Illustris galaxies than those observed in lens elliptical galaxies <cit.>. 
  These mismatches 
  point to either inadequacy in the theoretical model or systematic biases in the observational methods.
  
  Many of the observational constraints on f_ dm and γ come from strong lensing systems with elliptical galaxies as deflector. Strong gravitational lensing provides a robust probe of the total projected mass within the Einstein radius. Thus, combining the stellar kinematics with the lensing information can constrain the mass distribution in the deflector galaxy <cit.>. To constrain f_ dm, decoupling the baryonic and dark components in the total mass is necessary. In previous studies, the stellar mass was inferred from the spectral energy distribution to decouple the dark and baryonic components <cit.>. This stellar mass depends on the assumption of the stellar initial mass function (IMF) introducing an uncertainty by a factor of ∼3. However, the IMF can be constrained by making assumption on the mass profile <cit.>, or by decoupling the baryonic and dark components with other external constraints <cit.>.
    
  
  
  
  
  
  
  
  
  
  
  
  
  
  In this paper, we aim to constrain the γ distribution and f_ dm of elliptical galaxies from the lensing and kinematics data – independent of the SED-based stellar mass measurements – and to constrain the amount of adiabatic contraction (or expansion) in these galaxies. We model a sample of  galaxy–galaxy lenses to study their structural properties. These lenses are assembled from the SLACS survey <cit.>. 
  Previous SLACS analyses measured only the Einstein radius from the imaging data and then constrained the radially averaged logarithmic slope using stellar kinematics in combination with the imaging data.
  
  Due to several improvements in lens modelling techniques in the past decade, we can now constrain the logarithmic slope only from the imaging data, exploiting the richness of pixel-level information in the lensed arcs <cit.>. We model the SLACS lenses in our sample using state-of-the-art lens modelling techniques that simultaneously reconstruct the sources to extract the information contained in the lensed arcs. 
  Thus, we measure the local logarithmic slope at the Einstein radius only from the imaging data, independent of the stellar kinematics. We then combine the stellar kinematics with the lensing constraints to individually constrain the stellar and dark matter distributions, and infer the amount of adiabatic contraction in the dark matter distribution. Lensing-only measurement of the mass distribution is prone to the mass-sheet degeneracy <cit.>. We adopt MSD-invariant quantities as the lensing constraints in our joint lensing–dynamics analysis, and combining these constraints with the stellar kinematics allows us to constrain the MSD in the inferred mass distribution.
  
  
  
  This paper is organized as follows. In Section <ref>, we describe our lens sample and the imaging data. Next in Section <ref>, we describe our uniform modelling procedure for this sample. We report the structural properties of elliptical lens galaxies from the lens models in Section <ref>. We discuss our results in Section <ref> and summarize the paper in Section <ref>. Additionally in Appendix <ref>, we investigate the alignment between mass and light distributions. We adopt a flat Λ cold dark matter model as the fiducial cosmology with H_0=70 km s-1 Mpc-1, Ω_ m=0.3, and Ω_ b = 0.047. 
  The reported uncertainties are obtained from 16th and 84th percentiles of the corresponding posterior probability distributions. We use log x to express the natural logarithm and log_10x to express the common logarithm.
  
  
  
  § LENS SAMPLE
   
  
  Our lens sample consists of  galaxy–galaxy lenses from the SLACS survey. We first selected 50 galaxies from the full SLACS sample of 85 lenses by visually inspecting the lens images and selecting those: (i) without nearby satellite or line-of-sight galaxies, (ii) without highly complex source morphology, (iii) with HST imaging data in the F555W/F606W bands (hereafter V band), and (iv) not disc-like. Criteria (i) and (ii) are adopted so that we can uniformly apply our modelling procedure to the whole sample without needing to tweak the lens model settings on lens-by-lens basis. We adopt criterion (iii), because in the V band the deflector galaxy is relatively fainter in comparison with the lensed arcs than in the F814W band (hereafter I band), which makes it easier to decouple the deflector light from the lensed arcs during lens modelling. We provide the list of selected galaxies in Appendix <ref>.
  
  
  
   §.§ Imaging data
  
  Among the selected galaxies, some have imaging data from Advanced Camera for Surveys (ACS), and the rest from Wide Field and Planetary Camera 2 (WFPC2). The ACS images were taken with the F555W filter and the WFPC2 images were taken with the F606W filter. The images are obtained under the HST GO programs 10494 (PI: Koopmans), 10798 (PI: Bolton), 10886 (PI: Bolton), and 11202 (PI: Koopmans).
  
  The WFPC2 images were reduced for the original SLACS analysis <cit.>. We reduce the ACS images using the standard astrodrizzle software package <cit.>. The final pixel scale after drizzling is 0.05 arcsec.
  
  We obtain the point spread function (PSF) for each filter and camera combination using tinytim <cit.>.
  
  
  
   §.§ Stellar kinematics data
  
  We use the line-of-sight velocity dispersions of the lenses in our sample measured from the Sloan Digital Sky Survey (SDSS) spectra. The fibre radius is 1.5 arcsec and the typical seeing for the observations is 1.4 arcsec. <cit.> first measured the velocity dispersions from the SDSS reduction pipeline. <cit.> improved the measurements by updating the set of templates used to fit the spectra. We use this improved kinematics measurements in this study. These measurements are in good agreement with the Very Large Telescope (VLT) X-Shooter measurements of a subsample presented by <cit.>.
  
  <cit.> find a residual scatter in the joint lensing–dynamics analysis using the same kinematics data used in this study, which accounts for ∼6 per cent unaccounted systematic uncertainty in the measured kinematics. Therefore, we add 6 per cent uncertainty in quadrature to the measured uncertainties.
  
  
  
   §.§ Weak lensing data
  
  
  We incorporate weak lensing shear measurements of a sample of 33 SLACS lenses. The measurement pipeline is described by <cit.>, and the sample size of the analysed SLACS lenses is increased by <cit.>. <cit.> find that the shear measurements of these 33 lenses is consistent with the shear measurement of the Hyper-Suprime Cam (HSC) survey weak lensing measurements after weighting the HSC sample to match the stellar mass distribution of the SLACS lenses.
  
  Out of the 23 SLACS lenses in our sample, 11 have directly measured reduced shear γ̃_ shear≡γ_ shear/(1 - κ). For the remaining 12 lenses, we adopt the mean and scatter (which includes both the intrinsic scatter and the noise) of the measured reduced shears for the 33 lenses as the measured value and uncertainty, respectively. We adopt binned reduced shears only up to ∼100 kpc as the weak lensing constraint in our analysis. We do not use measurements beyond 100 kpc to avoid any potential bias from the 2-halo term as this term is not accounted for in our model. The centres of the adopted four bins are logarithmically spaced at 9.87 kpc,  17.78 kpc,  32.04 kpc, and  57.72 kpc.
  
  
  
  § LENS MODELLING
   
  
  We model the lenses using the lens modelling software lenstronomy, which is publicly available on GitHub[ <https://github.com/sibirrer/lenstronomy>] <cit.>. The robustness of lenstronomy in recovering lens model parameters has been verified through Time-Delay Lens Modelling Challenge <cit.>. Two separate participating teams used lenstronomy to successfully recover the hidden Hubble constant and lens model parameters with statistical consistency for TDLMC Rung 2. We package our modelling code into the dolphin pipeline[ <https://github.com/ajshajib/dolphin>], which is a wrapper for lenstronomy to uniformly model large lens samples. First in Section <ref>, we describe the components in our uniform lens model. Then in  Section <ref>, we describe the optimization procedure for the lens model and Bayesian inference of the model parameters. Next in Section <ref>, we assess the effect of the PSF on the measured logarithmic slopes.
  
  
  
   §.§ Model components
   
  
  We adopt the power-law ellipsoidal mass distribution (PEMD) for the deflector <cit.>. Although we aim to individually constrain the dark matter and stellar distributions from a joint lensing–dynamics analysis, first we only need to constrain MSD-invariant local lensing properties from our lens models, for which the PEMD model is sufficient. We combine these local lensing constraints with the stellar kinematics data in our joint lensing–dynamics analysis to individually constrain the dark matter and stellar distributions in Section <ref>. The convergence for the PEMD is given by
  
  
      κ (x, y) = 3 - γ/2( R_ E/√(q_ m x^2 + y^2 / q_ m))^γ-1,
  
  
  where R_ E is the Einstein radius, q is the axis ratio, and γ is the logarithmic slope for the mass distribution in 3D. For an isothermal profile, the logarithmic slope is γ=2. The on-sky coordinates (x, y) are rotated by position angle PAm from the (RA, DEC) coordinates to align the x–axis with the major axis of the projected mass distribution. We also adopt an external shear profile parametrized with the shear magnitude γ_ ext and the shear angle ϕ_ ext.
  
  We adopt a double Sérsic profile for the deflector's light distribution, as a single Sérsic profile leaves significant residual at the galaxy's centre <cit.>. The Sérsic profile is given by
  
  
      I (x, y) = I_ eexp[ -k {( √(x^2 + y^2/q_ L^2)/R_ eff)^1/n_ s - 1 }],
  	
  
  where R_ eff is the effective radius, I_ e is the surface brightness at R_ eff, q_ L is axis ratio, n_ s is the Sérsic index, and k is a normalizing constant so that R_ eff becomes the half-light radius <cit.>. The position angle for the light distribution is PAL. To constrain the degeneracy between the pairs of R_ eff and n_ s in the double Sérsic profile, we fix n_ s = 1 and n_ s = 4, i.e., the exponential profile and the de Vaucouleurs profile respectively <cit.>. For simplicity, we also join the ellipticity parameters q_ L and PA_ L between the two Sérsic profiles.
  
  We reconstruct the source galaxy's light distribution with a basis of shapelets and a Sérsic profile <cit.>. Similar shapelet-based source reconstructions have been successfully used to model a sample of 13 quadruply lensed quasars <cit.>, and to model time-delay lenses to measure the Hubble constant <cit.>. The order parameter n_ max determines the number of shapelets as N_ shapelets = (n_ max+1)(n_ max+2)/2. The scale size of the shapelets is ruled by the scaling parameter β.
  
  
  
   §.§ Optimization and inference
   
  
  We obtain the posterior probability distributions for our model parameters using the Markov chain Monte Carlo (MCMC) method. If the MCMC sampling is started from a point close to the maxima of the posterior, the chain can converge with relatively less computational time. Therefore, we first optimize the lens model to get a point close to the maxima of the posterior. We use the particle swarm optimization (PSO) method for this step <cit.>. To further make this optimization computationally efficient, we adopt the following optimization recipe:
  
  
  	
    * Join the deflector mass and light centroids and fix the logarithmic slope γ=2 and shear magnitude γ_ ext=0 for all the steps below.
  	
    * Create a mask for the lensed arcs. We provide the algorithm to automatically make the mask for the arcs in Appendix <ref>. Fix all the model parameters except for the deflector's light profile. Optimize the deflector light parameters masking the lensed arcs.
  	
    * To find a good starting point for the source light parameters, fix the lens model parameters and the deflector light parameters. Fix the Einstein radius R_ E and ellipticity parameters {q_ m, PA_ m} to the values measured by the SLACS analysis <cit.>. If such pre-determined values are not available, R_ E can be fixed to an approximate guess and the ellipticity parameters can be set to the values for the circular case. Fix the shapelet scale parameter β= 0.1 arcsec. Optimize only the remaining source light parameters to find an approximate position of the source on the source plane. Note, in this step the lensed arcs are not masked.
  	
    * Keep the deflector light parameters fixed and optimize for the source parameters and the PEMD parameters together. Keep β= 0.1 arcsec fixed in this step.
  	
    * Free β and optimize all the non-fixed parameters together.
  	
    * Repeat steps 2–5 with the current initial conditions from the previous step.
  
  
  We notice prominent residuals at the centre of the deflector after performing the above automated procedure. This is expected because the centres of elliptical galaxies are not perfectly described by Sérsic profiles and our signal-to-noise ratio is very high in the centre. To avoid bias in the model from this poor fitting of the deflector light profile at the centre, we mask out the central 0.4 arcsec and rerun the whole fitting procedure. The deflector light profile parameters for these systems are constrained from the light distribution that falls outside the central masked region. However for four systems – J0252+0039, J1112+0826, J1313+4615, and J1636+4707 – masking the deflector centre leads to even poorer quality fits as evaluated with the p-value of the χ^2 statistic. Therefore, we do not mask the deflector centres for these four systems.
  
  We set n_ max = 6 for most of the lens systems in our sample. However, for the following systems we have adopted the following n_ max values through trial-and-error with the above optimization procedure –  J0252+0039: 10, J0959+0410: 15, J1250+0523: 12, J1313+4615: 10, J1630+4520: 15.
  
  
  After the pre-sampling optimization, we then initiate the MCMC sampling from the  optimized lens model after step 6. During the sampling, we free the parameters γ and γ_ ext that were fixed in the optimization step. We also independently sample the centroids of the deflector mass and light distributions. We perform the MCMC sampling using emcee, which is an affine-invariant ensemble sampler <cit.>. We assure the convergence of the chain by checking that the median and standard deviation of the emcee walkers at each step have reached equilibrium.
  
  After our uniform modelling procedure for the 50 initially selected lenses, we vet for reliability of the lens models based on the following criteria: (i) absence of prominent model residuals indicative of poor source reconstruction, and (ii) the median of the inferred logarithmic slope has not converged to too low (γ≲1.4) or too high (γ≳ 2.8) values. The γ distribution can converge towards such extreme values if the lensed arcs are faint and thus the lensing information contained in the imaging data is not sufficient to constrain γ. When γ is too small (γ≲ 1.4), a central image will be produced in the lens models, which is not observed in the imaging data by looking at the color distribution. However, since we mask out the central region for most of our lens systems, a central image is not directly penalized in the likelihood term. Moreover, since the lensed arcs are relatively fainter, prominent residuals are not noticeable in the lens models even if the source reconstruction is poor. We treat such extreme values of γ as numerical artefacts due to weak constraining power of the imaging data and remove these systems from our sample. Note that criterion (ii) is practically a uniform prior median (γ) ∼𝒰(1.4, 2.8). 
  After this vetting procedure, we are left with  systems with reliable lens models. This subsample of 23 lenses is representative of the full sample of 85 SLACS lenses in terms of the stellar mass and velocity dispersion distributions, however our selection excludes the lens galaxies with relatively larger effective radii (R_ eff, V≳ 5 kpc; Figure <ref>). We perform the 1D Kolmogorov–Smirnov test for each of the log_10(M_⋆/M_⊙), σ_ los, and R_ eff, V distributions between the full SLACS sample and our subsample. The p-values are 0.37, 0.84, and 0.49, respectively. Therefore, the null hypothesis that our subsample is representative of the full SLACS sample cannot be rejected with high significance. For reference, a p-value of 0.05 would allow us to rule out the null hypothesis with 95 per cent confidence level. We show the lens images and the models in Figures <ref>, <ref>, and <ref>. We tabulate the marginalized posteriors of the lens model parameters in Table <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Effect of PSF on measured logarithmic slope
   
  
  We check the effect of the PSF choice on the measured logarithmic slopes. We adopt the following five PSF choices for this test:
  
  
  	
    * tinytim PSF with G2V star SED as our baseline PSF,
  	
    * tinytim PSF with Sc galaxy SED at redshift ⟨ z_ s⟩ =0.62, which is the average source redshift in our sample,
  	
    * tinytim PSF with elliptical galaxy SED at redshift ⟨ z_ d⟩=0.20, which is the average deflector redshift in our sample,
  	
    * a PSF created from re-centreing and stacking star cut-outs from the corresponding HST image, and
  	
    * effective PSF (ePSF) created from the stars from the corresponding HST image <cit.>.
  
  
  
  For PSF choices (iv) and (v), reliable PSFs could be extracted only for the ACS images. Therefore, we only check with PSF choice (iv) and (v) for nine systems with ACS images from our sample. We re-optimize our lens models for PSF types (ii) to (v). We compare the measured logarithmic slopes with these PSF choices in Figure <ref>. The mean deviation |Δ| ∼ 0.001–0.009 of the logarithmic slope distribution is negligible relative to the uncertainty of the individual logarithmic slope (typically σ_γ≳ 0.04) when the SED of the PSF is varied. 
  
  However, a scatter of approximately 0.02–0.03 is introduced in the logarithmic slope distribution from the choice of PSF.
  
  
  
  
  
  
  § STRUCTURAL PROPERTIES OF LENS GALAXIES
   
  In this section, we report our findings on the structural properties of the lens galaxies. 
  In Section <ref>, we present the distribution of the logarithmic slopes constrained from the lensing-only data.  
  Next in Section <ref>, we combine stellar kinematics with the lensing observables to infer the amount of contraction in our sample.
  
  
  
  
   §.§ The distribution of logarithmic slopes
   
  
  
  We plot the distribution of the SLACS lenses in our sample on the γ–R_ E/R_ eff plane in Figure <ref> left-hand panel. For comparison, we also plot the distribution of the 13 strongly lensed quasar systems from <cit.>. Strong lensing constrains the slope of the projected mass profile at the Einstein radius. Therefore, if we assume that elliptical galaxies are self-similar, then this plot illustrates the distribution of slopes of the projected mass density at different scale sizes. The median of the distribution is 2.09^+0.03_-0.04 and the intrinsic scatter of the distribution is 0.13 ± 0.02. We account for the uncertainty in individually measured γ by sampling 1000 sets of γ's assuming Gaussian uncertainty for each measurement and then obtaining the distribution of the medians and scatters from these 1000 sets.
  
  
  
  From this γ–R_ E/R_ eff distribution, we infer the dark matter fraction f_ dm≡ M_ dm / M_ total within R_ eff/2 using the approximate Bayesian computation (ABC) method <cit.>. We take R_ eff in I-band for the lenses in our sample from <cit.>. We model the mass distribution using the Navarro–Frenk–White (NFW) profile for the dark matter and the Hernquist profile for the luminous matter distribution <cit.>. The spherical NFW density profile is given by
  
  
      ρ_ NFW(r)  ≡ρ^ NFW_ s/(r/r^ NFW_ s)(1 + r/r^ NFW_ s)^2,
  
  
  where ρ_ s^ NFW is the density normalization and r_ s^ NFW is the scale radius. The Hernquist density profile is given by
  
  
      ρ_ Hern (r) ≡ρ^ Hern_ s/(r/r^ Hern_ s)(1 + r/r^ Hern_ s)^3,
  R_
  
  where ρ^ Hern_ s is the density normalization and r_ s^ Hern is the scale radius. The Hernquist scale radius in angular unit is R_ s^ Hern = R_ eff/1.8153. The NFW scale radius R_ s^ NFW is obtained from the empirical relation R_200≈ R_ eff / 0.015 given by <cit.> and a theoretical M–c relation, where c is the concentration c ≡ R_200 / R_ s^ NFW <cit.>. We fold in 0.2 dex scatter in the R_200–R_ eff relation and 0.16 dex scatter in the M–c relation. The redshift for the M–c relation are set at the SLACS lens redshifts. This model has two free parameters: f_ dm and its scatter σ_f_ dm. The dark matter fraction f_ dm sets the amplitude ratio between the NFW and Hernquist profiles. The posterior of p(f_ dm, σ_f_ dm| D) for the observables D is given by the Bayes' theorem as
  
  
      p(f_ dm, σ_f_ dm| D) ∝ p(D | f_ dm, σ_f_ dm)  p(f_ dm, σ_f_ dm).
  
  
  Here p(D | f_ dm, σ_f_ dm) is the likelihood and p(f_ dm, σ_f_ dm) is the prior. We take uniform priors f_ dm∼𝒰(0, 0.6) and σ_f_ dm∼𝒰(0, 0.4). 
  
  We perform a likelihood-free inference of p(f_ dm, σ_f_ dm| D) using the ABC method from the distribution of logarithmic slopes γ as observables. At first, we sample f_ dm and σ_f_ dm from the priors. For each sampled pair of f_ dm and scatter, we simulate 22 NFW+Hernquist convergence profiles. We get the local slope of these profiles by computing the quantity [1 + dlogκ(x) / dlog x]. The term [dlogκ(x) / dlog x] is the logarithmic slope of the convergence profile κ (x). We add 1 to this term to make it comparable to our parametrization γ for the logarithmic slope of the 3D density profile. We randomly sample one R_ E/R_ eff value from the Gaussian distribution corresponding to the uncertainty for each of the SLACS lens' R_ E/R_ eff values. In this way, the uncertainty in the measured R_ E and R_ eff is accounted for in our inference. Then, we obtain [1 + dlogκ(x) / dlog x] at these x=R_ E/R_ eff from the 22 simulated NFW+Hernquist profiles. At this point, we have a simulated distribution of logarithmic slopes γ_ sim corresponding to each sampled pair of f_ dm and σ_f_ dm. We adopt the Kullback–Leibler (KL) divergence as the distance between the observed γ distribution and the simulated γ_ sim distribution <cit.>. For computational simplicity, we assume that the scatter in γ is Gaussian distributed. The KL divergence between two Gaussian distributions p(x) = 𝒩(μ_1, σ_1) and q(x) = 𝒩(μ_2, σ_2) is given by
  
  
      D_ KL( p | q )    ≡∫_-∞^∞ p(x) log( p(x)/q(x)) dx 
         = log( σ_2/σ_1) + σ_1^2 + (μ_1 - μ_2)^2/σ_2^2 - 1/2.
  
  
  We take the γ_ sim distribution as p(x) and the observed γ distribution as q(x). We simulate for 2 × 10^5 pairs of f_ dm and scatter. From these simulations, we accept the ones with D_ KL < 0.001. The number of accepted simulations is ∼100. The distribution of f_ dm and scatter for these accepted simulations provides an approximation of their posterior p(f_ dm, σ_f_ dm| D) according to the ABC method (Figure <ref> right-hand panel). The 1D marginalized distributions give f_  dm = 0.23_-0.02^+0.03 and scatter =0.19_-0.04^+0.05 dex. In Figure <ref> left-hand panel, we show the 1σ and 2σ confidence intervals of [1 + dlogκ(x) / dlog x] corresponding to the median values of f_ dm and scatter. Inside the effective radius, the slope distribution of the Hernquist+NFW profile increases following the Hernquist profile as the stellar mass is dominant in this regime. However, the mass distribution transitions to be dark-matter dominated outside the effective radius, and thus the slope distribution shows a decreasing trend as the NFW profile is still at its inner r^-1 regime. For θ≫ R_ eff, the slope distribution would increase to have γ→ 3.
  
  
  
  
  We find no correlation between observed γ distribution and ratio f_ SIE≡σ_ los / σ_ SIE between the observed velocity dispersion and the SIE velocity dispersion (Figure <ref>). We calculate the SIE velocity dispersion is given by
  
  
      σ_ SIE = c √(R_ E D_ s/4  D_ ds),
  
  
  where c is the speed of light. If the underlying true mass profile follows a power law, then γ should positively correlate with f_ SIE [see equation (2.3) of <cit.>, cf. Figure 4 of <cit.>]. Therefore, there are two possible explanations for this lack of correlation: (i) the underlying true mass profile deviates from a power-law, and/or (ii) the noise is large enough that any potential correlation is washed away. We perform a linear regression for the γ–f_ SIE distribution with intrinsic scatter as a free parameter. We find that only 3.5 per cent or less intrinsic scatter is necessary with 95 per cent confidence to fit the data. This result indicates that the uncertainties in γ and f_ SIE are too large to detect the expected correlation for a power-law profile. This is a case of absence of evidence, not evidence of absence.
  
  
  
  
  
  
  
  
  
   §.§ Dark matter contraction
   
  
  When gas inside a dark matter halo cools and condensates at the centre as stars begin to form, the dark matter distribution also contracts in response. If the gas infall is slow and smooth, then the dark matter contraction is adiabatic <cit.>. Inversely, dark matter halo can also adiabatically expand, if there is gas outflow.
  
  We infer the contraction in the dark matter distribution in our lens galaxies by combining the lensing observables from our models with their measured stellar kinematics. For the dark matter distribution, we allow adiabatic contraction or expansion from the `pristine' NFW distribution. We follow the contraction formalism of <cit.>, which is originally based on the formalism of <cit.>. In this formalism, the initial radius r_ i of a dark matter particle and its final radius r_ f after contraction is related as
  
  
      r_ i M_ i (r_ i) = r_ f M_ f (r_ f),
  
  
  where M(r) is the total 3D mass within radius r. 
  Then, we have 
  
  
      M_ f (r_ f)    = M_ gal, f (r_ f) + M_ dm, f (r_ f) 
         = M_ gal, f (r_ f) + (1 -f_ gal) M_ i (r_ i),
  
  
  where f_ gal≡ M_ gal,i/M_ i = M_ gal,i/(M_ dm,i + M_ gal,i) is the fraction of the total mass that cools to form stars. We iteratively solve equations (<ref>) and (<ref>) to obtain the contraction factor Γ (r_ f) ≡ r_ f / r_ i. Following <cit.>, we furthermore adopt the halo response parameter ν to modify the relation between r_ i and r_ f as
  
  
      r_ i≡Γ^-ν (r_ f) r_ f.
  
  
  By varying ν in our models we can adjust the amount of contraction. For example, ν=0 corresponds to no contraction, ν=1 corresponds to full contraction as in <cit.>, and ν<0 corresponds to the case of adiabatic expansion. Cosmological hydrodynamical simulations have found less contraction than the model of <cit.>, e.g., the result of <cit.> corresponds to ν = 0.8 and the result of <cit.> corresponds to ν∼ 0.4 <cit.>. However, these simulations almost always find contracted dark matter halos to some extent at the masses of interest <cit.>.
  
  We adopt a composite mass distribution consisting of the contracted dark matter distribution and the stellar distribution to model the stellar kinematics of the lens galaxies in our sample. For computational simplicity, we use spherical cases of these profiles. This assumption is sufficient for the quality of the measured velocity dispersions <cit.>. The line-of-sight velocity dispersion can be expressed by solving the spherical Jeans equation as
  
  
      σ_ los^2 (R) = 2G/I(R)∫_R^∞𝒦_β(r/R) l(r) M(r)/r dr,
  
  
  where G is the gravitational constant, I(R) is surface brightness distribution, M(R) is the 3D enclosed mass, l(r) is the 3D luminosity density, and 𝒦_β is a function that depends on the parametrization of the anisotropy parameter β <cit.>. We adopt the Osipkov–Merritt anisotropy profile given by
  
  
      β(r) ≡ 1 - σ_ t^2/σ_ r^2 = r^2/r^2 + a_ ani^2 R_ eff^2	,
  
  
  where σ_ r^2 is the radial velocity dispersion, σ_ t^2 is the tangential velocity dispersion, and a_ ani is the anisotropy scaling factor <cit.>. For this parametrization, the function 𝒦_β has the form
  
  
      𝒦_β (u)    = u_ ani^2 + 1/2/(u_ ani+1)^3/2( u^2+u_ ani^2/u) tan^-1( √(u^2-1/u^2_ ani+1))  
                      - 1/2/u_ ani^2 + 1√(1 - 1/u^2),
  
  
  where u_ ani≡ a_ ani R_ eff/R <cit.>. The observed SDSS velocity dispersion is averaged over the aperture of 1.5 arcsec radius as
  
      σ_ ap^2 = ∫_ ap[I(R)σ_ los^2(R) ∘𝒮]  RdR dθ/∫_ ap[ I(R) ∘𝒮]  RdR dθ,
  
  where ∘𝒮 denotes convolution with the seeing.
  
  
  We adopt the NFW profile for the initial distribution of the dark matter, which is given by
  
  
      ρ_ NFW(r)  ≡ρ^ NFW_ s/(r/r^ NFW_ s)(1 + r/r^ NFW_ s)^2,
  
  
  where ρ_ s^ NFW is the density normalization and r_ s^ NFW is the scale radius. For the stellar mass distribution we first adopt a constant mass-to-light ratio (M/L) in Section <ref>, and then an M/L gradient in Section <ref>. We adopt a double Sérsic profile for the stellar light distribution similar to our lens models. To ensure robustness of the light profile fits, we fit the double Sérsic profile from 20×20 arcsec^2 cut-outs around the lens galaxies, which sufficiently contain the galaxy light that is above the background level. Before fitting the light profiles, we subtract the lensed arcs using our best fit lens models. We only obtain the best fit double Sérsic profile and adopt an uncertainty on the fit equivalent to 2 per cent uncertainty for the effective radius. This 2 per cent uncertainty is a conservative estimate, as the 2 per cent corresponds to 95 percentile of the uncertainty distribution of the effective radii in our lens models. We use the concentric Gaussian decomposition method from <cit.> to deproject the 2D stellar light or mass distribution into the corresponding 3D distribution. We take 30 Gaussian components with their standard deviations logarithmically spaced between 0.001 and 30 arcsec. This decomposition approximates the 2D stellar mass profile shapes within ∼0.2 per cent accuracy between 0.005 and 6 arcsec.
  
  
  
  
    §.§.§ Constant M/L for stellar mass distribution
   
  
  We first constrain the model parameters w ≡{ν, a_ ani, f_ gal, r_ s^ NFW, ρ_ s^ NFW} for individual lens systems. Here, r_ s^ NFW is the scale radius of the NFW profile and ρ_ s^ NFW is the normalization. 
  We impose a theoretical prior on the M_200–c_200 relation corresponding to ⟨ z_ d⟩ = 0.20 from <cit.> for the initial NFW halo. To incorporate the lensing constraints into this joint lensing–dynamics analysis, we fold in the posterior distributions of Einstein radius R_ E and the quantity R_ Eα^''_ E / (1 - κ_ E) from our lens models in Section <ref>. Here, α^''_ E is the double derivative of the deflection angle at R_ E and the κ_ E is the convergence at R_ E. The term R_ Eα^''_ E / (1 - κ_ E) is the mass-model-independent observable quantity from imaging data <cit.>. We also use the weak lensing measurements of the reduced shear, although weak lensing do not provide tight constraints due to large uncertainty in the measured shear of individual systems. We take uniform priors on ν∼𝒰(-0.5, 1), a_ ani∼𝒰(0.1, 5), and f_ gal∼𝒰(0, 0.25). Additionally we impose a prior on log_10 M_200 from <cit.> that depends on the measured stellar mass from the stellar population synthesis method assuming Chabrier IMF. We obtain the Chabrier IMF based stellar masses for the galaxies in our sample from <cit.>. The form of the prior is
  
  
      p(log_10 M_200|log_10 M_⋆^ Chab) 
             ≡𝒩(μ_ h + β_ h[log_10 M_⋆^ Chab - 11.3], σ_ h^2 ).
  
  
  The parameters in this prior are obtained by <cit.> by applying the SLACS selection function on SDSS galaxies that have HSC weak lensing data and assuming a model with M/L gradient for the stellar mass and an NFW profile for the dark matter. We take the uncertainty on these parameters to match with the dispersion observed in the SLACS galaxies by <cit.> for the same mass model to have μ_ h = 13.03 ± 0.25, β_ h = 1.42 ± 0.92, and σ_ h = 0.18 ± 0.16.
  
  
  
  
  We then infer the distribution of population-level parameters τ for our lens sample using hierarchical Bayesian inference. According to the Bayes' theorem, we can express the posterior of the population-level parameters as
  
  
      p(τ| D)    ∝ p(D |τ)  p(τ) 
         ∝ p(τ)  ∏_i^N∫dw_i p(D_i | w_i)  p(w_i |τ) 
         ∝ p(τ)  ∏_i^N∫dw_i p(w_i | D_i)/p(w_i) p(w_i |τ),
  
  
  Here, D ≡{ D_i} refers to the data set containing all the data D_i for individual lenses, p(D |τ) is the likelihood term on the population-level parameters τ, p(τ) is the prior, p(w_i |τ) governs the distribution of individual lens parameters given the population-level parameters, p(w_i | D_i) is the posterior of w_i for an individual system given its associated data D_i, and p(w) is the prior for the parameter set w as given above. The likelihood function p(D_i | w_i) for the ith lens system has the form
  
      p(D_i | w_i) ≡   𝒩(R_ E, i, σ_R_ E, i^2) ×𝒩(R_ Eα^''_ E / (1 - κ_ E), σ^2_[R_ Eα^''_ E / (1 - κ_ E)]) 
         ×𝒩(σ_ ap, σ_σ_ ap^2) ×∏_j 𝒩(γ̃_ shear, ij, σ^2_γ̃_ shear, ij),
  
  where σ_ ap is the observed velocity dispersion, and γ̃_ t, ij is the reduced weak-lensing shear in the jth radial bin. We assume that a_ ani, ν, and f_ gal are Gaussian-distributed for our lens sample with the Gaussians truncated where p(w) = 0. Thus, we take τ≡{μ_ν, σ_ν, μ_a_ ani, σ_a_ ani, μ_f_ gal, σ_f_ gal} as the population-level parameters.  Here, μ and σ refer to the mean and standard deviation of the population-level Gaussian distribution for the corresponding parameter. We take uniform priors μ_ν∼𝒰(-1, 1), μ_a_ ani∼𝒰(0.5, 5), and μ_f_ gal∼𝒰(0, 0.25). We take the same uniform prior 𝒰(-5, 1) for log_10σ_ν, log_10σ_a_ ani, and log_10σ_f_ gal. We perform the integral in equation (<ref>) using Monte Carlo integration. 
  In the Monte Carlo integration, we sample from p(w_i |τ) and sum over p(w_i | D_i) / p(w_i). We approximate p(w_i | D_i) with a Gaussian mixture model (GMM) of the corresponding samples drawn using MCMC. Before fitting the GMM, we smooth the distributions of the MCMC samples using the kernel density estimation (KDE) method with Silverman's rule for the bandwidth <cit.>. The number of components in the GMM for each individual lens posterior p(w_i | D_i) is selected using Bayesian information criterion (BIC) and it is typically within 5–10. We check the robustness of our hierarchical inference framework with simulated distributions of p(w_i | D_i) and find that the correct population-level parameter posterior p(τ| D) is recovered.
  
  
  We illustrate the posterior distributions of the population-level parameters in Figure <ref>, and tabulate their point estimates in Table <ref>.  In Figure <ref>, we highlight some of the important results from Figure <ref> along with their interpretations for galaxy properties. We infer the mean halo response parameter μ_ν = -0.06±0.04 and intrinsic scatter σ_ν≤ 0.092 (95 per cent upper limit), which is consistent within 1.5σ with no contraction on average for our lens sample.
  
  We compute the dark matter fraction f_ dm within R_ eff/2 directly from the decomposed dark and stellar mass distributions. We show the γ–f_ dm distribution in Figure <ref>. We find moderate anticorrelation between γ and f_ dm with the biweight mid-correlation r=-0.39 ± 0.15. The anticorrelation is expected as a shallower total mass profile slope will require higher contribution from the dark matter to lower the total mass profile's slope. The distribution of f_ dm in our sample has a mean of ⟨ f_ dm⟩ = 0.27±0.02 and a scatter of 0.14 ± 0.01.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    §.§.§ M/L gradient for stellar mass distribution
   
  
  Next, we incorporate a M/L gradient for the stellar mass distribution. We parameterize a power-law M/L gradient as 
  
  
      Υ (R) = Υ_ e( R/R_ eff)^-η,
  
  
  where Υ_ e is the M/L at R_ eff. However, this parametrization diverges for η > 0 as R → 0, the total stellar mass computed using the concentric Gaussian decomposition method stays finite. Moreover, the observed lensing properties and the kinematic properties are mostly sensitive to the mass distribution near the Einstein radius and the effective radius, respectively. As a result, the deviation from the exact relation near the centre (≲0.01 arcsec) in our Gaussian-decomposed approximation does not noticeably impact our analysis.
  
  We perform the same analysis from Section <ref>, but with the M/L gradient implemented in the stellar mass distribution. As a result, the model parameters for individual lenses are extended to w ≡{ν, a_ ani, f_ gal, r_ s^ NFW, ρ_ s^ NFW, η}. We only allow positive values for η as previous observations suggest that the M/L is higher at the centre than outer region of elliptical galaxies <cit.>. As we want to adopt an uninformative prior on the order of magnitude of η, we take a Jeffrey's prior on η as p(η) ∝ 1/η, which is equivalent to a uniform prior on logη. We set the bounds of the uniform prior on logη as p(logη) ∼𝒰(log 10^-3, 0). Furthermore, the population-level parameters for the Bayesian hierarchical inference are also extended to τ≡{μ_ν, σ_ν, μ_a_ ani, σ_a_ ani, μ_f_ gal, σ_f_ gal, μ_logη, σ_logη}.
  
  
  We infer the 95 per cent upper limit of exp(μ_logη) distribution to be 0.02. The mean halo response parameter for the sample is μ_ν = -0.03^+0.04_-0.05 with intrinsic scatter σ_ν≤ 0.074 (95 per cent upper limit), which is consistent with no contraction in the NFW halo (Figure <ref>). We show the distribution of the inferred f_ dm and the logarithmic slope γ in Figure <ref>. The distribution of f_ dm in our sample for the M/L gradient model has a mean ⟨ f_ dm⟩ = 0.28 ± 0.02 and an intrinsic scatter of 0.15 ± 0.02.
  
  To check for the impact of the M–c relation prior in this analysis, we perform the joint lensing–dynamics analysis without the M–c relation prior. We find that the uncertainties on the model parameters for each lens system expectedly increase without this prior, and the posteriors from with or without the prior are consistent within 1σ (Figure <ref>, Table <ref>). 
  
  
  We also check for the impact of our fiducial cosmology, in particular the Hubble constant, on our inference. We perform the joint lensing–dynamics analysis with H_0 = 67 km s^-1 Mpc^-1 and with H_0 = 74 km s^-1 Mpc^-1. We find all the posterior distributions to be consistent within 1σ. 
  Therefore, our inference is not sensitive to the choice of the Hubble constant.
  
  We show the surface density profiles and the deprojected 3D density profiles of the dark matter and stellar components of the total mass distribution for all the lens galaxies in Figures <ref> and <ref>, respectively. In Figure <ref>, we illustrate the median deviation of the NFW+stars model from the power law in the total surface density profile and the local logarithmic slope. The total density profile combining dark matter and stars deviates higher by ∼5 per cent (median value) at the Einstein radius with a median absolute deviation of ∼7 per cent within our sample. Here, the median absolute deviation is a proxy for the intrinsic scatter. The local slope of the density profile at the Einstein radius is shallower by ∼2 per cent (median value) with median absolute deviation of ∼2 per cent. These deviations are sensitive to our adopted prior on M_200. If we choose a prior with lower mean halo mass such as log_10 M_200∼𝒩 (12.9, 0.5) – which corresponds to the BOSS constant-mass (CMASS) galaxies with mean redshift z≈0.6 <cit.> – the resultant profiles would be within 0.8 per cent of the power laws at the Einstein radius (Figure <ref>).
  
  
  
  
  § DISCUSSION AND COMPARISON WITH PREVIOUS STUDIES
   
  
  In this section, we compare our results with previous studies based on observations and simulations, and interpret them in the context of elliptical galaxy evolution. First in Section <ref>, we discuss the central results of this study on the elliptical galaxy structure.
  
  We then place these results in the context of massive elliptical galaxy evolution in Section <ref>.
  We discuss our results on the M/L-gradient and its implication for the stellar IMF in Section <ref>. We discuss the implication of our results for time-delay cosmography in Section <ref>. We state the limitations of this study in Section <ref>.
  
  
  
   §.§ Elliptical galaxy structure
   
  
  
  The main results of this study on the elliptical galaxy structure are: (i) the dark matter distribution in massive elliptical galaxies at z∼0.2 are not contracted and close to an NFW profile on average (i.e., μ_ν∼ 0), and (ii) the total density profile combining the dark matter and stars is close to isothermal (i.e., ⟨γ⟩∼ 2). 
  
  
  
  
    §.§.§ Dark matter contraction
  
  
  Our result on the dark matter contraction agrees with <cit.>, who also find no contraction in SLACS galaxies using the lens models from <cit.>. Similarly, <cit.> and <cit.> find that the inner slope of the dark matter distribution is consistent with the NFW profile within the uncertainty for their samples of massive elliptical galaxy halos and group-scale halos, respectively. Therefore, these findings are consistent with our result. 
  
  However, our result contradicts the report of steeper central slopes (mean inner logarithmic slope ⟨β⟩ = 2.01^+0.19_-0.22) than that in the NFW profile (β = 1) by <cit.>. 
  Interestingly, <cit.> find their slope distribution to be bimodal, and the steep value given above corresponds to the mode with the larger mean. If <cit.> adopt a unimodal slope distribution, then the average inner slope is consistent with the vanilla NFW profile within 0.5σ for the M/L gradient model and within 1.3σ for the constant M/L model. Albeit, the unimodal distribution requires higher intrinsic scatter. Moreover, <cit.> note that their inferred halo mass is lower than the expectation of abundance matching studies due to the lack of a strong prior to constrain the halo scale size or its mass. Strong-lensing data is sensitive to the mass profile in the central region within the Einstein radius, which is typically ∼5–10 times smaller than the NFW scale radius. Thus, additional data or prior at scales larger than the NFW scale radius is necessary to robustly constrain the halo mass. We find that the degree of contraction depends on the halo mass prior in our analysis with a heavier prior on M_ 200 producing shallower inner slopes to fit the joint lensing–kinematics data. The sample of <cit.> has a similar stellar mass range as our sample and the difference between the mean redshifts of the samples Δ⟨ z ⟩ = 0.15 does not leave enough room to expand the halos from β∼2 to β∼1 within ∼1.44 Gyr. Therefore, we conclude that the differences between our result and that from <cit.> are largely caused by the difference in the adopted priors corresponding to the dark matter halo.
  
  
  
  
    §.§.§ Slope of the total density profile
  
  
  
  Our lensing-only models provide ⟨γ^ lensing⟩ = 2.08 ± 0.04 with a scatter of 0.13±0.02. From the joint lensing–dynamics analysis, we find the total density profile is shallower by approximately 5 per cent, which brings the sample mean of the logarithmic slope at the Einstein radius closer to the isothermal case (Figure <ref>). This near-isothermality of the total density profile agrees well with a multitude of pervious observations – e.g., based on strong-lensing only or jointly based on lensing and dynamics: <cit.>, and based on stellar dynamics: <cit.>.
  
  
  In Figure <ref>, we compare the distribution of the estimated logarithmic slopes γ constrained from the imaging data only in this study with those estimated by the SLACS analysis from combining stellar kinematics with the imaging data. We find no correlation between the estimated γ distributions from the two analyses with biweight mid-correlation r = 0.03 ± 0.17. However, the SLACS distribution has a mean of ⟨γ^ LD⟩ = 2.078 ± 0.027 and intrinsic scatter 0.16 ± 0.02 <cit.>. These values are consistent with our results within 1σ confidence level. For the 21 systems that have measured γ in <cit.>, we find Δ_γ^2 ≡∑_i^N=21(γ^ lensing_i,  this study - γ^ LD_i,  Auger+(10b))^2/σ_i,  total^2 = 21.76. The p-value assuming a χ^2-distribution for Δ_γ^2 with 21 degrees of freedom is 0.41. We can see in Figure <ref> that the two-component mass profile from lensing–dynamics can deviate from the lensing-only inference of the power-law profile towards either direction. However, the sample mean of the such deviations is smaller than ∼5 per cent near the Einstein radius, which explains the good agreement for ⟨γ⟩ between this study and <cit.>.  Thus, a correlation between the lensing-only local slope and the lensing–dynamics global slope is expected. However, such correlation cannot be detected within the noise. 
  
  
  
  We assume Osipkov–Merritt anisotropy profile as this profile follows the observation that the anisotropy profile in elliptical galaxies are isotropic at the centre and gradually becomes radial at the outskirts <cit.>. Our inferred anisotropy scaling factor a_ ani is consistent within 1σ with those reported by <cit.> for both a sample of time-delay lens galaxies and for a joint sample of time-delay lens galaxies and a subset of SLACS galaxies with integral-field unit spectroscopy (cf. table 6 in ). Note, there is significant overlap between our sample and the sample of <cit.>, but the adopted mass models for the kinematics analysis are different. The associated uncertainties reported by <cit.> are generally larger than our inference, which can be attributed to the mass model adopted by <cit.>, that maximally explores the mass-sheet degeneracy.
  
  
  
  
  
  
  In the next subsection, we place these two results in the context of elliptical galaxy formation and evolution.
  
  
  
  
  
  
   §.§ Evolution of massive elliptical galaxies
   
  
  The current paradigm for the formation of elliptical galaxies consists of two phases <cit.>. In the first phase up to z ≈ 2, gas condensates in a massive halo to form stars, and the dark matter distribution contracts as a result. In the second phase after z≈2, elliptical galaxies grow in size primarily through multiple major or minor mergers. Most of this later evolution is “dry” involving little gas and star formation, as evidenced by the almost uniformly old stellar populations in elliptical galaxies. If elliptical galaxy halos contract up to z≈2, our result then raises the question: how does the contracted halo at z ≈ 2 expands back to resemble the NFW profile at z∼0.2? Simulations have shown two principle mechanisms that can remove baryons from the central region and expand the dark matter halo – AGN-driven gas outflow, and dynamical heating through accretion of materials from the environment <cit.>. Although mergers play a crucial rule in the growth and evolution of elliptical galaxies, simulations find that dissipationless mergers do not change the shape of the dark matter profile <cit.>. Dissipational gas-rich mergers, in contrast, can make the profile steeper by transferring baryons towards the centre and thus further contracting the dark matter halo <cit.>. The close-to-isothermal nature of the total mass distribution can also be explained as the end result of rearranging the mass distributions through accretion of collision-less materials in gas-poor mergers <cit.>.
  
  
  To investigate the driving mechanism behind the halo expansion observed in our sample, we compare the structural properties of our galaxies with those from simulations that adopt varying baryonic physics (Figure <ref>). The simulation from <cit.> includes stellar and supernova feedback, and primordial abundance cooling. In contrast, the Magneticum and IllustrisTNG simulations includes AGN feedback in addition to stellar and supernova feedback, and the cooling mechanism includes metals <cit.>. For fair comparison with the simulations, we adopt their definition and compute the average logarithmic slope by fitting a power-law to the 3D density profile between 0.4R_ eff and 4R_ eff. We also compute the central dark matter fraction within a 3D aperture of radius R_ eff. Here, we treat the half-light radius R_ eff as the half-mass radius R_1/2, as the M/L gradient exponent in our sample is very small and consistent with zero. The compared IllustrisTNG simulation corresponds to z=0.2, which is the mean redshift of the galaxies in our sample. The fitting functions from <cit.> corresponding to Magneticum and Oser simulations – that are plotted in Figure <ref> – are shown to match with the simulated distributions over the redshift range z=0–2. In all of the panels of Figure <ref>, the general trend found in the simulations matches that in our observations. However, the distributions themselves do not overlap. The mismatch between our observation and IllustrisTNG simulation predominantly arises from the higher f_ dm^ 3D found in that simulation. However, the γ distribution of the IllustrisTNG simulation is generally consistent with our sample. Similarly, the γ^ lensing distribution obtained from fitting imaging-only data of simulated lenses in the SEAGLE simulation also reproduces the distribution from our observation (Figure <ref>). However, the f_ dm distribution of the SEAGLE simulation is not reported yet. In contrast, the Magneticum simulation produces higher radially averaged γ distribution than our observed one, but reproduces a f_ dm distribution that is generally consistent with our observation (cf. figure 8 of ). Note, the simulations adopted different IMFs – Chabrier IMF in the IllustrisTNG, EAGLE, and Magneticum simulations, and Salpeter IMF in the Oser simulation. Adopting the Salpeter IMF – which is consistent with the galaxies in our sample (Section <ref>) – would lead to lower f_ dm in the corresponding simulations, thus bringing the f_ dm closer to agreement with our observation. However, if the IMF mismatch cannot fully explain the observed differences, changes in the baryonic prescriptions in the simulations need to be fine tuned. The slope of the correlation in our observed f_ dm–γ distribution falls between the ones from the Magneticum and Oser simulations, which is consistent with the fact that the AGN feedback is a necessary driver in the evolution of these galaxies. To match the slope of the Magneticum simulation with our observation, weaker AGN and stellar feedback would be necessary. However, that would also push up the simulated γ distribution and increase the discrepancy with our observation. Note that the lower values of f_ dm than simulations is not unique to our lensing-based analysis, as previous lensing-based studies also find similarly lower f_ dm distributions <cit.>. Moreover, the Magneticum and IllustrisTNG simulations themselves have offsets between their predicted f_ dm distributions due to differences in their implementations of baryonic physics despite the similarities of the adopted feedback mechanisms. 
  The offset between simulations highlights that the effect of AGN feedback depends critically on the specific implementation. Thus our observations cannot test AGN feedback in general, but they are a powerful empirical test of specific implementations.
  
  
  Similarly, the trend in our observed distribution of f_ dm^ 3D–R_ eff matches with the ones from the simulations, however the distributions themselves do not match (Figure <ref>). The moderate correlation (r=0.46±0.07) between the size of the galaxies and the central dark matter fraction supports that these galaxies have grown predominantly through minor mergers. As minor dissipataion-less mergers mostly deposit material at the outer region of the halo, the galaxy's half-mass radius increases as a result. Thus, the dark matter fraction within the half-mass radius also increases and becomes correlated with the half-mass radius. 
  
  
  The velocity dispersion distribution of the IllustrisTNG simulation is smaller than the measured distribution for our sample. <cit.> note that this difference between observation and simulation in velocity dispersion, along with the mismatch in the f_ dm distribution, point to potentially insufficient implementation of baryonic physics in the IllustrisTNG simulation.
  
  
  
    In summary, our constraints on the elliptical galaxy structure in the context of previous observations and simulations are consistent with the following formation and evolution scenario for massive elliptical galaxies. The first stage of the formation of elliptical galaxies is through dissipational processes at z≳ 2, when most of their present day stars are formed. The dissipation leads to contraction in the dark matter halo and the resultant total density profile observed in cosmological numerical simulations is steeper than the isothermal case <cit.>. After z ≈ 2, the growth of the elliptical galaxies is dominated by gas-poor dissipationless mergers, explaining their growth in size without the addition of younger stellar populations <cit.>. This growth mechanisms is consistent with the R_ eff–f_ dm correlation observed in our sample. Multiple dissipationless mergers decrease the total density profile of the galaxies to bring it close to isothermal, and increase the half-mass radius and the central dark matter fraction with decreasing redshift <cit.>. Furthermore, AGN feedback expands back the contracted dark matter halos <cit.>, which is supported by the slope of our observed γ^ LD_[0.4, 4]R_ eff–f_ dm^ 3D distribution. Dynamical heating from accretion may also play a role in expanding the dark matter halos in addition to the AGN feedback, however we do not find any indication either in favor or against the presence of dynamical heating in our sample. Gas-rich dissipational mergers can also happen at some point of the galaxies' growth history. However, the small intrinsic scatter in the halo response parameter (σ_ν≲ 0.1) around the mean μ_ν∼ 0 in our sample, and the ages of the stellar populations <cit.>, indicate that such gas-rich mergers are relatively rare, even though their contribution cannot be ruled out or confirmed conclusively given the present precision of numerical simulations and observations <cit.>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Stellar IMF and M/L gradient
   
  
  
  
  
  Our inferred stellar mass M_⋆^ LD distribution from the joint lensing–dynamics analysis support a heavier stellar IMF such as the Salpeter IMF. In Figure <ref>, we illustrate the distribution of the IMF mismatch parameter as defined by <cit.> as
  
  
      α^ Salp≡M_⋆^ LD/M_⋆^ SPS, Salp,
  
  
  where M_⋆^ SPS, Salp is the SPS-based measurement of the stellar mass assuming the Salpeter IMF. If the true underlying IMF follows the Salpeter function, then α^ Salp = 1 is expected. For the Chabrier IMF, in comparison, the IMF mismatch parameter should be α^ Salp≈ 0.56 or log_10α^ Salp≈ -0.25. The mean IMF mismatch parameter of our sample is ⟨log_10α^ Salp⟩_ median = 0.00 ± 0.03 with an intrinsic scatter of 0.11^+0.03_-0.02. Our observed distribution of α^ Salp–σ_ e is in good agreement with the fitting function provided by <cit.> 
  
  
      log_10α^ Salp = (0.38 ± 0.04) ×log_10( σ_ e/200  km s^-1) + (-0.06 ± 0.01 ),
  
  
  which is obtained from fitting SLACS <cit.> and ATLAS^ 3D <cit.> galaxies. This relation has a root-mean-square (rms) scatter of 0.12 dex.
  
  
  Several previous strong-lensing studies also found evidence for heavier IMF such as the Salpeter <cit.>. Multiple non-lensing studies similarly found evidence for a heavy IMF in elliptical galaxies – e.g., from dynamics of local elliptical galaxies <cit.>, and based on single stellar population (SSP) models <cit.>.
  
  
  Whereas the observed IMF in the Milky Way stellar populations is closer to a Chabrier IMF <cit.>, our observation of a Salpeter-like IMF in massive elliptical galaxies indicates that stars formed in a different environment in these elliptical galaxies. Star formation in a turbulent environment with high gas density – that can result from gar-rich mergers and accretion events at the early phase of elliptical galaxy formation – can lead to a heavy IMF such as the Salpeter IMF <cit.>.
  
   
  
  
  We find that the M/L-gradient exponent is very small and consistent with zero, as the 95 per cent upper limit for the sample mean of the M/L-gradient exponent η is 0.02. This small value is consistent with <cit.>, who found that most of their lens galaxies individually favor η < 0.1. Our result, however, is in tension with <cit.>, who inferred μ_η=0.24 ± 0.04 and σ_η=0.09±0.02 from a larger SLACS subsample and HSC weak-lensing information. Some fraction of the difference can be explained by the difference in the adopted light profile. Whereas <cit.> adopt a de Vaucouleurs profile for I band, we adopt a double Sérsic profile in V-band. The double Sérsic profile fits the light distribution better than the de Vaucouleurs profile. We observe that our double Sérsic profile fits are generally steeper than the de Vaucouleurs fit of the SLACS galaxies from <cit.>, which explains the steeper mass profile than the light profile reported by <cit.>. Such steeper mass profile than the light profile can not arise from a color gradient between I and V bands, as observationally I-band light profile is steeper than the V band as the color gets bluer radially outward <cit.>. Furthermore, <cit.> adopt an isotropic profile for orbital anisotropy in their dynamical model, whereas we adopt an Osipkov–Merritt anisotropy profile, which can potentially account for some part of the observed discrepancy by trading radial anisotropy for M/L gradient to reproduce the observed velocity dispersion.
  
  
  
  Several recent studies based on SPS modelling have pointed to a radially varying IMF in local elliptical galaxies, where the IMF in the central ≲2 kpc region is bottom-heavy – even super-Salpeter – and it gets more bottom-light radially outward to match with the Chabrier IMF in the outer region <cit.>. Such a radially varying IMF is consistent with the scenario where the central stellar population is dominated by in situ stars that formed following the Salpeter IMF in a highly turbulent environment at high redshift (z ≳ 2), and the stellar population in the outer region mostly comes from merged satellites that had a different star-forming environment closer to the Milky Way and thus the Chabrier IMF. However, this radially varying IMF is interpreted from a radially varying M/L ratio based on the SPS models. Thus, at face value our result of very small M/L gradient appears to be in tension with the SPS-based studies mentioned above. However, the M/L gradient observed in the SPS-based studies is dominated by a prominent gradient at <0.1R_ eff, whereas our joint lensing–dynamics analysis is sensitive near 0.5R_ eff–R_ eff where the SPS-based M/L flattens out. Thus, the apparent difference between our result and SPS-based M/L gradient results can be reconciled. Indeed, <cit.> point out that the observed Salpeter IMF in joint lensing–dynamics studies is consistent with the SPS-based result of a radially varying IMF, if the measurement scales are taken into account [cf. Figure 17 of <cit.>].
  
  
  Our inferred star formation efficiency is larger by approximately a factor of 2 than those reported by studies using abundance matching at the relevant redshift and halo mass assuming the Chabrier IMF <cit.>. Adopting the Salpeter IMF would make the star formation efficiency derived from abundance matching consistent with our result.
  
  
  
  
  
  
   §.§ Implication for time-delay cosmography
   
  
  In time-delay cosmography, the measurement of the Hubble constant H_0 requires an accurate constraint on the surface density profile near the Einstein radius. Similar to this study, analyses by the Time-Delay COSMOgraphy (TDCOSMO) have adopted a power-law mass distribution – among other choices – for lens modelling <cit.>. A sample of seven strong-lensing systems with measured time delays provide a ∼ 2 per cent measurement of the Hubble constant <cit.>. However, relaxing the assumption on the power-law profile and thus allowing the full range of mass-sheet degeneracy in the lens models inflates the uncertainty to ∼8 per cent <cit.>. As a result, precise knowledge of the internal structure of elliptical galaxies is currently the limiting factor for the time-delay measurement of H_0. Therefore, one option to shrink back the uncertainty on H_0 is to apply an external prior on the internal structure of time-delay lens galaxies, and thus narrow down the allowed range of the MSD. Alternatively, <cit.> constrain the MSD by incorporating the structural information from the power-law lens models presented in this study and their measured kinematics under the assumption that the TDCOSMO and SLACS lenses are drawn from the same parent population. The additional hypothesis and the addition of external information reduces the uncertainty on H_0 to 5 per cent. However, including the structural information from the SLACS galaxies also increases the mean surface density of the time-delay lens galaxies near the Einstein radius by ∼9.5 per cent from the power-law estimate, which also translates to a ∼9.5 per cent decrease in the inferred H_0. This increased estimate of the surface density is consistent within the uncertainty with our joint lensing–dynamics analysis, as we find on average ∼5 per cent higher surface density near the Einstein radius for ours stars+NFW model than the power-law model (see Figure <ref>).
  
  
  TDCOSMO analyses have also adopted an NFW+stars model in addition to the power-law model and found that the two models are in very good agreement <cit.>. In contrast, we observe that the NFW+stars model deviate from the power law on average by ∼5 per cent at the Einstein radius. We show that the deviation of the NFW+stars model from the power-law model in our analysis is largely driven by the adopted prior on the halo mass. As strong-lensing and kinematics data mostly constrain the mass distribution in the central region, the total halo mass or the dark matter profile is poorly constrained from lensing–dynamics analysis without a prior on the halo mass, or on the mass distribution at scales larger than the NFW scale radius. Our empirical prior on the halo mass M_200 is obtained from the HSC weak-lensing measurements of the SDSS galaxies, which is the parent sample of the SLACS galaxies. Moreover, the selection function for the SLACS galaxies is accounted for in this prior. If instead, we had adopted a prior on M_200 with a smaller mean by ∼0.3 dex, the resultant NFW+stars density profile would have matched the power-law model very well on average (Figure <ref>). A higher M_200 translates to a higher normalization of the NFW profile, thus the total density profile becomes dark-matter-dominated at a relatively smaller radius and the slope of the total density profile gets shallower. We provide two possible explanations for why our NFW+stars model deviates from the power-law whereas the ones in the TDCOSMO analyses do not.
  
  
  The first possible explanation is that there may be a bias in the NFW profile normalization in the TDCOSMO models due to weak prior constraints on the mass distribution at large scales. TDCOSMO analyses have generally adopted a prior only on the NFW scale radius, which only has a weak constraint on M_200 as the M–c relation is not imposed. As a result, the resultant halo mass may be biased low as lens models tend to produce lower halo masses without strong priors on the mass distribution at scales larger than the NFW scale radius <cit.>. A smaller halo mass by ∼0.3 dex may not be noticeable within the measurement uncertainty when looking at individual lens galaxies, but it may be identifiable by investigating the M_⋆–M_ 200 relation for the TDCOSMO lenses. However, due to the typically large intrinsic scatter in M_200, a sample size of seven may not be sufficient to identify such a bias.
  
  
  The second possible explanation is that the SLACS lens galaxies and the TDCOSMO lens galaxies do not belong to the same galaxy population. Even if the predecessors of the SLACS galaxies at z=0.6 represent the time-delay lens galaxy population, they may have evolved since z=0.6 to have a different internal structure at z=0.2 <cit.>. A M_200 prior with larger mean leads to a larger dark matter fraction f_ dm. As f_ dm is observed to increase from higher redshift to lower redshift <cit.>, it is possible that a lower f_ dm in the time-delay lenses with mean deflector redshift ⟨ z_ d⟩=0.6 makes the power-law model and the NFW+stars model to be in good agreement. 
  
  We note that an increasing f_ dm with decreasing redshift in the TDCOSMO lenses would also be consistent at least qualitatively with the weak trend observed in the individual H_0 measurements from these lenses <cit.>. As increasing f_ dm can create larger positive deviation between the NFW+stars model and the power-law model, the resultant H_0 from the power-law model would shift towards higher values with decreasing redshift. However, this scenario requires a combination of both of the above explanations where high-redshift galaxies are less vulnerable to an NFW normalization bias due to their lower underlying f_ dm, and lower redshift galaxies are more vulnerable to the normalization bias due to their higher underlying f_ dm.
  
  
  
  <cit.> suggest that a hypothetical 10 per cent positive bias in the Hubble constant from time-delay cosmography would point to the presence of a large core in the dark matter distribution. We show that such a shift in the Hubble constant can also be explained by a relatively higher normalization in the vanilla NFW profile, if the SLACS lens galaxies and the TDCOSMO lens galaxies are structurally self-similar.
  
  
  
  
  
  
  
  
  
  
  
   §.§ Limitations of this study
   
  
  Although we do not account for external convergence from the line-of-sight structures in our lens models, this simplification has no impact on our results, as we discuss in this paragraph. <cit.> find the mean external convergence for a subsample of 33 SLACS lenses to be ⟨κ_ ext⟩ = 0.005_-0.009^+0.022. However, this subsample is curated to only select systems with small local overdensities [see <cit.> for details]. Our sample of 23 lenses has a large overlap with this sample of 33 lenses from <cit.>, thus the overlapping fraction of our sample has negligible mean external convergence. The external convergences for the remaining systems are likely to be not negligible and our joint lensing–dynamics analysis can be suspect of bias due to not correcting for the external convergence in our joint lensing–dynamics analysis. However, we only adopt MSD-invariant quantities that are not sensitive to the external convergence – such as the Einstein radius R_ E and the local differential term R_ Eα^''_ E/(1 - κ_ E) – as the lensing constraints in our joint lensing–dynamics analysis. Thus, our joint lensing–dynamics analysis is insensitive to the external convergence.
  
  
  
  
  In our adiabatic contraction model, we assume that the stellar mass distribution initially resembled the NFW profile following the adiabatic contraction model of <cit.>. This is certainly not true as the stellar mass delivered to the elliptical galaxies through mergers did not resemble the NFW profile at the time of the merger. However, cosmological hydrodynamical simulations have shown that such a simple model can reproduce the contracted profile of the dark matter even after gas-rich mergers, albeit with a modification in the adiabatic invariant quantity rM(r) <cit.>. We adopted the modification of <cit.>, which can approximate the contracted profiles of both <cit.> and <cit.>. However, since cosmological hydrodynamical simulations have so far been unable to match all the properties of their observational counterparts, our simple model for adiabatic contraction may not be fully justified to truthfully represent the interplay between the baryonic and the dark matter distributions. Such inconsistencies may be identified by an alternative analysis of the same sample that adopts the generalized NFW profile for the dark matter distribution, and compare the results with the ones from our the adiabatic contraction model. We leave such explorations for future studies.
  
  
  
  
  § SUMMARY
   
  
  We uniformly modelled a sample of  SLACS lenses and constrained their structural properties from a joint lensing–dynamics analysis. The lens modelling in this study is different from the original SLACS analysis, in which the lens images were modelled by fixing the logarithmic slope to γ=2 and then the logarithmic slopes were inferred from the stellar kinematics. In contrast, in this study we first estimate the logarithmic slopes only from the lensing observables, i.e., the lens image. We then combine the stellar kinematics to constrain the amount of contraction in the dark matter distribution for two models of stellar mass distribution: (i) with constant M/L and (ii) with M/L gradient. We summarize the main results of this paper below.
  
  
  	
    * From the combination of lensing and kinematic observables, we constrain the average halo response parameter μ_ν = -0.06 ± 0.04 with intrinsic scatter σ_ν≤ 0.092 (95 per cent upper limit) for a constant stellar M/L model. For a stellar M/L gradient model, we find μ_ν = -0.03^+0.04_-0.05 and σ_ν≤ 0.074. Our results are consistent with a dark matter halo described by an NFW profile with no contraction nor expansion. For comparison, the <cit.> model corresponds to ν = 1, the contraction in <cit.> simulations correspond to ν = 0.8, and the contraction in <cit.> simulations correspond to ν∼ 0.4 – which are all ruled out by our result. Our results are consistent with a scenario in which elliptical galaxies grow by dissipational processes at z≳2, steepening their dark matter halos. At later times, AGN feedback – with potential additional contributions from dynamical heating through accretion events – expand the dark matter halos back to an NFW profile, on average.
  	
    * The distribution of logarithmic slopes γ for the power-law model constrained from the imaging-only data has a median ⟨γ^ lensing_ PL⟩ = 2.08±0.03 and an intrinsic scatter 0.13 ± 0.02. This is consistent with the slope distribution from the SLACS analysis ⟨γ^ LD_ PL⟩ = 2.078 ± 0.027 with an intrinsic scatter of 0.16±0.02. We find that the NFW+stars profile constrained from our joint lensing–dynamics analysis only deviates by ≲ 5 per cent on average near the Einstein radius (R_ E∼ R_ eff/2 for our sample), with even smaller deviation at smaller scales. The small deviation in the mean compared to the intrinsic scatter explains the good agreement between the average local logarithmic slope from lensing-only data and the radially averaged logarithmic slopes from <cit.>. 
  	
  
  	
    * For the stellar M/L gradient model, we find that most galaxies do not require a significant M/L gradient. The 95 per cent upper limit for the sample mean of the exponent η in our M/L-gradient model Υ∝ (R/R_ eff)^-η is 0.02, which corresponds to ≲5 per cent decrease in M/L between 0.1R_ eff and R_ eff. Moreover, the inferred stellar masses from joint lensing–dynamics analysis for the galaxies in our sample is consistent with the Salpeter IMF, with the IMF mismatch parameter ⟨log_10α^ Salp⟩_ median = 0.00 ± 0.03 with an intrinsic scatter 0.11^+0.03_-0.02. Such a heavy IMF in the central regions of the elliptical galaxies can be explained by a different star-forming environment than in the Milky Way, e.g., the presence of turbulence in high gas density.
  	
  	
  
  
  
  In the future, larger samples of galaxy–galaxy lenses at different redshifts will be able to further constrain the evolutionary tracks of elliptical galaxies. Such larger samples can be assembled from past surveys with available high-resolution imaging and ancillary data, e.g., the full SLACS sample <cit.>, the Strong Lensing Legacy Survey (SL2S) sample <cit.>, and the SLACS for the MASSES (S4TM) sample <cit.>. Current surveys such as the Dark Energy Survey (DES) and the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys are producing new galaxy–galaxy lens candidates for confirmation and follow-up on the order of hundreds <cit.>.
  Future surveys, e.g., the Vera Rubin Legacy Survey for Space and Time, the Nancy Grace Roman Space Telescope, and Euclid, will increase the number of newly discovered galaxy–galaxy lenses to thousands <cit.>. 
   An automated and uniform modelling pipeline will be essential to study such large samples of lenses and this paper has taken the initial steps towards such an automated pipeline.
  
  
  
  § ACKNOWLEDGEMENTS
  
  We thank the anonymous referee for many useful comments that helped us improve this manuscript. We thank Matthew Auger for sharing the SLACS weak lensing measurements from <cit.> in digital form. We express gratitude to Adriano Agnello, Elizabeth Buckley-Geer, Thomas Collett, Frederic Courbin, Xuheng Ding, Aymeric Galan, Martin Millon, Veronica Motta, Sampath Mukherjee, Dominique Sluse, and Chiara Spiniello for providing suggestions that improved this analysis and manuscript. We additionally thank Matthew Auger, Chris Fassnacht, and Leon Koopmans for helpful discussions. We also thank the SLACS team for collecting the wonderful data used in this paper. AJS and TT were supported by the National Aeronautics and Space Administration (NASA) through the Space Telescope Science Institute (STScI) grant HST-GO-15320. AJS was additionally supported by a Dissertation Year Fellowship from UCLA Graduate Division. This research was supported by the U.S. Department of Energy (DOE) Office of Science Distinguished Scientist Fellow Program. TT acknowledges support by the Packard Foundation through a Packard Research fellowship and by the National Science Foundation through NSF grants AST-1714953 and AST-1906976. The development of dolphin is supported by NASA through the STScI grant HST-AR-16149.
  
  This work used computational and storage services associated with the Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research and Education’s Research Technology Group. AJS thanks Smadar Naoz for providing access to additional computing nodes on the Hoffman2 Shared Cluster.
  
  This research made use of lenstronomy <cit.>, dolphin (<https://github.com/ajshajib/dolphin>), fastell <cit.>, numpy <cit.>, scipy <cit.>, astropy <cit.>, jupyter <cit.>, matplotlib <cit.>, seaborn <cit.>, sextractor <cit.>, emcee <cit.>, colossus <cit.>, pomegranate <cit.>, and chainConsumer (<https://github.com/Samreay/ChainConsumer>).
  
  
  
  § DATA AVAILABILITY
  
  
  The HST imaging data used in this paper are publicly available from Mikulski Archive for Space Telescopes (MAST). The other ancillary measurements – i.e., the stellar kinematics and the weak lensing measurements – are obtained from previous studies and we provide references to these previous studies. These measurements can be obtained either directly from the paper or by requesting the corresponding author of the related paper. The lens modelling codes lenstronomy and dolphin used in this paper are publicly available on GitHub.
  
  
  
  
  
  
  
  
  
  
  mnras
  
  
  
  
  
  
  
  
  
  
  
  § MASS AND LIGHT ALIGNMENTS
   
  
  Alignment between the dark and baryonic components of mass distribution can be used to validate the predictions from simulations. In elliptical lens galaxies, the misalignment between the mass and light is observed to be within ±10 in the absence of large external shear effects (≳0.1), whereas larger misalignments are usually accompanied with large external shear <cit.>. This observation agrees well with the Illustris simulation <cit.>. However, there have been mismatching reports in the literature on the correlation between the light and mass ellipticities or axis ratios, which can most likely be attributed to different selection functions. The differences in the selection functions can arise from different lens-finding methods and from different population of lenses – e.g., quads or doubles, galaxy–galaxy lenses or lensed quasars. For example, <cit.> report strong correlation between mass and light ellipticities. In comparison, <cit.> find weak to no correlation between the mass and light ellipticities. 
  
  In this appendix, we investigate the alignment between mass and ellipticities and their correlation with some other model parameters. We present the results in Appendix <ref> and discuss them in Appendix <ref>.
  
  
  
   §.§ Results
   
  
  
  
  
  
  We illustrate the distributions of these quantities in Figure <ref>. In the following subsections, we present the offset between the centroids and the misalignments between ellipticity magnitudes and orientations.
  
  
  
    §.§.§ Centroid
  
  
  We find an rms scatter of 67 ± 5 mas in the centroid offset. For our fiducial cosmology, the 68 per cent and 95 per cent upper limits of the absolute offsets are 218 ± 19 pc and 426 ± 68 pc, respectively (Figure <ref>, first panel). 
  
  
  
  
    §.§.§ Ellipticity
  
  We find a weak correlation between the axis ratios of the mass and light distributions with a biweight mid-correlation of r=0.52± 0.10 (Figure <ref>, third panel). Biweight mid-correlation is similar to the Pearson's r coefficient. However, it depends on the median instead of the mean, thus it is more robust against outliers <cit.>. We get the uncertainty in the mid-correlation by sampling the axis ratios of all the lenses 1000 times from their posterior distributions.
  
  We check if the lens systems are double-image or quadruple-image systems by checking the number of images produced by the lens models for a hypothetical point source located at the source centres. We find 19 out of 23 lenses to be double-image systems and the remaining 4 to be quadruple-image systems. The majority of doubles explain the lower ellipticity and shear distributions observed in the SLACS sample compared to that of the quadruply lensed quasars (Figure <ref>).
  
  
  
    §.§.§ Position angle
  
  
  While comparing the position angles of the mass and light, we ignore 4 lenses with q_ m > 0.9 and q_ L > 0.9, as the estimate of the position angle in low-ellipticity cases can be unreliable. Out of the remaining 19 lenses, we find 16 to have misalignment angle ΔPA≤12 deg. The three systems with ΔPA>12 deg have relatively higher external shear γ_ ext≳ 0.05 (Figure <ref>, right-hand panel).
  
  
  
   §.§ Discussion
   
  We find moderate correlation between axis ratios of the mass and light (biweight mid-correlation r=0.55±0.09). Although <cit.> and <cit.> both report strong correlation between mass and light ellipticities, we compute the biweight mid-correlation from the values reported by these authors to find r=0.26 ± 0.07 (weak correlation) from the galaxy–galaxy lenses of <cit.> and r=0.91 ± 0.02 (very strong correlation) from the lensed quasars of <cit.>. Additionally, <cit.> find strong correlation (r = 0.77 ± 0.16) from a sample of 9 galaxy–galaxy lenses from the Cambridge And Sloan Survey Of Wide ARcs in the skY (CASSOWARY). In contrast, <cit.> find very weak to no correlation between mass and light for a sample of 13 quadruply lensed quasars (r=0.13 ± 0.10). Similarly, <cit.>, <cit.>, and <cit.> also find no correlation between the mass and light ellipticities. Some of the differences between these studies, including ours, can be due to the differences in the data quality, analysis techniques, and adopted models. Some of these differences can also be due to the selection function, for example between the galaxy–galaxy lenses and lensed quasars. SLACS is deflector-selected and mostly comprises doubly lensed objects, whereas quadruply lensed quasars are expected to be strongly selected to favor big inner caustics <cit.>. Therefore, they favor high ellipticity and shear. Moreover, the SLACS sample has a smaller average redshift than the sample of quasar lenses, e.g., from <cit.>. Thus, the line-of-sight effect is much less important for SLACS than for the lensed quasars. As <cit.> model their quasar lens sample using HST imaging of comparable quality, similar models, and the same modelling software as in this paper, the observed differences between the SLACS sample and the quasar lens sample in Figure <ref> have to arise from the differences in the selection functions.
  
  We find that in most galaxies with moderate ellipticity (q ≤ 0.9), the major axes of the mass and light distributions are well aligned within ± 12 deg. The lenses with large misalignment (ΔPA>12 deg) also have relatively larger (γ_ ext≳ 0.05) external shear within the sample. This result is consistent with previous studies  <cit.>. The absence of systems with large misalignment angle and low external shear is consistent with the prediction of galaxy formation models that highly misaligned orbits in isolated galaxies are unstable and thus rare <cit.>. Only in blue starburst galaxies – unlike the galaxies in our sample – constant gas-flow can sustain highly misaligned orbits <cit.>.
  
  
  
  
  § ALGORITHM TO MASK THE LENSED ARCS
   
  
  We illustrate the algorithm for creating the mask for the lensed arcs in Figure <ref>.
  
  
  
  
  
  
  
  
  § REMAINING LENS MODELS
   
  
  In this appendix, we show the rest of the lens models  in Figures <ref> and <ref> continuing from Figure <ref>.
  
  
  
  
  
  
  
  
  
  § SELECTED SLACS GALAXIES FOR MODELLING
   
  
  We provide the list of selected galaxies in Table <ref>.
  
  
  
  
  
  
  
  
  
  