Two-dimensional kinematics and dynamical modelling of the ‘Jackpot’ gravitational lens from deep MUSE observations

We present results from the first spatially resolved kinematic and dynamical modelling analysis of the unique SDSSJ0946+1006 (‘Jackpot’) triple-source lens system, where a single massive foreground 𝑧 = 0 . 222 galaxy multiple-images three background sources at different redshifts. Deep IFU spectroscopic data were obtained using the MUSE instrument on the VLT, which, compared to previous single-slit observations, provides full azimuthal area coverage, high sensitivity (5 hour integration) and high angular resolution (0.5 arcsec FWHM). To account for the strong continuum contributions from the 𝑧 = 0 . 609 source, a multiple-component stellar template fitting technique is adopted to fit to the spectra of both the lens galaxy and the bright lensed background arc simultaneously. Through this, we robustly measure the first and second moments of the two-dimensional stellar kinematics out to about 10 kpc from the centre of the lens, as well as resolving the inner profile inwards to ∼ 1 kpc. The two-dimensional kinematic maps show a steep velocity dispersion gradient and a clear rotational component. We constrain the characteristic properties of the stellar and dark matter (DM) mass components with a sufficiently flexible parameterised dynamical model and an imposed lensing mass and find a DM density slope of 𝛾 = 1 . 73 + 0 . 17 − 0 . 26 , i.e. significantly steeper than an unmodified NFW profile ( 𝛾 = 1) and consistent with a contracted DM halo. Our fitted models have a lensing-equivalent density slope of 𝜂 = 0 . 96 ± 0 . 02, and thus we confirm most pure lensing results in finding a near isothermal profile for this galaxy.


INTRODUCTION
The structure of the most massive elliptical galaxies provides a window onto the early history of galaxy formation.The currently-popular "two-phase" formation history of such galaxies (e.g.Oser et al. 2010) -with an early starburst-driven phase and subsequent accretion through mergers -offers a vital framework for understanding the assembly and evolution of galaxies, and the distribution of stars and dark matter (DM) within them.
Lacking ordered gas dynamical probes, stellar kinematics represent the main method to study the mass distribution in massive elliptical galaxies (Cappellari 2016).However, in the absence of detailed observations from high resolution data, analyses at any significant redshift are limited to simple dynamical Jeans type models (e.g.Cappellari 2008), as opposed to more general triaxial orbit based models (e.g.Schwarzschild 1979;van den Bosch et al. 2008).Additionally, due to the dependence on luminous mass tracers, of which there is an absence of at large radii, and the sensitivity to many model dependant degeneracies (e.g.Bender et al. 1994;Carollo et al. 1995;Gerhard & Binney 1996;Romanowsky & Kochanek 1997;van den Bosch 1997;Weĳmans et al. 2009;Oldham & Auger 2016), galaxy dynamics alone do not provide sufficient information to disentangle the stellar and dark mass contributions, and hence place robust constraints on the physics of the DM particle.
★ Contact e-mail:hannah.c.turner@durham.ac.uk In the rare cases of galaxies that gravitationally lens more distant sources, additional constraints are made available.At its simplest, strong gravitational lensing provides a robust single estimate of the mass projected within the Einstein radius (Treu 2010).With a more sophisticated pixel-based analysis (e.g.following Dye & Warren 2005), lensing can also probe the slope and shape of the mass density profile in the vicinity of the lensed images, at large radii where dynamical studies are not as sensitive (e.g.Ritondale et al. 2019;Shajib et al. 2021;Etherington et al. 2022).However, lensingonly studies can be susceptible to degeneracies in the lens modelling and inherently lack sensitivity to the distribution of matter far from the lensed arcs.These degeneracies are independent of those affecting kinematic studies (Courteau et al. 2014).
The unification of mass and structure constraints from kinematic and dynamical modelling with the larger scale mass information from strong lensing allows further insight into the intrinsic properties of galaxies and the nature of DM.Access to a diverse range of spatial scales makes combined lensing and dynamical studies powerful in disentangling the stellar and the DM mass distributions of lens galaxies at their characteristic radii, and thus breaking the degeneracies between these two components (e.g.Treu & Koopmans 2004;Koopmans et al. 2006;Barnabè et al. 2009;Auger et al. 2010a;Treu et al. 2010;Oldham & Auger 2018;Shajib et al. 2021).For example, measurements of the Einstein radius from lensing studies can provide constraints on a lens galaxy's total enclosed mass, and by combining stellar dynamics with lensing studies, galaxy-scale strong lensing can provide robust measurements of the stellar IMF (e.g.Treu et al. 2010;Auger et al. 2010a;Smith et al. 2015).The sensitivity of joint lensing and dynamical studies to different mass scales also gives direct insight into the apparent and surprising near-isothermality of mass in early-type galaxies (ETGs) (e.g.Koopmans et al. 2009;Auger et al. 2010b;Li et al. 2018).More pertinently, this 'bulge-halo conspiracy' describes how the total mass distribution of ETGs can be described by a power law, but the mass profile of neither the baryonic nor dark matter components can be described by a power law on their own (Treu & Koopmans 2004;Treu et al. 2006;Humphrey & Buote 2010).For an up-to-date review of galaxy-scale strong lensing, including the application of lensing-plus-dynamics studies, see Shajib et al. (2022).
Double-source-plane lenses (DSPLs), or compound lenses, are a rare and valuable type of gravitational lens system, occurring when a single foreground lens galaxy simultaneously multiply-images two background source galaxies at different redshifts.The best studied example of a DSPL to date is the 'Jackpot'  = 0.222 lens (SDSSJ0946+1006), discovered serendipitously by Gavazzi et al. (2008) as part of the Sloan Lens ACS survey (Bolton et al. 2006).The Jackpot system consists of a bright ring at  = 0.609 and a further ring at a greater radius, indicating a more distant second source ( spec = 2.035) from which constraints on the cosmological parameters can be obtained (Collett & Auger 2014;Smith & Collett 2021).A further multiply-imaged source at  ≈ 6 has been reported by Collett & Smith (2020), making Jackpot a triple-source-plane lens system.
The Jackpot system hosts one of only a few cases of a dark substructure detected through lensing perturbations (see Vegetti et al. 2010).The substructure is cited as having a mass high enough that one would expect it to host a luminous galaxy, as well as having a surprisingly high central density.The primary lens has also been reported to have a steep density slope by Minor et al. (2021), a claim that sits in contention with that of earlier findings in Collett & Auger (2014).The authors note that the peculiar inferred properties of the subhalo could be due to a deviation from the CDM paradigm with respect to the particle physics of DM, such as dark matter self-interactions (see Colín et al. 2002;Vogelsberger et al. 2012;Zavala et al. 2019;Turner et al. 2021), in which case one would expect to detect many more highly concentrated substructures in future surveys.They also note, however, that the substructure properties could be affected by the lack of generality in their model, stating that a more flexible host galaxy combined with spatially resolved kinematics could provide stronger constraints for the subhalo concentration.
As a result of its astrophysical significance and applications in addressing open cosmological questions, the lensing properties of the Jackpot have been intensively studied.In contrast to this, and despite the additional advantages that kinematic data can bring to breaking degeneracies, only relatively limited kinematic measurements exist for the lens system (e.g.Auger et al. 2009;Sonnenfeld et al. 2012;Spiniello et al. 2015).These studies used single-slit observations, such as the 1 arcsec-wide slit used in Sonnenfeld et al. (2012), yielding measurements to ∼ 1 arcsec from the centre of the lens, and the 2 × 1.5 arcsec slit of Spiniello et al. (2015), and thus cover only a moderately small radial distance.As a result, this limits the area from which measurements can be obtained to within the bright Einstein ring, effectively excluding the radial ranges at which DM contributions become more significant and resulting in restricted spatial information.
In this work, we present a combined kinematic and dynamical analysis of the Jackpot lens galaxy.We apply template fitting methods to deep integral field unit (IFU) spectroscopic data from the MUSE instrument, and employ anisotropic dynamical Jeans models that are robustly constrained by the lensing mass at the Einstein radius, to measure the total 2D-projected density profile slope of the lens galaxy.This paper is organised as follows.In Section 2, we describe our IFU data and highlight how they differ from and improve on those of previous studies.In Section 3, we detail our use of multiple-component template fitting to the observed galaxy spectra, in order to extract the stellar kinematics of the foreground lens galaxy.Section 4 describes a set of kinematic predictions obtained through dynamical modelling of a gNFW + stars model with varying parametrisations.The observed stellar kinematics are then compared with the dynamical model predictions to recover the best-fitting DM density slope and mass, and the total projected logarithmic density slope for a mass profile with our best-fitting parameters.In Section 5, we discuss our findings and place them in the context of previous studies of the Jackpot lens.We also consider the robustness of our assumptions and the fundamental limitations of our modelling approach.Findings are summarised in Section 6.

