Galaxy structure with strong gravitational lensing: decomposing the internal mass distribution of massive elliptical galaxies

ABSTRACT

We investigate how strong gravitational lensing can test contemporary models of massive elliptical (ME) galaxy formation, by combining a traditional decomposition of their visible stellar distribution with a lensing analysis of their mass distribution. As a proof of concept, we study a sample of three ME lenses, observing that all are composed of two distinct baryonic structures, a ‘red’ central bulge surrounded by an extended envelope of stellar material. Whilst these two components look photometrically similar, their distinct lensing effects permit a clean decomposition of their mass structure. This allows us to infer two key pieces of information about each lens galaxy: (i) the stellar mass distribution (without invoking stellar populations models) and (ii) the inner dark matter halo mass. We argue that these two measurements are crucial to testing models of ME formation, as the stellar mass profile provides a diagnostic of baryonic accretion and feedback whilst the dark matter mass places each galaxy in the context of LCDM large-scale structure formation. We also detect large rotational offsets between the two stellar components and a lopsidedness in their outer mass distributions, which hold further information on the evolution of each ME. Finally, we discuss how this approach can be extended to galaxies of all Hubble types and what implication our results have for studies of strong gravitational lensing.

galaxies: formation, galaxies: evolution; gravitational lensing: strong

1 INTRODUCTION

The work of Hubble (1926) famously classified galaxies into three groups: ellipticals, spirals, and irregulars. Today, using samples of hundreds of thousands of galaxies, these classifications have been broadly established to hold across all galaxies, in the local Universe (Hoyos et al. 2011; Vika et al. 2013; Vulcani et al. 2014) and at high redshift (Buitrago et al. 2008; Bruce et al. 2012; Chevance et al. 2012; Van Der Wel et al. 2012). The bulk structure of a galaxy can be quantified by its one-dimensional projected surface brightness profile. The Sersic function, |$\mu ({R}) \propto R^{\frac{1}{n}}$| (Sersic 1968; de Vaucouleurs 2005), has proven extremely useful for this purpose across the entire Hubble sequence. The bulk of the light in a typical elliptical galaxy is well described by n ∼ 4 and an exponential disc structure by n ∼ 1. Although most galaxies can be labelled with a specific morphology on the Hubble diagram without much ambiguity, they can also exhibit subdominant structures with other morphologies (see Graham 2013). For example, spirals have bulges (Thomas & Davies 2006; MacArthur, González & Courteau 2009), ellipticals have discs (Thomas et al. 2002; Bluck et al. 2014), and irregulars may show signs of both bulges and discs (Margalef-Bentabol et al. 2016). This has motivated descriptions that use multiple Sersic profiles (Lackner & Gunn 2012; Bruce et al. 2014; Vika et al. 2014; Kennedy et al. 2016) to decompose galaxies into their constituent physical structures.

Fitting the light distribution of a galaxy with a superposition of light profiles is a challenging and highly degenerate problem (e.g. Häußler et al. 2013). For example, the extended wings of a central bulge can be difficult to separate from structures further out (e.g. a disc), as they blend with one another and the background sky emission. Bulges, discs, and bars may all appear as compact central structures in a galaxy, which are hard to decouple from one another (Kormendy 1995). Absorption of light by dust can also impact the fit, leading to more ambiguities in the interpretation of a galaxy’s structure (Rest et al. 2001; Lauer et al. 2005).

Using integral field spectroscopy (IFS), the SAURON (Bacon et al. 2001), and ATLAS3D (Cappellari et al. 2011) surveys demonstrated the importance, when trying to infer a galaxy’s structure, of having data sensitive to its mass. For example, Emsellem et al. (2004, 2007, 2011) found that |$\sim 90{{\ \rm per\ cent}}$| of galaxies in a volume-limited sample showed some level of ordered rotation in their kinematics despite showing no signs of a disc in their light (Krajnović et al. 2011, 2013). IFS data also reveal that biases in the inferred galaxy structure may arise due to the 3D inclination (Devour & Bell 2017) and the triaxiality (Méndez-Abreu et al. 2010) of the different structures within galaxies, which limits inferences based on their 2D projected morphologies (Cappellari 2008; Emsellem et al. 2011; Weijmans et al. 2014).

In this work, we propose that strong gravitational lensing can both mitigate the degeneracies that arise when fitting a galaxy’s 2D light distribution and provide key insight on its underlying physical mass structure. Strong lensing is the deflection of light from a background source around a foreground lens galaxy, giving rise to multiple images of the source with characteristic distortions (see Kochanek 2004, for an overview). These distortions encode information on the foreground lens’s mass distribution and make it easier to separate the different galaxy components in comparison to using photometry alone. Conveniently, lensing data can be extracted from the same CCD imaging as the photometry, in contrast to IFS data that must be obtained separately: at significant cost in observing time, and typically at much coarser spatial resolution.

For this proof-of-concept study, we use our open-source lens modelling software pyautolens¹ (Nightingale & Dye 2015; Nightingale, Dye & Massey 2018) to fit the projected optical luminosity distribution of three isolated massive elliptical (ME) galaxies with velocity dispersions of 160–250 km s⁻¹, whilst simultaneously constraining their underlying mass structure via a strong lensing analysis. We recover the distribution of both light and mass in projection along the line of sight, enabling these two inferences of galaxy structure to be compared directly and circumventing ambiguities due to dust. Crucially, this approach allows us to confirm whether or not features apparent in the light distribution correspond to genuine physical structures.

This question is central to tests of ME galaxy formation in the LCDM cosmology, which predicts that they assemble their stellar mass both by dissipative star formation and by mergers with other galaxies. Each of these processes should give rise to physically distinct components in the phase-space structure of MEs that may have observable signatures in their stellar mass surface density profiles. In particular, active galactic nuclei (AGNs) activity at the centre of haloes more massive than |$\sim 10^{11}\, \mathrm{M_{\odot }}$| is thought to act as a ‘thermostat’ that shuts down in situ growth of their central galaxy by suppressing radiative cooling of fresh cold gas from a hot circumgalactic reservoir (Croton et al. 2006; Bower et al. 2006, 2017). The progenitors of present-day MEs form rapidly in overdense regions at z ≳ 2–4, hence the stars they form in situ (i.e. before AGN suppression) are characterized by high-metallicity and high-phase-space density (Baugh, Cole & Frenk 1996; Hopkins et al. 2005; De Lucia et al. 2006). By the present day, these stellar populations have old ages and red colours.

Although in situ star formation is suppressed in MEs, DM mass growth continues in an approximately self-similar fashion (e.g. Guo & White 2008). The most massive haloes at the present day, which are much more massive than the ‘threshold’ mass for AGN suppression, coalesce relatively recently. Their immediate progenitor haloes are themselves likely to be above the threshold mass, and hence also to host gas-poor galaxies with ‘red and dead’ metal-rich stellar populations. This leads to a picture in which the bulk of stellar mass in MEs is assembled through dissipationless mergers of several equally important fragments. These fragments have a broad range of DM halo masses but a narrow range of stellar masses (De Lucia & Blaizot 2007).

The structure of the resulting MEs (represented to lowest order by the scale and shape of their surface brightness profile) is therefore dominated by the ‘initial conditions’ and gravitational dynamics of these mergers, in which DM plays a dominant role (e.g. Cole et al. 2000; Boylan-Kolchin, Ma & Quataert 2006; Naab, Johansson & Ostriker 2010; Hilz, Naab & Ostriker 2013; Laporte et al. 2013). Broadly speaking, forward models of ME formation predict a composite structure arising from the superposition of phase-mixed stellar debris from each progenitor, each with a characteristic profile and scale (similar to the assembly of stellar haloes around less massive galaxies described by Amorisco 2015). For typical ME assembly histories, the dominant accreted components are expected to have similar amplitude, shape, and radial extent (e.g. Cooper et al. 2013, 2015). Consequently, the aggregate stellar density profile is (on average) unlikely to show strong inflections that correspond to transitions between regions of the galaxy dominated by debris from different progenitors. Indeed, observations (e.g. Kormendy et al. 2009) show that to good approximation, the overall structure may be described as a single component with Sersic n ≳ 4 (Schombert 2015).

The strongest observational features are therefore most likely to arise between the more extended aggregate accreted component and the centrally concentrated component, corresponding to an in situ, high-redshift ‘nugget’ (with much higher phase-space density as the result of dissipative formation). Sufficiently deep imaging of local MEs observes this (Huang et al. 2013, 2018b; Spavone et al. 2017), including components with different ellipticities, orientations (isophotal twists) (Oh, Greene & Lackner 2016a,), and centres (Goullaud et al. 2018). Spectroscopic studies of MEs are also consistent with this picture, whereby the stellar ages and alpha-element abundances of their central regions are both observed to increase with velocity dispersion (Greene et al. 2013, 2015; Raskutti, Greene & Murphy 2014), trends that appear washed out at larger radii where stars are anticipated to have been more recently accreted.²

Observations such as these have been interpreted as evidence for a model in which an ‘outer envelope’ arises from prolonged and significant accretion on to a high-density stellar core. This model, often referred to as ‘two-phase’ assembly (e.g. Oser et al. 2010), is broadly equivalent to the predictions of the LCDM models described above in the ME regime. Examples of massive compact ellipticals, thought to correspond to these dissipative cores before significant accretion, have been observed at high redshift but are extremely rare at the present day (e.g. Trujillo et al. 2006; Szomoru et al. 2011; Oldham et al. 2017). Overall, observations appear to be consistent with the evolution of the galaxy population as a whole predicted by LCDM models, once selection effects are taken into account (Laporte et al. 2013; Xie et al. 2015; Furlong et al. 2017; Roy et al. 2018). Furthermore, indirect tracers of dark matter mass (e.g. environment or size) detect correlations consistent with the above picture in populations of elliptical galaxies (Lani et al. 2013; Sonnenfeld et al. 2015; Huang et al. 2018a).

Nevertheless, decomposing MEs into multiple components remains very challenging, because their two (or more) superposed elliptical components (with similar stellar populations) are often well described by a single elliptical profile. This makes them the ideal case study for investigating how strong lensing can aid the decomposition of a galaxy’s structure. We will demonstrate that strong lens modelling: (i) allows us to confirm that the structures we see visually correspond to genuine mass components; (ii) provides direct access to the stellar mass distribution of each component (without stellar population modelling); and (iii) infers the central (≈10 kpc) dark matter halo mass of each galaxy. With a sample of just three objects, our conclusions primarily focus on discussing the utility of this method. In the future, existing samples of hundreds of lenses will enable a direct test of LCDM expectations for the assembly of MEs, in particular the fundamental relationship between halo mass, stellar mass distribution, and galaxy history. Future samples of 100 000 strong lenses (Collett 2015) will allow such an analysis to be generalized over the entire Hubble sequence.

This paper is structured as follows. In Section 2, we describe our sample of three ME galaxies selected for detailed study. In Section 3, we describe the pyautolens method for simultaneous photometric and strong lensing analysis. In Section 4, we investigate the mass structure of our lenses. In Section 5, we discuss the implications of our measurements, and we give a summary in Section 6. We assume a Planck 2015 cosmology throughout (Planck Collaboration XIII 2015).

2 DATA REDUCTION AND LENS SAMPLE

2.1 HST data reduction

We adopt a modified data-reduction pipeline that uses a combination of in-house tools and those from the standard Space Telescope Science Institute library. We first correct images for charge transfer inefficiency using the arctic software (Massey et al. 2010, 2014). We then bias subtract, flat-field and coadd multiple exposures using the calacs and astrodrizzle packages to create a final data product. Following Rhodes et al. (2007), we combine the images using a square drizzle kernel and a pixel fraction of 0.8, to a final pixel scale of 0.03 arcsec. To determine the point spread function (PSF) of this final image, we measure the focus of the telescope during each exposure, via the quadrupole moments of stars throughout that exposure (cf. Harvey et al. 2015), then coadd (appropriately rotated) tinytim models of each exposure’s PSF (Krist, Hook & Stoehr 2011). We store a single, 2D model of the combined PSF at the location of the lens, with a matching pixel scale of 0.03 arcsec, with the peak of the PSF the centre of a pixel.

This procedure produces a sky-subtracted, stacked image of each lens, which is used for the analysis. The variance in each pixel is estimated as |$\sqrt{\sigma ^2_{\mathrm{BG}} + d_{\mathrm{j}}}$|⁠, where |$\sigma ^2_{\mathrm{BG}}$| is the estimated variance in the background and d_j is the total counts detected in pixel j (accounting for the variation in exposure time due to dithering). The variances and background sky level are included as part of the modelling procedure (see N18).

2.2 Lens sample

The preliminary sample of three lenses used in this work is taken from the Sloan Lens ACS Survey (SLACS; e.g. Bolton et al. 2008), where galaxies were initially selected from the Sloan Digital Sky Survey (SDSS) main (Richards et al. 2002) and luminous red galaxy (LRG; Eisenstein et al. 2001) catalogues, as containing spectral emission lines from more than one redshift in a single fibre. The SLACS lens galaxies therefore primarily consist of ME galaxies with no known differences to other similar galaxies in the main or LRG parent samples (e.g. Bolton et al. 2006; Treu et al. 2009), other than a bias towards higher total mass.

Imaging of SLACS lenses is available across a range of wavelengths and spatial resolutions. We analyse monochromatic images, and have two primary considerations in selecting a waveband: (i) high spatial resolution to better constrain the lens model and (ii) longer wavelength coverage to be more sensitive to the lens galaxy’s old stellar populations (which contribute most significantly to its mass). We choose HST/Advanced Camera for Surveys (ACS) imaging taken with the F814W filter, which provides high-resolution imaging at rest-frame near-infrared (NIR) wavelengths. We also require that this imaging is available with an exposure time of at least one Hubble orbit. This work is focused on showing what can be measured for individual systems, so selection effects are not important for this work’s conclusions.

