Project Dinos I: A joint lensing–dynamics constraint on the deviation from the power law in the mass profile of massive ellipticals

The mass distribution in massive elliptical galaxies encodes their evolutionary history, thus providing an avenue to constrain the baryonic astrophysics in their evolution. The power-law assumption for the radial mass profile in ellipticals has been sufficient to describe several observables to the noise level, including strong lensing and stellar dynamics. In this paper, we quantitatively constrained any deviation, or the lack thereof, from the power-law mass profile in massive ellipticals through joint lensing– dynamics analysis of a large statistical sample with 77 galaxy–galaxy lens systems. We performed an improved and uniform lens modelling of these systems from archival Hubble Space Telescope imaging using the automated lens modelling pipeline dolphin. We combined the lens model posteriors with the stellar dynamics to constrain the deviation from the power law after accounting for the line-of-sight lensing effects, a first for analyses on galaxy–galaxy lenses. We find that the Sloan Lens ACS Survey (SLACS) lens galaxies with a mean redshift of 0.2 are consistent with the power-law profile within 1.1 𝜎 (2.8 𝜎 ) and the Strong Lensing Legacy Survey (SL2S) lens galaxies with a mean redshift of 0.6 are consistent within 0.8 𝜎 (2.1 𝜎 ), for a spatially constant (Osipkov–Merritt) stellar anisotropy profile. We adopted the spatially constant anisotropy profile as our baseline choice based on previous dynamical observables of local ellipticals. However, spatially resolved stellar kinematics of lens galaxies are necessary to differentiate between the two anisotropy models. Future studies will use our lens models to constrain the mass distribution individually in the dark matter and baryonic components.


INTRODUCTION
Strong gravitational lensing is the phenomenon where multiple images of a background object appear due to the deflection of light by the gravitational effect of a massive foreground object.The foreground object is typically a massive galaxy or a galaxy cluster.Strong lensing at the galaxy scale is a valuable probe of the mass distribution of galaxies yielding numerous applications in astrophysics and cosmology (for a review, see Shajib et al. 2022b).
For galaxy-scale deflectors, a background galaxy can be imaged into multiple lensed arcs or a full Einstein ring.Furthermore, multiple images of a point-like source -e.g., of a supernova (Goobar et al. 2017) or an active galactic nucleus (AGN; Walsh et al. 1979) -can also appear.Galaxy-galaxy lens systems, i.e., systems that do not contain a prominent resolvable point source in the background galaxy, are more commonly discovered than multiply-imaged point source systems, as supernovae and AGNs are rarer than galaxies.Some of the largest samples of galaxy-galaxy lens systems to date are the Sloan Lens ACS 1 (SLACS) survey (Bolton et al. 2006), the Strong Lensing Legacy Survey (SL2S; Gavazzi et al. 2012), the SLACS for the Masses (S4TM; Shu et al. 2015), the BOSS 2 Emission-Line Lens Survey (BELLS; Brownstein et al. 2011), and the BELLS for the GALaxy-Ly EmitteR sYstems (BELLS GALLERY; Shu et al. 2016).These samples of galaxy-galaxy strong lenses have been used to study the internal structure of elliptical galaxies by modelling the imaging data from the Hubble Space Telescope (HST) and ground-1 Advanced Camera for Surveys 2 Baryon Oscillation Spectroscopic Survey based telescopes (e.g., Auger et al. 2009Auger et al. , 2010a;;Treu et al. 2010;Sonnenfeld et al. 2013b;Sonnenfeld et al. 2015).

Radial form of the mass profile in massive ellipticals
In the current paradigm, massive elliptical galaxies are considered to be the end product of hierarchical mergers of dark matter haloes along with the stars and gas contained within them.In this formation scenario, the elliptical galaxies undergo two processes, which can sometimes co-occur.
For one process, the baryonic gas dissipatively cools within the dark matter halo, starting from an initial Navarro-Frenk-White (NFW) mass distribution at  ≫ 2 (Navarro et al. 1996(Navarro et al. , 1997)).The dark matter halo also contracts in response to the baryonic infall, becoming steeper than the () ∝  −1 profile at the inner region (Blumenthal et al. 1986;Gnedin et al. 2004;Abadi et al. 2010).For the other process, galaxies undergo multiple mergers and accretion events to eventually reach a dynamically pressure-supported state.These events cause the individual mass density profiles of the dark matter and the baryons, and thus the total density profile, to change in various ways.While dissipationless mergers do not lead to a change in the dark matter density profile (Gnedin et al. 2004;Ma & Boylan-Kolchin 2004), dissipational mergers lead to a steeper slope at the inner region (Sonnenfeld et al. 2014).Furthermore, baryonic feedback mechanisms -for example, feedback from stellar wind, supernovae, and AGN -can cause an outflow of gas, making the dark matter halo expand in response, which leads to a shallower total density profile (Naab et al. 2007;Duffy et al. 2010;Johansson et al. 2012;Dubois et al. 2013).These feedback mechanisms can also counteract the initial dark matter contraction in the first place (Governato et al. 2012;Pontzen & Governato 2012).Indeed, recent cosmological hydrodynamical simulations, such as the Magneticum and the Illus-trisTNG, predict a shallowing trend in the density profile slope with decreasing redshift at  < 1, thus positing gas-poor mergers being the dominant channel for the growth of massive ellipticals (Remus et al. 2017;Wang et al. 2019).Therefore, observationally constraining the total density profile of elliptical galaxies at various redshifts would be a valuable approach to studying their formation history and mechanisms.By fine-tuning the baryonic physics prescriptions used in the cosmological hydrodynamical simulations to reproduce the density profile trends of these galaxies, we can gain insights into the underlying physical processes that drive their evolution.
The total density profile in elliptical galaxies is found to be well approximated by a power-law mass distribution (i.e., () ∝  − ) using observations from the combination of strong lensing and stellar dynamics (e.g., Treu et al. 2006;Shajib et al. 2021), from stellar dynamics only (e.g., Cappellari 2016;Derkenne et al. 2021), from X-ray observations (e.g., Humphrey & Buote 2010), and from a combination of strong lensing, weak lensing, and dynamics (e.g., Gavazzi et al. 2007).Although neither the dark matter nor the baryons individually follow the power law, the total density profile following the power law is called the 'bulge-halo conspiracy' (Dutton & Treu 2014;Treu et al. 2006).The total density profile in elliptical galaxies is typically observed to be slightly steeper than isothermal (i.e.,  ∼ 2.1; Auger et al. 2010b;Ritondale et al. 2019;Shajib et al. 2021).
More recently, however, Etherington et al. (2023b) suggest that the actual mass distribution may potentially deviate from the power-law profile based on their observation that the surface mass density correlates differently with the logarithmic slope  of the same systems measured from lensing-only and joint lensing-dynamics analyses.However, such a potential deviation from the power-law profile is yet to be quantitatively constrained with a high confidence level, and unaccounted systematic uncertainties pertaining to the lens or dynamical modelling are not yet ruled out as the source of the discrepancy mentioned above.Nonetheless, it is essential to use an accurate model of the mass distribution in the lensing galaxies before extracting observed structural properties to compare with the predictions from simulations.Otherwise, modelling systematics can bias the observed properties of the lensing galaxies, invalidating their use to confirm or rule out simulation predictions.
Furthermore, any significant deviation from the power-law mass profile has important implications for time-delay cosmography (for recent reviews, see Treu et al. 2022;Birrer et al. 2022).The 2.4 per cent measurement of the Hubble constant from the  0 Lenses in the COSMOGRAIL's Wellspring (H0LiCOW) collaboration assumes power-law mass profile in the lensing galaxies (e.g., Suyu et al. 2013;Wong et al. 2020).If the power-law mass profile is accurate, then this measurement rules out systematics in the local  0 measurement based on the cosmic distance ladder of type Ia supernova calibrated with Cepheids (e.g., Riess et al. 2022) and confirms the existence of new physics beyond the Λ cold dark matter (ΛCDM) cosmology (for a review, see Di Valentino et al. 2021).However, any potential deviation from the power-law profile in elliptical galaxies can be a source of major systematic in the time-delay cosmography.Therefore, the Time-Delay Cosmography (TDCOSMO) collaboration simultaneously constrains the Hubble constant and the galaxy mass profile by combining stellar kinematics with the lensing information (Birrer et al. 2020;Shajib et al. 2023).