DATA
Deep integral field unit (IFU) adaptive-optics-assisted spectroscopic data (project ID 0102.A-0950 as described in Collett & Smith 2020;Smith & Collett 2021) were obtained from a 5.2 hour total integration time with the MUSE instrument on the VLT.When compared to previous single-slit observations, this data provides full azimuthal area coverage, high sensitivity and a high angular resolution (0.5 arcsec FWHM).As a result, we can measure the two-dimensional stellar kinematic properties out to ∼ 10 kpc (∼ 2.7 arcsec) from the centre of the lens, as well as resolving the inner profile inwards to ∼1 kpc.Measuring the kinematics out to larger radii with greater precision than previous studies were able to achieve helps to break the degeneracies between the stellar and DM mass components by probing the radii where DM contributions become more significant.Thus, allowing us to span the relevant range for projected mass slope measurements at the Einstein radius that allow us to place our results into the context of previous pure lensing studies.
Fig. 1 shows a collapsed MUSE image of the lens and its environment.Denoted are the extent of the stellar kinematic measurements and the locations of the foreground lens galaxy and the first Einstein ring.Some of the structure visible at low surface brightness is caused by MUSE sensitivity variations, but the overall asymmetry, with an extended plume to the north, is reproduced in other imaging, e.g.Fig. 3 of Sonnenfeld et al. (2012).
The physical scale resolution of our observations (1.86 kpc for a 0.5 arcsec FWHM at  = 0.222) is an order of magnitude coarser than that of the typical dynamical analyses of nearby early-type galaxies.For example, the ATLAS 3D survey (Cappellari et al. 2011) which, with an angular FWHM of 1.5 arcsec and galaxy redshift of  ≲ 0.01, has a physical spatial resolution typically ∼0.15 kpc.However, the modelling techniques that we use here are also routinely applied to galaxies with comparable physical resolutions to our data.For example, the MaNGA Survey galaxies (Law et al. 2016), which have a poorer median spatial resolution of 2.54 arcsec FWHM for a median redshift of  = 0.037, leading to a physical FWHM of 1.8 kpc.

STELLAR KINEMATICS
In order to construct dynamical models and constrain the 2Dprojected total mass profile slope of the lens, robust measurements of the spatially resolved stellar kinematic properties must be obtained.To achieve this, stellar template fitting is employed, and measurements of the first and second order moments of the lens galaxy's 2D stellar kinematics are acquired.To account for the strong continuum contributions from the higher-redshift source, a multiple-component fitting technique must be adopted.The implementation of spatial binning methods is necessary to achieve a high enough signal-to-noise (S/N) ratio for precise kinematic measurements.

Kinematic Template Fitting
Prior to fitting spectral templates to our measured galaxy spectra, we implement the VorBin Cappellari & Copin (2003) two-dimensional adaptive spatial binning method.Fundamentally, this works by evaluating pixels that are near to each other and grouping them together in order to achieve an approximately constant S/N ratio.The data described in this work have been binned into 53 bins.
To measure the kinematics of the lens galaxy, we make use of the stellar template fitting software ppxf (Cappellari & Emsellem 2004;Cappellari 2017), which implements penalised pixel-fitting to extract the moments of the line of sight velocity distribution from galaxy spectra.This method works by searching a library of template spectra over a range of metallicities and ages, in this case an ensemble of simple stellar populations (SSPs), to fit to the observed spectra from each bin and thus take kinematic measurements.With this approach, we consider a wavelength range of 4699 Å to 7408 Å and apply the penalised pixel-fitting method to each bin individually.The template library is taken from the E-MILES stellar population models (Vazdekis et al. 2016), chosen for their broad spectral range (1680 to 50000 Å), good resolution (FWHM = 2.5 Å from 3540 Å to 8950 Å) and age/metallicity coverage (-1.79 < [M/H] < +0.26 and ages above 30 Myr).Unlike other SSP models that do not extend far enough in to the blue end of the spectrum for us to recover the relatively distant z = 0.609 source, this broad spectral range allows us to fit to the younger spectral features of the source galaxy, as well as the older features of the lens elliptical galaxy.This basic approach works well on the very central regions of the system where there exists little contamination from the source light.
In order to place robust constraints on the mass model, it is desirable to map the kinematic measurements in a range of radii that span the Einstein radius, given our mass constraint at this distance and the high area coverage made available by our MUSE data.However, there are a number of bins in these outer regions where strong continuum contributions from both the source and lens galaxies are present.This demands more complexity in our template fitting than the standard ppxf treatment provides and makes it necessary to fit to both lens and source simultaneously, as whilst the source emission lines can be simply masked out, the bright source continuum can not sufficiently be masked without a significant loss of spectra.This is demonstrated in the first panel of Fig. 2, which shows strong source emission lines and prominent Balmer absorption lines in our observed data that are not being fit by the total fit template.Whilst these lines are mostly present in the spectra of the bluer source galaxy, not accounting for them will affect our ability to recover the kinematics of the lens galaxy as ppxf otherwise compromises by using templates from an older but higher-sigma population to account for the source.This motivates the necessity for a multiple-component fitting approach, whereby this problem is alleviated through the addition of both a second set of stellar templates to model the source galaxy, and a set of gas emission templates.The improvement of the total fit to the observed spectra obtained through the addition of these additional components can be seen in the second and third panel of Fig. 2 respectively, and is described in Section 3.2.