These criteria yield a sample of ∼40 ME lenses, which SDSS spectroscopy separates into three bins of velocity dispersion (uncorrected for aperture effects): low (230 < σ_SDSS < 270 km s⁻¹), intermediate (310 < σ_SDSS < 350 km s⁻¹), and high (σ_SDSS > 360 km s⁻¹). We select one lens randomly from each bin, resulting in the lenses SDSSJ0252+0039 (low σ_SDSS, hereafter SLACS1), SDSSJ1250+0523 (intermediate σ_SDSS, hereafter SLACS2), and SDSS1430+4105 (high σ_SDSS, hereafter SLACS3). None are cD or brightest cluster galaxies, but SLACS3 is in the field of a nearby cluster (Treu et al. 2009). SLACS1 and SLACS3 were both observed in HST programme GO-10886 and SLACS2 in programme GO-10494. Table 1 gives a summary of each lens’s observed properties (Bolton et al. 2008), including σ_SDSS, and the source and lens redshifts, z_src and z_lens. Table 2 shows the results of previous SLACS analyses – in particular photometry, stellar population modelling and stellar dynamics (Auger et al. 2010,), and mass modelling (Bolton et al. 2008).

Table 1.

Properties of the three SLACS lenses studied in this paper. Columns 2 and 3 give their Right Ascension and Declination in J2000 coordinates. Columns 4, 5, and 6 give the redshifts of the lens and source, and the velocity dispersion of the lens (uncorrected for aperture effects) (Bolton et al. 2008). Column 7 notes whether the lens is near a known cluster (Treu et al. 2009).

Target name	RA	Dec	zlens	zsrc	σSDSS	Near a cluster?
					(km s−1)
SDSSJ0252+0039 (SLACS1)	02h 52′45.21″	+00°39′58″	0.2803	0.9818	164 ± 12	–
SDSSJ1250+0523 (SLACS2)	12h 50′28.26″	+05°23′49″	0.2318	0.7953	252 ± 14	–
SDSSJ1430+4105 (SLACS3)	14h 30′04.10″	+41°05′57″	0.2850	0.5753	322 ± 32	MaxBCGJ217.49493+41.10435 (z = 0.270)

Table 2.

Previous measurements of the stellar and mass distributions of the SLACS lenses, taken from the literature. Columns 2–5 give the F814W-band magnitude (non-dust corrected), effective radius, axial ratio, and position angle of the lens galaxy, taken from Bolton et al. (2008) and computed by fitting each image with a de Vaucouleurs light profile. Columns 6–9 give the Einstein radius, axial ratio, position angle, and the number of source components, taken from Bolton et al. (2008) and inferred by fitting each lens with an SIE mass profile and parametric source. Columns 10 and 11 give the stellar mass for a Chabrier and Salpeter IMF and are taken from Auger et al. (2010). Columns 12–13 give the dark matter fractions for a Chabrier IMF and Salpeter IMF from Auger et al. (2010). Column 14 gives the power-law density slope inferred via a joint lensing and stellar dynamics analysis from Auger et al. (2010).

Target	I814	re (’)	q (b/a)	θ (°)	|$\mathrm{ log}\left[\frac{M_{\mathrm{ Ein}}}{\mathrm{M}_{\odot }}\right]$|	bSIE	q (SIE)	θ (SIE)	Nsrc	|$\mathrm{ log}\left[\frac{M_{\rm *}^{\mathrm{ Chab}}}{\mathrm{M}_{\odot }}\right]$|	|$\mathrm{ log}\left[\frac{M_{\rm *}^{\mathrm{ Sal}}}{\mathrm{M}_{\odot }}\right]$|	fDM (Chab)	fDM (Sal)	αSD
name	(obs.)
SDSSJ0252+0039 (SLACS1)	18.04	1.39	0.94	97.2	11.25	1.04	0.93	106.2	3	11.21 ± 0.13	11.46 ± 0.13	0.42 ± 0.18	−0.03 ± 0.32	1.57 ± 0.12
SDSSJ1250+0523 (SLACS2)	16.70	1.81	0.97	114.8	11.26	1.13	0.96	130.8	5	11.53 ± 0.07	11.77 ± 0.07	0.31 ± 0.11	−0.22 ± 0.20	2.30 ± 0.12
SDSSJ1430+4105 (SLACS3)	16.87	2.55	0.79	120.7	11.73	1.52	0.68	111.7	6	11.68 ± 0.12	11.93 ± 0.11	0.63 ± 0.10	0.35 ± 0.18	2.06 ± 0.18

2.3 Sample description

Based on previous SLACS results (Tables 1 and 2) and the reduced images in Fig. 1, we give a brief description of each lens in our sample. For clarity, Fig. 1 also shows an image of the background galaxies after subtracting the lens light (middle column), and de-lensed reconstructions of that source (right column), which are computed using the highest likelihood lens model found by pyautolens (which is described next). Neither the lens subtraction nor source reconstruction shows large differences when plotted using other high-likelihood lens models, thus these images are indicative of the analysis for any good lens model.

Figure 1.

The observed images (left column), lens-subtracted images (middle column), and source reconstructions (right column) for the three lenses in our sample, SLACSJ0252+0039 (top row), SLACSJ1250+0523 (middle row), and SLACSJ1430+1405 (bottom row). The lens subtractions and source reconstructions are computed using the highest likelihood model found for each image; however, their visual appearance does not change significantly for other high-likelihood models. The source reconstructions are computed using an adaptive pixellation and are represented as a Voronoi tessellation. The kpc per arcsec scale of each lens is 4.281, 3.729, and 4.335 kpc arcsec⁻¹, respectively.

2.3.1 SLACS1

SLACS1 is pictured in the top row of Fig. 1. The lens galaxy has an elliptical visual morphology extending smoothly from its centre to ∼3 arcsec. The lensed source shows arcs above and below the lens galaxy, which the source reconstruction reveals are a compact knot of light just outside the fold caustic. This lens is therefore a doubly imaged system with minimal extended structure in the source’s surface brightness.

The lens is the lowest mass galaxy in our sample, with an Einstein mass |$M_{\rm Ein} = 11.25 \, \mathrm{M}_{\odot }$| within its Einstein radius R_Ein = 4.40 kpc, a stellar mass (Chabrier IMF) |$\, \mathrm{M}_{\odot }^{\rm Chab} = 11.21\pm 0.13 \, \mathrm{M}_{\odot }$|⁠, and a velocity dispersion σ_SDSS = 164 ± 12 km s⁻¹. This translates to a stellar mass fraction within R_Ein of |$f^{\rm Chab}_{\rm *,Ein} = 0.40\pm 0.12$| for a Charbrier IMF. Notably, it has the lowest density slope (inferred via lensing and stellar dynamics) in the whole SLACS sample, with a value α_SD = 1.57 ± 0.12.

The HST imaging of SLACS1 had three faint galaxies not associated with the lens or source, which overlapped the masked region within which the analysis is performed. These were subtracted using a linear light profile fitting routine and their pixel variances were increased to infinity such that the analysis ignored them. The object at (x, y) = (3.0 arcsec, 1.5 arcsec) was included in the lens model as a singular isothermal sphere fixed to these coordinates. However, omitting it was found to have no impact on the results discussed in this work.

2.3.2 SLACS2

SLACS2 is pictured in the middle row of Fig. 1. The lens galaxy again visually shows one smooth extended component, but extending further, to ∼3.5 arcsec. The lensed source is complex, with a near-circular ring of light, two radial arcs protruding into the lens’s centre, and six distinct flux peaks. The source reconstruction reveals two galaxies on either side of the inner caustic, with features indicative of merging visible in their extended surface brightness profiles, in particular a disrupted tidal tail trailing each source. These tails form the inner ring of light, radial arcs, and extended arcs outside the lens in the image-plane.

Although this lens was chosen from the intermediate-σ_SDSS bin with a value σ_SDSS = 252 ± 14 km s⁻¹, it turns out to be the highest mass object in our sample with |$M_{\rm Ein} = 11.26 \, \mathrm{M}_{\odot }$| within R_Ein = 4.18 kpc and |$\, \mathrm{M}_{\odot }^{\rm Chab} = 11.77\pm 0.07 \, \mathrm{M}_{\odot }$|⁠. Its stellar mass fraction is also the highest in our sample with |$f^{\rm Chab}_{\rm *,Ein} = 0.68\pm 0.11$|⁠. Finally, its density slope α_SD = 2.30 ± 0.12 makes it one of the steeper slopes in SLACS and the steepest in our sample.

2.3.3 SLACS3

SLACS3 is pictured in the bottom row column of Fig. 1. The lens again shows one smooth extended component. The source shows two bright knots of light in a double image configuration, which the source reconstruction shows is the source galaxy’s central bulge. However, also visible is a wealth of additional extended structure surrounding this bulge, which in the image-plane forms multiple giant arcs around the lens. The source may be a face-on spiral galaxy, where the arms on the opposite side of the caustic are not visible due to the reduced magnification. Alternatively, it is a complex merging system.

This object has |$M_{\rm Ein} = 11.73 \, \mathrm{M}_{\odot }$| within R_Ein = 6.53 kpc, |$\, \mathrm{M}_{\odot }^{\rm Chab} = 11.68\pm 0.12 \, \mathrm{M}_{\odot }$|⁠, and σ_SDSS = 322 ± 32 km s⁻¹. It has a stellar mass fraction |$f^{\rm Chab}_{\rm *,Ein} = 0.33\pm 0.09$| and density slope α_SD = 2.06 ± 0.18. Unlike the other lenses in our sample, SLACS3 is near a cluster (Treu et al. 2009).

3 METHOD

When an extended source is gravitationally lensed, light rays emanating from different regions of the source galaxy trace different paths through the lens galaxy. This provides a projected and extended view of the lens’s gravitational potential, information that can be recovered with knowledge of the source’s unlensed light distribution. Thus, we are exploiting the lensed source’s surface brightness profile, in contrast to the previous photometric studies of SLACS lenses that used only the source’s position to measure the lens galaxy’s Einstein mass, M_Ein (e.g. Auger et al. 2010). It is this exploitation of the source’s extended information that allows us to constrain the lens’s underlying mass structure.

To perform this analysis, we use our new lens modelling software pyautolens,³ which is described in Nightingale et al. (2018, N18 hereafter), building on the works of Warren & Dye (2003, WD03 hereafter), Suyu et al. (2006, S06 hereafter), and Nightingale & Dye (2015, N15 hereafter). We refer readers to these works for a full description of pyautolens. Key points to note are:

• The lens galaxy’s light and mass distributions are fitted simultaneously;

• The source’s surface brightness distribution is reconstructed on an adaptive pixel grid (see the right column of Fig. 1);

• The Bayesian framework of S06 is used to objectively determine the most probable source reconstruction, and the complexity of the lens model is also chosen objectively via Bayesian model comparison (using multinest; Feroz, Hobson & Bridges 2009; Park et al. 2018);

• The method is fully automated and requires no user intervention (after a brief initial set-up) for the analysis presented in this work.

Table 3 summarizes the lens models that are available to fit the lenses in our sample, alongside the free parameters associated with each. The equations and a description of each profile are given in Appendix  A. In brief, a light profile is used to generate an intensity map I, which fits and subtracts the lens’s light. A mass profile is used to generate a surface density profile, κ(x), which is integrated to compute a deflection angle map |$\vec{{\alpha }}_{\rm x,y}$|⁠, allowing one to ray-trace image pixels to the source-plane for the source reconstruction. When multiple light or mass profiles are modelled simultaneously, their individual intensities or deflection angles are summed and their parameters are given numerical subscripts (e.g. x_l1, n_l2, etc.).

Table 3.

The light and mass profiles used to model each strong lens. Column 1 gives the model name, and column 2 whether it models the lens’s light and mass or just its mass. Column 3 states which mass component it represents. Column 4 gives its associated parameters, where all orientation angles θ are defined counterclockwise from the positive x-axis.

Model	Component(s)	Represents	Parameters
Sersic	Light+	Stellar Matter	(xl, yl) – Profile centre	ql – Axial ratio	θl – Orientation
	Mass		Il – Intensity	Rl – Effective radius	nl – Sersic index
			Ψl – Mass-to-light ratio	Γl – Radial gradient
Exponential (Exp)	’	’	Identical to Sersic with fixed nl = 1
Singular isothermal	Mass	Total (Stellar	(x, y) – Profile centre	q – Axial ratio	θ – Orientation
Ellipsoid (SIE)		+Dark Matter)	θE – Einstein radius	–	–
Singular isothermal sphere (SIS)	’	’	Identical to SIE with fixed q = 1 and θ omitted
Spherical NFW (NFWSph)	’	Dark Matter	(xd, yd) – Profile centre	–	–
			κd – Halo normalization	rs – Scale radius (fixed to 30 kpc)
Elliptical NFW (NFWEll)	’	’	(xd, yd) – Profile centre	qd − profileaxis − ratio
			κd – Halo normalization	rs – Scale radius (fixed to 30 kpc)
Generalized spherical NFW (gNFWSph)	’	’	(xd, yd) – Profile centre	γd – Inner slope
			κd – Halo normalization	rs – Scale radius (fixed to 30 kpc)
Shear	’	Line of sight	γsh – Magnitude	θsh – Orientation