Project Dinos
In this paper, we introduce Project Dinos (PI: Shajib)3 , whose primary goal is to study elliptical galaxy evolution using strong-lensing galaxies.This project aims to constrain the mass distribution in the elliptical galaxies at various redshifts, both in the total mass and individually in dark matter and baryons, and compare them with cosmological hydrodynamical simulations to understand the impact of baryonic feedback and mergers on the elliptical galaxy evolution at  ≲ 1.
In this first paper of the Project Dinos series, we quantified any deviation from the power-law profile in elliptical galaxies.We constrained the logarithmic slope of the mass profiles from a large statistical sample collated from multiple archival samples: SLACS, SL2S, and BELLS.The principal criteria to adopt these specific surveys are large sample size, diversity in the covered redshift range, and availability of the HST archival imaging and stellar kinematic measurements.
In most of the previous analyses of these lens samples, the imaging data were only used to measure the Einstein radius by modelling the mass profile of the system with a singular isothermal ellipsoid (SIE).The logarithmic slope of the mass distribution was then constrained by combining the stellar dynamics with the lensing-inferred enclosed projected mass.However, if the data quality is sufficient, it is also possible to constrain the local logarithmic slope of the mass profile directly from the imaging data (Birrer 2021).Then, combining the stellar dynamics with the lensing constraints can provide us with more detailed information on the mass structure, for example, by breaking the degeneracy between the normalizations of the stellar and dark matter components in the lensing galaxy (e.g., Shajib et al. 2021).Furthermore, the power-law mass models for a large sample of galaxy-galaxy lens systems with available kinematics are required to provide a prior on the mass profile shape to improve the Hubble constant measurement from time-delay cosmography (Birrer et al. 2020(Birrer et al. , 2021)).In this paper, we modelled HST archival images of galaxy-galaxy lens systems with more flexible models than the simple SIE model adopted by previous studies on these same datasets.In this way, we extracted more lensing information from the imaging data beyond only constraining the Einstein radius.Furthermore, our sample is the largest to date with state-of-the-art power-law models of strong lenses.Future papers in this series will investigate the dark matter and baryonic distribution in elliptical galaxies and their evolution across redshift using the sample presented in this paper.
This paper is organized as follows.In Section 2, we present our lens sample, the HST imaging data used for lens modelling, and the available ground-based stellar velocity dispersion measurements.In Section 3, we describe the lens modelling procedure using an automated and uniform modelling framework.We then compare our lens modelling results with previous studies in Section 4. In Section 5, we describe the hierarchical Bayesian framework for combining the lens model posteriors with the stellar dynamics to quantify any deviation from the power-law mass model.We discuss the implications of our results in Section 6 and conclude the paper in Section 7. Throughout the paper, we adopt the flat ΛCDM model as our fiducial cosmology with  0 = 70 km s −1 Mpc −1 and Ω m = 0.3.

LENS SAMPLES AND DATA
We chose an initial lens sample that consists of lens systems from the SLACS, SL2S, and BELLS samples.We used all the available archival HST images of these systems in the ultraviolet and visible (UVIS) bands.The specific HST filters used to model each lens system are provided in Appendix A. We also obtained the redshift and velocity dispersion measurements for these systems from the literature.

Lens sample selection
In this section, we describe the parent lens survey for each of our chosen subsamples and how we selected these initial subsamples for lens modelling.The lenses in the SLACS survey were found by selecting targets with higher-redshift emission lines and lowerredshift absorption features measured using the spectrum obtained from the Sloan Digital Sky Survey (SDSS; Bolton et al. 2006).The sample is then further refined by selecting systems showing clear evidence of multiple imaging of the background source by the lensing galaxy.Following Shajib et al. (2021), we chose a subsample of 50 systems to model out of the full sample of 85 'Grade A' lenses (described in Auger et al. 2009).This subsample was selected by visually inspecting the lens images and excluding systems where (i) there are nearby satellites or line-of-sight galaxies, (ii) there is highly complex source morphology, (iii) the deflector galaxy is disc-like, or (iv) there is no archival HST imaging data for the system in the visible band (F555W or F606W).Requirements (i) and (ii) were imposed due to simpler lens systems being more suitable for uniform modelling with a higher probability of success, whereas requirement (iii) was set so that we only have elliptical galaxies in the sample.We imposed requirement (iv) so that the lens light profile of the galaxies can be measured uniformly across the SLACS sample.
Another sample analysed in this paper is the SL2S sample.The sample was identified using images from the Canada-France-Hawaii Telescope Legacy Survey (Gavazzi et al. 2012(Gavazzi et al. , 2014;;More et al. 2012).The candidates were identified by finding blue features around elliptical galaxies consistent with the presence of gravitationally lensed arcs.The samples were confirmed using follow-up HST imaging and Very Large Telescope (VLT), Keck, or Gemini spectroscopy.Out of the 56 lenses in the sample, we selected 31 lenses primarily based on the availability of visible band HST images for that particular system.However, we also excluded two systems (SL2SJ0213−0743, SL2SJ1405+5502), which have -band HST images, but the nearby foreground objects are too complex to mask effectively.The last sample used in the analysis is the BELLS sample.The sample consists of candidates discovered by higher redshift emission lines within the BOSS galaxies (Brownstein et al. 2011) in a similar method as the SLACS lenses.Our sample consists of the 22 'Grade A' systems for the BELLS sample.'Grade B' lenses and below were excluded as they either lack the counter images or have an insufficient signal-to-noise ratio (/) for modelling.

Processing of the HST imaging data
The imaging data for the selected lenses are from HST's Advanced Camera for Surveys (ACS), Wide Field and Planetary Camera 2 (WFPC2), and Wide Field Camera 3 (WFC3).Table 1 provides the HST program IDs, instruments, and filters for the imaging data of our lens sample.
We combined images with multiple exposures using the astrodrizzle software package (Avila et al. 2015).The drizzling process with astrodrizzle also allowed us to remove cosmic rays from the single exposure image.The drizzling pixel scale was chosen based on the camera used: 0.04 arcsec for WFC3, 0.05 arcsec for ACS, and 0.10 arcsec for WFPC2.An exception was made for the SLACS WFPC2 images, where we used the data reduced by Auger et al. (2009) with the pixel scale set to 0.05 arcsec.Auger et al. (2009) employed a custom data-reduction pipeline where the sub-pixel shifts between exposures were used to allow the drizzling pixel scale to be significantly smaller than the detector pixel scale.To achieve a higher image quality after drizzling, we left the orientation of the final output image with respect to the RA-Dec coordinate system to be the same as the individual raw exposures.We then used Source-Extractor to estimate the mean background light and subtract it from each image (Bertin & Arnouts 1996).
We find up to a 0.5 arcsec global offset in the world coordinate system (WCS) between images from different bands due to differences in their absolute astrometric reference.This often occurs when comparing images from different HST cameras.Thus, to ensure that the coordinate systems of the images are correctly aligned across all filters, we constrained the centroid of the deflector galaxy in each filter individually by fitting a single Sérsic profile to the light distribution within the central 0.6 arcsec radius of the deflector galaxy using lenstronomy (Birrer & Amara 2018;Birrer et al. 2021).We then corrected the relative offsets between the WCSs using the difference between the fitted centroids.We find this method sufficient in reducing the modelling residuals caused by image misalignment between the filters.As shown in Table 1, the HST images for some SLACS lenses only have a single exposure.Therefore, astrodrizzle cannot remove cosmic rays from these images.We used the deepCR software program for these images to automatically detect and mask the cosmic rays (Zhang & Bloom 2020).deepCR uses a neural network trained on HST images to detect and mask cosmic rays.Figure 1 shows an example of the mask produced by deepCR for the system SDSSJ0819+4534.
Since images from different cameras and filters have different point spread functions (PSF), we simulated the PSFs for each combination of camera and filter with tinytim (Krist et al. 2011) using a G2V star's spectral energy distribution.When producing the PSF, we also accounted for zero-point offsets produced by the time-dependent aberration changes on the mirror (Krist et al. 2011).The pixelized PSFs were set to have the same pixel sizes as the corresponding drizzled HST images.Shajib et al. (2021) find that different tinytim configurations and various empirical methods to estimate the PSF (such as ePSF from stars within the corresponding image; Anderson & King 2000) produce negligible differences in sample mean of the model parameters.

Measured stellar velocity dispersions
We obtained the measured redshift of the source and deflector galaxies and the single-aperture stellar velocity dispersion of the deflector galaxy from the literature.For the SLACS systems, we used the same SDSS velocity dispersion measurements as Birrer et al. (2020).The circular fibre used by the SDSS program had a diameter of 3 arcsec, and the nominal full width at half-maximum (FWHM) of the seeing was 1.4 arcsec in the -band (Bolton et al. 2008;Birrer et al. 2020).The velocity dispersion for the SLACS systems was then measured using templates described in Shu et al. (2015).Birrer et al. (2020) find that the uncertainty in the SDSS velocity dispersion measurements was underestimated with an unaccounted systematic uncertainty of 6 per cent.We also accounted for such underestimated uncertainties as described in Section 5.3.For the SL2S lenses, we used the velocity dispersions measured by previous studies using long-slit spectra obtained from the Keck Observatory, the VLT, and the Gemini North Telescope (Sonnenfeld et al. 2013b;Sonnenfeld et al. 2015, see therein for the slit dimensions and seeing).The velocity dispersions for the BELLS systems were obtained from BOSS spectroscopy as part of the SDSS-III program.The SDSS-III program used an upgraded optical spectrograph having a smaller fibre diameter of 2 arcsec and the observations had a 1.3 arcsec nominal seeing in the -band (Brownstein et al. 2011).

LENS MODELLING
We used the automated lens modelling pipeline dolphin4 to perform uniform lens modelling for all the lenses in our sample (Shajib et al. 2021).dolphin uses the publicly available software package lenstronomy5 as its modelling engine (Birrer & Amara 2018;Birrer et al. 2021).The effectiveness of the lenstronomy package was demonstrated with the Time-Delay Lens Modelling Challenge, where multiple participating teams used the software package to successfully recover the hidden lens model parameters (Ding et al. 2021).lenstronomy has also been used in many other applications related to strong gravitational lenses, such as in time-delay cosmography (Birrer et al. 2019;Shajib et al. 2020Shajib et al. , 2022a) ) and in studying dark matter substructures (Gilman et al. 2020).Whereas lenstronomy provides a powerful and feature-rich API to allow fine-tuning of the lens model settings for individual systems, dolphin vastly simplifies the modelling procedure for the user through a semi-automated decision-making process without requiring a large amount of investigator time spent for fine-tuning.
In Section 3.1, we briefly review the theory of strong gravitational lensing and introduce the relevant terminologies.Then, in Section 3.2, we discuss the model components used to describe the deflector and source galaxy's mass and light distributions.Section 3.3 presents the methods used to obtain the model posteriors from the data.Finally, in Section 3.4, we present our final lens model sample and discuss our quality assurance procedure for the lens models.