Multiple Component Fitting
Fig. 3 shows examples of both a central bin almost entirely dominated by lens light and a highly contaminated bin at the Einstein radius.To each bin spectrum, we mask out prominent subtraction residuals from the brightest sky lines and fit stellar kinematic components corresponding to the lens and the source galaxies and also a small number of gas components.To do this we define a total of four components to be fit: • a combination of SSP templates at redshift  ≈ 0.222 to represent the lens; • a combination of SSP templates at redshift  ≈ 0.609 to represent the bright source; • a gas template at the source redshift, corresponding to the Balmer series emission lines from H to H and with fixed Case-B recombination flux ratios; • a further gas template at redshift z = 0.609 for the [O ii] doublet1 .
Figure 2. A zoomed-in view of the total fit to the galaxy spectrum for the 1 star, 2 star and 2 star + gas template fitting methods.The 1 star fit highlights how the standard usage of ppxf gives a poor fit to the composite spectra, i.e. pixels contaminated by light from the Einstein ring.We can see strong emission lines in the observed data that are not being matched by the total fit template derived from the 2 star fitting method.The 2 star + gas panel demonstrates the significantly improved fit of the total model to the observed spectrum.Also indicated are the absorption features at the lens redshift, and the emission peaks at the source redshift.
Despite the main focus of this fitting being the extraction of the lens kinematics, each component has its own velocity and velocity dispersion in order to recover the lens component without bias 2 .
A combination of the best fitting SSP templates is determined and averaged to make a template model for each of the components and the fractional contribution of each model to the observed spectrum is optimized in order to construct the total fit.The top panel of Fig. 3 shows that for bins dominated by lens light, the total fit is almost entirely constructed from SSP templates at redshift  = 0.222.In contrast to this, bins such as the one shown in the bottom panel are composed of a relative contribution of the lens, source and gas templates.
The addition of the two gas components to account for the Balmer 2 Recall that the data were binned according to the continuum S/N ratio to optimise the recovery of the lens stellar kinematics, and therefore the binning is poorly configured for spatially resolving the source kinematics.
emission lines and the oxygen doublet allows us to obtain a further improvement to the fit with significantly reduced residuals.This can be seen in the second and third panels of Fig. 2 when compared with the first panel.The emission features that are not well fit by the 1 star model, such as the deep absorption wings on either side of the H- emission line (at ∼6600 Å in the observed frame), are now present in the total fit of the 2 star + gas model.

Kinematic Results
Fig. 4 shows the derived velocity map for an ∼ 8×10 arcsec 2 area and clearly displays the axis of rotation of the lens galaxy, giving maximum velocities of ∼ ±100 km s −1 about a kinematic axis with a misalignment with the photometric axis of the order 10°, dependent on the radius at which they are measured3 .This rotation has been hinted at in previous works, with Sonnenfeld et al. (2012) noting evidence for some rotation in their analysis; however, the paper states that the stellar kinematics of the lens galaxy are dominated by pressure support, rather than rotation.Our data confirms that the kinematics are indeed dispersion dominated, with  2 rot ≪  2 .Fig. 4 demonstrates the way in which our high-areacoverage MUSE data and multiple-component fitting allow us to fully map the 2D-kinematic properties of the Jackpot lens galaxy, out to a much greater radius than previous single-slit studies, allowing us to now fully characterise the rotation proposed by Sonnenfeld et al. (2012).
Fig. 4 also shows the velocity dispersion map of the lens galaxy and illustrates the way in which  falls from ∼ 280 km s −1 at the centre of the lens to ∼ 230 km s −1 at a radius of ∼ 2 arcsec, and highlights a discernible velocity dispersion gradient in the inner region.
Fig. 5 shows the radial velocity dispersion profile derived from this study, as compared with that of Spiniello et al. (2015), Sonnenfeld et al. (2012) and Auger et al. (2009).Our measurements agree with those of previous studies in the regions with overlapping coverage (i.e. the inner ∼2 arcsec).In the regions at larger radii than this, our measured profile exhibits the "quadrupole" structure seen in the velocity dispersion map in Fig. 4, with velocity dispersions that are lower on the major axis and higher on the minor.This behaviour is still found if we adopt very large bins in the outer regions; we further discuss the implications of this result in Section 5.2.

DYNAMICAL MODELLING WITH JAM
In this section, we model the spatially-binned kinematic measurements using the anisotropic Jeans model approach, as implemented in jam (Cappellari 2008(Cappellari , 2012)).We impose a robustly constrained aperture lensing mass to reduce the freedom in the models.By optimising the model parameters as described in this section, we maximise the likelihood of the observed root-mean-square velocities,  rms , obtained from the observed velocities and velocity dispersions in Section 3.3 and given as  rms = √  2 +  2 .As seen by the extended envelope in Fig. 1, and further reflected in Fig. 4, we measure our stellar kinematics out to ∼ 10 kpc (∼ 2.7 arcsec) from the centre of the lens.However, preliminary tests of our dynamical models showed indications that the Jackpot is not dynamically simple at greater radii, and we note here the possible evidence for tidal interactions found in Sonnenfeld et al. (2012).This presents a fundamental limitation of how well we can model the lens kinematics in the outer region, as this component cannot be accommodated by our simple dynamical models, which are limited to oblate axisymmetric mass and luminosity distributions.If indeed just tidal debris, we do not expect this to be very dominant in mass, but including the tracers in this region in the modelling could bias the recovery of the mass components of interest.
As a result, we determine that our Jeans models are not appropriate to model the lens kinematics at greater radii, and as such, all dynamical models described herein are fit to the exclusion of our nine outermost bins.We thus opt to restrict the radius for dynamical model predictions to 1.95 arcsec, which is still sufficient to allow projected mass slope measurements at the Einstein radius.We further discuss this choice and its implications in Section 5.2.The measured kinematics for this restricted region are shown in Fig. 6, the final panel of which shows the fractional increase of the second velocity moment by ordered rotation, which further demonstrates that the galaxy is dispersion dominated within the fitted radius.