pyautolens fits the lens galaxy’s light whilst simultaneously modelling its mass and reconstructing the source. Therefore, before each source reconstruction, a light profile for the lens is computed, convolved with the instrumental PSF and subtracted from the observed image, with the source analysis subsequently performed on this residual image. The light and mass models are therefore sampled from the same non-linear parameter space, meaning that the Sersic profile(s) used to fit the lens’s light can be translated to stellar density profiles (assuming a mass-to-light profile) and incorporated into the mass model. Thus, the light profiles may be constrained by both their fit to the lens galaxy’s light and the strongly lensed image of the source galaxy. We can therefore take two independent approaches to lens modelling that make different assumptions about the lens’s total mass distribution:

• Total-mass model: The model represents all of the lens galaxy’s mass (i.e. both stellar and dark matter), and the lens galaxy’s light profile is not used to constrain the mass model. Thus, the lens’s light profile is constrained only by a fit to the surface photometry of the galaxy (and not by the strong lensing analysis), as in the previous studies of galaxy structure referenced in the introduction. The total-mass model used in this work is a singular isothermal ellipsoid (SIE).

• Decomposed-mass model: The model splits the lens galaxy’s mass into stellar and dark matter components, with the lens galaxy’s light profile included in the former. Thus, the lens’s light profile is constrained by both its fit to the lens galaxy’s light and by the strong lensing analysis – thereby requiring that each component of the light profile corresponds to a genuine mass structure. The decomposed-mass models used in this work assume either a single elliptical Sersic or two-component Sersic (with Sersic index free) + Exponential (Exp, with Sersic index fixed to 1) profile for the stellar mass. A spherical Navarro–Frenk–White (NFWSph) profile is used for the dark matter. Our choice of an Exp in the two-component model is not an attempt to represent a specific galaxy structure (e.g. a disc). Instead, it is motivated by model comparison, which will be shown in the next section.

The ‘Results’ section of this work will compare both of the approaches above, to quantify the impact of including strong lensing information in measurements of galaxy structure. Both of the mass models above include an external shear term, which accounts for the contribution of line-of-sight structures to the strong lensing signal.

We stress that we are not comparing the total-mass model and decomposed-mass models so as to try and determine which ‘better’ fits the lensing data, nor do we compare their inferred density profiles in this work. Such a comparison would require us to fit a more complex total-mass profile (e.g. with a variable inner density slope) and to consider well-known degeneracies like the source-position transformation (Schneider & Sluse 2013). Our use of a total-mass model is motivated only to fit the lens galaxy’s light in a way that is independent of the mass model, as shown in N18.

When multiple light profiles are fitted (e.g. Sersic + Exp), one has to make a decision as to whether certain parameters that are shared by both components are independent of one another, like their centres (x_l and y_l), orientations (θ_l), and mass-to-light ratios (Ψ_l). Rather than making ad hoc assumptions about these parameters, pyautolens uses Bayesian model comparison to objectively determine the lens model complexity. At different phases of the analysis, different parametrizations of the light and mass models are fitted using multiple independent non-linear searches, computing the Bayesian evidence of each model via multinest. We favour a more complex model when its Bayes factor (ratio of Bayesian pieces of evidence) is 20 or above (considered ‘strong’ evidence in Bayesian statistics). Thus, the choice of lens model incorporates the principle of Occam’s razor, where a more complex model not only has to fit the data better, but better enough to justify its extra parameters and complexity. N18 demonstrated that this is an effective way to choose the appropriate model complexity, in simulations where the correct answer is known. The specifics of each model comparison phase are described next in their corresponding ‘results’ section.

When quoting the best-fitting lens models, the ‘most-probable’ sets of parameters are quoted, corresponding to the median of their one-dimensional marginalized posterior probability distribution. Errors are quoted at 3σ, corresponding to the 0.3 and 99.7 percentiles of this posterior.

4 RESULTS

We now describe how Bayesian model comparison informs the model of each lens galaxy. We quote the increase in evidence between two models, Δln ϵ = ln ϵ_model1 − ln ϵ_model2, referred to as the Bayes factor. In Bayesian statistics, a Bayes factor of 10 would roughly correspond to a ‘substantial’ preference of model 1, whereas our criterion of 20 is considered ‘decisive’.

The overall normalization of the evidence depends on how pyautolens’s adaptive analysis is initialized (see N18), which is performed independently for the total-mass and decomposed-mass profile analyses. This means that the evidence values lnϵ for these two profiles cannot be compared directly (which, as discussed in Section 3, is not something we perform in this work anyway). Instead, a fair comparison must compare Δln ϵ between different models of the total-mass profile and decomposed-mass profiles. For example, in the next section, we compare Δln ϵ for the Sersic and Sersic + Exp models for both the total-mass and decomposed profiles, but we do not compare Δln ϵ values between the total-mass profile and decomposed-mass profiles themselves. The inferred parameters for individual lens models are discussed throughout this section and the tables containing their values are at the end of the script in Appendix  B.

4.1 Light profile components

The number of light profile components is chosen by comparing two models, a single Sersic profile and a two-component Sersic + Exp model where the two components are centrally and rotationally aligned (x_l1 = x_l2, y_l1 = y_l2, and θ_l1 = θ_l2). We perform two independent runs of this phase, for each image, using: (i) a total-mass profile (SIE); (ii) a decomposed-mass profile (Light + NFWSph) with shared mass-to-light ratio (Ψ_l1 = Ψ_l2) between the multiple light components. The results are presented in Table 4 and address the following question:

• Are the lenses single-component or multicomponent systems?

  For all three images, the single Sersic model gives a significantly lower Bayesian evidence than the multicomponent model. This occurs regardless of whether we use the total-mass or decomposed-mass model. For the total-mass profile, the Bayes factors preferring the two-component models are Δln ϵ = 568 for SLACS1, Δln ϵ = 224 for SLACS2, and Δln ϵ = 297 for SLACS3. For the decomposed models, they are Δln ϵ = 539 for SLACS1, Δln ϵ = 106 for SLACS2, and Δln ϵ = 535 for SLACS3. Thus, the inclusion of lensing does not change the number of physical components we associate with each galaxy. However, it leads to lower Bayes factors in SLACS1 (−29) and SLACS2 (−118) and a higher factor in SLACS3 (238).

Table 4.

The Bayesian evidence of each model used in the light profile model comparison phase, comparing the Sersic and Sersic + Exp models for the lens galaxy light profile. For each image, this phase computes two models, assuming: (i) a total-mass SIE profile (rows 1, 3, and 5); (ii) a decomposed Light + NFWSph profile (rows 2, 4, and 6). Column 1 gives each image’s target name, column 2 the type of mass profile, column 3 the mass-to-light ratio assumption, column 4 the profile of the second light component (all exponential in this table), and columns 5 and 6 the Bayesian evidence values of each model, where bolded values show those chosen via model comparison.

Target name	Mass profile	Mass-to-light	Second light	Single Sersic	Sersic+Halo
		ratio	component profile		(Shared geometry)
|$\mathbf {SLACSJ0252\!+\!0039}$|	Total	N/A	Exponential	212 290.8	212 858.4
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Shared	Exponential	212 273.4	212 812.6
|$\mathbf {SLACSJ1250\!+\!0523}$|	Total	N/A	Exponential	216 990.6	217 214.5
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Shared	Exponential	217 074.0	217 180.3
|$\mathbf {SLACSJ1430\!+\!4105}$|	Total	N/A	Exponential	217 423.3	217 720.7
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Shared	Exponential	217 110.8	217 645.2

The parameters of the light profiles fitted in this section are given in Tables B1 and B2. For each lens, the single Sersic model goes to values of n_l1 > 3.5 and non-circular axial ratios (q_l1 around 0.9), both typical of massive elliptical galaxies. For the two-component model, the first component is centrally compact (lower R_l1), concentrated (higher n_l1), and round (q_l1 near 1.0), whereas the second component is more extended and flatter.

By comparing the parameters of the total-mass and decomposed-mass models given in Tables B1 and B2, we can determine whether the use of lensing changes the lens galaxy’s inferred light profile and therefore structure. For both the Sersic and Sersic + Exp models, consistent fits are obtained either by a photometric-only analysis, or one that includes lensing. Thus, our use of lensing information has (in the analysis so far) not affected the light profile we infer (as might be expected given Δln ϵ in SLACS1 and SLACS2 was lower for the decomposed-mass profile than the total-mass profile). There are two ways one can interpret this. On the one hand, it could indicate that the light profile fits from photometry alone are a reliable tracer of the lens’s underlying mass structure. On the other hand, if the assumed light profile is not a representative description of the lens’s underlying mass structure, lensing would be unable to change the model in a way that alters the inferred light profile anyway.

In Table 5, we compare a Sersic + Exp and Sersic + Sersic model, using both the total-mass and decomposed-mass profiles. For all images, we find that the Sersic + Exp model provides a higher Bayesian evidence than the Sersic + Sersic model. Furthermore, the Sersic indexes of the second Sersic components, n_l2, are all consistent with n_l2 = 1.0. Therefore, the Sersic + Exp model is used hereafter.

Table 5.

The Bayesian evidence of each model used in the light profile model comparison phase, comparing the Sersic + Exp and Sersic + Sersic models for the lens galaxy light profile. For each image, this phase computes two models, assuming: (i) a total-mass SIE profile (rows 1, 3, and 5); (ii) a decomposed Light + NFWSph profile (rows 2, 4, and 6). Column 1 gives each image’s target name, column 2 the type of mass profile, column 3 the mass-to-light ratio assumption, and columns 4 and 5 the Bayesian evidence values of each model, where bolded values show those chosen via model comparison. The 6th column shows the Sersic index inferred for the second component of the Sersic + Sersic model.

Target name	Mass profile	Mass-to-light	Second component	Second component	Sersic index
		ratio	exponential	Sersic
|$\mathbf {SLACSJ0252\!+\!0039}$|	Total	N/A	212 858.4	212 855.5	|$n_{\rm l2} = 1.16^{+0.23}_{-0.32}$|
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Shared	212 812.6	212 791.4	|$n_{\rm l2} = 1.15^{+0.27}_{-0.29}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	Total	N/A	217 214.5	217 189.2	|$n_{\rm l2} = 1.23^{+0.24}_{-0.26}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Shared	217 180.3	217 172.2	|$n_{\rm l2} = 1.26^{+0.24}_{-0.29}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	Total	N/A	217 720.7	217 721.3	|$n_{\rm l2} = 0.86^{+0.27}_{-0.30}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Shared	217 645.2	217 640.2	|$n_{\rm l2} = 1.03^{+0.24}_{-0.36}$|

4.2 Light profile geometry

Given a model with two components of stellar material is preferred for all three lenses, we next investigate whether these two components are geometrically aligned or offset. This is performed by comparing four Sersic + Exp models where their geometry is: (i) aligned (the results of the previous phase); (ii) offset rotationally (θ_l1 ≠ θ_l2); (iii) offset centrally (x_l1 ≠ x_l2 and y_l1 ≠ y_l2); and (iv) offset in both. We perform three independent runs of this phase, using a total-mass profile (SIE) or decomposed-mass profile (Light + NFWSph) with mass-to-light ratio(s) between the two light components shared (Ψ_l1 = Ψ_l2) or independent (Ψ_l1 ≠ Ψ_l2). These results are presented in Table 6 and address the following two questions:

• Are the two components rotationally offset from one another?

  For each lens in Table 6, there are two models that include a rotational offset: one where the centres are aligned (column 6) and one where they are offset (column 8). We quote the largest Bayes factors from either column. For all three lenses, rotational offsets are detected for both the total-mass model (Δ lnϵ = 54, 35, and 120) and decomposed model (Δ lnϵ = 53, 57, and 200). The inclusion of lensing therefore increases the significance of two detections.

• Are the two components centrally offset from one another?

  As for the rotational offsets, we quote here the largest Bayes factors available from columns 7 and 8 of Table 6. For the total-mass models, we detect a centroid offset in SLACS2, with a Bayes factor Δ lnϵ = 21. This increases to Δ lnϵ = 24 for the decomposed model with independent mass-to-light ratios, which also gives a detection in SLACS3 with a Bayes factor Δ lnϵ = 53. This is the first case where the decomposed-mass model has deviated from the total-mass model.

Table 6.

The Bayesian evidence of each model used in the light profile geometry model comparison phase, comparing different light profile geometries for the lens galaxy light profile. For each image, the phase is run three times using a model that: (i) assumes a total-mass SIE profile (rows 1, 4, and 7); (ii) assumes a decomposed Light + NFWSph mass profile where light components share their mass-to-light ratio (rows 2, 5, and 8); and (iii) where these mass-to-light ratios are instead independent (rows 3, 6, and 9). Column 1 gives each image’s target name, column 2 the type of mass profile, column 3 the mass-to-light ratio assumption, column 4 the profile of the second light component (all exponential in this table), and columns 5–8 the Bayesian evidence values of each light profile geometry, where bold-face values show those chosen via model comparison.

Target name	Mass profile	Mass-to-light	Second light	Sersic+Halo	Sersic+Halo	Sersic+Halo	Sersic+Halo
		ratio	component	(Shared	(Independent θl)	(Independent xl /	(Independent
			profile	geometry)		yl)	geometry)