Introduction to strong lensing formalism
In the single-lens-plane regime, the effect of gravitational lensing is described with the lens equation, which maps the source plane vector  to the image plane vector  as where () is the deflection angle.The lensing potential  relates to the deflection angle as For a 2D projected surface mass density Σ(), we can define the dimensionless convergence quantity as where Σ crit is the lensing critical surface mass density, defined as where  is the speed of light,  is the gravitational constant, and  l ,  s ,  ls are the angular diameter distances from the observer to the lens, from the observer to the source, and from the lens to the source, respectively.The lensing potential is then connected to the convergence via the Poisson equation Therefore, the distortion of the source galaxy's light in a lensing system, described by (), can be used to probe the surface mass density Σ().

Definitions of model components
The lens model consists of three main components: the mass profile of the deflector galaxy, the light profile of the deflector galaxy, and the light profile of the source galaxy.We modelled the mass profile of the deflector galaxy with a power-law ellipsoidal mass distribution (PEMD; Barkana 1998).The convergence of this profile is given by where  E is the Einstein Radius,  is the logarithmic slope of the mass distribution, and  m is the axis ratio.The coordinates ( 1 ,  2 ) are defined by rotating the (RA, Dec) coordinate system by the major axis's position angle  m .We also added a 'residual shear' field (commonly referred to as the 'external shear' in the literature) in the lens model parametrized by a shear magnitude  shear and angle  shear .This residual shear field accounts for additional shear contributed by the line-of-sight structures and any unaccounted angular structure in the central deflector (e.g., Etherington et al. 2023a).
For the light profile of the deflector galaxy, we used the elliptical Sérsic profile given by where  eff is the effective radius of the profile,  s is the Sérsic index,  L is the ratio between semi-major and semi-minor axis,   = 1.999  − 0.327 is the normalizing factor to make  eff the half-light radius, and  0 is the amplitude (Sérsic 1968).Similar to the mass profile, the coordinates ( 1 ,  2 ) are defined by rotating the (RA, Dec) coordinate system by the major axis's position angle  L .
As a single Sérsic profile often cannot accurately fit the light distribution at the galaxy's centre (Claeskens et al. 2006;Suyu et al. 2013), we used two Sérsic profiles to allow more radial freedom in fitting the galaxy's light distribution.To reduce the degeneracies between the two Sérsic profiles, we fixed the Sérsic indices of the profiles as  s = 4 (i.e., the de Vaucouleurs profile) and  s = 1 (i.e., the exponential profile), respectively.To minimize the number of unnecessary free parameters in the model, we joined the ellipticity parameters  L and  L between the two Sérsic profiles.
We find that the two Sérsic profiles are still inadequate in some cases to reproduce the light profile of the deflector at the very centre (within ∼0.4 arcsec).Thus, to avoid a large residual in the centre introducing any bias, we adopt a circular mask at the centre of the deflectors to exclude the residuals within the masked region from the computed likelihood.In our default setup, we set the radius of the mask to be 0.4 arcsec.However, for systems with lensed arcs falling within the central 0.4 arcsec, we chose a smaller radius to avoid masking out the arcs.Nearby galaxies along the line of sight exist in some lens systems, for example, SDSSJ0944−0147 or SDSSJ1215+0047.For these systems, we manually introduce an appropriately sized circular mask to block out the light from the lineof-sight galaxies.The mask settings for each lens system are given in Appendix A.
For the source galaxy's light distribution, we used the combination of an elliptical Sérsic profile and a basis set of multiple shapelet components (i.e., weighted 2D Gauss-Hermite polynomials; Refregier 2003;Birrer et al. 2015).The extent of all the shapelet components simultaneously scale with a single parameter .The order parameter  max determines the number of shapelet components  shapelets = ( max + 1)( max + 2)/2.Initially, we set  max = 6 as the default setting for all systems.However, for systems with complicated source light profiles, we adopted a higher value of shapelet order parameter,  max , chosen by trial and error to improve the model's fit.The specific choices made for  max in each system are provided in Appendix A.
We simultaneously fit all the available HST filters for a given system.In doing so, we allowed the amplitudes and scale radii of all the light profiles to be free across different filters while imposing the same deflector mass profile for all the filters.
The lens model discussed in this section was designed to increase the success rate in uniformly modelling the galaxy-galaxy lenses in our sample.However, several systems contain complex featureseither in the source light or in the deflector environment -that must be manually modelled on a case-by-case basis.For this paper, we consider this type of lens system as a 'failure case' and exclude those from our final lens model sample.A more detailed discussion of how we selected our final lens model sample is given in Section 3.4.

Modelling procedure
We sampled from the posterior probability distribution function (PDF) of the model parameters with the Markov chain Monte Carlo (MCMC) method using emcee (Goodman & Weare 2010;Foreman-Mackey et al. 2013).To achieve faster convergence in the MCMC sampling, we obtained a point close to the maxima of the posterior PDF using the particle swarm optimization (PSO) method (Kennedy & Eberhart 1995).We then initiated the random walkers for the MCMC sampling within a small region around the best solution from the PSO step.We visually inspected for MCMC convergence by ensuring that the median and standard deviation of the parallel walker positions had reached equilibrium for at least 1000 steps.
We improved the efficiency of finding the posterior probability distribution by following the same optimization sequence or 'recipe' as Shajib et al. (2021).The optimization recipe first adopted a SIE mass model for the deflector with no residual shear (i.e.,  shear = 0).The Einstein radius of the system,  E , was also fixed to a value obtained from previous analyses (i.e., Auger et al. 2009;Sonnenfeld et al. 2013a;Brownstein et al. 2011), while the ellipticity parameters were set to the values for the circular case ( m = 1).For the first step of the optimization sequence, the pipeline optimized the deflector galaxy's light profile parameters with all the other parameters fixed.While optimizing for the deflector galaxy's light profile, the pipeline also created a mask to cover the lensed arcs.Then, the source light parameters were optimized after fixing the lens light and deflector mass parameters.The later steps involved freeing the previously fixed parameters and finding the best-fitting values for all the parameters simultaneously.
Since we do not know the order of magnitude of the shapelet scale parameter , we adopted a log-uniform prior or Jeffreys prior on , that is, () ∝ 1/.We also adopted an empirically motivated uniform prior such that the deflector mass profile's position angle  m and the light profile's position angle  L are aligned within 15 deg.To achieve faster convergence in the optimization, we set the edges of this uniform prior to fall as a one-sided Gaussian function with a small standard deviation (0.03 deg) instead of having a sharp cut.This choice was motivated by the finding of previous studies that the misalignment between the deflector mass and light profiles is less than 10 deg when the residual shear magnitude is small ( ext < 0.1) (e.g., Treu et al. 2009;Sluse et al. 2012;Shajib et al. 2019).The close alignment between the deflector mass and light profiles is also supported by the Illustris simulation (Xu et al. 2017).

Quality assurance for lens models
Once the MCMC chains had converged, we excluded the lens systems that fell into any of these scenarios: (i) the reconstructed lens image missing prominent lens features observed in the HST image, (ii) poor source reconstruction indicated by prominent model residuals, (iii) the reconstructed source retaining the distorted shape of the lensed arc, or (iv) the logarithmic slope  converging to extremely atypical values ( < 1.4 or  > 2.8).The  values usually diverge in such a way when the imaging data do not contain enough lensing information to constrain it (Shajib et al. 2021).Examples include when the lensed arcs are too faint or when no counter-image can be detected in the image.For lens systems without any detectable counter-image above the noise level, our pipeline tends to reconstruct the source to retain the distorted shape of the single arc visible in the system.
Using the modelling procedure and selection process described above, we successfully modelled 77 lens systems from an initial sample of 105 lens systems.We break down the number of successfully modelled lens systems for each survey in Table 2. Figures 2 and 3 illustrate the best-fitting models for two systems, SDSSJ1627−0053 and SL2SJ1401−5544, as examples.All the other models are illustrated in Appendix B.
We provide the point-estimates of the lens model parameters and 1 uncertainties for the SLACS, SL2S, and BELLS lenses in Tables 3, 4, and 5, respectively.In Figure 4, we show the distribution of the redshift, central velocity dispersion, and other model parameters for the SLACS, SL2S, and BELLS lenses.In the next section, we compare our lens model parameters with those from previous studies and estimate the modelling systematics in the logarithmic slope  based on this comparison.

MODEL COMPARISON WITH PREVIOUS STUDIES
In this section, we compare our lens model posteriors with previous analyses.We compare our parameter posteriors to those previously obtained from imaging data only with lens models similar to ours in Section 4.1 and to those obtained from previous joint lensingdynamics analysis in Section 4.2.