Mass Model
To construct our dynamical mass models, the mass distribution of the lens galaxy is represented as a combination of the stellar mass deprojected from the observed light (with constant stellar mass-to-light ratio), a spherical dark matter halo, and an excess central mass component, as further detailed in this section.The model components are parameterised in terms of their fractional contribution to the lensing mass inside the Einstein radius.Therefore, all models explored are consistent with the lensing configuration.

Stellar Mass
We follow the common approach of de-projecting the luminosity density using a multi-gaussian expansion (MGE) (as per Emsellem et al. 1994;Cappellari 2002) fit to high-resolution imaging from HST, where the MGE projected surface brightness is given as where ( ′ ,  ′ ) are the polar coordinates on the plane of the sky ( ′ ,  ′ ). is the number of Gaussian components with total luminosity given by   , observed axial ratio of 0 ≤  ′  ≤ 1 and width along the major axis of   , as per Cappellari (2002).
As this method would be unreliable in the presence of the bright Einstein arc, we fit our MGE to an image from which the arc has been subtracted using a lens reconstruction model, as is described in Section 4 of Etherington et al. (2022).Specifically, we use an image from HST (F814W-band) from which a PyAutoLens-fitted (Nightingale et al. 2018(Nightingale et al. , 2021) ) parameterised source model has been subtracted.This image, with a rest-frame wavelength of ∼6700 Å, is not expected to be sensitive to any modest variations in age and metalicity, and therefore its luminosity in this band is assumed to trace the stellar mass surface density reasonably well (in shape, but not normalisation) in the absence of any IMF gradients.The parameterised model assumes the source galaxy to be well-described intrinsically by a smooth Sérsic-profile galaxy.Although this method yields a less precise source-subtraction than the alternative pixelised-source approach, it is preferred here to avoid overfitting of the lens galaxy light.
As implemented in the MgeFit code, the luminosity profile is measured in a number of elliptically defined sectors of this residual image of Jackpot and the MGE fits a series of Gaussians to the profile, describing the intensity and shape of the total surface brightness.The PSF of the HST image was approximated with a single Gaussian of 0.1 arcsec FWHM; however, little difference was found in the recovered parameters when this value was varied within reasonable limits.The projected mass density in stars is then assumed to be  proportional to the projected luminosity density, i.e. with a constant stellar mass-to-light ratio.Table 1 presents the stellar mass density MGE model for the best fitting model parameterisation of the Jackpot lens, as described in Section 4.3.Fig. 7 shows the F814W-band HST image, the HST image with the smooth Sérsic-profile source model subtracted, and the residuals obtained from subtracting the MGE luminosity profile from the source-subtracted model.In the absence of the source, the lens is parameterised simply by an ellipse with a Sérsic profile.We note Table 1.MGE ★ for the Jackpot galaxy.The columns represent, left to right, the projected surface mass density multiplied by the best fitting stellar fraction parameter in the free  jam models, the MGE width, and axis ratio.that whilst the source-subtracted model fails to reproduce all of the observed features in the arcs that reflect real structures in the source, this treatment is sufficient to allow the natural robustness of MgeFit to follow the true luminosity distribution of the lens galaxy.Fig. 8 shows an isophote plot of both the source-subtracted image of Jackpot and the MGE model.The MGE model surface brightness is in good agreement with that of the galaxy in the inner ∼ 2 arcsec (i.e.where our kinematics are measured), but poorly reproduces the shape of the observed isophotes beyond this, where the light from the diffuse outer envelope becomes significant, as also seen in Fig. 1.We note the presence of this same effect in Fig. 1 of Posacki et al. (2015), but find that in comparison, our treatment of the HST image works in reducing the contamination from the bright arcs and improves the model fit in the inner regions.

Dark Matter Halo
A similar process to the one described in Section 4.1.1 is followed to obtain a second series of Gaussians describing the galaxy DM surface density.Here, the DM halo is assumed to be spherical and have density well-described by a generalized NFW profile (gNFW; Zhao 1996) of the form where  is the physical radius and with  = 1 corresponding to the original NFW slope.Here  s is the scale radius, fixed at 100 arcsec based on reasonable assumptions for the virial radius of the order 600 kpc (as per Gavazzi et al. 2007), a halo concentration parameter of c vir ≈ 6 (from a mass-concentration relationship from Macciò et al. 2008) and given the NFW halo density profile.It was found that varying this value slightly (i.e. by 20%) had very little effect on the recovered kinematic estimates, as is expected since the observational constraints are well inside  s .
In our use of the gNFW profile, we are not assuming any specific physical origin of any difference with respect to the NFW profile, but a slope of  > 1 could represent, for example, a contraction in response to the baryonic mass (e.g Blumenthal et al. 1986;Gnedin et al. 2004).

Excess Central Mass
We include an additional 'excess' central mass component,  cen , through the mechanism that jam uses to model a central black hole.In our implementation, this is understood to subsume any real point mass, i.e. a black hole, as well as any centrally concentrated mass in excess of a constant stellar mass-to-light ratio.This excess mass is modelled as an additional, very small Gaussian component and is explored in the range 0 ≤  cen ≤ 7×10 10 M ⊙ , giving the mass model a further degree of flexibility at smaller radii.We further discuss the simplification of a mass component of ∼zero radius in Section 5.3.

Model Normalisation
In order to derive the total-mass surface density, the Gaussians resulting from the MGE fits to the luminous and dark components, as described in Sections 4.1.1 and 4.1.2respectively, are combined with the excess central mass described in Section 4.1.3.A projected luminous fraction (i.e. the stellar mass as a fraction of the Einstein mass from lensing) is explored in the range 0 ≤  ≤ 1.We reduce the freedom in the dynamical models by rigidly enforcing a robustlyconstrained lensing mass  E at the Einstein radius, normalising the MGE such that where  ★ (<  E ),    (<  E ) and  cen are the projected mass contributions from stars, dark matter and the central excess, to the total Einstein-aperture lensing mass.Given the redshifts of the source and lens in the Jackpot system, the measured  E = 1.397 arcsec (Collett & Auger 2014) yields  E = 3.08 × 10 11  ⊙ .