|$\mathbf {SLACSJ0252\!+\!0039}$|	Total	N/A	Exponential	212 858.4	212 912.7	212 870.2	212 890.7
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Shared	Exponential	212 812.6	212 865.5	212 804.1	212 871.3
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Independent	Exponential	212 819.4	212 867.3	212 861.7	212 862.5
|$\mathbf {SLACSJ1250\!+\!0523}$|	Total	N/A	Exponential	217 214.5	217 228.1	217 223.5	217 249.0
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Shared	Exponential	217 180.3	217 194.0	217 210.9	217 219.7
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Independent	Exponential	217 183.2	217 216.6	217 205.2	217 240.4
|$\mathbf {SLACSJ1430\!+\!4105}$|	Total	N/A	Exponential	217 720.7	217 840.3	217 753.2	217 852.4
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Shared	Exponential	217 645.2	217 780.8	217 663.0	217 804.7
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Independent	Exponential	217 676.5	217 823.4	217 694.5	217 876.0

The inferred model parameters are given in Table B3. The detected rotational offsets are of the order of Δθ ≈ 70° for SLACS1, Δθ ≈ 80° for SLACS2, and Δθ ≈ 23° for SLACS3. SLACS2 and SLACS3 showed centroid offsets of the order of Δr ≈ 0.09 and 0.344 kpc, respectively. Table B3 confirms that the same geometric offsets are inferred between the two components, regardless of whether a total-mass or decomposed-mass model is assumed (even for SLACS3, where the centroid offset’s Bayes factor was below 20 for the total-mass model). Thus, the inclusion of lensing is not changing the model that we infer, but it is reaffirming that what we see in lens galaxy’s light is physically genuine.

Fig. 2 shows the projected geometry of the Sersic + Exp model in each lens. In SLACS1 and SLACS3, the Exp component can be seen to extend to much larger radii than the Sersic component, whereas in SLACS2 the inner component is of comparable size.

Figure 2.

A schematic of the mass distribution in each lens galaxy, using the geometrically offset Sersic + Exp model chosen for each lens in Table 6. The red represents the Sersic component and blue the Exp component and each ellipse is drawn at the half-light radius R_l of each component. These figures use the parameters of the decomposed-mass profile, but as shown in Table B3 the total-mass profile gives consistent values.

We now compare the geometry of the total-mass model’s SIE profile to that of the lens’s two-component light profiles, using the parameters given in Table B3. For SLACS1 and SLACS3, the SIE fit assumes an orientation close to that of the extended Exp component. Next, we will show that the Exp component dominates their inner mass distribution, thus the SIE profile is simply tracing the lens’s overall mass. For SLACS2 (where the Exp component contributes less to M_Ein), the SIE model assumes an orientation θ = 106°, halfway between the two light profile components (θ_l1 = 71° and θ_l2 = 146°). It also appears more spherical than these two components (q = 0.96, compared to q_l1 = 0.88 and q_l2 = 0.87), suggesting that the best-fittingSIE geometry is a compromise between the two components of the decomposed-mass model. The large centroid offset found in SLACS3 shows the same behaviour, whereby the total-mass profile’s centre (x = 0.013 arcsec and y = 0.018 arcsec) is half-way between that of the two-component light model (x_l1 = 0.03 arcsec, y_l1 = 0.03 arcsec, x_l2 = −0.01 arcsec, and y_l2 = −0.03 arcsec).

4.3 Mass-to-light profile

Next, we establish whether the two components of each lens share the same mass-to-light ratio and whether this ratio varies with radius. The results of using one or two mass-to-light ratios were presented in Section 4.2 and are shown in Table 6. This addresses the following question:

• Do the two components have different mass-to-light ratios?

  By comparing rows 2 and 3, 5, and 6, and 8 and 9 of Table 6, we can infer whether the mass-to-light ratios of the two components are the same. For SLACS1, we find no evidence for different mass-to-light ratios, whereas for SLACS2 and SLACS3 we prefer the model with independent Ψ_l’s with Bayes factors of Δ lnϵ = 21 and 72.

  Next, we run a new model comparison phase to investigate the radial variation, comparing two models (which assume the geometric offsets and mass-to-light ratio parametrizations preferred in the previous phases): (i) a model introducing a radial gradient Γ_l in the mass-to-light ratio of the Sersic component; (ii) a model introducing a radial gradient in the Exp component. The results for this phase (alongside the no radial gradient cases) are given in Table 7, addressing the following question:

• Does the mass in either component vary radially compared to its light?

  Table 7 shows that support for a radial gradient was not found in SLACS1 or SLACS3, but was in SLACS2 in the Sersic component with an evidence of Δ lnϵ = 30 and the Exp component with an evidence Δ lnϵ = 54. These detections signify a decrease of mass relative to light in the lens’s central regions, leading to a shallower inner mass profile. We therefore assume they are indicating the same effect and choose the highest evidence detection that is in the Exp component.

Table 7.

The Bayesian evidence of each model used to test for a radial variation in the mass-to-light profile of each lens. For each image, the phase is run three times using a model that: (i) assumes the same light profile geometry and mass-to-light ratios as determined in the previous phase (column 6); (ii) assumes this model again, but with a radial variation in the mass-to-light profile Γ of the Sersic component (column 7); and (iii) in the Exp component (column 8). Column 1 gives each image’s target name, column 2 the type of mass profile, column 3 the mass-to-light ratio assumption, column 4 if there is a rotational offset between the Sersic and Exp components, and column 5 if there is a centroid offset. Columns 6–8 give the Bayesian evidence of each model, where bolded values show the model chosen.

Target name	Mass profile	Independent Ψl?	Rotationally	Centrally offset?	No radial	Radial variation Γl	Radial variation Γl
			offset?		variation	(Sersic)	(Exp)
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Yes	Yes	No	212 867.3	212 867.4	212 847.4
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Yes	Yes	Yes	217 246.6	217 276.0	217 300.9
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Yes	Yes	Yes	217 922.6	217 918.8	217 931.1

The mass-to-light ratios of the first component Ψ_l1 and second component Ψ_l2 of each decomposed-mass model are given in Table B4. For SLACS1, the inferred values of Ψ_l1 and Ψ_l2 are consistent with one another, whereas for SLACS3 they are not, as expected from the results of model comparison. In SLACS3, when only one shared value of Ψ_l is assumed, Ψ_l is forced to a value that is not consistent with the inferred value of Ψ_l2. This negatively impacts the inferred mass model and highlights that lensing can distinguish the mass-to-light ratios of individual galaxy components.

For SLACS2, independent mass-to-light ratios were preferred by model comparison, but Table B4 shows they are consistent with one another. This is because, when this model was preferred, the radial gradient Γ_l was not included in the mass model. Thus, the omission of the radial gradient leads us to falsely prefer a model in which Ψ_l1 and Ψ_l2 are different. The negative value of Γ_l2 indicates that the mass model is placing less mass in the centre of SLACS2. This may explain why, for SLACS2, the increases in evidence throughout the previous phases were typically much lower than SLACS1 or SLACS3, or resulted in decreases of Δln ϵ. Decomposed-mass models assuming a constant mass-to-light profile were unable to obtain a sufficiently shallow density profile in the galaxy’s central regions, resulting in a worse fit compared to the total-mass profile. SLACS2 is the only lens whose source galaxy light extends all the way to its very centre (<1 kpc) in the form of two radial arcs. This may be why we are sensitive to its inner mass profile and thus able to detect a radial gradient, but find no constraints for SLACS1 or SLACS3.

4.4 Dark matter profile

Next, we run a model comparison phase assuming more complex models for the dark matter halo of each lens galaxy. This will tell us whether we can go beyond the spherical NFWSph model and inform us what impact the assumption of a relatively simple dark matter model has on our inferred stellar mass density profile. Both model comparisons use the mass-to-light profile chosen in the previous section; however, the mass-to-light gradient chosen for SLACS2 is omitted as the complexity of this parameter space was too difficult to sample robustly. The model comparison is composed of four phases: (i) a NFWSph (the result of the previous phase); (ii) an elliptical NFW profile, NFWEll, whose centre is fixed to the Sersic component but whose rotational angle is free to vary; (iii) a NFWSph where the centre is free to vary (and rotational angle fixed to the Sersic components); and (iv) a spherical generalized NFW profile, gNFWSph. These results are presented in Table 8 and address the following two questions:

• What do we learn about the dark matter geometry?

  Table 8 shows that none of the models using an elliptical NFW profile produce a Bayes factor above 20, meaning that there was insufficient information in the lens data to constrain the ellipticity of the dark matter halo or detect if it is geometrically offset from the baryonic mass. The models using a NFWSph with offset centres also failed to produce a detection.

• Can we constrain the inner dark matter slope?

  Table 8 shows that the spherical gNFW model did not produce a Bayes factor above 20 for any lens, again meaning this model could not be constrained using the lens data available.

Table 8.

The Bayesian evidence of models testing the use of alternative dark matter NFW profiles, for each lens. All phases assume the same light profile geometry and mass-to-light ratios as determined in the previous phase (except the mass-to-light gradient is omitted for SLACS2). For each image, the phase is run four times assuming a: (i) spherical NFWSph profile; (ii) an elliptical NFWEll profile with rotational angle free to vary and centre fixed to the Sersic component; (iii) a spherical NFWSph with centre free to vary; and (iv) a spherical generalized NFW gNFWSph profile. Column 1 gives each image’s target name and column 2 the type of mass profile. Columns 3–6 give the Bayesian evidence of each model, where bolded values show the model chosen.

Target name	Mass profile	NFWSph	NFWEll	NFWEll	gNFWSph
			(Free θd)	(Free xd and yd)
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	212 854.8	212 860.1	212 853.7	212 852.8
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	217 246.7	217 259.5	217 264.2	217 255.9
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	217 922.0	217 925.5	217 928.6	217 930.1

The results of model comparison demonstrate that dark matter profiles more complex than a NFWSph cannot be constrained reliably using the lens imaging in this work. Nevertheless, in Table B5 the inferred parameters of the dark matter distributions for these models are shown. The axial ratios for the elliptical NFW profiles are all around 0.95 and therefore close to spherical. It is possible therefore that model comparison prefers a spherical NFW simply because the inner dark matter of each lens is close to spherical. The inner slopes of the gNFW profiles tend to values below the NFW value of 1, but are consistent with it within errors, again raising the possibility that these Bayes factors are low because the true profiles are close to 1. Given the low Bayes factors, we advise caution in the interpretation of these values and will seek to better understand how reliable our dark matter measurements are in future work.

Table B5 also shows the inferred parameters that control the stellar mass distribution for models that assume a NFWSph, NFWEll, and gNFWSph for the dark matter profile. For all three lenses, the inferred parameters governing the stellar mass profile (Ψ_l, R_l, n_l, I_l) are consistent with one another regardless of the assumed dark matter profile. The one exception is Ψ_l2 for SLACS3, which shows lower values for the NFWEll (6.43) and gNFWSph (6.78) models than the NFWSph (7.71). These differences are small enough to have no significant impact on integrated quantities of the mass profile (e.g. masses, mass fractions, etc.).

We thus conclude that our lens models are sensitive to the inner dark matter halo masses, but cannot measure their geometry or inner density slope. Our inferred stellar mass profiles are insensitive to assumptions about dark matter. This was expected, given that the light profiles we infer for the lens galaxy are consistent for both the total-mass and decomposed-mass profiles; if dark matter played an important role, we would expect to see the decomposed model adjust its stellar mass distribution to compensate for the missing complexity in the dark matter component. Nevertheless, the lack of degeneracy between stellar and dark matter is somewhat surprising and will be investigated further in the future. We believe this is because of the differences between the stellar and dark matter distribution’s density profiles and projected orientations.

4.5 Masses and surface density profiles

The highest Bayesian evidence model for the decomposed-mass model fit to each lens galaxy is as follows:

• SLACS1 – A Sersic + Exp + NFWSph model, where the Sersic and Exp components are rotationally offset (≈70°), centrally aligned, and share the same mass-to-light ratio.

• SLACS2 – A Sersic + Exp + NFWSph model, where the Sersic and Exp components are rotationally offset (≈80°), centrally offset (≈0.1 kpc), and have independent mass-to-light ratios including a radial gradient in Exp component.

• SLACS3 – A Sersic + Exp + NFWSph model, where the Sersic and Exp components are rotationally offset (≈20°), centrally offset (≈0.3 kpc), and have independent mass-to-light ratios.

Table 9 gives mass estimates for each lens galaxy, using its highest Bayesian evidence decomposed model. The masses of the different components within radii of 10 and 500 kpc, using circular apertures, are shown. Mass estimates are broadly consistent with those found in previous SLACS analysis, especially given the uncertainties that arise due to varying the initial mass function. Dark matter halo masses of the order of |$10^{12-13}\, \mathrm{M}_{\odot }$| are consistent with the expectation that we are studying massive elliptical galaxies. However, it should be noted that our lens model only constrains the inner dark matter halo profile directly. The extrapolation of this profile to large radii and thus mass estimates at ≳50 kpc are dependent on what assumptions we make about their scale radii, which for this work were fixed to |$10\, R_{\rm l}$|⁠. Our dark matter mass estimates at large radii may be subject to change if our assumptions about the scale radius are not reliable.

Table 9.

The mass measurements of each galaxy using the final model chosen by Bayesian model comparison, estimated within circular apertures of radii 10 and 500 kpc. For each galaxy, the mass (column 3), stellar mass (column 4), Sersic mass (column 5), Exp mass (column 6), and dark matter mass (column 7) are given.

Target name	Radius	Total mass	Stellar mass	Sersic mass	Exponential mass	Dark matter mass
	(kpc)	(⁠|$10^{11} \, \mathrm{M}_{\odot }$|⁠)	(⁠|$10^{11} \, \mathrm{M}_{\odot }$|⁠)	(⁠|$10^{11} \, \mathrm{M}_{\odot }$|⁠)	(⁠|$10^{11} \, \mathrm{M}_{\odot }$|⁠)	(⁠|$10^{11} \, \mathrm{M}_{\odot }$|⁠)