Comparison with previous lensing-only analyses
We compare our lens model posteriors with those obtained by Shajib et al. (2021) and Etherington et al. (2022), who also adopted the PEMD for the lens mass model.Shajib et al. (2021) modelled 23 SLACS lenses, all of which are included in our analysis.Out of the 43 SLACS lenses modelled by Etherington et al. (2022), 21 overlap with our analysis, and 14 overlap with the analysis by Shajib et al. (2021).
For our measured Einstein radii, we find a root-mean-square (rms) difference of 1.5 per cent with those from Shajib et al. (2021) and 5.6 per cent with those from Etherington et al. (2022) (Figure 6).The high rms difference in the Einstein radii between our analysis and Etherington et al. (2022) can be attributed to the outliers SDSSJ0912+0029, SDSSJ1213+6708, SDSSJ1218+0830, for which the lensed arcs have relatively low / (see the deflectorlight-subtracted images illustrating the arcs in the fourth columns of Figures B1 and B2).These three outliers differ in the Einstein radius by more than 10 per cent between our analysis and Etherington et al. (2022).Interestingly, all of these three outliers were considered failure cases in the automated modelling done by Shajib et al. (2021), further highlighting the low / in the lensed arcs as an impediment to robust lens modelling.Excluding these outliers, the rms difference between our analysis and Etherington et al. (2022) falls down to 2 per cent, consistent with the expected systematic level empirically estimated by Bolton et al. (2008), Sonnenfeld et al. (2013a), and Etherington et al. (2022).We also exclude these three outliers in other comparisons for the remainder of this section, but we include them in our hierarchical analysis in Section 5.
We find a significant difference in the logarithmic slope  between different studies.The median absolute difference in  is 8.2 per cent between our lens models and those from Shajib et al. (2021), and 9.7 per cent between our lens models and those from Etherington et al. (2022).These differences are much larger than the median 3 per cent statistical uncertainty we estimated from the MCMC chains (Figure 6).Our measured logarithmic slopes correlate only weakly with those from Shajib et al. (2021) (Pearson correlation coefficient  = 0.30 ± 0.11) and Etherington et al. (2022) ( = 0.17 ± 0.16).The weak correlation between our  measurements and Shajib et al. ( 2021)'s is noteworthy, as both analyses used the dolphin modelling pipeline with very similar lens model configurations.A major difference, however, is that we use multi-band modelling, whereas Shajib et al. (2021) only modelled a single HST filter.The apparent differences in the logarithmic slope of the mass profile between the studies suggest that there might be an additional modelling systematic uncertainty not accounted for by the statistical uncertainty.
Since Etherington et al. (2022) used a different source reconstruction method (adaptive brightness-based pixelisation and regularisation grid) implemented in a different modelling pipeline PyAu-toLens (Nightingale et al. 2021), we estimated the modelling systematic uncertainty in  by comparing our model predicted values with theirs, excluding the three outlier systems mentioned above.We did not use the comparison with the values from Shajib et al. ( 2021) since the difference is smaller.Thus, the comparison with Table 3. Lens model parameters for the SLACS systems.Here,  E is the Einstein radius,  is the logarithmic slope of the mass profile,  m is the mass axis ratio,  m is the mass position angle,  shear is the residual shear magnitude,  shear is the residual shear angle,  L is the light axis ratio,  L is the light position angle, and  eff is the effective or half-light radius of the light profile.The angles  m ,  shear , and  L are defined north of east.The point estimates are the medians of the 1D marginalised posteriors, and the 1 uncertainties represent the 16 th and 84 th percentiles.For the logarithmic slope of the mass profile , we provide an additional modelling systematic uncertainty.Additional details about estimating the systematic uncertainty are provided in Section 4.1.The lens light parameters correspond to the -band filter (F606W, F600LP, or F555W).The effective radii  eff, were obtained from separately fitting larger -band HST image cutouts -that contain the full extent of the deflector galaxies -after subtracting the lens-model-predicted arcs (see Section 5.2 for more details).This table can be easily retrieved at https://www.projectdinos.com/dinos-i.Etherington et al. (2022) provides the most conservative estimate of the systematic uncertainty pertaining to lens modelling.We assumed that systems with higher 'lensing information' would allow less modelling systematics to manifest in the  measurement.We first quantified this lensing information using a weighted / quantity as

Lens system
where the summation index  goes over all the pixels on the lensed arcs,   is the flux in the lensed arcs after subtracting off the deflector light, and   is the total noise level.We included the pixels on the arcs only with fluxes at least three times the background noise level.Among these pixels, however, we excluded the ones that belong to blobs with less than five contiguous pixels to disregard spurious blobs that are less likely to be lensed features.The ad hoc weight   for    the th pixel is given by where   is the radial distance of the pixel from the centre of the deflector galaxy,   is the azimuthal angle of the pixel, and  ref is a reference angle which we chose to correspond to the brightest pixel on the arc.The exponents  and  tune the impact of this weighting scheme.For  = 0, the definition in equation ( 8) becomes the standard total / for the pixels over the lensed arcs.The rationale for the choice of weighting terms in equation ( 9) is as follows.The constraint on  comes from the radial stretching of the source, which is the next leading-order lensing constraint on the radial mass profile after the Einstein radius.This radial stretching is most effectively  is adjusted from the measured values to correspond to a circular aperture with radius  eff /2 (Jorgensen et al. 1995).As the distribution of key properties, such as  lens and  eff / E , are different among the samples, selection effects may potentially affect the accuracy of our analysis.Therefore, in our hierarchical analysis in Section 5.3, we consider the dependence of the potential deviation from the power law on  lens and  eff / E .The kernel density estimates of the 1D distributions (in the panels along the diagonal) are obtained using only the point estimates of the corresponding parameters.
constrained by the differential radial thickness of the arcs and then by the tangential stretch for a symmetrical (e.g., elliptical) mass profile (Birrer 2021).Therefore, the radial and angular terms in the equation ( 9) provide more weights to pixels away from the Einstein radius (corresponding to more radial stretch) and away from the reference angle (corresponding to more tangential stretch), respectively.This definition for the lensing information I will be higher for a system with lensed arcs than a system with only point images with the same total /.The angular weight term is taken as a multiplicative factor on the fractional radial distance so that the total information is simply the total / for the case of a circular lens that creates a perfect Einstein ring from a point source.We obtain the values  ≈ 16 and  ≈ 0 by maximizing the anti-correlation between log   and log I for all the lensing systems in our sample.We maximized the anti-correlation since we expect lower  uncertainty for systems with higher I.We optimized  and  as real numbers but approximated the optimized values to the nearest integers.For the optimized values of  ≈ 16 and  ≈ 0, log I and log   have a Pearson correlation coefficient of  = −0.71,whereas the 'unweighted' total / of the arcs (i.e., corresponding to  = 0) leads to a lower anti-correlation with  = −0.65.We adopted the form of modelling systematics dependent on the lensing information I to be  ,sys =  max ,sys tanh (I scale /I) , (10) , respectively).For the three outlier systems in the comparison between this paper and Etherington et al. ( 2022) (marked with red points; SDSSJ0912+0029, SDSSJ1213+6708, and SDSSJ1218+0830), the lensed arcs have low / and the Einstein radii differ by more than 10 per cent.If we exclude these outliers, we find the rms difference between the different analyses to be <2 per cent.2022) is  = 0.17 ± 0.16.We reillustrate the left-hand panel in Figure 7 after adding the estimated systematic uncertainty levels to the error bars.
where we set  max ,sys = 1/3 and I scale is the scaling parameter that determines how fast  ,sys changes with I in our sample.The total uncertainty on  for an individual lens system is then given by for both our values and those from Etherington et al. (2022).We obtain I scale ≃ 290 by minimizing the following penalty function which computes the absolute difference of a reduced  2 quantity from 1 for  lens = 18 systems that are common between our sample and that from Etherington et al. (2022), excluding the three outliers marked in Figure 5.We illustrate the estimated modelling systematic uncertainty  ,sys as a function of I in Figure 7.In the hierarchical analysis done in Section 5, we considered the total uncertainty for all the  measurements after adding in quadrature the statistical uncertainties with the systematic ones, that is,  ,sys estimated using equation (10).Despite the differences observed between the studies for the individual lens model posteriors, the sample means from these studies agree well.Our estimated sample median of the logarithmic slope  = 2.04 ± 0.04 for the SLACS lenses is in agreement within 1 with those from Shajib et