Anisotropic Modelling and Parameter Search
The main goal of this study is to measure the slope of the total 2Dprojected mass profile.To minimise any bias in this quantity, we adopt a model that is sufficiently flexible to reproduce the observed kinematics.
The normalised MGE descriptions of the surface brightness and the total-mass surface density, along with the variable parameters described in this section, are used to calculate a prediction of the projected  rms field for an anisotropic axisymmetric galaxy model.The prediction includes smoothing by the MUSE PSF approximated Free  1.73 +0.17 as a single Gaussian with FWHM of 0.5 arcsec4 .The predicted second moments are calculated at the luminosity-weighted Voronoi bin centres for comparison with the observed data.We explore parameter space for two distinct model sets with velocity dispersion ellipsoids aligned with the cylindrical (, ) polar coordinate system (the cylindrical jam method, Cappellari 2008Cappellari , 2012)); a set of models with the DM density slope as a free parameter, and a further set of models with the DM density fixed as a NFW profile.Within the cylindrically-aligned paradigm, the orbital anisotropy parameter is defined as   = 1− 2  / 2  (Cappellari 2008).We explore the parameter space using a Markov chain Monte Carlo (MCMC) Ensemble sampler, as described by Foreman-Mackey et al. ( 2013), and summarise the free parameters of our models described above as follows: • , the DM density power law slope.Values in the range 0.5 ≤  ≤ 3 are explored, with  = 1 describing the standard NFW slope; • , the orbital anisotropy parameter which, in the cylindricallyaligned case, describes the ratio of the radial velocity dispersion to the vertical component.Here we consider values of −0.6 ≤  ≤ 0.6, where a negative value of beta indicates a relatively larger vertical velocity dispersion.
• , the galaxy inclination, with a lower limit of 35°imposed by the minimum observed axial ratio of the MGE Gaussians describing the distribution of the kinematic-tracer population.This is defined such that an inclination of 90°corresponds to the edge-on case.
We define the prior probability density function (PDF) for all of our parameters with flat priors in , ,  ★ ,  cen and , and impose physically motivated constraints on the extreme values as described above.

Results
Using the likelihood derived from the predicted and observed  rms , and given the priors above, we sample the posterior PDF for the five model parameters.Fig. 9 shows the marginalised parameter constraints, which are given in Table 2.The results of our modelling can be summarised as follows: • The preferred excess central mass is well-constrained to be ∼ 8 × 10 9 M ⊙ for the free  models, and ∼ 7 × 10 9 M ⊙ for the NFW models.If this component indeed represents only a central black hole, this would make it somewhat over-massive given the galaxy properties.For a galaxy such as the Jackpot lens with a central velocity dispersion of ∼ 280 km s −1 , from the standard black hole mass vs. sigma relation (van den Bosch 2016), one would expect a central black hole with mass ∼ 1.6 × 10 9 M ⊙ (with a scatter of 0.49 ± 0.03).Given this discrepancy, it seems unlikely that all excess central mass is contributed by a black hole.We also note the potential for a systematic overestimation of the central mass component, as described in Appendix A and further discuss the implications of this in Section 5.3.
• Our inferred DM density slope, taken from the model set with  as a free parameter, is 1.73 +0.17 −0.26 , significantly steeper than an unmodified NFW profile (i.e. = 1).Such a slope could represent a baryon-contracted halo appropriate to a massive galaxy (Sonnenfeld & Cautun 2021).We find that models with a slope flatter than NFW are strongly disfavoured, but not disallowed.
• In the free  models, there is a strong degeneracy between  ★ and , which prevents us from obtaining closed inference on the stellar fraction.Instead, a range of 0.15 <  ★ < 0.57 is somewhat weakly favoured.Thus, the inferred mass budget inside the Einstein radius for the models with a free DM density slope is  ★ :  DM :  cen = 0.38 : 0.59 : 0.03, albeit with substantial uncertainty.For the models with a fixed NFW-like DM density, this parameter is much more tightly constrained and the mass budget is 0.71 : 0.27 : 0.02.The expected stellar mass fraction under the assumption of a Chabrier IMF is  Chab ★ = 0.26 ± 0.07, whilst with a Salpeter IMF we expect  Salp ★ = 0.46 ± 0.13 (Auger et al. 2009).Our preferred stellar mass fraction for the free  models is broadly consistent with that of a Salpeter IMF.For the NFW models, we find a stellar mass fraction heavier than the predictions of a Salpeter IMF and that is inconsistent with a Chabrier IMF.
• In the free  case, a DM mass fraction of  DM = 0.59 +0.24 −0.19 was inferred.This sits in good agreement with that of Gavazzi et al. (2008), who found a surprisingly high DM mass fraction inside the effective radius (2.0 arcsec for the Jackpot) of  DM (<  eff ) = 0.73 ± 0.09.
• The velocity ellipsoid is strongly constrained to be nearly isotropic in both model sets, with orbital anisotropy parameters of  = −0.03± 0.03 and  = −0.01 ± 0.04 for the free  and NFW cases respectively, which is consistent with the low values typically found, e.g. for  > 200 km s −1 galaxies in the ATLAS 3D survey (Cappellari et al. 2011).The model likelihood is rather insensitive to the orbital anisotropy, although this is perhaps to be expected given the context of the assumed cylindrically-aligned coordinate system.
• The stellar mass as a fraction of the total lensing mass is highly degenerate with the free DM density power law slope; this is to be expected as, to first order, either reducing the stellar fraction at the expense of the DM fraction, or flattening the DM profile, act to effectively predict larger root-mean-square velocities at larger radius.Despite  being free to explore values up to a slope of 3, we see that our models do not exceed a slope of ∼ 2 due to the hard prior imposed on the lower limit of the stellar fraction, and the strong anti-correlation between these two parameters.This behaviour is further demonstrated in Fig. 10, which shows that the models with the steepest DM profile are also the models with the lowest stellar mass and smallest excess central mass.The apparent cut-off at  ∼ 2 corresponds to the case where the DM component accounts for the totality of the dynamical mass; there is simply no further flexibility in the mass budget for a steeper  to be explored.
• In the model set with a fixed NFW-like DM density slope, the stellar fraction is anti-correlated with  cen .We see that the models with a high stellar fraction prefer smaller excess central masses as, in the absence of flexibility from a free DM slope parameter, both components act to account for any compact additional central mass, so are free to compensate one another.The model set with a free  does not demonstrate this behaviour, as there is instead the freedom for interplay between the three mass components in the models, constrained by the total lensing mass.In these cases, the greater concentration of DM in the central regions reduces the necessity for such a large excess central mass component, and suppresses the stellar contribution.
• We found the galaxy inclination to be completely unconstrained in the free  models, and only somewhat constrained in the NFW-like models, with an inferred inclination value of 51 +21 −10 .More 'edge-on' inclinations up to 90°were not strongly excluded, but were instead disfavoured.In both cases, the inclination shows no significant covariance with the parameters of interest.Fig. 11 shows the 2D-projected mass profiles for each of the individual mass components described in Section 4.1, for both the free  and NFW cases.
The large range of profiles and relative contributions in each model demonstrate the way in which the free  models have a large degree of flexibility to effectively 'trade off' mass components at different characteristic radii.These models prefer a steep ( > 1) DM density slope, and hence their range of DM projected mass profiles here are shallower than seen in the NFW models.The median contribution of the DM mass component is significantly greater for  ≲ 3 arcsec than in the NFW models, as this model set has the freedom to allow for a greater concentration of DM mass in the central regions whilst still being constrained by the total lensing mass, thus dominating the total mass in the centre of the galaxy.
Conversely, the models with a DM density slope fixed at NFWlike show DM mass contributions that do not dominate the total projected mass until well outside the Einstein radius ( ≳ 4 arcsec), and instead prefer a more substantial stellar mass fraction.Indeed, Fig. 11 shows that, in the absence of a steep enough DM halo, the NFW-like models are almost solely dominated by the stellar contribution in the central regions.Notwithstanding uncertainty in the relative component contributions, the projected mass profile at the Einstein radius is tightly constrained by the data, and the flexibility afforded to the model parameters results in almost indistinguishable recovered total mass slopes.
Figs. 12 and 13 show the  rms predictions from our best model, which reproduces the data well over the full range of measurements, compared to predictions of a model without a DM halo, that attributes the gravitational potential solely to luminous matter and any additional central mass, and of a further model with no excess central mass component.As we only constrain our total mass around the Einstein radius, our three contributing mass parameters (stellar, DM and central) are free to compensate one another, e.g. in the absence of the excess central mass required to reproduce the observed higher velocities in the innermost regions (where the  cen component would dominate), consequent adjustment between the DM and stellar components are necessary to produce the steep inner  rms profile.The kinematic predictions from the best fitting dynamical model exhibit a clear gradient in  rms , with values ranging from ∼ 280 km s −1 at the centre of the lens to ∼ 230 km s −1 in the outer regions, in close agreement with the observed kinematics, with relatively small residuals.The models that do not include all three mass components, however, unsurprisingly fail to predict the observed kinematics whilst simultaneously constraining model parameters that are realistic, as the models are not sufficiently flexible.Fig. 13 shows that the angular structure of the model without an excess central mass component is similar to that of the best model, while the no-DM model demonstrates prominent high  rms lobes along the major axis.The model without a DM halo substantially overestimates the  rms in a range 0.2 arcsec ≲  ≲ 1.2 arcsec.This is a result of the total mass distribution of the lens being solely constrained by the centrallyconcentrated luminous MGE in this case, and lacking an extended mass component.The model without an excess central mass component successfully reproduces the observed kinematics at all radii, but does so at the expense of requiring an unrealistic stellar mass fraction (  ★ = 0.09) and a very steep DM density profile ( = 1.97).This is consistent with the - ★ panel in Fig. 9 and the dark blue data points in Fig. 10.While formally consistent with the lensing and dynamics,  relative to the tabulated values of Auger et al. (2009), a model with such a low stellar mass fraction would imply a stellar IMF that is a factor of 2-3 lighter even than the Chabrier IMF.