|$\mathbf {SLACSJ0252\!+\!0039}$|	10.0	3.12 ± 0.18	2.01 ± 0.29	0.43 ± 0.05	1.58 ± 0.10	1.11 ± 0.30
|$\mathbf {SLACSJ0252\!+\!0039}$|	500.0	33.53 ± 8.08	2.09 ± 8.18	0.43 ± 0.05	1.66 ± 1.16	31.14 ± 8.21
|$\mathbf {SLACSJ1250\!+\!0523}$|	10.0	4.78 ± 0.03	1.65 ± 0.48	1.127 ± 0.25	0.38 ± 0.13	3.12 ± 0.72
|$\mathbf {SLACSJ1250\!+\!0523}$|	500.0	91.12 ± 20.00	2.36 ± 20.00	1.67 ± 0.34	0.69 ± 0.28	88.66 ± 20.39
|$\mathbf {SLACSJ1430\!+\!4105}$|	10.0	8.17 ± 0.04	6.19 ± 0.16	2.63 ± 0.07	3.56 ± 0.15	1.97 ± 0.18
|$\mathbf {SLACSJ1430\!+\!4105}$|	500.0	67.90 ± 4.74	12.04 ± 5.05	2.98 ± 0.11	9.06 ± 0.45	55.96 ± 5.07

The surface density profiles of each lens’s preferred model are displayed in Fig. 3. The Sersic component is displayed in red, the Exp component in green, and the NFW component in black. The vertical blue lines indicate the region in which light from the lensed source is observed (i.e. in which the lensing analysis directly constrains the mass profile). These plots are derived using circular apertures centred on each component of the mass model, therefore omitting any geometric offsets. For each lens, the first component dominates the inner (<2 kpc) regions, with the second component beginning to dominate around the Einstein radius at 5–10 kpc. At larger radii (>15 kpc), where our analysis becomes an extrapolation, the dark matter begins to dominate. Of the total stellar mass, the Exp component is the dominant component in SLACS1 and SLACS3, contributing |$80{{\ \rm per\ cent}}$| (⁠|$80{{\ \rm per\ cent}}$|⁠) and |$57{{\ \rm per\ cent}}$| (⁠|$75{{\ \rm per\ cent}}$|⁠) within 10 kpc (500 kpc), respectively. For SLACS2, it makes up just |$35{{\ \rm per\ cent}}$| (⁠|$21{{\ \rm per\ cent}}$|⁠) within 10 kpc (500 kpc).

Figure 3.

One-dimensional surface density profiles of SLACS1 (left), SLACS2 (middle), and SLACS3 (right), showing the stellar mass density profile of the Sersic component (red), Exp component (green), and the NFW dark matter profile (black). Densities are estimated using circular apertures centred on each lens galaxy. The dashed vertical blue lines show the Einstein radius of each lens and the thick blue lines the radial extent of where the lensed source’s light is visible. These indicate where the lensing data directly constrain the mass model. The density profile of the most probable model is given by the solid lines. To compute the models at 3σ errors, given by the dashed lines, the density profiles of all models sampled by multinest are weighted by their sampling probability and averaged so as to include the covariance between different models.

4.6 The role of lensing

Fig. 4 shows the residuals of the highest likelihood Sersic model and chosen Sersic + Exp model, for the decomposed-mass profiles. For all three lenses, the chosen two-component models can be seen to improve the fit in the central regions of the lens galaxy, contributing to the increase in Bayesian evidence. However, in SLACS1 and SLACS3, improvements can also be seen visually in the residuals of the lensed source reconstruction, indicating that the strong lensing analysis is also contributing to the choice of the two-component model. Although improvements cannot be seen ‘by eye’ in SLACS2, we will now show they are also present in the fit.

Figure 4.

A comparison of the residuals of the Sersic and final chosen Sersic + Exp models using the decomposed-mass profile for all three images. The colourmap of the residuals is plotted using a symmetric log normalization and is trimmed to a small range in values to emphasize the difference in residuals between the two models. For all three lenses, the lens light fit can be seen to improve for the Sersic + Exp model, and in SLACS1 and SLACS3 the source residuals also improve.

We now seek a more quantitative measure of how much of this signal can be attributed to lensing. We define |$\chi ^2_{\mathrm{ Src}}$|⁠, the sum of χ² values in pixels where the lensed source is located (source contribution values >0.2, see N18), as a measure of the goodness of fit for the source. This comprises a different number of image pixels for the total-mass and decomposed-mass profile fits, therefore we renormalize each |$\chi ^2_{\mathrm{ Src}}$| value to that of the total-mass profile by multiplying by the ratio of number of pixels in each fit. We convert this to |$\mathrm{ln}\, \epsilon _{\mathrm{ Src}} = -0.5\chi ^2_{\mathrm{ Src}}$|⁠, the likelihood contribution of the source fit. The values of |$\mathrm{ln}\, \epsilon _{\mathrm{ Src}}$| quoted correspond to the highest likelihood model at the end of each analysis, although there is little variation between models.

pyautolens’s Bayesian framework also measures the evidence of a model via a set of source-plane regularization terms, which quantify the complexity of the source reconstruction (such that worse mass models correspond to a more complex and lower evidence solution). χ² is therefore not the only term that should be inspected to address how much the lensing analysis has improved. However, we find that the regularization terms are driven by the quality of the lens subtraction, as opposed to the mass model, thus we omit them. Therefore, we define |$\Delta \mathrm{ln}\, \epsilon _{\mathrm{ Src}} = \mathrm{ln}\, \epsilon _{\mathrm{ Src},\mathrm{ f}} - \mathrm{ln}\, \epsilon _{\mathrm{ Src},\mathrm{ i}}$| (i.e. final preferred model − initial Sersic model), to describe how much the lensing data contribute to the increase in evidence for a more complex model. Finally, we compare this value for the decomposed- and total-mass profiles, to ensure their different light subtractions do not impact our comparison.

For the total-mass profiles, we obtain |$\Delta \mathrm{ln}\, \epsilon _{\mathrm{ Src}} = -7$| for SLACS1, |$\Delta \mathrm{ln}\, \epsilon _{\mathrm{ Src}} = 9$| for SLACS2, and |$\Delta \mathrm{ln}\, \epsilon _{\mathrm{ Src}} = 50$| for SLACS3. Since we use the same SIE mass profile throughout, the source reconstruction only improves as a consequence of an improved lens subtraction leaving fewer residuals in the regions of the image the source is located, which is the case for SLACS3. For SLACS1 and SLACS2, the single Sersic profile produced a sufficiently good lens subtraction over these regions. For the decomposed models, we obtain |$\Delta \mathrm{ln}\, \epsilon _{\mathrm{ Src}} = 116$|⁠, |$\Delta \mathrm{ln}\, \epsilon _{\mathrm{ Src}} = 61$|⁠, and |$\Delta \mathrm{ln}\, \epsilon _{\rm{Src}} = 152$| for SLACS1, SLACS2, and SLACS3 respectively. Thus, more complex light profiles are clearly necessary to better fit the lensing data.

4.7 Physical plausibility and Bayesian evidence

The Bayesian evidence is a powerful diagnostic in determining the lens model complexity. However, caution must be exercised in relying on just the increase in evidence as a means to determine the physical plausibility of a lens model. A model may fit the data better, and therefore give a large Bayes factor, but the actual reason for this may be a different scenario physically. Comparing different models, with different assumptions, is a useful way to test for this. For instance, offsets in the rotation angles and centroids of the Sersic and Exp components were seen for both mass model assumptions (total and decomposed), increasing our confidence that this is a physically genuine component of the lens model.

There are models where a large Bayes factor does not definitely rule out another physically plausible scenario. An example of this is our detection of a radial gradient in the stellar mass distribution. This model was chosen because the data required a more shallow mass distribution in the central regions of the lens galaxy of SLACS2, in a way our assumed dark matter profiles (including a gNFW) could not produce. It is plausible a more flexible dark matter model capable of removing more mass from the centre (e.g. a cored profile) would produce the same Bayes factor. It would be difficult to distinguish between these possibilities using just our lensing data and demonstrates a limitation of the method.

5 DISCUSSION

5.1 The formation and structure of massive elliptical galaxies

We quantified the structure of three massive elliptical (ME) strong lenses using two different approaches. In the first approach, we fit a light profile to each galaxy that was neither constrained by nor used to constrain the mass model of the same galaxy that we obtained from a strong lensing analysis. Two-component light profiles were preferred over single-component light profiles with this approach for all the three MEs, with Bayes factors of the order of 200–500 (translating to a ≥10σ detection). The overall light profiles of the one- and two-component fits were visibly similar, signifying that the improvement in the fit associated with the two-component model was spread over many thousands of image pixels. Nevertheless, detailed inspection of the image residuals revealed that improvements from the two-component model could be seen in the fit to the lens galaxy.

The Sersic component of each lens was more compact (lower R_l), rounder (higher q_l), and more concentrated (higher n_l) than the exponential component. Photometric decompositions of (non-lensing) MEs have inferred the same two-component (or more) structure of MEs (Huang et al. 2013; Oh et al. 2016). They argue that this supports the forward models of ME formation (e.g. two-phase formation) outlined in the introduction, where the compact Sersic component is an old ‘red nugget’ core of each ME formed in a rapid and dissipative event at high redshift (z > ∼3). The second component is a surrounding envelope of stellar material built via passive growth of mergers thereafter. However, with photometry alone, detecting this decomposition relies on an inflexion between the surface photometry of the two components being observable (e.g. Oh et al. 2016) or a χ² improvement in the fit (e.g. D’Souza et al. 2014). Due to projection effects, it is not clear such an inflexion is always observable and χ² thresholds are challenging to quantify robustly due to the potential of overfitting.

By incorporating the lens’s light profile into the mass model, such that it is simultaneously constrained by the lensing analysis, we are able to independently confirm whether these two components are physically distinct. Incorrect fits due to projection effects or overfitting must lead to a degradation of the lensing analysis. For all three lenses, the two-component light model (now incorporated into the mass model) improved the lensing analysis compared to the single-component fit. Therefore, lensing confirms that our three MEs are comprised of (at least) two physically distinct components that are consistent with forward models of ME formation (e.g. Cooper et al. 2013, 2015).

As outlined in the introduction, the outer envelope that we model as an individual exponential profile is anticipated to consist of many separate components, each corresponding to an accretion event in the MEs past. Therefore, what we model as an exponential should be viewed as the superposition of these subcomponents where its inferred properties represent some average of all of the (most significant) accretion events in the MEs’ history. Attempting to decompose this component further, into two or more components, is beyond the scope of this work, but could necessitate our lensing-based approach, given that the blurring of these components in the photometric imaging likely makes them inseparable via just light profile fitting.

Our limited sample of just three objects means that a more rigorous comparison between observations and theory is not yet warranted. We therefore turn our attention to what the inclusion of strong lensing provides and argue that, once we have analysed a sufficiently large sample, lensing will provide information vital to distinguishing models of ME formation.

5.2 What does lensing observe?

Lensing reveals two crucial pieces of additional information about each ME that are not measured when using photometry alone. The first is the mass-to-light profile of each stellar component of the galaxy. This does not require stellar population models, circumventing the debated form of the initial mass function (Treu et al. 2009) and could of course act an independent test of the IMF alongside studies such as Auger et al. (2010,) and Sonnenfeld et al. (2015). In one lens, the two components had different mass-to-light ratios (SLACS3) and in another we observed a radial gradient (SLACS2 – see Sonnenfeld et al. 2018 for radial gradients in other SLACS lenses). The mass-to-light profiles for SLACS2 and SLACS3 removed mass from the central component relative to its outer envelope. This could indicate that internal baryonic processes (e.g. AGN feedback) can play a role in reshaping the stellar mass profile, or be a signature of supermassive black hole scouring (Faber et al. 1996; Trujillo et al. 2004; Merritt 2006). Larger samples are necessary to investigate this further.

The second, and perhaps most crucial piece of information lensing measures, is the central mass of each ME’s dark matter halo. Having a direct measurement of each galaxy’s halo mass means that (using large samples) we can tie each galaxy (statistically) to its initial formation in the early Universe (e.g. via halo merger trees). This can circumvent selections bias (e.g. progenitor bias) and offers a direct means by which to test ME formation in cosmological simulations. Combined with our inferences on the stellar density profile, we can distinguish what role an ME’s environment (and thus merger history) plays in shaping its stellar distribution compared to internal baryonic processes.

5.3 Structural geometry

The ‘outer envelope’ of each ME was more elliptical than the central component (q_l = 0.65–0.8) and rotationally offset in projection, with two of these offsets exceeding 70°. Isophotal twists are seen in other works (Huang et al. 2013; Oh et al. 2016), but rarely at such large angles. This is likely because most studies use lower mass galaxies; larger rotational offsets are indeed seen by Goullaud et al. (2018) in ellipticals of similar masses to ours. These offsets are consistent with models of ME formation, provided that: (i) the outer envelope aligns with the direction of preferential merger accretion in its past (i.e. its surrounding environment) and (ii) the central nugget’s orientation can be decoupled from this direction during its formation. This decoupling is observed in the EAGLE simulations (Velliscig et al. 2015) and explained by a highly energetic event at z > 3 (e.g. a major merger) decoupling the galaxies’ central baryonic structure from its surrounding environment, leaving the baryonic material accreted afterwards to trace it. More detailed theoretical studies are warranted, but large rotational offsets may hold an imprint of an ME’s high-redshift formation, for example the merger ratio, direction, and abundance of cold gas (Hopkins et al. 2009).