Comparison with previous lensing-dynamics analyses
The SLACS and SL2S samples were previously modelled with the SIE model (Bolton et al. 2008;Auger et al. 2009;Brownstein et al. 2011;Sonnenfeld et al. 2013a).The total density slope of the profile was then obtained using the stellar velocity dispersion.Comparing the Einstein radius measurements constrained by our lensing-only analysis to those from these previous joint lensing-dynamics studies, we find a rms difference of 3.3 per cent for the SLACS lenses, 5.6 per cent for the SL2S lenses, and 5.5 per cent for BELLS lenses (Figure 8).
The rms difference between our measured Einstein radii and those in the literature is higher than the expected uncertainty of 2 per cent that is empirically estimated by Bolton et al. (2008) and Sonnenfeld et al. (2013a).We note that the Einstein radii obtained by the previous studies are from lens modelling only and do not depend on the stellar kinematics.This discrepancy thus can arise from the difference be- tween the SIE model adopted by Bolton et al. (2008) and Auger et al. (2009) and the PEMD adopted by us.Indeed, by comparing our Einstein radii values to those obtained by Shajib et al. (2021) using the same PEMD for a SLACS subsample, we find a rms difference of 1.5 per cent, consistent within the measurement uncertainty.Other modelling differences could also contribute to the discrepancy mentioned above.For example, previous lensing-dynamics analyses modelled and subtracted the deflector galaxy's light from the imaging data before lens modelling, where Bolton et al. (2008) and Sonnenfeld et al. (2013a) used the de Vaucouleurs' profile, and Brownstein et al.
(2011) used B-spline fits.In our analysis, however, we simultaneously fit the deflector light profile (using two Sérsic profiles) during the lens model optimization.However, the absolute differences of the sample mean between our Einstein radii and those from the literature (Bolton et al. 2008;Auger et al. 2009;Sonnenfeld et al. 2013a;Brownstein et al. 2011) are 0.01 arcsec, 0.01 arcsec, and 0.02 arcsec for the SLACS, SL2S, and BELLS samples, respectively.Therefore, the systematic bias in the Einstein radius at the sample level is less than 2 per cent between the SIE model and the PEMD.
When comparing our lensing-only analysis to previous joint lensing-dynamics analyses from Auger et al. (2010b, SLACS) and Sonnenfeld et al. (2013b, SL2S), we find the median absolute difference in the logarithmic slope  to be 11.3 per cent for the SLACS systems and 6.1 per cent for the SL2S systems (Figure 8).Alternatively, the rms difference in  is 16.8 per cent for the SLACS systems and 14 per cent for the SL2S systems.We find only a weak correlation between the  measurements, where the Pearson correlation coefficient is  = 0.28 ± 0.10 for the SLACS lenses and  = 0.26 ± 0.16 for the SL2S lenses.However, our sample mean of  agrees within 1 with those from previous lensing-dynamics analyses.For the SLACS lenses, the sample median for  from our analysis is 2.04 ± 0.04 and that from the previous lensing-dynamics analysis is 2.08 ± 0.03 (Auger et al. 2010a).Similarly for the SL2S lenses, the sample median for  from our analysis is 2.08 ± 0.06 and that from the previous lensing-dynamics analysis is 2.01 ± 0.04 (Sonnenfeld et al. 2013b).
As discussed in Section 2.3, Birrer et al. (2020) find a previously unaccounted systematic uncertainty of 6 per cent in the SDSS velocity dispersion measurements for the SLACS systems, which corresponds to an additional 12 per cent uncertainty on the measured  from their lensing-dynamics analysis.After accounting for this additional systematic uncertainty, we find for the SLACS systems Assuming a  2 distribution for Δ 2 ,SLACS with 28 degrees of freedom, we obtain a -value of 0.67.For the SL2S systems -where we did not include any additional systematic uncertainty to the measured  from the lensing-dynamics analysis -we find Δ 2 ,SL2S = 27.8.For a  2 distribution with 15 degrees of freedom, this Δ 2 ,SL2S provides a -value of 0.02.However, the Δ 2 ,SL2S quantity is dominated by SL2SJ0904−0059, which has an atypically low  LD = 1.48 ± 0.11.When we exclude this system as an outlier, then Δ 2 ,SL2S = 18.2 provides a -value of 0.2 for 14 degrees of freedom.Therefore, we cannot rule out the hypothesis that our lensing-only  measurements are consistent with those from the previous lensing-dynamics analyses for both the SLACS and SL2S systems, given the uncertainties.Thus, if the power-law model assumption is correct, then the expected strong correlation in the individual  measurements between our lensing-only and previous lensing-dynamics analyses could have been smeared by the large uncertainties.We can not, therefore, conclude on a deviation from the power law based on our finding of the weak correlation.

DEVIATION OF THE MASS PROFILE FROM A POWER LAW
In this section, we combine the lens model posteriors with the velocity dispersions measured from single-aperture spectra to constrain any deviation from the power law in the total density profile.We parametrize the deviation from the power-law form using the masssheet transformation (MST).We briefly introduce the MST in Section 5.1, and describe the formalism and assumptions used to predict the velocity dispersion from our lens models in Section 5.2.Then, in Section 5.3, we describe the hierarchical Bayesian analysis that we used to constrain the deviation from the power law at the population level.
1  7, the grey error bars show the total uncertainty by adding the statistical and systematic uncertainties in quadrature.The 12 per cent systematic uncertainty added to the  LD measurements from the SLACS is based on the 6 per cent systematic uncertainty on SDSS velocity dispersion measurements reported by Birrer et al. (2020).The rms difference between our  and  E measurements with those from the literature is 15.9 per cent and 4.5 per cent, respectively.Using the -values of the corresponding  2 statistics, we cannot rule out the consistency between our  measurements and those from the previous lensing-dynamics analysis, given the uncertainties.We also do not detect any significant dependence of the  E difference on the mass axis ratio  m .

Mass-sheet transformation
The MST (Falco et al. 1985) is a multiplicative transformation of the lens equation (i.e., equation 1).This transformation modifies the lens equation as The convergence field is then transformed as The MST has the property where image positions remain invariant while the source position scales with  such that As the MST changes the convergence field and thus the mass distribution of the lensing galaxy, the model-predicted stellar velocity dispersion changes with the transformation as well.The line-of-sight stellar velocity dispersion of the lensing galaxy integrated within an aperture scales approximately as 16) (Birrer et al. 2020).We can therefore estimate the value of the MST parameter  that transforms the adopted power-law model into the 'true' mass distribution by comparing the measured velocity dispersion  ap,meas and the predicted velocity dispersion  ap,PL by the power-law model as In addition to lensing caused by the convergence from the lensing galaxy (denoted with  gal ), the overdensity or underdensity of the structures along the line of sight with respect to the average matter density of the Universe can also produce additional magnification or demagnification.We denote the 'external' convergence from these line-of-sight structures with  ext .Therefore, the effective total lens mass distribution of the system, , is then given by The MST can be either associated with transforming the mass profile of the main lensing galaxy with the 'internal MST' parameter  int , ignoring the line-of-sight structures in the lens modelling, or a combination of both.The total transformation affecting the stellar kinematics is thus given by  total = (1 −  ext ) int (19) (Shajib et al. 2023).To break this degeneracy, we adopted the estimated  ext from Birrer et al. (2020) for the SLACS systems and from Wells et al. (in preparation) for the SL2S systems.For systems without individual  ext measurements in these two samples, we use the overall sample distribution of  ext .The  ext estimates for the BELLS systems are not available from the literature.We leave their  ext estimation for a future study and exclude this sample for the following analysis in this paper.

Modelling of the stellar dynamics
The aperture-integrated, line-of-sight velocity dispersion  ap provides a lensing-independent method to estimate the mass profile of the lensing galaxy.The stellar velocity dispersion can be obtained by solving the Jeans equation, which describes the phase space distribution of a galaxy with a 3D stellar density profile  * () located within gravitational potential Φ().In the case of spherical symmetry, the velocity dispersion can be decomposed into its radial component  r () and its tangential component  t ().The stellar anisotropy of the system is parametrized as For a relaxed system, the spherical Jeans equation can be expressed as (Binney & Tremaine 2008;Birrer et al. 2020).We can solve this spherical Jeans equation to obtain where  () is the 3D total enclosed mass profile.The line-of-sight velocity dispersion is then given by where Σ * () is the 2D projected stellar density profile,  is the 2D projected radius (Binney & Mamon 1982).The luminosity-weighted line-of-sight velocity dispersion  ap within an aperture is given by where * S denotes PSF convolution with the seeing (Treu & Koopmans 2004;Suyu et al. 2010).
We chose a spatially constant anisotropy profile as our baseline model, which is consistent with dynamical observables obtained from both long-slit and integral field spectra of local elliptical galaxies (Gerhard et al. 2001;Cappellari et al. 2007;Cappellari 2008).In this model, we treat  ani as a free parameter.We also test with an alternative choice of the anisotropy profile given by the Osipkov-Merritt parameterization where  ani is a unitless scale factor (Osipkov 1979;Merritt 1985).
To compute the model-predicted stellar velocity dispersion, the 3D mass profile () was obtained corresponding to the 2D convergence profile in equation ( 6).Assuming a constant mass-to-light ratio, we obtained the 2D projected stellar density profile Σ * () from the surface brightness profile  () constrained from large (∼16×16 square-arcsec) cutouts of the HST images.We fit the surface brightness profile with either a single or double Sérsic profile on the HST images.These cutouts are larger than those used in lens modelling to more robustly estimate the surface brightness profile of the lensing galaxy up to its full extent.The assumed value of the mass-to-light ratio does not matter since it cancels out in equation ( 24).To obtain the 3D stellar mass density  * , we decomposed the projected stellar density profile Σ * () into multi-Gaussian expansion (MGE; Emsellem et al. 1994;Cappellari 2002) and then deprojected the Gaussian profiles into 3D.