Total Mass Profile Slope
To compare our result to previous (non-dynamical) lensing studies of the Jackpot, we must relate our composite profile results to the power-law total density profile used in the lensing literature, (e.g.Collett & Auger 2014;Minor et al. 2021;Etherington et al. 2022) 5 .The sum of a stellar component and a dark matter halo does not yield a power-law total density profile, yet making this approximation has historically been acceptable, due to the so-called bulge-halo conspiracy (Dutton & Treu 2014): for lensing ellipticals, the sum of stellar and dark matter components is approximately isothermal over the length-scales probed by strong-lensing constraints.
To define an equivalent power law index from our composite model, we use all of the Gaussian components for the profile (constructed as described in Section 4.1) to compute the 2D-projected mass profile.For our free  and NFW model sets, the local slope of this profile, at the Einstein radius, is 1.03 ± 0.03 and 1.07 ± 0.04 respectively.As this is a profile of integrated mass, a value larger than unity corresponds to a shallower-than-isothermal density profile.This quantity is not directly comparable to the lensing literature, as lensing does not strictly measure the local density slope at the Einstein ring.We instead calculate the power-law profile with the equivalent lensing effect as our inferred composite model.Collett (2014) and Kochanek (2020) showed that the slope inferred from lens modelling is sensitive to the radial derivative of the deflection angles at the location of the lensed images.The parameter we need to calculate is the dimensionless quantity , which is well-constrained by lensing data as per Kochanek (2020) and written as where   is the mean convergence at the Einstein radius and  ′′ is the second derivative of the deflection profile at   .For models that use a power law relation with surface mass density ∝  −  , this quantity is given by  = 2( − 1). (5) We therefore use our composite mass model to calculate  from Equations 4, and convert this to a lensing-equivalent power-law slope  using Equation 5.This gives us the value of  that should be Figure 13.Kinematic maps for the root-mean-square velocities,  rms , obtained from the data, alongside the best fitting models for each scenario described in Fig. 12.Also shown are the model residuals.The best fitting model exhibits a clear  rms gradient and the residuals are ∼ ±40 km s −1 .The no-DM model exhibits higher residuals than the best fitting model and completely fails to predict the observed root-mean-square velocities obtained from the observed data.measured from a lensing-only study for a mass profile with our best-fitting parameters.As shown in Fig. 14, sampling from the posterior distribution, we find a projected logarithmic density slope of  = 0.96 ± 0.02 for our free  models and  = 0.93 ± 0.02 for our NFW models.This indicates a density profile that is marginally shallower than the isothermal case, where  = 1.

Comparisons with the literature
The near isothermal density slope that we have measured is reasonable given previous studies of populations of massive elliptical galaxies.Koopmans et al. (2006); Grillo et al. (2008); Duffy et al. (2010) all found that for the ensemble of lens galaxies, the average total density slope is approximately isothermal.
For the Jackpot lens specifically, our lensing-equivalent total density-profile slope of  = 0.96±0.02(for the free  model) is in good agreement with previous measurements of 1.00±0.03,1.1±0.1,1.03 ± 0.02, 1.06 ± 0.03 and 1.01 ± 0.18 from Gavazzi et al. (2008), Sonnenfeld et al. (2012), Collett & Auger (2014), Etherington et al. (2022) (lensing only) and Etherington et al. (2023) (lensing + dynamics), respectively6 .Our model also does an excellent job of predicting the lensing deflection at the location of the second ring, with a deflection angle of 1.89 ± 0.03 arcseconds.This is ∼0.2 arcseconds less than the Einstein radius of the second ring, but a slight underestimate is to be expected since our model neglects the presence of mass in the first source.Collett & Auger (2014) inferred an SIS Einstein radius of 0.16 ± 0.02 arcsec for the first source.Adding this to our deflection angle yields a second Einstein radius of 2.05±0.04arcsec, entirely consistent with the 2.07 ± 0.02 arcsec measured by Gavazzi et al. (2008).Additionally, we infer a DM density slope of 1.73 +0.17 −0.26 which is significantly steeper than an unmodified NFW profile (i.e. = 1), but sits in good agreement with the lensing + dynamics DM profile slope of 1.7 ± 0.2 found for this galaxy by Sonnenfeld et al. (2012).
Conversely, our derived total density slope is inconsistent with that from the Minor et al. (2021) study, who find a surprisingly steep density slope of  = 1.32 ± 0.04 from lensing alone.The paper reports a correlation between the derived density slope and the inferred subhalo mass; thus, a shallower lens slope would have implications for the claim of an unusually dense and massive halo.Their footnote 2 cites however a fairly modest decrease of 25% in subhalo mass if the slope was confirmed to be approximately isothermal.
When drawing such comparisons, it is important to note the contrast in (and limitations of) the methods used in dynamical and pure lensing studies.In dynamical analyses, restrictive assumptions are often made on the lens galaxy axisymmetry and orbital structure; pure lensing studies typically assume simple power laws and a simplified linear external shear.As shown in this study, and also found in Etherington et al. (2023), lensing-only studies appear to predict marginally steeper projected density slopes than lensing + dynamics studies do.If the observed discrepancies stem from a lack of complexity in the dynamical modelling, one would expect that as you relax the simplifications and introduce spatially resolved structure, the recovered slopes from the two methods should be in closer agreement.As we have shown, this is not the case, with our reported projected density slope showing less consistency with pure lensing slopes than that of the lensing + dynamics measurement from Etherington et al. ( 2023)'s one-aperture kinematics.Instead, we consider the possibility that the disparity arises from either deficiencies in the lens modelling, or more subtle limitations in the kinematics that can only be solved through more sophisticated dynamical modelling techniques such as Schwarzschild models (Schwarzschild 1979).