In two lenses, the centre of the two stellar components was observed to be offset, by ∼100 pc (SLACS2) and ∼300 pc (SLACS3). Whilst these offsets were seen in the light profile fit, without lensing it would have been unclear whether or not they were due to dust in the central regions of the galaxy (Goullaud et al. 2018). These offsets are not two bimodal peaks of material in the central regions of the galaxy, but instead caused by lopsidedness in the stellar mass distribution at large radii (e.g. >10 kpc; Gomer & Williams 2018). Thus, the offsets are tracing how recently material was accreted on to the ME’s outskirts, given that material accreted more recently will not have had sufficient time to settle on to orbits centred on its overall potential well. Centroid offsets could in principle be used to estimate the time since a previous (major) merger; however, the details are likely very complicated and depend on the merger in-fall velocities, directions, and the galaxy’s mass distribution. Unfortunately, we cannot be guided by simulations due to their softening lengths being above 300 pc (e.g. Schaller et al. 2015).

5.4 Implications for strong lensing

Whilst our focus has been on showing how lensing can benefit photometric studies of galaxies, the reverse is also true. That is, for lensing systems in which the source does not provide a large amount of constraining information (e.g. it is only doubly imaged, compact, or point-like), fitting the lens’s light will better constrain its mass profile (given adequate assumptions about the mass-to-light ratio profile). Furthermore, the freedom of well-known lensing degeneracies like the mass-sheet transformation (Falco, Gorenstein & Shapiro 2002; Schneider & Sluse 2013; Schneider 2014) and source-position transformation (Wertz & Orthen 2018) will be restricted, as the lens’s stellar mass profile must coincide with the photometric data. These degeneracies also assume axisymmetry and may be alleviated by the geometric offsets detected in this work.

There are a number of techniques aiming to use strong lensing as a precise probe of cosmology (e.g. Oguri et al. 2012; Suyu et al. 2013; Collett & Auger 2014). These methods require an accurate mass model for the lens. The mass profiles typically used in these studies do not capture the rotational and centroid offsets observed in this work (although see Birrer et al. 2019), which have a unique impact on ray tracing and could potentially impact the cosmological inference. Once we have analysed larger samples, we will investigate what impact, if any, these assumptions have on the inferred cosmology. Simultaneously fitting the lens’s light and mass may also improve attempts to detect dark matter substructures via strong lensing (e.g. Vegetti et al. 2014), where baryonic structures like an inclined disc mimic the substructure signal (Hsueh et al. 2016, 2017, 2018).

Investigating the ‘Fundamental Surface of Quads’, Gomer & Williams (2018) showed that the positions of quadruply imaged quasars provide evidence for rotational and centroid offsets in the mass distributions of lens galaxies, consistent with the models we have presented here. However, the authors only investigated offsets between the stellar and dark matter components in galaxies. Based on our work, we would suggest that their inferences may instead be due to offsets between the different baryonic components of their lens galaxies. However, the selection function of strongly lensed quasars leads to a different lens galaxy population and further study is necessary to confirm this.

5.5 Galaxy structure with strong lensing

Strong lensing promises to be a powerful technique to study galaxy structure. Euclid, LSST, and other surveys are set to find of the order of 100 000 strong lenses in the next decade (Collett 2015), spanning the entire Hubble diagram. Galaxies of later-type morphology will produce more distinct and thus measurable lensing effects than MEs, due to their distinct bulge–disc structure (Orban De Xivry & Marshall 2009; Hsueh et al. 2017). Whereas this study focused on each galaxy’s two most dominant stellar components, future studies can investigate how sensitive lensing is to smaller galaxy structures (e.g. bars and pseudo-bulges), which are difficult to separate photometrically from the bulk of the galaxy’s light. The light from structures that dominate well outside the Einstein radius of a lens, such as the stellar halo (D’Souza et al. 2014; Spavone et al. 2017), will also be observable in future surveys and, in conjunction with weak lensing, will provide further constraints on the distribution of stellar and dark mass.

Our technique complements currently ongoing large IFS surveys like the Mapping Nearby Galaxies at Apache Point Observatory (Yan et al. 2016). These surveys are taking resolved spectra for 10 000 galaxies, enabling the comparison of their kinematic mass structure and photometry. However, these IFS surveys are restricted to the local Universe (e.g. z < 0.02), whereas lensing will span a broad range in lens galaxy redshifts (from z ∼ 0.05 to z > 3). Thus, evolutionary studies of galaxy structure are only possible with lensing.

6 SUMMARY

We fitted three strong gravitational lenses with the lens modelling software pyautolens. For each system, we compared two mass models: (i) a total-mass model, where the lens galaxy’s structure is inferred independently from strong lensing and (ii) a decomposed-mass model, where the lens galaxy’s structure is folded into the lensing analysis and thus constrained by it. This allowed us to compare whether the structures we see in each lens galaxy’s light translate to physically genuine mass structures that are necessary to fit the lensing data.

All three lens galaxies were massive ellipticals (ME). For each, we inferred that they are structurally composed of two components (a double Sersic model) as opposed to one component (a single Sersic model). The two-component model consists of a compact and round central nugget, surrounded by a flatter and more extended outer envelope. The two-component model was preferred, irrespective of whether the lensing constraints were used. However, including lensing increased the significance of each detection and thus confirmed these structures are physically genuine (as opposed due to overfitting, dust, or projection effects). Furthermore, the inclusion of lensing enabled us to measure two quantities that a photometric-only analysis cannot: (i) the stellar mass distribution (without stellar population models) and (ii) the inner dark matter halo mass of each galaxy.

We interpreted this two-component model in the context of forward models of ME galaxies (Cooper et al. 2013, 2015). The central components are consistent with a central ‘red nugget’ (e.g. Trujillo et al. 2006; Oldham et al. 2017) formed at high redshift (z > 2) due to a highly dissipative event, as evidenced by their structural parameters (high Sersic index, low effective radius). The outer components are extended envelopes of material (low Sersic index, high effective radius) accreted via mergers from redshift 2 onwards. With just three objects, a rigorous test of ME formation was not warranted. However, by directly measuring dark matter masses, in the future we will be able to statistically tie a large sample of ME’s to their counterparts in simulations and use their inferred stellar density profiles to decouple what role baryonic feedback processes and accretion play in their evolution. Thus, strong lensing provides a new test of ME formation, especially given recent theoretical works highlighting correlations between halo mass and stellar density profile.

We investigated the geometry of these two components. We detect that they are rotationally offset from one another in projection, in two cases by over 70°. This suggests that the dissipative event that forms the central nugget changes its orientation relative to its surrounding local environment, which the outer envelope (being aligned with the direction of preferential accretion) traces. Such offsets are consistent with the EAGLE simulations (Velliscig et al. 2015). Centroid offsets between each component were also detected (∼200 kpc) and interpreted as a lopsidedness in the mass distribution of each galaxy (Gomer & Williams 2018). This provides information on when the last episode of significant accretion took place, but may be challenging to extract because these offsets are below the softening length of modern cosmology simulations (Schaller et al. 2015).

Over the next decade, surveys such as Euclid and LSST are set to amass samples in excess of 100 000 strong lenses, spanning the entire Hubble diagram from redshifts 0 to 3 and beyond. The analysis demonstrated in this work offers a completely new view of a galaxy’s structure, which is not limited to the ME galaxies investigated here. Direct measurements of a galaxy’s host dark matter halo are crucial to understanding each galaxy’s place in hierarchical structure formation, whilst their stellar density profiles provide insight into the role baryonic physics has in shaping each galaxy. In the future, we will seek to understand how lensing can further our understanding of the formation and evolution of galaxies of any Hubble type.

ACKNOWLEDGEMENTS

We thank the anonymous referee whose suggestions improved the quality of the work. The authors thank Ian Smail, Russell Smith, Simon Dye, and Andrew Robertson for helpful conversations. JWN is supported by the UK Space Agency, through grant ST/N001494/1. RJM is supported by a Royal Society University Research Fellowship. APC is supported by the UK Science and Technology Facilities Council, through grant ST/P000541/1.

SLACS1 and SLACS3 were both observed in HST programme GO-10886 and SLACS2 in programme GO-10494. This work used the DiRAC Data Centric system at Durham University, which is operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National e-Infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the UK National e-Infrastructure.

Footnotes

1

https://github.com/Jammy2211/PyAutoLens

2

These expectations for accretion on to MEs are significantly different to those for Milky Way-like galaxies, which correspond to the regime of in situ-dominated growth in lower mass haloes largely unaffected by AGN feedback. In that regime, the analogous accreted stellar component (the stellar halo) is expected to be dominated by only one or two accretion events of high mass ratio, and to comprise only a small fraction of the total mass. Stellar halo progenitors arrive on more circular orbits (giving rise to streams and other unrelaxed structures) and consist of stellar populations very different to those of the central galaxy (Bullock & Johnston 2005; Amorisco 2015; De Las Aguas Robustillo Cortés, González & Verdugo 2017).

3

The pyautolens software is open-source and available from https://github.com/Jammy2211/PyAutoLens. The results in this work were computed using an earlier Fortran build of autolens.

APPENDIX A: LENS MODELS

The lens light profile intensities I and mass model deflection angles α_x,y (via integration of the convergence κ) are computed using the elliptical coordinate system |$\kappa (\vec{\xi })$|⁠, where |$\xi = \sqrt{{x}^2 + y^2/q^2}$|⁠. Parameters associated with the lens’s total-mass profile have no subscript, light matter profile the subscript ‘l’, dark matter profile the subscript ‘d’, and external shear the subscript ‘sh’. The model centres, axial ratios, and rotation angles are, in most cases, treated as free parameters. Multiple component models are computed by summing each individual component’s intensities (for a light profile) or deflection angles (for a mass profile) and their parameters are given matching numerical subscripts (e.g. x₁, x_l2, etc.). The light profile intensities are computed using an adaptive oversampling routine and mass model deflection angles an adaptive numerical integrator and bi-linear interpolation, which are described in N18. The lens model changes throughout the analysis pipeline and Section 4 describes how the final lens model is chosen.

The lens’s light profile is modelled using one or more elliptical Sersic profiles

$$\begin{eqnarray*} I_{\rm Ser} (\xi _{\rm l}) = I_{\rm l} \exp \bigg \lbrace -k_{\rm l} \bigg [ \bigg (\frac{\xi _{\rm l}}{R_{\rm l}} \bigg)^{\frac{1}{n_{\rm l}}} - 1 \bigg ] \bigg \rbrace , \end{eqnarray*}$$

(A1)

which have up to seven free parameters: (x_l, y_l), the light centre, q_l, the axial ratio, θ_l, the orientation angle, I_l, the intensity at the effective radius R_l, and n_l, the Sersic index. k_l is a function of n_l. The exponential light profile I_Exp(ξ_l) corresponds to n_l = 1.

For total-mass modelling, a softened power-law ellipsoid (SPLE) density profile of the form

$$\begin{eqnarray*} \kappa _{\rm pl} (\xi) = \frac{(3 - \alpha)}{1 + q} \bigg (\frac{\theta _{\rm E}}{\xi } \bigg)^{\alpha - 1} \end{eqnarray*}$$

(A2)

is used, where θ_E is the model Einstein radius in arc seconds. The power-law density slope is α, and setting α = 2 gives the SIE model. The inclusion of an external shear field is supported, which introduces two additional parameters, the shear magnitude γ_sh and orientation of the semimajor axis measured counterclockwise from north, θ_sh.

For decomposed-mass modelling, the Sersic profile given by equation (A1) is used to give the light matter surface density profile

$$\begin{eqnarray*} \kappa _{\mathrm{ Ser}} (\xi _l) = \Psi _l \bigg [\frac{q_{\rm l} \xi _{\rm l}}{R_{\rm l}}\bigg ]^{\Gamma _{\rm l}} I_{\mathrm{ Ser}} (\xi _l) \, \, , \end{eqnarray*}$$

(A3)

where Ψ_l gives the mass-to-light ratio in electrons per second and Γ_l can fold a radial dependence into the conversion of mass to light. The assumption used in N18 of a constant mass-to-light ratio is given for Γ_l = 0 and models that assume this value or include Γ_l as a free parameter are investigated. If there are multiple light profile components, these can either assume a single shared value of Ψ_l or two independent values. The exponential convergence profile κ_Exp(ξ_l) corresponds to n_l = 1.

The dark matter component is given by a generalized Navarro–Frenk–White (NFW) profile, representing the universal density profile predicted for dark matter haloes by cosmological N-body simulations (Navarro, Frenk & White 1996; Zhao 1996; Moore et al. 1998) and with a volume mass density given by

$$\begin{eqnarray*} \rho = \frac{\rho _\mathrm{ s}}{(r/r_\mathrm{ s})^\gamma (1 + r/r_\mathrm{ s})^{3-\gamma }}. \end{eqnarray*}$$

(A4)

The halo normalization is given by ρ_s, and r_s, the scale radius, is fixed to the value r_s = 10R_l (Bullock et al. 2001) and the NFW profile corresponds to the solution γ_d = 1. Coordinates for the NFW profile are scaled by r_s, giving the scaled elliptical coordinate η_d = ξ_d/r_s. Deflection angles for spherical NFW models use the analytic solutions given by Golse & Kneib (2002) and N18, whereas elliptical models integrate using the formalism of Keeton, Gaudi & Petters (2005).

APPENDIX B: LENS MODELLING RESULTS