Bayesian hierarchical framework
We performed a Bayesian hierarchical analysis to constrain the mean deviation from the power law at the population level while also adopting a prior on the anisotropy profile parameter at the same population level.We performed this using hierArc6 , a software program for Bayesian hierarchical analysis of strong lensing systems (Birrer et al. 2020).
Here, we provide a brief description of the hierarchical framework.Let   be the data for the  th individual lens and  ≡ {  }  be the dataset containing all the data from  number of lenses.For this analysis, the lensing data for each dataset   is summarized with the lens model posteriors for the logarithmic slope   and the Einstein radius  E, .The dataset   additionally includes the measured velocity dispersion  ap,meas, , the estimated external convergence  ext, , and the single or double Sérsic profile parameters fitted to the lens galaxy's surface brightness from Section 5.2.We add a 2 per cent systematic uncertainty (Shajib et al. 2021) in quadrature to the measured  eff uncertainty provided in Tables 3 and 4. We find from a test that our final results are robust against a more conservative choice of a uniform uncertainty level of 5 per cent in  eff .For the  th lens system, the individual-level parameter set   includes the internal MST parameter  int, and the stellar anisotropy parameter - ani, for the constant anisotropy model or  ani, for the Osipkov-Merritt model.In the hierarchical framework, the population distribution of individual-level parameters   are described with population-level parameters Ξ. Bayes' theorem gives the posterior (Ξ | ) of the population-level parameters as where ( | Ξ) is the likelihood function and (Ξ) is the prior of the population-level parameters.While (  |   ) is the likelihood function for the  th lens system and (  | Ξ) describes the distribution of individual lens parameters given the population-level parameters.Now, we explain all the population-level model parameters Ξ used in this analysis.Following Birrer et al. (2020), we described the population distribution of the internal MST parameter  int using a Gaussian distribution with mean   int and scatter   int .We incorporated linear dependencies of   int on the redshift  and the scale ratio  eff / E as Here,  ref is the reference value of   int for  = 0.3 and  eff / E = 1, and   int and   int are linear slope parameters.
To describe the population distribution of  ani for the constant anisotropy model, we used a Gaussian distribution with mean   ani and scatter   ani .For the Osipkov-Merritt model, we described the population distribution of the anisotropy scale factor  ani as a Gaussian distribution with mean   ani and scatter   ani   ani .
As mentioned in Section 2.3, Birrer et al. (2020) find that the uncertainty in the velocity dispersion measurements for the SLACS lenses from the SDSS spectra was underestimated.Therefore, we included a systematic uncertainty term set by the fractional uncertainty  SDSS ,sys .The total uncertainty   for the measured velocity dispersion  ap of an individual SLACS lens is therefore given by where  ,stat is the estimated statistical uncertainty for the velocity dispersion measurement.Thus, for the constant anisotropy model, the population-level parameters sampled by hierArc are For the Osipkov-Merritt model, the population-level anisotropy parameters are replaced with   ani and   ani .We tabulate the adopted priors on these population-level parameters in Table 6.
The posterior distribution of the population-level parameters obtained using hierArc is shown in Figure 9, and the corresponding point estimates are provided in Table 7. Figure 10 shows the goodness of fit for the SLACS and SL2S model-predicted velocity dispersions using the maximum likelihood estimators of the population-level parameters.The model-predicted velocity dispersions are derived for both the power-law profile and the mass-sheet transformed profile assuming the constant or the Osipkov-Merrit anisotropy profile.
For the constant anisotropy model, the mean reference MST parameter is  ref −0.05 deviates from the power-law profile with a higher confidence level of 3.2.However,  ref  int is consistent within 0.9 between the constant anisotropy and Osipkov-Merritt models.Although the likelihood ratio favours the Osipkov-Merritt model by a factor of ∼250 (i.e., Δ log L ≈ 5.5) 7 , this is due to the Osipkov-Merritt model being more capable of fitting the population scatter with a Gaussian distribution.Single-aperture velocity dispersionsas used in this analysis -have no information to differentiate between different choices of anisotropy model.Therefore, spatially resolved kinematics are necessary to differentiate between the two anisotropy models.In this paper, we adopted the constant anisotropy model as the baseline, given its consistency with dynamical observables of local massive ellipticals, whereas the Osipkov-Merritt model does not allow tangential orbits that are often observed in the local ellipticals (e.g., Gerhard et al. 2001;Cappellari et al. 2007;Cappellari 2008).
Accounting for the distribution of redshift and  eff / E in the cor- 7 Here, The likelihood ratio is equivalent to the Bayesian information criterion as both models have the same number of free parameters.We extracted the maximum log-likelihood value for each model from the MCMC chain.
We checked that the second-highest log-likelihood value in each chain is much closer to the maximum log-likelihood value relative to the difference in the log-likelihood between the two models.Thus, we ensured that the maxima of the likelihood functions were well sampled, giving robustness to our numerical estimation of the maximum-likelihood ratio.
responding samples, the SLACS galaxies are consistent with the power-law profile (i.e.,  ref

𝜆 int
= 1) within 1.1, and so are the SL2S galaxies within 0.8, for the constant anisotropy model.For the Osipkov-Merritt model, the SLACS and SL2S galaxies are consistent within 2.8 and 2.1, respectively, with the power-law profile.

DISCUSSION
In this section, we first discuss the impact of the anisotropy profile choice on joint lensing-dynamics analysis in Section 6.1.Then, in Section 6.2, we discuss the effect of the internal MST on the local logarithmic slope of the density profile.Finally, we discuss the redshift evolution of the logarithmic slope obtained from the different modelling methods in Section 6.3.