Robustness of Assumptions
As noted in Section 4, there exists a possible degeneracy between the models that maximise the likelihood in the inner and outer regions of the galaxy.There seems present an orthogonal angular dependence, such that we see a relatively high velocity dispersion along the minor axis at large radii, but conversely along the major axis at small radii.This is particularly evident in the velocity dispersion panel of Fig. 4, where we suspect that this behaviour may be related to a greater influence from DM at this radii.
Given the apparent non-axisymmetry implied by the low surfacebrightness envelope at large radii, and the possible signatures for past interactions (as speculated by Sonnenfeld et al. 2012), we excluded measurements from the nine Voronoi bins at  ≳ 2 arcsec from our preferred modelling.If we instead fit to the full, unrestricted range of data, we recover a resulting projected density slope of  = 1.042 ± 0.02, which sits in slightly closer agreement with the lensing only studies, but the model now provides a poorer fit to the kinematics, especially in the outer regions.
The kinematics in non-axisymmetric mass distributions can in principle be tackled using more general dynamical models, such as the orbit-based approach of Schwarzschild (1979) (e.g.see Poci & Smith 2022).This approach would perhaps mitigate the limitations imposed by our simple model, but would in turn demand much more stringent requirements on the data with a necessity for a very high signal-to-noise ratio that is unfeasible with the present observations.This is especially true when dealing with the pervasive contamination of source light at the first Einstein radius.

Central Mass in Excess of a Constant Mass-to-Light Ratio
In the construction of our dynamical models we allowed for an additional compact mass, that is not described by the luminosity distribution or the NFW profile, and used this to describe any excess central mass.We find that the preferred excess central mass of our free  models is well constrained to be ∼ 8.23 × 10 9 M ⊙ , which if attributed to a central black hole only, would be an outlier relative to black hole scaling relations (e.g.Gebhardt et al. 2000;Tremaine et al. 2002;Thomas et al. 2016;van den Bosch 2016).Given the  BH −  relation derived by van den Bosch (2016), it is expected that the true black hole contribution to the central mass component will be less than ∼ 10 9 M ⊙ .Although the form of the scaling relations at the highest masses is still uncertain (e.g.Thomas et al. 2016), it is unlikely that the Jackpot lens galaxy truly harbours a ∼ 10 10  ⊙ black hole.
A more plausible explanation is that in our analysis we assume a constant / ★ but, whilst we do not expect variations in age in galaxies of this type, there may be a metallicity gradient present, leading to a modest M/L ★ gradient such as described by Tortora et al. (2010).Moreover, if there is a radial gradient in the stellar initial mass function as reported by Martín-Navarro et al. (2015); La Barbera et al. (2017); van Dokkum et al. (2017) (but see Alton et al. 2017;Vaughan et al. 2018), then a much more substantial / ★ gradient may be present, and indeed Collett et al. (2018) saw exactly this for a nearby lens.These works would suggest a steep increase in mass within ≲1 kpc.Whilst we expect that the I-band image would be a faithful tracer of the stellar mass profile for a modest age and metalicity variation, this might not be true in the case of radial IMF variation; however, such a gradient should be absorbed into our central mass component, to first order, at the resolution of the present data.
In Appendix A we show that tests with synthetic MUSE data suggest a potential for overestimation of the central mass.In the absence of a central mass in the input data, a mass comparable to that predicted by the standard  BH −  relation was recovered.However, this is an order of magnitude smaller than the excess central mass recovered from the real data.

CONCLUSIONS
We have presented results from a kinematic and dynamical analysis of the Jackpot lens galaxy using new data obtained from a 5 hour MUSE integration, in order to constrain the 2D-projected total mass profile slope.To account for contamination from the source galaxy light, we implemented a multiple component fitting technique adapted from the ppxf code that extracts the lens galaxy kinematics to first and second order.Simple gNFW + stars dynamical models were constructed with parameterised orbital anisotropies, DM density power-law slopes, stellar mass fractions and excess central mass components, and a robustly constrained aperture lensing mass was imposed.The posterior PDF for the model parameters was sampled, and a chi-squared likelihood maximised to derive the projected total density slope.This is the first 2D spatially resolved kinematics study for this system, and confirms the significant signature of rotation detected in previous studies.We measure rotation about the minor axis of  ≈ ±100 km s −1 and a steep decrease in velocity dispersion from a central value of  ≈ 290 km s −1 to  ≈ 200 km s −1 in the outer regions.Notwithstanding the strong presence of rotation, the galaxy is dispersion-dominated at all radii.The kinematic measurements are consistent with those of previous single-slit studies of the Jackpot lens (i.e.central velocity dispersions of  sonn12 = 287 ± 11 km s −1 and  spin15 = 300 ± 22 km s −1 ), but are now fully mapped out in two-dimensions.
From the jam dynamical modelling, we infer a mass budget inside the Einstein radius that is dominated by stars (∼70 %) if the halo slope is fixed to the NFW shape.For modified halos, we infer a larger DM fraction, and a DM density slope of  = 1.73 +0.17 −0.26 , which is significantly steeper than NFW (i.e. = 1).This is in agreement with previous results for the Jackpot itself (Sonnenfeld et al. 2012) and, at face value, supports the scenario in which DM haloes contract in response to the presence of a massive baryonic component.Indeed, our measured value of  is consistent with the gNFW slope of  = 1.57used by (Sonnenfeld & Cautun 2021) to model contracted haloes for massive lens galaxies.While similar conclusions have been reached by some ensemble studies of lens galaxies (e.g.Grillo 2012), others find that unmodified NFW haloes are preferred by the data (e.g.Shajib et al. 2021).
Our fitted models yield a 2D-projected total mass profile slope for the Jackpot lens of 1.03 ± 0.03, and a lensing-equivalent projected logarithmic density profile slope of  = 0.96 ± 0.02.Thus we confirm most-pure lensing results in finding a near isothermal profile (e.g.Collett & Auger 2014;Etherington et al. 2022).Our profile is inconsistent with the surprisingly steep slope measurement of Minor et al. (2021).
The main goal of this work, and of ongoing extensions relating to the Jackpot, is to suppress the remaining systematic errors and degeneracies, so as to fully exploit the cosmological potential offered by this unique lens system.An improved analysis of the lensing properties, exploiting multi-band imaging for all three sources, is presented by Ballard et al. (in preparation).Future extensions to our work will incorporate the measured kinematics simultaneously with the lensing information, to create an even more detailed picture of the Jackpot system.
Additionally, having obtained spatially-resolved, high sensitivity, high resolution MUSE data for a larger sample of lens galaxies, we have begun to incorporate the techniques described in this paper to obtain robust ensemble total and DM density profile slope measurements.This will furthermore allow us to place improved constraints on the stellar initial mass function, the distribution of mass within galaxies and the structure of DM haloes, thus ultimately furthering our understanding of the intrinsic properties of galaxies and the nature of DM itself.