The most probable lens models discussed in Section 4 and their 3σ errors are given in Tables B1–B4.

Table B1.

The inferred most probable parameters (with 3σ errors) of the Sersic profile, using a total-mass Sersic + SIE model (rows 1, 3, and 5) and decomposed Sersic + NFWSph model (rows 2, 4, and 6). The inferred parameters of the Sersic profile are consistent for both the models, for all three lenses.

Target name	Mass profile	Mass-to-light ratio	Parameters
|$\mathbf {SLACSJ0252\!+\!0039}$|	Total	N/A	|$R_{\mathrm{l}} = 1.92^{+0.57}_{-0.54}$|	|$n_{\mathrm{l}} = 3.90^{+0.47}_{-0.47}$|	|$q_{\mathrm{l}} = 0.940^{+0.028}_{-0.029}$|	|$\theta _{\mathrm{l}} = 92^{+20}_{-20}$|
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	N/A	|$R_{\mathrm{l}} = 1.78^{+0.42}_{-0.28}$|	|$n_{\mathrm{l}} = 3.60^{+0.40}_{-0.32}$|	|$q_{\mathrm{l}} = 0.917^{+0.025}_{-0.021}$|	|$\theta _{\mathrm{l}} = 94.4^{+8.9}_{-7.6}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	Total	N/A	|$R_{\mathrm{l}} = 2.10^{+0.50}_{-0.50}$|	|$n_{\mathrm{l}} = 5.31^{+0.49}_{-0.54}$|	|$q_{\mathrm{l}} = 0.883^{+0.016}_{-0.017}$|	|$\theta _{\mathrm{l}} = 71.2^{+5.1}_{-4.3}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	N/A	|$R_{\mathrm{l}} = 2.11^{+0.61}_{-0.55}$|	|$n_{\mathrm{l}} = 5.31^{+0.58}_{-0.52}$|	|$q_{\mathrm{l}} = 0.882^{+0.020}_{-0.020}$|	|$\theta _{\mathrm{l}} = 70.9^{+4.8}_{-4.4}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	Total	N/A	|$R_{\mathrm{l}} = 1.93^{+0.26}_{-0.27}$|	|$n_{\mathrm{l}} = 4.73^{+0.28}_{-0.29}$|	|$q_{\mathrm{l}} = 0.924^{+0.012}_{-0.011}$|	|$\theta _{\mathrm{l}} = 87.0^{+4.5}_{-4.7}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	N/A	|$R_{\mathrm{l}} = 2.41^{+0.29}_{-0.35}$|	|$n_{\mathrm{l}} = 5.13^{+0.25}_{-0.30}$|	|$q_{\mathrm{l}} = 0.9019^{+0.0084}_{-0.0078}$|	|$\theta _{\mathrm{l}} = 85.6^{+2.4}_{-2.5}$|

Table B2.

The inferred most probable parameters (with 3σ errors) of the Sersic + Exp profiles with geometric parameters aligned between the two components. The results are for a total-mass Sersic + Exp + SIE model (rows 1, 3, and 5) and decomposed Sersic + Exp + NFWSph model (rows 2, 4, and 6) where the mass-to-light ratios of the two components are the same. The inferred parameters of the Sersic and Exp profiles are consistent for both the models, for all three lenses.

Target name	Mass profile	Mass-to-light ratio	Parameters
|$\mathbf {SLACSJ0252\!+\!0039}$|	Total	N/A	|$R_{\mathrm{l1}} = 0.199^{+0.060}_{-0.046}$|	|$n_{\mathrm{l1}} = 2.81^{+0.87}_{-0.76}$|	|$q_{\mathrm{l1}} = 0.976^{+0.032}_{-0.043}$|	|$\theta _{\mathrm{l1}} = 93^{+10}_{-10}$|
–	–	–	|$R_{\mathrm{l2}} = 1.01^{+0.12}_{-0.11}$|	|$q_{\mathrm{l2}} = 0.863^{+0.060}_{-0.056}$|	–	–
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Single	|$R_{\mathrm{l1}} = 0.183^{+0.050}_{-0.045}$|	|$n_{\mathrm{l1}} = 2.68^{+0.73}_{-0.70}$|	|$q_{\mathrm{l1}} = 0.984^{+0.020}_{-0.026}$|	|$\theta _{\mathrm{l1}} = 97.5^{+7.9}_{-7.9}$|
–	–	–	|$R_{\mathrm{l2}} = 0.96^{+0.11}_{-0.11}$|	|$q_{\mathrm{l2}} = 0.841^{+0.035}_{-0.048}$|	–	–
|$\mathbf {SLACSJ1250\!+\!0523}$|	Total	N/A	|$R_{\mathrm{l1}} = 0.60^{+0.18}_{-0.16}$|	|$n_{\mathrm{l1}} = 3.61^{+0.57}_{-0.58}$|	|$q_{\mathrm{l1}} = 0.875^{+0.017}_{-0.018}$|	|$\theta _{\mathrm{l1}} = 71.0^{+4.2}_{-4.0}$|
–	–	–	|$R_{\mathrm{l2}} = 2.01^{+0.24}_{-0.21}$|	|$q_{\mathrm{l2}} = 0.977^{+0.032}_{-0.047}$|	–	–
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Single	|$R_{\mathrm{l1}} = 0.57^{+0.15}_{-0.14}$|	|$n_{\mathrm{l1}} = 3.50^{+0.73}_{-0.71}$|	|$q_{\mathrm{l1}} = 0.875^{+0.022}_{-0.024}$|	|$\theta _{\mathrm{l1}} = 71.6^{+5.0}_{-4.9}$|
–	–	–	|$R_{\mathrm{l2}} = 2.02^{+0.30}_{-0.26}$|	|$q_{\mathrm{l2}} = 0.978^{+0.031}_{-0.049}$|	–	–
|$\mathbf {SLACSJ1430\!+\!4105}$|	Total	N/A	|$R_{\mathrm{l1}} = 0.608^{+0.065}_{-0.063}$|	|$n_{\mathrm{l1}} = 2.95^{+0.21}_{-0.21}$|	|$q_{\mathrm{l1}} = 0.9329^{+0.0084}_{-0.0082}$|	|$\theta _{\mathrm{l1}} = 81.5^{+2.5}_{-2.6}$|
–	–	–	|$R_{\mathrm{l2}} = 2.87^{+0.22}_{-0.22}$|	|$q_{\mathrm{l2}} = 0.710^{+0.038}_{-0.043}$|	–	–
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Single	|$R_{\mathrm{l1}} = 0.571^{+0.055}_{-0.068}$|	|$n_{\mathrm{l1}} = 2.727^{+0.19}_{-0.20}$|	|$q_{\mathrm{l1}} = 0.928^{+0.007}_{-0.008}$|	|$\theta _{\mathrm{l1}} = 81.8^{+2.4}_{-2.7}$|
–	–	–	|$R_{\mathrm{l2}} = 2.65^{+0.19}_{-0.24}$|	|$q_{\mathrm{l2}} = 0.674^{+0.033}_{-0.041}$|	–	–

Table B3.

The inferred most probable parameters (with 3σ errors) of the Sersic + Exp profiles with the alignment of geometric parameters corresponding to the highest evidence model chosen in the light profile geometry phase (Table 6). The results are for a total-mass Sersic + Exp + SIE model (rows 1, 4, and 7), decomposed Sersic + Exp + NFWSph model (rows 2, 5, and 8) with shared mass-to-light ratio and independent mass-to-light ratios (rows 3, 6, and 9) are given.

Target name	Mass profile	Mass-to-light ratio	Parameters
|$\mathbf {SLACSJ0252\!+\!0039}$|	Total	N/A	|$R_{\mathrm{l1}} = 0.144^{+0.033}_{-0.028}$|	|$n_{\mathrm{l1}} = 2.26^{+0.71}_{-0.65}$|	|$q_{\mathrm{l1}} = 0.907^{+0.053}_{-0.053}$|	|$\theta _{\mathrm{l1}} = 41^{+10}_{-20}$|
–	–	–	|$R_{\mathrm{l2}} = 0.879^{+0.11}_{-0.095}$|	|$q_{\mathrm{l2}} = 0.830^{+0.049}_{-0.052}$|	|$\theta _{\mathrm{l2}} = 105.0^{+9.2}_{-9.4}$|	|$\theta = 99^{+10}_{-10}$|
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Single	|$R_{\mathrm{l1}} = 0.141^{+0.036}_{-0.030}$|	|$n_{\mathrm{l1}} = 2.36^{+0.76}_{-0.64}$|	|$q_{\mathrm{l1}} = 0.889^{+0.059}_{-0.057}$|	|$\theta _{\mathrm{l1}} = 38^{+20}_{-10}$|
–	–	–	|$R_{\mathrm{l2}} = 0.811^{+0.12}_{-0.093}$|	|$q_{\mathrm{l2}} = 0.824^{+0.042}_{-0.044}$|	|$\theta _{\mathrm{l2}} = 106.2^{+8.0}_{-8.3}$|	–
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Independent	|$R_{\mathrm{l1}} = 0.165^{+0.048}_{-0.042}$|	|$n_{\mathrm{l1}} = 2.90^{+0.93}_{-0.87}$|	|$q_{\mathrm{l1}} = 0.894^{+0.062}_{-0.070}$|	|$\theta _{\mathrm{l1}} = 36^{+10}_{-10}$|
–	–	–	|$R_{\mathrm{l2}} = 0.81^{+0.13}_{-0.11}$|	|$q_{\mathrm{l2}} = 0.818^{+0.047}_{-0.051}$|	|$\theta _{\mathrm{l2}} = 105.9^{+8.5}_{-8.4}$|	–
|$\mathbf {SLACSJ1250\!+\!0523}$|	Total	N/A	|$R_{\mathrm{l1}} = 0.99^{+0.29}_{-0.28}$|	|$n_{\mathrm{l1}} = 4.48^{+0.60}_{-0.53}$|	|$q_{\mathrm{l1}} = 0.882^{+0.020}_{-0.020}$|	|$\theta _{\mathrm{l1}} = 70.9^{+4.7}_{-4.4}$|
–	–	–	|$R_{\mathrm{l2}} = 2.11^{+0.31}_{-0.28}$|	|$q_{\mathrm{l2}} = 0.75^{+0.18}_{-0.20}$|	|$\theta _{\mathrm{l2}} = 163^{+18}_{-14}$|	|$\theta = 104^{+20}_{-10}$|
–	–	–	|$x_{\mathrm{l1}} = 0.0202^{+0.0028}_{-0.0024}$|	|$y_{\mathrm{l1}} = 0.0098^{+0.0030}_{-0.0027}$|	|$x_{\mathrm{l2}} = 0.049^{+0.025}_{-0.024}$|	|$y_{\mathrm{l2}} = 0.022^{+0.027}_{-0.027}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Single	|$R_{\mathrm{l1}} = 0.61^{+0.23}_{-0.21}$|	|$n_{\mathrm{l1}} = 3.62^{+0.50}_{-0.44}$|	|$q_{\mathrm{l1}} = 0.874^{+0.019}_{-0.020}$|	|$\theta _{\mathrm{l1}} = 72.0^{+4.3}_{-4.6}$|
–	–	–	|$R_{\mathrm{l2}} = 2.06^{+0.33}_{-0.24}$|	|$q_{\mathrm{l2}} = 0.920^{+0.14}_{-0.15}$|	|$\theta _{\mathrm{l2}} = 149^{+16}_{-17}$|	–
–	–	–	|$x_{\mathrm{l1}} = 0.021^{+0.0021}_{-0.0025}$|	|$y_{\mathrm{l1}} = 0.010^{+0.0029}_{-0.0027}$|	|$x_{\mathrm{l2}} =0.0361^{+0.0024}_{-0.0025}$|	|$y_{\mathrm{l2}} = 0.0277^{+0.0031}_{-0.0028}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Independent	|$R_{\mathrm{l1}} = 1.07^{+0.28}_{-0.26}$|	|$n_{\mathrm{l1}} = 4.58^{+0.56}_{-0.62}$|	|$q_{\mathrm{l1}} = 0.880^{+0.020}_{-0.018}$|	|$\theta _{\mathrm{l1}} = 71.2^{+4.2}_{-4.4}$|
–	–	–	|$R_{\mathrm{l2}} = 2.11^{+0.44}_{-0.32}$|	|$q_{\mathrm{l2}} = 0.70^{+0.11}_{-0.10}$|	|$\theta _{\mathrm{l2}} = 162^{+16}_{-12}$|	–
–	–	–	|$x_{\mathrm{l1}} = 0.0201^{+0.0025}_{-0.0024}$|	|$y_{\mathrm{l1}} = 0.0091^{+0.0025}_{-0.0024}$|	|$x_{\mathrm{l2}} = 0.056^{+0.028}_{-0.032}$|	|$y_{\mathrm{l2}} = 0.036^{+0.022}_{-0.026}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	Total	N/A	|$R_{\mathrm{l1}} = 0.571^{+0.058}_{-0.043}$|	|$n_{\mathrm{l1}} = 2.873^{+0.612}_{-0.530}$|	|$q_{\mathrm{l1}} = 0.924^{+0.021}_{-0.019}$|	|$\theta _{\mathrm{l1}} = 93.0^{+4.3}_{-4.6}$|
–	–	–	|$R_{\mathrm{l2}} = 2.83^{+0.24}_{-0.25}$|	|$q_{\mathrm{l2}} = 0.71^{+0.016}_{-0.021}$|	|$\theta _{\mathrm{l2}} = 69.0^{+2.0}_{-1.8}$|	|$\theta = 80.2^{+1.3}_{-0.9}$|
–	–	–	|$x_{\mathrm{l1}} = 0.028^{+0.0024}_{-0.0021}$|	|$y_{\mathrm{l1}} = 0.031^{+0.0033}_{-0.0027}$|	|$x_{\mathrm{l2}} = 0.003^{+0.021}_{-0.020}$|	|$y_{\mathrm{l2}} = -0.045^{+0.021}_{-0.022}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Single	|$R_{\mathrm{l1}} = 0.598^{+0.049}_{-0.063}$|	|$n_{\mathrm{l1}} = 2.89^{+0.16}_{-0.18}$|	|$q_{\mathrm{l1}} = 0.9214^{+0.0078}_{-0.0077}$|	|$\theta _{\mathrm{l1}} = 92.8^{+2.5}_{-2.8}$|
–	–	–	|$R_{\mathrm{l2}} = 2.88^{+0.16}_{-0.20}$|	|$q_{\mathrm{l2}} = 0.681^{+0.031}_{-0.028}$|	|$\theta _{\mathrm{l2}} = 70.3^{+1.9}_{-2.4}$|	–
–	–	–	|$x_{\mathrm{l1}} = 0.0283^{+0.0010}_{-0.0011}$|	|$y_{\mathrm{l1}} = 0.0307^{+0.0010}_{-0.0011}$|	|$x_{\mathrm{l2}} = -0.017^{+0.013}_{-0.016}$|	|$y_{\mathrm{l2}} = -0.0095^{+0.014}_{-0.018}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Independent	|$R_{\mathrm{l1}} = 0.587^{+0.042}_{-0.053}$|	|$n_{\mathrm{l1}} = 2.86^{+0.14}_{-0.15}$|	|$q_{\mathrm{l1}} = 0.922^{+0.0076}_{-0.0074}$|	|$\theta _{\mathrm{l1}} = 92.0^{+2.2}_{-2.3}$|
–	–	–	|$R_{\mathrm{l2}} = 2.80^{+0.13}_{-0.17}$|	|$q_{\mathrm{l2}} = 0.683^{+0.029}_{-0.029}$|	|$\theta _{\mathrm{l2}} = 70.0^{+1.7}_{-2.2}$|	–
–	–	–	|$x_{\mathrm{l1}} = 0.0280^{+0.0011}_{-0.0010}$|	|$y_{\mathrm{l1}} = 0.0306^{+0.0010}_{-0.0011}$|	|$x_{\mathrm{l2}} = -0.013^{+0.011}_{-0.015}$|	|$y_{\mathrm{l2}} = -0.002^{+0.013}_{-0.017}$|