Impact of the anisotropy profile on joint lensing-dynamics analysis
Our result highlights the importance of the stellar anisotropy profile choice in the dynamical modelling when constraining the mean deviation from the power law.By comparing our baseline choice of the constant anisotropy model with the Osipkov-Merritt model, we find that the latter provides a higher significant deviation (3.2) from the power-law profile than the former one (1.1),although both choices are consistent with each other at 0.9.Spatially resolved kinematics of lens galaxies are essential to differentiate between these anisotropy profile choices to mitigate the associated systematic uncertainty in constraining any potential deviation from the power law, and we leave this investigation for a future study.
In this paper, we directly constrained the deviation from the power law, or the lack thereof, through a joint lensing-dynamics analysis that has several improvements in the dynamical modelling over previous such analyses (e.g., Auger et al. 2010b;Sonnenfeld et al. 2013a).First, we used a stellar anisotropy model with one degree of freedom instead of the isotropic model that has no degree of freedom.The isotropic model is a special case of our constant anisotropy model with  ani = 0, which is 1.2 away from our population mean for for the stellar anisotropy.The population distribution of the corresponding anisotropy parameter was modelled using a Gaussian distribution with mean   ani or   ani and scatter of   ani   ani or   ani for the constant or Osipkov-Merritt model, respectively.We accounted for a systematic uncertainty in the SDSS velocity dispersion measurements for the SLACS lenses using the factor  SDSS ,sys .The point estimates of these parameters are provided in Table 7.The mean reference parameter  ref  int is consistent between the two anisotropy models within 0.9.The constant anisotropy model is consistent with the power-law profile (i.e.,  ref  int = 1) within 0.9.However, the Osipkov-Merritt model deviates from the power-law model with a larger confidence level of 3.2.We adopted the constant anisotropy model as our baseline choice, as this model is consistent with the dynamical observables in local elliptical galaxies (e.g., Gerhard et al. 2001;Cappellari 2008).Our result highlights the impact of the stellar anisotropy model choice in constraining the deviation from the power-law profile (e.g., as suggested by Etherington et al. 2023b), which can be informed by spatially resolved kinematics in the future.
ani .Second, we used in the dynamical modelling a more accurate surface brightness profile obtained from the full extent of the lens galaxy instead of the Hernquist or De Vaucouleurs's profile (Hernquist 1990;de Vaucouleurs 1948), which can often be inadequate to describe lens galaxy's surface brightness over its full extent and thus potentially induce systematics.Our state-of-the-art analysis provides a direct approach to quantify any deviation from the power law, which is an improvement over previous studies investigating such deviation based on inconsistent correlations between lensing-only and lensingdynamics analyses within the scope of the power-law model.We find consistency with the power-law profile at 1.1 confidence level for our baseline choice of the constant anisotropy model.Therefore, = 1 for the pure power-law case -given the individual-level parameters ,  E , and  eff for that system.Thus, the uncertainty term  ,model incorporates the uncertainties from ,  E , and  eff .The right-hand panel shows the kernel density estimation plot of the normalised difference between  model and  meas .With the constant anisotropy profile, the best-fit mass-sheet transformed profiles predict similar velocity dispersion values as the power-law model, as expected from the consistency of  ref  int with the case with no MST (i.e.,  int = 1) within 1 ( ref  int = 0.91 +0.10 −0.09 ).However, using the Osipkov-Merritt anisotropy profile, we find that for a majority of the systems, the model-predicted velocity dispersions are higher than the measured values, and mass-sheet transformed profile is required to fit the measured values better ( ref the observed discrepancy by Etherington et al. (2023b) between the correlation of the mass surface density with lensing-only  measurements and that with previous lensing-dynamics  measurements will likely be alleviated with a more flexible and robust model for stellar dynamics, without requiring a significant departure from the power law.The discrepancy may also be alleviated to some further extent if additional modelling systematic uncertainties on  are accounted for, as done in this analysis (Section 4.1).

Impact of the internal MST on the logarithmic slope
The internal MST changes the shape of the physical mass distribution in the deflector galaxy.Thus, it impacts the local logarithmic slope.
In this subsection, we illustrate, with an example, the change in the logarithmic slope when an internal MST is performed on the mass profile.The internal MST does not add mass at a distance far away ( ≫  eff ) from the deflector galaxy.Therefore, to ensure that the transformed mass profile goes to zero at a large radius, that is, to have lim →∞  ′ () = 0, we use an 'approximate' MST given by Blum et al. (2020).Here,  s is the scale radius where the mass sheet truncates.We can deproject this transformed convergence profile from 2D into the 3D density profile To illustrate an example, we apply the internal MST on a mass profile with  = 2.05 and  E = 1.2 arcsec.These fiducial values are adopted to match the median values from our sample of SLACS and SL2S lenses.We ensure that the approximate MST recreates the same lensing effect observed in the imaging data as the pure MST by setting  s = 10 E as the lensing image is only sensitive to the inner regions ( ≲ 2 E ) of the mass profile of the lensing galaxy (Schneider & Sluse 2013;Birrer et al. 2020).For  int = 0.91, the mass profile slope at the Einstein radius decreases from ( E ) = 2.05 to  MST ( E ) = 1.99, which represents a deviation by 3 per cent.The change in the shape of the total mass profile -with a decrease in the density inside  E and an increase outside -due to the MST is also illustrated in Figure 11.
We note that different choices of  s would lead to slightly different shapes of the mass-sheet transformed profiles, as shown in Figure 11.Thus, we do not have sufficient information to obtain the complete shape of the density profile at all radii.Additional information, such as a strong prior on  s , or performing the joint lensing-dynamics analysis with a composite model of dark matter and baryonic distributions with empirical priors on their shapes (e.g., Shajib et al. 2021), would be necessary to obtain an accurate density profile shape at all radii.These investigations are left for future papers in the Project Dinos series.

Evolution of total density profile in elliptical galaxies
Cosmological hydrodynamical simulations, such as the Magneticum and IllustrisTNG, predict the 'global' logarithmic slope averaged over a radius range containing  eff .This global logarithmic slope decreases with decreasing redshift at  < 1, as mainly gas-poor mergers make the density profile in massive ellipticals shallower over time (Remus et al. 2017;Wang et al. 2019;Wang et al. 2020).Although this prediction should be tested using observations in principle, we refrain from doing so for several reasons.
First, we do not detect a significant evolution in the mean MST parameter   int , as the large uncertainty on the   int parameter makes The 3D density profile is normalised such that  PL (0.1 eff ) = 1, where the subscript 'PL' denotes the power-law case.Different colours represent different core radius  s used for the approximate MST (equation 29).The initial mass profile of the system is a power-law profile with  = 2.05 and  E = 1.2 arcsec.The vertical grey band represents the 1 distribution of  E / eff for our lens sample, illustrating the scale range containing the local lensing information.At  =  E with  s = 10  E , the local logarithmic slope deviates only by ∼3 ± 3 per cent from the power-law profile, and the transformed profile is consistent with the power-law profile within 1 over the illustrated range of radius.
the redshift dependence consistent with zero (Table 7).This large   int uncertainty can be attributed to the large systematic uncertainties added to the measured logarithmic slopes , to a larger extent for the SL2S systems (Table 4) owing to their short-exposure images having low / on the lensed arcs.A future study will model the SL2S lenses with new, deeper HST imaging to constrain the redshift evolution more tightly (Sheu et al., in preparation).
Second, we cannot uniquely constrain a global (i.e., radiallyaveraged) logarithmic slope using our analysis framework to make a fair comparison with the predictions from simulations.Although several previous studies have compared the logarithmic slope with those from the simulations (e.g., Bolton et al. 2012;Sonnenfeld et al. 2013b), those comparisons were performed with the assumption of a power-law mass model, where the local logarithmic slope matches with the global logarithmic slope by definition.A lensing-only analysis provides the local logarithmic slope (Birrer & Treu 2019).Thus, they cannot be compared to the global logarithmic slope obtained from simulation.Although joint lensing-dynamics analysis can be used to obtain the global logarithmic slope, the adopted model for such analysis requires the flexibility to produce the true mass distribution in galaxies.Our framework of using the MST to constrain any deviation from the power-law lens model only applies a global scaling on the power-law distribution when fitting dynamical observables without actually changing the shape of the mass profile (Birrer et al. 2020).However, applying an approximate internal MST to the power-law lens model to uniquely recover the true mass distribution is also not possible due to the impact of  s choice on the transformed profile, as discussed in Section 6.2.Future studies will use an empirically motivated composite mass model to constrain the actual mass profile shape from joint lensing-dynamics analysis to perform a fair test of the simulation predictions (Sheu et al., in preparation).
Third, the logarithmic slope can depend on other galaxy properties, such as the velocity dispersion and the stellar mass density (Auger et al. 2010a;Sonnenfeld et al. 2013a).Thus, if these properties are not taken as control variables, galaxies at the different redshifts within the same sample (or, 'super'-sample as assembled in this study) may be subject to different selection functions.Such selection functions, if unaccounted for, can lead to erroneous detection of a trend in the logarithmic slope across redshift (see Sonnenfeld et al. 2013b).
Fourth, the logarithmic slope is only one structural parameter among several other structural properties predicted or assumed by the simulations.For example, the dark matter fraction within a given aperture and the stellar initial mass function (IMF).To comprehensively test the predictions from the simulations to learn about the baryonic astrophysics that has shaped the massive ellipticals, all these structural parameters must be compared simultaneously, not just the logarithmic slope.However, to obtain these additional properties from observations, a composite mass model will be required to individually describe the dark matter and baryonic distributions.
Based on the above discussion, we recommend future studies testing the predictions from simulations with strong lensing to adopt these elements: • perform joint lensing-dynamics analysis with a composite mass profile, where dark matter and baryonic matter distribution in the deflectors are individually described, • control for potential selection effects or covariances with other galaxy properties, such as the central velocity dispersion or the stellar surface density, and • test simulation predictions for multiple structure properties simultaneously, for example, the inner dark matter fraction, the stellar IMF, etc., in addition to the global logarithmic slope.

CONCLUSION
In this paper, we present uniform lens models of 77 galaxy-galaxy lens systems.These 77 systems were collated from the SLACS, SL2S, and BELLS samples.To uniformly model these systems, we used the automated modelling pipeline dolphin, which uses lenstronomy as the modelling engine.This dolphin pipeline successfully modelled these 77 systems out of 103 attempted ones, which is a high success rate for an automated lens modelling pipeline.The remaining systems could be modelled by carefully tuning the models by hand, as previous studies did.However, we excluded those 'failure cases' from modelling in this paper, as the minor increase in the statistical power from a larger sample of 103 systems is not worth the additional requirement for investigator time to hand-tune the lens models on a case-by-case basis for the remaining systems.
We combined our lens model posteriors with models to constrain any potential deviation from the power-law profile in the mass distribution.We parametrized the deviation from the power-law mass model using the internal MST parameter  int , where  int = 1 implies no deviation.We performed a hierarchical Bayesian analysis using the SLACS and SL2S systems from our modelled sample.We incorporated the lensing effect of the line-of-sight structures in this analysis using estimates from the literature, a first for analyses on galaxy-galaxy lenses.We excluded the BELLS systems in this step, as the estimates for their line-of-sight lensing effect are currently unavailable.The key results from this study are as follows: • We constrain the sample mean of the internal MST parameter with dependence on redshift and Einstein radius as   int = 0.91 +0.10 −0.09 − 0.12 +0.31 −0.26 ( − 0.3) − 0.04 +0.04 −0.04 ( eff / E − 1), for the constant anisotropy model.Therefore, at  = 0.3 and  eff / E = 1, the population mean for massive ellipticals is consistent with the power-law profile (i.e.,   int = 1) within 0.9.
• For an alternative choice of the stellar anisotropy profile using the Osipkov-Merrit parametrization, the deviation from the power law has a higher significance of 3.2 (i.e.,  ref  int = 0.81 +0.06  −0.05 ).However,  ref  int is consistent between the two models within 0.9.This result highlights the potential impact of the stellar anisotropy model in constraining any deviation from the power-law profile by combining lensing and dynamical observables.The constant anisotropy profile is consistent with the dynamical observables in local massive ellipticals, therefore, we chose this model as our baseline.However, spatially resolved kinematics of lens galaxies will be necessary to differentiate between these stellar anisotropy models.
• We find our logarithmic slope measurements have a median absolute deviation of 8-10 per cent with those from Etherington et al. (2022), who used the same PEMD for the lens model but a different source reconstruction algorithm and software package.This difference is larger than the median 3 per cent statistical uncertainty obtained from our posteriors and points to potential residual systematics in the lens modelling.We find this difference to be larger for lens systems with lower-/ arcs.Therefore, we estimated the modelling systematic uncertainty dependent on a 'weighted' total / quantity of the arc (Section 4.1).
• We find the Einstein radii constrained by different studies using the PEMD for the lens model agree very well within 2 per cent, which is the expected level of consistency empirically estimated by Bolton et al. (2008).However, the Einstein radii given by the studies using the SIE mass model can differ by 3-6 per cent, falling beyond the level of expected consistency.
This study presents the largest sample to date of galaxy-galaxy lenses with uniform, power-law models.However, this decade will see larger samples with similarly exquisite imaging data from the HST and velocity dispersion measured using large, ground-based telescopes (e.g., Tran et al. 2022).The sizes of recently discovered samples from the Dark Energy Survey, the DECam Local Volume Exploration Survey, and the Dark Energy Spectroscopic Survey are O (10 3 ) (Jacobs et al. 2019a,b;Huang et al. 2020;Storfer et al. 2022;Zaborowski et al. 2023).In the near future, the sample size of newly discovered systems will rapidly grow to O (10 4 )-O (10 5 ) thanks to current or future surveys carried out by the Vera Rubin Observatory, the Nancy Grace Roman Space Telescope, and the Euclid (Oguri & Marshall 2010;Collett 2015).As employed in this paper, an automated framework will be essential to tackle the computational challenge these very large samples will bring forth.Additionally, these large samples will provide a much higher statistical constraining power to precisely track the evolutionary trends in the structural properties of massive ellipticals and shed light on the precise nature of the baryonic astrophysics at play through comparison with simulations.
Table A1.HST imaging filters and other settings used in modelling the SLACS and SL2S systems.For each system, we manually optimised the shapelet order  max for the source-light model and the radius  inner of the inner mask used to avoid residuals in the deflector galaxy's centre.The 'foreground mask' column indicates whether we included an additional mask for the line-of-sight contaminants in the associated filter.Table A2.HST imaging filters and other settings used in modelling the BELLS systems.For each system, we manually optimised the shapelet order  max for the source-light model and the radius  inner of the inner mask used to avoid residuals in the deflector galaxy's centre.The 'foreground mask' column indicates whether we included an additional mask for the line-of-sight contaminants in the associated filter.

Lens system
Survey Filter Shapelets  max  inner (arcsec) Foreground mask The lens models for the other imaging filters are of similar quality (as shown in Figure 3).The lens models for the remaining SL2S systems are shown in Figures B6 and B7.The lens systems SL2SJ0226-0406 (in Figure B5) and SL2SJ1427+5516 (in Figure B7) have a disc-like morphology and are not included in the analysis in Section 5. 1"

Figure 1 .
Figure 1.An example demonstration of deepCR in masking cosmic rays for a single-exposure image (in F435W band) of SDSSJ0819+4534.The left-hand panel shows the unmodified HST image, and the right-hand panel shows the HST image with the masked pixels shown in purple

Figure 2 .
Figure 2. Best-fitting lens model for SDSSJ1627−0053.The columns (from left to right) represent the observed HST image, the model-based reconstruction, the normalised residual, and the reconstructed source light distribution.The top row corresponds to the F555W filter, and the bottom row corresponds to the F814W filter.The flux values are in units of e − s −1 .All the scale bars are in units of arcsec.

Figure 3 .
Figure 3. Best-fitting lens model for SL2SJ1401+5544.The columns (from left to right) represent the observed HST image, the model-based reconstruction, the normalised residual, and the reconstructed source light distribution.The top row corresponds to the F475X filter, and the bottom row corresponds to the F600LP filter.The flux values are in units of e − s −1 .All the scale bars are in units of arcsec.

Figure 6 .
Figure 6.Comparison of the measured logarithmic slopes between our analysis (denoted by  ours ) and other lensing-only analyses from Shajib et al. (2021) and Etherington et al. (2022) (denoted by  Shajib+21 and  Etherington+22 , respectively).The three outlier systems from Figure5are similarly marked with red.We find only a weak correlation between the measured logarithmic slope  from different analyses.For example, the Pearson correlation coefficient between the  values from our analysis andEtherington et al. (2022) is  = 0.17 ± 0.16.We reillustrate the left-hand panel in Figure7after adding the estimated systematic uncertainty levels to the error bars.

Figure 7 .
Figure 7.Estimated modelling systematic uncertainty on the logarithmic slope .The left-hand panel shows the estimated systematic uncertainty (black line) as a function of the lensing information I defined in equation (8).The modelling systematic is estimated using the difference in the  values (blue points with blue error bars) of individual systems between the models from this analysis and those fromEtherington et al. (2022), excluding the three outliers marked in Figures5 and 6.The grey error bars show the total uncertainty by adding the estimated systematic uncertainty with the statistical one in quadrature.The middle panel shows the total uncertainties (grey error bars) for the measured  values in our full sample (blue points, with the blue error bars showing statistical uncertainties).The lens systems with higher lensing information I have less systematic uncertainty, which was estimated using the model shown by the black line in the left-hand panel.The right-hand panel shows a comparison of the measured logarithmic slopes between our analysisEtherington et al. (2022) with the systematic uncertainty of the lenses added to the grey error bars.
09 , consistent with the power-law profile (i.e.,  ref  int = 1) within 0.9.The point estimates of the slope parameters   int = −0.12+0.31 −0.26 and   int = −0.04+0.04 −0.04 are consistent with having no dependence on the redshift and  eff / E , respectively.For the Oskipkov-Merritt model, the mean reference MST parameter  ref  int = 0.81 +0.06

Figure 9 .
Figure 9.Posteriors of the population-level parameters from our hierarchical Bayesian analysis.The population distribution of the internal MST parameter  int was described with a Gaussian distribution using mean  int (,  eff / E ) =  ref  int +   int ( − 0.3) +   int ( eff / E − 1) and scatter   int .Here, we compare the posterior between the spatially constant anisotropy model (orange, parametrized with  ani ) and the Osipkov-Merritt model (green, parametrized with  ani ) for the stellar anisotropy.The population distribution of the corresponding anisotropy parameter was modelled using a Gaussian distribution with mean   ani or   ani and scatter of   ani   ani or   ani for the constant or Osipkov-Merritt model, respectively.We accounted for a systematic uncertainty in the SDSS

Figure 11 .
Figure11.The deviation of the 3D density profile (left-hand panel) and the local logarithmic slope  ( ) ≡  (log )/ (log  ) (right-hand panel) for the MST of a power-law profile.The 3D density profile is normalised such that  PL (0.1 eff ) = 1, where the subscript 'PL' denotes the power-law case.Different colours represent different core radius  s used for the approximate MST (equation 29).The initial mass profile of the system is a power-law profile with  = 2.05 and  E = 1.2 arcsec.The vertical grey band represents the 1 distribution of  E / eff for our lens sample, illustrating the scale range containing the local lensing information.At  =  E with  s = 10  E , the local logarithmic slope deviates only by ∼3 ± 3 per cent from the power-law profile, and the transformed profile is consistent with the power-law profile within 1 over the illustrated range of radius.

Figure B1 .
Figure B1.Best-fitting lens models for the first 10 of the 34 SLACS lenses.The first column shows the HST images of the lenses in the visible band (F555W or F606W).The second column shows the model-based reconstruction and the third column shows the normalised residual.We visualise the lensed arcs in the fourth column by subtracting the deflector light profile from the HST image.The fifth column shows the reconstructed source light and the sixth column shows the magnification map.The lens models for the other imaging filters are of similar quality (as shown in Figure 2).The lens models for the remaining SLACS systems are shown in Figures B2,B3 and B4.

Figure B2 .
Figure B2.Lens images and models for the next 10 out of the 34 SLACS lenses.This continues from Figure B1. 1

Figure B3 .
Figure B3.Lens images and models for another 10 out of the 34 SLACS lenses.This continues from Figure B1 and B2.

Figure B4 .Figure B5 .
Figure B4.Lens images and models for the last 4 out of the 34 SLACS lenses.This continues from Figure B1, B2 and B3.

Figure B6 .
Figure B6.Lens images and models for the next 10 out of the 24 SL2S lenses.This continues from Figure B5.

Figure B7 .Figure B8 .
Figure B7.Lens images and models for the last 4 of the 24 SL2S lenses SLACS.This continues from Figure B5 and B6.

Figure B9 .
Figure B9.Lens images and models for the next 9 out of the 19 BELLS lenses.This continues from Figure B8.

Table 1 .
Description of the HST programs from which we assemble our archival lens sample.Programs marked with an asterisk are snapshot programs that acquired only single-exposure images.

Table 2 .
Number of systems with available HST imaging, systems we selected to attempt modelling, and systems successfully modelled for the lens samples considered in this paper.

Table 4 .
Lens model parameters for the SL2S systems.See the caption of Table3for the column descriptions.

Table 5 .
Lens models parameters for the BELLS systems.See the caption of Table3for the column descriptions.Due to the lack of -band imaging for the BELLS systems, the lens light parameters were obtained from HST images in the -band (F814W).
Distributions of deflector redshift  lens , central velocity dispersion   eff /2 , Einstein radius  E , logarithmic slope , residual shear magnitude  shear , and the ratio of half-light radius to Einstein radius  eff / E for the SLACS (purple), SL2S (orange), and BELLS (green) systems.The velocity dispersion   eff /2 Comparison of the measured Einstein radii between our analysis (denoted by  ours E ) and other lensing-only analyses from Shajib et al. (2021) and Etherington et al. (2022) (denoted by (Auger et al. 2009;Brownstein et al. 2011;Sonnenfeld et al. 2013and the Einstein radius  E (right-hand panel) between our analysis, which used the PEMD, and previous studies that used the SIE model to obtain  E and joint lensing-dynamics (LD) analysis to obtain (Auger et al. 2009;Brownstein et al. 2011;Sonnenfeld et al. 2013a).The SLACS, SL2S, and BELLS systems are represented with purple, orange, and green colours, respectively.As with Figure

Table 6 .
Description and adopted priors of the population-level parameters in our hierarchical Bayesian analysis.Population mean of internal MST parameter at  = 0.3 and  eff / E = 1 (see equation 27)  ref   int Linear dependency of  int on redshift  (see equation 27)   int ∼ U (−1, 1)   int Linear dependency of  int on  eff / E (see equation 27)   int ∼ U (−1, 1)   int 1 Gaussian scatter in  int   int ∼ Loguniform(0, 0.5)   ani Population mean of the anisotropy parameter  ani in the constant anisotropy model   ani ∼ U (−0.49, 1)   ani 1 Gaussian scatter of  ani in the constant anisotropy model   ani ∼ Loguniform(0.01,0.50)   ani Population mean of anisotropy scale factor in the Osipkov-Merritt model (equation 25)   ani ∼ Loguniform(0.1,5)   ani 1 Gaussian scatter of  ani in the Osipkov-Merritt model is given by   ani   ani   ani ∼ Loguniform(0, 1)

Table 7 .
Point estimates of the population-level parameters from our hierarchical Bayesian analysis.The parameters are described in Table6.
Goodness of fit of the dynamical observables between the power-law model and the power-law model with the MST.The goodness of fit is illustrated with the normalised difference between the model-predicted velocity dispersion  model and the measured velocity dispersion  meas , which is computed as(  model −  meas )/ √︃  2 ,model +  2 ,meas .Here,  ,model and  ,meas are the uncertainties for  model and  meas , respectively.The model-predicted velocity dispersion  model for each individual system was computed using the most likely values of the population-level parameters -except for taking  ref