Figure 1 .
Figure 1.Collapsed MUSE image of the Jackpot triple-source lens system and its environment.The red dashed line denotes the extent of the stellar kinematic measurements, and the white contours highlight the location of the foreground lens galaxy and the first Einstein ring.Inset is a F814W image from HST of the Jackpot lens and the first and second Einstein rings.

Figure 3 .
Figure 3.A demonstration of the multiple components used in the total fit to the galaxy spectrum, as demonstrated for both a central bin almost entirely dominated by lens light and a bin near the Einstein radius that is highly contaminated from the lensed source.The stellar kinematic components correspond to the lens galaxy at  = 0.222 and the source galaxy at  = 0.609.Also present are the two gas components accounting for the [O II ] doublet and Balmer emission lines.

Figure 4 .
Figure 4. Velocity and velocity dispersion maps for an ∼ 8×10 arcsec 2 region of the Jackpot, overlaid on contours of the galaxy flux.The upper panel shows a clear signature of rotation along the major axis and the lower panel shows a falling velocity dispersion gradient.

Figure 5 .
Figure 5.The velocity dispersion of each Voronoi bin as a function of radius from this work.For comparison, we show measurements from Spiniello et al. (2015), Sonnenfeld et al. (2012), and SDSS(Auger et al. 2009).Our data are broadly consistent with previous measurements in the inner regions, but show a clear decline in velocity dispersion towards a larger radius, where the precision of earlier datasets is lower.

Figure 6 .
Figure 6.The measured kinematics from the restricted data, corresponding to the region used for dynamical modelling.The final panel shows the fractional increase of the second velocity moment caused by ordered rotation.

Figure 7 .
Figure 7.The F814W-band HST image, the HST image with the parameterised source model subtracted and the residuals from the MGE model subtraction.Each panel corresponds to a field of view of 9 arcsec on a side.

Figure 8 .
Figure 8. Wide-field and zoom-in isophote plots of the source-subtracted image of Jackpot (black).The contours of the MGE model surface brightness are overlaid in red.The MGE is a good fit in the central regions, but deviations are evident at large radii due to the outer stellar envelope.

Figure 9 .
Figure 9.The posterior PDF for the free  and NFW model parameters.The contours show the 68 and 95% confidence regions.The parameters explored are: the inner slope of the DM density profile, ; the orbital anisotropy parameter, ; the stellar mass as a fraction of the total lensing mass,  ★ ; any central mass in excess of a constant stellar mass-to-light ratio,  cen ; the galaxy inclination, .The diagonal plots show the marginalised posterior densities for each parameter.

Figure 10 .
Figure 10.The stellar mass as a fraction of the total lensing mass for all free  models, as a function of the DM density slope.Coloured points represent the preferred excess central mass and the black dashed line denotes contours for the corresponding DM mass fraction.

Figure 11 .
Figure 11.The 2D-projected mass profiles for each individual mass component: the stellar (green), DM (pink) and excess central mass (blue) components, as well as the total projected mass profile (black), for both the model with a free DM slope parameter (left) and the model with NFW-like DM (right).To enable comparisons, the slopes derived by Collett & Auger (2014) and Minor et al. (2021) have also been plotted.The vertical and horizontal dotted lines denote the Einstein radius and mass, respectively.

Figure 12 .
Figure 12.The mean and azimuthal range of the root-mean-square velocity predictions from our best model (green) compared to the best possible model without an excess central mass (purple), and the best model predictions without accounting for the mass in DM (orange).Also shown are the measured  rms from each Voronoi bin as a function of radius (black).The horizontal error bars of the observed data represent the width of the Voronoi bins.The model set without a DM halo fails to predict the observed kinematics, while the model set without an excess central mass can successfully reproduce the kinematics, although in this case the recovered stellar mass and DM halo profile are unreasonable (see text).

Figure 14 .
Figure 14.The 2D-projected total logarithmic density slope obtained from each of our two model sets; the free  model and the NFW-like model.The slopes derived by Sonnenfeld et al. (2012), Collett & Auger (2014), Minor et al. (2021) and Etherington et al. (2022) and their respective errors have been included for comparison.

Figure A1 .
Figure A1.Kinematic maps showing the  rms fields of the synthetic data.We show two models, one created from the best-fitting free  model and one representing a 'vanilla' scenario with no excess central mass and an NFW-like halo.The left hand panel shows the high resolution model, the middle panel is the noise-free spatially binned model and the third is the synthetic model after adding observational noise.

Figure A2 .
Figure A2.The marginalised posterior densities for each parameter for the real data and the two sets of mock data, generated from the 'best' and 'vanilla' models.The dashed line represents the 'truth' value used to generate the mock data.The parameters explored are: the inner slope of the DM density profile, ; the orbital anisotropy parameter, ; the stellar mass as a fraction of the total lensing mass,  ★ ; any central mass in excess of a constant stellar mass-to-light ratio,  cen ; the galaxy inclination, .

Table 2 .
The median and 68% confidence bounds for the model parameters from both of our model sets.

Table A1 .
The median and 68% confidence bounds for the recovered model parameters from our mock data.