Table B4.

The inferred most probable (with 3σ errors) inferred mass-to-light ratios Ψ and radial gradient Γ the Sersic + Exp profiles with the alignment of geometric parameters corresponding to the highest evidence model chosen in the light profile geometry phase (Table 6). The results are for decomposed Sersic + Exp + NFWSph model.

Target name	Mass profile	Mass-to-light ratio	Parameters
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Single	|$\Psi _{\mathrm{l}} = 5.26^{+0.083}_{-0.076}$|	–	–
|$\mathbf {SLACSJ0252\!+\!0039}$|	Decomposed	Independent	|$\Psi _{\mathrm{l1}} = 4.17^{+0.368}_{-0.332}$|	|$\Psi _{\mathrm{l2}} = 5.64^{+0.231}_{-0.246}$|	–
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Single	|$\Psi _{\mathrm{l}} = 1.27^{+0.062}_{-0.0085}$|	–	–
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Independent	|$\Psi _{\mathrm{l1}} = 1.66^{+0.266}_{-0.245}$|	|$\Psi _{\mathrm{l2}} = 1.00^{+0.197}_{-0.187}$|	–
|$\mathbf {SLACSJ1250\!+\!0523}$|	Decomposed	Radial Gradient	|$\Psi _{\mathrm{l1}} = 1.56^{+0.301}_{-0.311}$|	|$\Psi _{\mathrm{l2}} = 1.20^{+0.247}_{-0.228}$|	|$\Gamma _{\mathrm{l2}} = -0.24^{0.174}_{-0.154}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Single	|$\Psi _{\mathrm{l}} = 3.46^{+0.082}_{-0.096}$|	–	–
|$\mathbf {SLACSJ1430\!+\!4105}$|	Decomposed	Independent	|$\Psi _{\mathrm{l1}} = 3.37^{+0.570}_{-0.69}$|	|$\Psi _{\mathrm{l2}} = 7.71^{+0.386}_{-0.355}$|	–

Table B5.

The inferred most probable (with 3σ errors) inferred dark matter lensing masses κ_d, dark matter axial ratios q_d, dark matter inner density slopes γ_d, mass-to-light ratios Ψ, effective radii R_l, Sersic indexes n_l, and intensities I_l. These models use the Sersic + Exp profiles with the alignment of geometric parameters corresponding to the highest evidence model chosen in the light profile geometry phase (Table 6). The results are shown for three dark matter profiles, a spherical NFW NFWSph, an elliptical NFW NFWEll, and spherical generalized NFW gNFWSph.

Target name	Dark matter profile	Parameters
|$\mathbf {SLACSJ0252\!+\!0039}$|	NFWSph	|$\kappa _{\mathrm{d}} = 0.042^{+0.014}_{-0.011}$|	–	–	–
|$\mathbf {SLACSJ0252\!+\!0039}$|	–	|$\Psi _{\mathrm{l}} = 5.26^{+0.083}_{-0.076}$|	|$R_{\mathrm{l1}} = 0.141^{+0.036}_{-0.030}$|	|$n_{\mathrm{l1}} = 2.36^{+0.76}_{-0.64}$|	|$I_{\mathrm{l1}} = 0.448^{+0.125}_{-0.143}$|
|$\mathbf {SLACSJ0252\!+\!0039}$|	–	–	|$R_{\mathrm{l2}} = 0.811^{+0.12}_{-0.093}$|	–	|$I_{\mathrm{l2}} = 0.075^{+0.021}_{-0.029}$|
|$\mathbf {SLACSJ0252\!+\!0039}$|	NFWEll	|$\kappa _{\mathrm{d}} = 0.044^{+0.016}_{-0.017}$|	|$q_{\mathrm{d}} = 0.891^{+0.031}_{-0.034}$|	–	–
|$\mathbf {SLACSJ0252\!+\!0039}$|	–	|$\Psi _{\mathrm{l}} = 5.11^{+0.475}_{-0.412}$|	|$R_{\mathrm{l1}} = 0.186^{+0.065}_{-0.057}$|	|$n_{\mathrm{l1}} = 3.07^{+0.98}_{-0.75}$|	|$I_{\mathrm{l1}} = 0.282^{+0.062}_{-0.063}$|
|$\mathbf {SLACSJ0252\!+\!0039}$|	–	–	|$R_{\mathrm{l2}} = 0.885^{+0.19}_{-0.12}$|	–	|$I_{\mathrm{l2}} = 0.064^{+0.020}_{-0.024}$|
|$\mathbf {SLACSJ0252\!+\!0039}$|	gNFWSph	|$\kappa _{\mathrm{d}} = 0.041^{+0.025}_{-0.031}$|	|$\gamma _{\mathrm{d}} = 0.688^{+0.542}_{-0.578}$|	–	–
|$\mathbf {SLACSJ0252\!+\!0039}$|	–	|$\Psi _{\mathrm{l}} = 5.71^{+0.653}_{-0.736}$|	|$R_{\mathrm{l1}} = 0.166^{+0.042}_{-0.045}$|	|$n_{\mathrm{l1}} = 2.83^{+0.78}_{-0.70}$|	|$I_{\mathrm{l1}} = 0.333^{+0.072}_{-0.069}$|
|$\mathbf {SLACSJ0252\!+\!0039}$|	–	–	|$R_{\mathrm{l2}} = 0.820^{+0.22}_{-0.18}$|	–	|$I_{\mathrm{l2}} = 0.070^{+0.013}_{-0.009}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	NFWSph	|$\kappa _{\mathrm{d}} = 0.104^{+0.019}_{-0.020}$|	–	–	–
|$\mathbf {SLACSJ1250\!+\!0523}$|	–	|$\Psi _{\mathrm{l1}} = 1.66^{+0.266}_{-0.245}$|	|$R_{\mathrm{l1}} = 1.07^{+0.28}_{-0.26}$|	|$n_{\mathrm{l1}} = 4.58^{+0.56}_{-0.62}$|	|$I_{\mathrm{l1}} = 0.106^{+0.014}_{-0.016}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	–	|$\Psi _{\mathrm{l2}} = 1.00^{+0.197}_{-0.187}$|	|$R_{\mathrm{l2}} = 2.11^{+0.44}_{-0.32}$|	–	|$I_{\mathrm{l2}} = 0.025^{+0.008}_{-0.07}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	NFWEll	|$\kappa _{\mathrm{d}} = 0.090^{+0.022}_{-0.021}$|	|$q_{\mathrm{d}} = 0.976^{+0.021}_{-0.021}$|	–	|$I_{\mathrm{l2}} = 0.0^{+0.0}_{-0.0}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	–	|$\Psi _{\mathrm{l1}} = 1.72^{+0.376}_{-0.311}$|	|$R_{\mathrm{l1}} = 0.91^{+0.42}_{-0.39}$|	|$n_{\mathrm{l1}} = 4.24^{+0.76}_{-0.65}$|	|$I_{\mathrm{l1}} = 0.115^{+0.013}_{-0.017}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	–	|$\Psi _{\mathrm{l2}} = 1.11^{+0.326}_{-0.276}$|	|$R_{\mathrm{l2}} = 2.17^{+0.49}_{-0.33}$|	–	|$I_{\mathrm{l2}} = 0.025^{+0.007}_{-0.008}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	gNFWSph	|$\kappa _{\mathrm{d}} = 0.105^{+0.033}_{-0.031}$|	|$\gamma _{\mathrm{d}} = 0.846^{+0.540}_{-0.631}$|	–	–
|$\mathbf {SLACSJ1250\!+\!0523}$|	–	|$\Psi _{\mathrm{l1}} = 1.86^{+0.322}_{-0.365}$|	|$R_{\mathrm{l1}} = 0.91^{+0.44}_{-0.40}$|	|$n_{\mathrm{l1}} = 4.28^{+0.71}_{-0.66}$|	|$I_{\mathrm{l1}} = 0.114^{+0.016}_{-0.011}$|
|$\mathbf {SLACSJ1250\!+\!0523}$|	–	|$\Psi _{\mathrm{l2}} = 0.97^{+0.254}_{-0.289}$|	|$R_{\mathrm{l2}} = 2.21^{+0.48}_{-0.39}$|	–	|$I_{\mathrm{l2}} = 0.025^{+0.009}_{-0.008}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	NFWSph	|$\kappa _{\mathrm{d}} = 0.054^{+0.020}_{-0.019}$|	–	–	–
|$\mathbf {SLACSJ1430\!+\!4105}$|	–	|$\Psi _{\mathrm{l1}} = 3.37^{+0.570}_{-0.691}$|	|$R_{\mathrm{l1}} = 0.587^{+0.042}_{-0.053}$|	|$n_{\mathrm{l1}} = 2.86^{+0.14}_{-0.15}$|	|$I_{\mathrm{l1}} = 0.166^{+0.053}_{-0.047}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	–	|$\Psi _{\mathrm{l2}} = 7.71^{+0.386}_{-0.355}$|	|$R_{\mathrm{l2}} = 2.80^{+0.13}_{-0.17}$|	–	|$I_{\mathrm{l2}} = 0.016^{+0.006}_{-0.006}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	NFWEll	|$\kappa _{\mathrm{d}} = 0.077^{+0.025}_{-0.026}$|	|$q_{\mathrm{d}} = 0.960^{+0.021}_{-0.024}$|	–	–
|$\mathbf {SLACSJ1430\!+\!4105}$|	–	|$\Psi _{\mathrm{l1}} = 3.30^{+0.612}_{-0.710}$|	|$R_{\mathrm{l1}} = 0.632^{+0.051}_{-0.054}$|	|$n_{\mathrm{l1}} = 2.93^{+0.17}_{-0.15}$|	|$I_{\mathrm{l1}} = 0.153^{+0.038}_{-0.042}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	–	|$\Psi _{\mathrm{l2}} = 6.43^{+0.431}_{-0.452}$|	|$R_{\mathrm{l2}} = 3.04^{+0.19}_{-0.22}$|	–	|$I_{\mathrm{l2}} = 0.014^{+0.005}_{-0.006}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	gNFWSph	|$\kappa _{\mathrm{d}} = 0.099^{+0.036}_{-0.042}$|	|$\gamma _{\mathrm{d}} = 0.774^{+0.431}_{-0.545}$|	–	–
|$\mathbf {SLACSJ1430\!+\!4105}$|	–	|$\Psi _{\mathrm{l1}} = 3.31^{+0.742}_{-0.821}$|	|$R_{\mathrm{l1}} = 0.59^{+0.061}_{-0.059}$|	|$n_{\mathrm{l1}} = 2.84^{+0.17}_{-0.18}$|	|$I_{\mathrm{l1}} = 0.170^{+0.042}_{-0.041}$|
|$\mathbf {SLACSJ1430\!+\!4105}$|	–	|$\Psi _{\mathrm{l2}} = 6.78^{+0.568}_{-0.611}$|	|$R_{\mathrm{l2}} = 2.93^{+0.23}_{-0.24}$|	–	|$I_{\mathrm{l2}} = 0.015^{+0.005}_{-0.007}$|