# Strong gravitational lensing's 'external shear' is not shear

Amy Etherington, James W. Nightingale, Richard Massey, Sut-Ieng Tam, Xiaoyue Cao, Anna Niemiec, Qiuhan He, Andrew Robertson, Ran Li, Aristeidis Amvrosiadis, et al.

► To cite this version:

Amy Etherington, James W. Nightingale, Richard Massey, Sut-Ieng Tam, Xiaoyue Cao, et al. Strong gravitational lensing's 'external shear' is not shear. Monthly Notices of the Royal Astronomical Society, 2024, 531 (3), pp.3684-3697. 10.1093/mnras/stae1375. hal-03982577

HAL Id: hal-03982577

https://hal.science/hal-03982577v1

Submitted on 26 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Strong gravitational lensing’s ‘external shear’ is not shear

Amy Etherington<sup>1,2</sup> James W. Nightingale<sup>1,2★</sup> Richard Massey<sup>1,2</sup> Sut-Ieng Tam<sup>1</sup>,<sup>3</sup>

XiaoYue Cao<sup>1,4,5</sup> Anna Niemiec<sup>1</sup> Qiuhan He<sup>1,2</sup> Andrew Robertson<sup>1,6</sup> Ran Li<sup>4,5</sup>

Aristeidis Amvrosiadis<sup>1,2</sup> Shaun Cole<sup>1,2</sup> Jose M. Diego<sup>1,7</sup> Carlos S. Frenk<sup>1,2</sup> Brenda L. Frye<sup>1,8</sup>

David Harvey<sup>1,9</sup> Mathilde Jauzac<sup>1,2,10,11</sup> Anton M. Koekemoer<sup>1,10</sup> David J. Lagattuta<sup>1</sup>

Samuel Lange<sup>1,2</sup> Marceau Limousin<sup>1,11</sup> Guillaume Mahler<sup>1,1</sup> Ellen Sirks<sup>1,12</sup> and Charles

L. Steinhardt<sup>13</sup>

Accepted 2024 May 29. Received 2024 May 2; in original form 2023 January 26

## ABSTRACT

The distribution of mass in galaxy-scale strong gravitational lenses is often modelled as an elliptical power-law plus ‘external shear’, which notionally accounts for neighbouring galaxies and cosmic shear along our line of sight. A small amount of external shear could come from these sources, but we show that the vast majority does not. Except in a handful of rare systems, the best-fitting values do not correlate with independent measurements of line-of-sight shear: from weak lensing in 45 *Hubble Space Telescope* images, or in 50 mock images of lenses with complex distributions of mass. Instead, the best-fit external shear is aligned with the major or minor axis of 88 per cent of lens galaxies; and the amplitude of the external shear increases if that galaxy is discy. We conclude that ‘external shear’ attached to a power-law model is not physically meaningful, but a fudge to compensate for lack of model complexity. Since it biases other model parameters that *are* interpreted as physically meaningful in several science analyses (e.g. measuring galaxy evolution, dark matter physics or cosmological parameters), we recommend that future studies of galaxy-scale strong lensing should employ more flexible mass models.

**Key words:** gravitational lensing: strong – galaxies: structure.

## 1 INTRODUCTION

Gravitational lensing is the deflection of light rays by nearby concentrations of matter and their associated gravitational fields. If the light ray should pass straight through an object as massive as a galaxy, it can be deflected along multiple routes around the galaxy, and appear distorted into arcs or an ‘Einstein ring’. Such galaxy-scale strong lensing has been used to infer the distribution of mass in massive elliptical galaxies (Gavazzi et al. 2007; Koopmans et al. 2009; Auger et al. 2010; Bolton et al. 2012; Sonnenfeld et al. 2013; Etherington et al. 2022), to infer their dark matter content, stellar mass-to-light ratios, and inner structure (Massey, Kitching & Richard

2010; Sonnenfeld et al. 2012; Shu et al. 2015, 2016a; Oldham & Auger 2018; Nightingale et al. 2019). If the background source is variable, measurements of time delays between multiple images can be used to measure cosmological parameters (Kilbinger 2015) or the Hubble constant (Suyu et al. 2017; Wong et al. 2019; Harvey 2020; Gomer & Williams 2021). If the lens galaxy contains small substructures, which would be a smoking gun of the ‘clumpy’ cold dark matter (CDM) model, they would also perturb the multiple images (Natarajan & Springel 2004; Vegetti et al. 2010, 2014; Hezaveh et al. 2016; Li et al. 2016, 2017; Despali et al. 2019; Ritondale et al. 2019; Amorisco et al. 2022; He et al. 2022a, b; Nightingale et al. 2024).

When modelling the distribution of mass to fit strong lensing data, two additional free parameters are frequently included. The (amplitude and angle of) ‘external shear’ is intended to represent the cumulative deflection of light by all other gravitational potentials

© 2024 The Author(s).

<footnote n="1">¹Department of Physics, Centre for Extragalactic Astronomy, Durham University, South Rd, Durham, DH1 3LE, UK</footnote>

<footnote n="2">²Department of Physics, Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="3">³Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan</footnote>

<footnote n="4">⁴School of Astronomy and Space Science, University of Chinese Academy of Sciences, Beijing 100049, China</footnote>

<footnote n="5">⁵National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100012, China</footnote>

<footnote n="6">⁶Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="7">⁷Instituto de Física de Cantabria (CSIC-UC), Avda. Los Castros s/n, E-39005 Santander, Spain</footnote>

<footnote n="8">⁸Department of Astronomy/Steward Observatory, University of Arizona, 933 N Cherry Ave., Tucson, AZ 85721, USA</footnote>

<footnote n="9">⁹Laboratoire d'Astrophysique, EPFL, Observatoire de Sauverny, CH-1290 Versoix, Switzerland</footnote>

<footnote n="10">¹⁰Aix Marseille University, CNRS, CNES, LAM, F-13388 Marseille, France</footnote>

<footnote n="11">¹¹Space Telescope Science Institute, 3700 San Martin Dr, Baltimore, MD 21218, USA</footnote>

<footnote n="12">¹²Sydney Consortium for Particle Physics and Cosmology, School of Physics, The University of Sydney, NSW 2006, Australia</footnote>

<footnote n="13">¹³Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, København DK-Ø2100, Denmark</footnote>

<footnote>*E-mail: james.w.nightingale@durham.ac.uk</footnote>
along the line of sight. Indeed, the best-fitting value of external shear matched a model of the line of sight for three of six galaxy lenses studied by Wong et al. (2011). Subsequent papers have even proposed using external shear as relatively high signal-to-noise measurements of the ‘cosmic’ shear along individual lines of sight (Birrer et al. 2017; Desprez et al. 2018; Kuhn et al. 2020; Fleury, Larena & Uzan 2021; Hogg et al. 2022). However, the best-fitting values of external shear did not match the lines of sight to the three other galaxies studied by Wong et al. (2011) – and, in general, best-fitting values are much larger than expected for both galaxy lenses (Keeton, Kochanek & Seljak 1997; Witt & Mao 1997; Hilbert et al. 2007; Suyu et al. 2010; Barrera et al. 2021) and cluster lenses (Robertson et al. 2020; Limousin, Beauchesne & Jullo 2022). Using mock observations, Cao et al. (2022, hereafter C22) demonstrated that the best-fit shear can be incorrect if the model of the mass distribution is missing complexity.

In this paper, we compare the external shear measured from Hubble Space Telescope (HST) imaging – of strong lensing galaxies from the Sloan Lens ACS (SLACS) Survey, (Bolton et al. 2008a), BOSS Emission Line Lens GALaxy-Ly α EmitteR sYstems (GALLERY) survey (Shu et al. 2016b), and four lenses in clusters – against independent measurements of the shear along the same line-of-sight, observed as weak lensing of adjacent galaxies. To gain understanding, we also analyse C22’s mock HST imaging, generated without external shear, but fitted with shear as a free parameter.

Comparing independent measurements of shear will test the hypothesis that strong lensing external shear is a real, physical quantity. Strong and weak lensing measurements average over different spatial scales and are obtained at different redshifts, so they might not be identical; but they should be strongly correlated. Analysing three lenses in the Cosmic Evolution Survey (COSMOS) field (Faure et al. 2008) with sophisticated statistics, Kuhn et al. (2020) measured a smaller covariance between strong and weak lensing shears than the difference between individual systems, indicating that more data were required. In this work, with a much larger sample of galaxies, we simply aim to detect a correlation between the two probes.

We define relevant concepts from lensing theory in Section 2. We then describe our observed and mock data in Section 3 and our analysis methods in Section 4. We present our results in Section 5 and investigate possible causes in Section 6. We interpret these in a wider context and conclude in Section 7. Throughout this paper, we assume a Planck 2015 cosmological model (Ade et al. 2016).

## 2 GRAVITATIONAL LENSING THEORY

Gravitational lensing describes the deflection of light rays from distant sources, by matter along its path to an observer, through angle $\alpha$. This maps 2D coordinates of light on the distant source plane $\beta$ to coordinates where they are observed on the foreground image plane $\theta = (\theta_1, \theta_2)$, via the lens equation

$$ \beta = \theta - \alpha(\theta). \qquad (1) $$

If the gravitational lens is much thinner than its angular diameter distance from the observer $D_1$, its distribution of mass can be treated as a 2D surface density projected along the line-of-sight $\Sigma(\theta) = \int \rho(D_1\theta_1, D_1\theta_2, z) dz$. Following the notation of Narayan & Bartelmann (1996), the deflection angle $\alpha = \nabla\psi$ is the vector gradient of a 2D lensing potential

$$ \psi(\theta) = 2 \frac{D_{ls}}{D_1 D_s} \frac{2G}{c^2} \int \Sigma(\theta') \ln|\theta - \theta'| d^2\theta', \qquad (2) $$

$$ = 2 \frac{D_{ls}}{D_1 D_s} \int \Phi(D_1 \theta, z) dz, \qquad (3) $$

where $\Phi$ is the 3D Newtonian potential, and the pre-factor (which involves angular diameter distances to the lens, to the source, and from the lens to the source) reflects the geometrical efficiency of a lens: peaking if it is half-way between the source and the observer.

The Jacobian of the lens equation (1) is thus

$$ A_{ij} \equiv \frac{\partial \beta_i}{\partial \theta_j} = \delta_{ij} - \frac{\partial \alpha_i}{\partial \theta_j} = \delta_{ij} - \frac{\partial^2 \psi}{\partial \theta_i \partial \theta_j}. \qquad (4) $$

### 2.1 Weak lensing

If the source is much smaller than the scale of local variations in the gravitational field, the Jacobian can be approximated as constant

$$ A_{ij} \approx \begin{pmatrix} 1 - \kappa - \gamma_1 & -\gamma_2 \\ -\gamma_2 & 1 - \kappa + \gamma_1 \end{pmatrix} \qquad (5) $$

$$ = (1 - \kappa) \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} - \gamma \begin{pmatrix} \cos 2\phi & \sin 2\phi \\ \sin 2\phi & \cos 2\phi \end{pmatrix}, \qquad (6) $$

where the convergence

$$ \kappa = \frac{1}{2} \nabla^2 \psi, \qquad (7) $$

and where the two components of shear

$$ (\gamma_1, \gamma_2) = \left( \frac{1}{2} \left( \frac{\partial^2 \psi}{\partial \theta_1^2} - \frac{\partial^2 \psi}{\partial \theta_2^2} \right), \frac{\partial^2 \psi}{\partial \theta_1 \partial \theta_2} \right) \qquad (8) $$

can also be expressed in terms of shear amplitude $\gamma^2 = \gamma_1^2 + \gamma_2^2$ and angle $\phi = \frac{1}{2} \arctan(\gamma_2/\gamma_1)$. The convergence magnifies a source, and the shear changes its shape. Strictly, these quantities are only observable in combination as ‘reduced shear’ $g_i = \gamma_i/(1 - \kappa)$. However, in the weak lensing regime, $\kappa \ll 1$, so $g_i \approx \gamma_i$.

### 2.2 Strong lensing

If light from one side of a source is deflected differently to light from the other side, it can appear distorted in the image plane as an arc; it is also possible to see multiple images of a single source, if more than one solution exists with different $\alpha$ and $\theta$. To reconstruct $\alpha(\theta)$ it is usual to note that

$$ \kappa(\theta) = \frac{\Sigma(\theta)}{\Sigma_{\text{crit}}} \quad \text{with constant } \Sigma_{\text{crit}} = \frac{c^2}{4\pi G} \frac{D_s}{D_1 D_{ls}}, \qquad (9) $$

which is equal to the mean surface mass density within the Einstein radius, $R_{\text{Ein}}$. For axisymmetric lenses the Einstein radius is uniquely defined by the radius of the circular tangential critical curve that is produced where the magnification diverges in the lens plane. This occurs where the tangential eigenvalue of the Jacobian (equation 4) $\lambda_t = 1 - \kappa - \gamma$ is equal to zero. For asymmetric lenses, the definition of Einstein radius must be generalized; we choose to use the effective Einstein radius

$$ R_{\text{Ein,eff}} = \sqrt{\frac{A}{\pi}}, \qquad (10) $$

where A is the area enclosed by the tangential critical curve.

When considering (typically early-type) galaxy-scale lenses, it is common practise to parametrize the surface mass distribution as an elliptical power-law (Suyu 2012)

$$ \Sigma(\theta_1, \theta_2; b, q, \gamma) = \frac{3 - \gamma}{1 + q} \left( \frac{b}{\sqrt{\theta_1^2 + \theta_2^2 / q^2}} \right)^{\gamma-1}, \qquad (11) $$

where $b \ge 0$ is the angular scale length (referred to in some papers as the Einstein radius, but distinct from the more robust effective
Einstein radius in equation 10), $0 < q \le 1$ is the projected minor to major axis ratio of the elliptical isodensity contours, and (confusingly denoted) $\gamma$ is the logarithmic slope of the mass distribution in 3D (for an 'isothermal' distribution, $\gamma = 2$). If we also allow the mass to be translated to central coordinates ($\theta_1^s, \theta_2^s$) and reoriented to position angle $\phi^{\text{PL}}$, which we measure counterclockwise from the positive $\theta_1$-axis, the model has six free parameters.

The primary lens may not be the only source of shear. Any 'external' component due to other galaxies or clusters near the lens or along the ray path, and constant on scales larger than $b$ (rather than the size of the source) is modelled as two more free parameters

$$ (\gamma_1^{\text{ext}}, \gamma_2^{\text{ext}}) = \gamma^{\text{ext}} (\cos(2\phi^{\text{ext}}), \sin(2\phi^{\text{ext}})), \quad (12) $$

where $\gamma^{\text{ext}}$ is the amplitude and $\phi^{\text{ext}}$ is the angle of the shear measured counterclockwise from the $\theta_1$-axis. This is applied as an additional component of $\alpha(\theta)$. It does not change $\kappa(\theta)$.

# 3 DATA

## 3.1 Mock lens galaxies

We analyse a set of 50 mock lens images, representative of data from the SLACS survey. They were generated by C22 for an investigation into the systematic errors induced by the elliptical power-law model. We summarize the simulation procedure below; a detailed description can be found in section 2.4 of that paper.

The surface mass density of the lens galaxy comprises two components: a dark matter halo, parametrized by the spherical generalized Navarro, Frenk & White (gNFW) profile (Zhao 1996; Navarro, Frenk & White 1997; Cappellari et al. 2013), plus visible stellar matter, parametrized by a multiple Gaussian expansion (MGE; Cappellari 2002). The model parameters of the gNFW and MGE profiles of each lens galaxy are set to the best-fitting parameters from fits of these distributions to Sloan Digital Sky Survey-Mapping Nearby Galaxies at APO (SDSS-MANGA) stellar dynamics data, derived by Li et al. (2019) using the Jeans anisotropic model (JAM) method. The position angle of each Gaussian component in the MGE is fixed, however their axis ratios are free to vary, allowing for elliptical gradients in the mass distribution.

The light distribution of the source galaxy is modelled by a single Sérsic profile (Graham & Driver 2005) with effective radius $R_{\text{eff}} = 0.15$ arcsec, Sérsic index $n = 1$, and axis ratio $q = 0.7$. The position in the source plane ($x_s, y_s$) is drawn from a Gaussian distribution with mean 0 arcsec and standard deviation 0.1 arcsec, and the position angle is uniformly selected between 0° and 180°. The light from the source galaxy is ray-traced from the source plane at $z = 0.6$ to the image plane through the lens equation (equation 1), to simulate its lensed appearance. Further, to mimic observational effects, the image is convolved with a Gaussian PSF with 0.05 arcsec standard deviation, and sampled by 0.05 arcsec square pixels. A flat background sky of 0.1 electrons per second is assumed, and an exposure time of 840 s is used to add Poisson noise from the source and background sky. The signal-to-noise ratio of the brightest pixel in the synthetic images is set to ~50, by adjusting the intensity of the Sérsic source accordingly. No external shear was simulated in the mock data.

## 3.2 Observed lens galaxies

We analyse three sets of galaxy–galaxy strong lenses. These include 42 lenses from the SLACS survey (Bolton et al. 2008a) that were

<table><caption>**Table 1.** Table of parameters for the four galaxy–galaxy lenses in the outskirts of clusters, the lens name refers to the cluster in which the lens resides.</caption><thead><tr><th>Lens name</th><th>RA</th><th>Dec.</th><th>z<sub>1</sub></th><th>z<sub>s</sub></th></tr></thead><tbody><tr><td>MACS1149-GGL18</td><td>177.410247</td><td>22.352017</td><td>0.544</td><td>—</td></tr><tr><td>Abell370-GGL19</td><td>39.963013</td><td>-1.534783</td><td>0.375</td><td>2.371</td></tr><tr><td>MACS1149-GGL20</td><td>177.402816</td><td>22.436607</td><td>0.544</td><td>1.806</td></tr><tr><td>RX J2129-GGL1 (snail)</td><td>322.428780</td><td>0.108071</td><td>0.235</td><td>1.610</td></tr></tbody></table>

fitted without significant residuals by Etherington et al. (2022, hereafter E22)'s automated pipeline. Most are isolated field galaxies. They were found by searching for high-redshift emission lines in the spectra of low-redshift galaxies obtained through a 3 arcsec fibre. They were then imaged by the HST Advanced Camera for Surveys (ACS) in the F814W band, and processed into stacked images with 0.05 arcsec pixels. We also reprocessed these to measure weak lensing, following the procedure described by Tam et al. (2020), which supersamples the pixels to 0.03 arcsec. We exclude lenses J1143–0144 and J1420+6019, for which only one exposure was obtained.

We analyse 15 lenses from the GALLERY survey (Shu et al. 2016b) that were modelled by E23. These are also field galaxies, found by searching for compact Lyman $\alpha$-emitting source galaxies in spectra with a 2 arcsec fibre. They were imaged with the HST Wide-Field Camera 3 (WFC3) in the F606W band and processed, following Shu et al. (2016a), into 0.04 arcsec pixels. We do not attempt to measure weak lensing shear in these data.

We analyse four galaxy–galaxy strong lenses in the outskirts of galaxy clusters, where we expect a 5–15 per cent true external shear. Before beginning any analysis, we searched archival HST F814W imaging, and selected lenses with multiple imaging of sources that are extended a similar amount as the field lenses. Positions and redshifts of the selected lenses are given in Table 1. For MACS1149-GGL18 no source redshift has been recorded, where necessary we assume a source redshift of $z_s = 1.5$ and test that the results do not change significantly when we change this assumption over a range of source redshifts from 0.5 to 2. Two of our selected lenses had been previously modelled by Desprez et al. (2018), although constrained using only the positions of multiple images, rather than all the pixels. We analysed the HST data similarly to the SLACS lenses, except for the 'cosmic snail'. For that lens alone, we do not measure weak lensing, but use the independent estimate of shear from Desprez et al. (2018)'s model IV of the galaxy cluster, constrained by cluster-scale strong lensing, and shown by Desprez et al. (2018) to be consistent with ground-based measurements of weak lensing.

# 4 METHODS

## 4.1 Weak lensing analysis

We identified galaxies on lines of sight adjacent to strong lenses using SExTRACTOR (Bertin & Arnouts 1996), and measured their shapes using the PYRRG (Harvey et al. 2019) implementation of the Rhodes, Refregier & Groth (2000) shear measurement method. This estimates the mean reduced shear in a patch of sky, by averaging galaxies' apparent shapes

$$ \epsilon_i = \epsilon_i^{\text{int}} + G\gamma_i, \quad (13) $$

which have been transformed by weak lensing (Section 2.1) from an unknown intrinsic shape $\epsilon^{\text{int}}$. The 'shear responsivity' $G$ varies as a
function of galaxy flux, and its overall scaling has been calibrated on simulated images with known shear (Leauthaud et al. 2007). Under the assumption that galaxies’ intrinsic shapes are randomly oriented, i.e. $\langle \epsilon_i^{\text{int}} \rangle \approx 0$,

$$ \langle \gamma_i \rangle = \left\langle \frac{\epsilon_i - \epsilon_i^{\text{int}}}{G} \right\rangle \approx \left\langle \frac{\epsilon_i}{G} \right\rangle - \left\langle \frac{\epsilon_i^{\text{int}}}{G} \right\rangle \approx \frac{\langle \epsilon_i \rangle}{\langle G \rangle}. \quad (14) $$

Following Massey et al. (2007), we assume that the median redshift of the lensed galaxies is $z \sim 1.26$. Thereafter, following Smail, Ellis & Fitchett (1994), we treat them all as being at this effective redshift. None of our results change significantly if we adjust this value.

We average weak lensing shear measurements from the ~140 galaxies within 60 arcsec of the strong lens galaxy (no weights are applied to the galaxies that are averaged). The precision of this measurement is limited by the randomness in the distribution of the intrinsic shapes

$$ \sigma_{\text{int}}^2 = \left\langle \left( \frac{\epsilon_i^{\text{int}}}{G} \right)^2 \right\rangle. \quad (15) $$

We measure $\sigma_{\text{int}} \sim 0.3$, consistent with Leauthaud et al. (2007), and hence uncertainty $\sigma_{\text{int}}/\sqrt{140} \approx 0.02$ on each component of mean shear. This is similar to uncertainty on our strong lensing measurements. None of our results change significantly if we use a 45 or 90 arcsec aperture instead.

Although the line-of-sight directly through each galaxy-scale lens is not in the weak lensing regime, we assume that $\langle g_i \rangle \approx \langle \gamma_i \rangle$ still holds since the vast majority of adjacent lines of sights will be only weakly lensed. Nor do we attempt to model and subtract the weak shear due to the galaxy-scale lens itself. Doing so would mix the weak lensing and strong lensing analyses; and it is unnecessary at our achieved level of precision because the near-circular symmetry of most lenses means that the lens contributes negligibly to $\langle \gamma_i \rangle$ inside a 60 arcsec circular aperture.

## 4.2 Strong lensing analysis

We analysed all data using the automated strong lens modelling software PYAUTOLENS<fnref n="1" /> (Nightingale, Dye & Massey 2018; Nightingale et al. 2021b). This fits parameters of the lens model using all of the pixels in an image (not just e.g. locations of the centre of light, as in previous works).

The pipelines used to fit the mock and observed data are described fully in C22 and E23, respectively. Briefly, we model the distribution of mass in both mock and real data using an elliptical power-law (equation 11) plus external shear (equation 12). We then repeat the fit, fixing external shear $\gamma^{\text{ext}} = 0$. We model the distribution of light in real lens galaxies using a double Sérsic profile with a centre that is free to vary independently to that of the mass distribution, and for the source galaxy using an adaptive Voronoi mesh of pixels. For the mock data, we use C22’s fit in which the lens light is perfectly subtracted and the source light is modelled as an elliptical Sérsic. C22 also performed fits using a Voronoi mesh for the source light. However, since the mock data were created assuming a Sérsic source, the model we chose can perfectly describe the source, so any systematics we observe will be solely due to mismatch between the model and truth of the mass distribution, which is the point of interest in this study.

In Appendix A, we perform two tests of our shear measurements: (i) we refit every strong lens model including one or two nearby

bright galaxies explicitly as singular isothermal sphere and; (ii) use a different source analysis which assumes a Delaunay mesh. We recover consistent shear magnitudes and position angles across the lens sample. Results presented in the main paper therefore cannot be explained as due to lens models missing line-of-sight galaxies or systematics in the source analysis.

## 4.3 Multipole perturbations of an ellipse

In Section 6.1, we shall investigate whether the strong lensing external shears depend on deviations from an elliptically symmetric distribution of mass. Specifically, we shall quantify the multipole deviations of two types of contour: the iso-convergence contour at $\kappa = 1$ of the gNFW + MGE distributions used to create mock data, and the critical curves of both the mock and the observed galaxies. Note, that this definition follows the use of isophote shapes in studies of galaxy morphology (Ryden 1992; Chaware et al. 2014) and sometimes in studies of lens mass morphology (O’Riordan & Vegetti 2024), but is different from the perturbations produced by additional mass with $m = 4$ moments in Chu et al. (2013).

To implement this measurement, we store contours as a 2D array of points in polar coordinates [$\phi_{\text{contour}}$, $R_{\text{contour}}$], then calculate perpendicular deviations of each point from the true ellipse

$$ R_{\text{el}}(\phi) = \frac{a\sqrt{1 - \epsilon^2}}{\sqrt{1 - \epsilon^2 \cos^2(\phi - \phi_{\text{el}})}}, \quad (16) $$

where $a$ is the major axis, $\phi_{\text{el}}$ is the major axis orientation, and $\epsilon$ is the eccentricity (defined as $\epsilon^2 = 1 - b^2/a^2$ where $b$ is the minor axis). The deviations are then parametrized using multipoles

$$ \delta R_m(\phi; a_m, b_m) = \sum a_m \cos(m(\phi - \phi_{\text{el}})) + b_m \sin(m(\phi - \phi_{\text{el}})), \quad (17) $$

where $m$ is the order of the multipole, and $a_m$ and $b_m$ are the magnitude of the deviations with symmetry along or at 45° to the major and minor axes, respectively. We then perform a non-linear search to fit the model

$$ R(\phi; a_4, b_4) = R_{\text{el}}(\phi) + \delta R_4(\phi; a_4, b_4), \quad (18) $$

to the radial values of the contour. We assume uniform priors on the free parameters in the fit over a reasonable range and fit for them using the nested sampling algorithm DYNESTY. We assume the residual errors can be described by a Gaussian distribution and maximize the likelihood

$$ \mathcal{L}(R | R_i, \sigma) = \prod_i \left[ \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(R(\phi_i) - R_i)^2}{2\sigma^2}\right) \right], \quad (19) $$

where $R_i$ are the radial values of the contour and $R(\phi_i; a_4, b_4)$ are the model-predicted values from equation (18) at each angular coordinate in the contour $\phi_i$. Curves with best-fitting values of $a_4 > 0$ are ‘discy’; those with $a_4 < 0$ are ‘boxy’ (see Fig. 1).

# 5 RESULTS

## 5.1 SL shears do not correlate with WL shears

Strong lensing measurements of shear $\gamma^{\text{SL}}$ (obtained as the best-fit external shear $\gamma^{\text{ext}}$) typically have amplitudes up to an order of magnitude larger than weak lensing measurements $\gamma^{\text{WL}}$ along the same line of sight. The mock lenses have mean best-fit $|\langle \gamma^{\text{SL}} \rangle| = 0.019 \pm 0.002$, despite the true values all being $\gamma^{\text{ext}} = 0$. Our measurement using a PL + ext mass model is consistent with C22’s

<footnote n="1">¹The PYAUTOLENS software is open source and available from https://github.com/Jammy2211/PyAutoLens</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Examples of boxy ($a_4 = -0.1$) and discy ($a_4 = +0.1$) perturbations to an ellipse. The perturbations shown are ~10 times larger than those typically observed. In both cases the perturbation at 45°, $b_4 = 0$ (see equation 17).</figcaption></figure>

value $|\langle\gamma^{\text{SL}}\rangle| = 0.015$ using pixel-based source reconstructions. The real lenses have mean best-fit shear $|\langle\gamma^{\text{SL}}\rangle| = 0.098 \pm 0.011$, which is much larger than both the measured weak lensing shear $|\langle\gamma^{\text{WL}}\rangle| = 0.028 \pm 0.002$.

Strong lensing measurements of shear do not correlate with weak lensing measurements (Fig. 2). To make this comparison (for the real lenses only), we first define rotated coordinate systems such that $\gamma_1^{\text{SL}} = \gamma_1^{\text{SL}} \ge 0$ and $\gamma_2^{\text{SL}} = 0$. Thus, we need plot only three of the four components of shear. Second, we compensate for the different redshifts of the strongly lensed and weakly lensed sources by rescaling values of $\gamma_1^{\text{SL}}$ by $(D_s/D_{ls})z_i^* = z_i(D_{ls}/D_s)z_i^* = 1.26$, i.e. the effective value at the redshift of the weakly lensed galaxies (see equation 2). This scaling is exact only if the external shear is both real and dominated by neighbouring structures at the same redshift as the lens (Wong et al. 2011 found that 5/8 of the shear is from neighbouring structures). In any case, the rescaling is by a factor with mean of only 1.26 and rms 1.06, and our conclusions do not change if the rescaling is omitted or normalized to a different redshift. If strong and weak lensing measure the same quantity, we then expect $\gamma_1^{\text{WL}}$ to correlate with $\gamma_1^{\text{SL}}$, and $\gamma_2^{\text{WL}}$ to scatter around zero. We find that $\langle\gamma_2^{\text{WL}}\rangle = -0.004 \pm 0.003$ is on average below zero, and its scatter (0.02) is consistent with uncertainties calculated from the distribution of weak lensing shears. The best-fit slope $\gamma_1^{\text{WL}} = (-0.06 \pm 0.04)\gamma_1^{\text{SL}}$ actually infers a negative correlation, however this does not take into account the uncertainty on the strong lensing shears and so the uncertainty is likely underestimated. The Pearson correlation coefficient $-0.19 \pm 0.22$ implies that, if there is a correlation, there is also a large amount of scatter.

There are eight lenses for which $\gamma^{\text{ext}} \approx \gamma^{\text{WL}}$, including two of the four lenses which reside in the outskirts of clusters (see Section 5.4 for further discussion). However, there does not appear to be anything unique about these lenses that would make the shear possible to measure in these cases.

## 5.2 SL shears are (suspiciously) aligned with the mass

For both mock and real data, the best-fit ‘external’ shear is usually aligned with the major axis of the lens mass ($\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{ext}}^{\text{PL+ext}} \le 30^\circ$) or with its minor axis ($\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{ext}}^{\text{PL+ext}} \ge 60^\circ$): see Fig. 3. If external shears were truly measuring external perturbations, their orientations would be random (modulo intrinsic alignments between the shape of a galaxy and its surrounding tidal field, but these are much smaller than our achieved measurement precision; Zhang et al. 2023). The preference for aligning with the mass distribution again

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Values of the shear along the lines of sight to 39 galaxy–galaxy lenses, independently measured using strong lensing ‘external’ shear $\gamma^{\text{SL}}$ and weak lensing $\gamma^{\text{WL}}$. Shears are oriented such that $\gamma_2^{\text{SL}} = 0$, and rescaled to be at the same effective redshift. If strong and weak lensing shears were identical, all points would lie on the dashed lines. We instead find that external shears inferred from strong lensing are consistently larger than those measured by weak lensing, and not aligned.</figcaption></figure>

suggests an ‘internal’ shear that compensates for the inability of a power-law model to represent the more complex true distribution of mass. Furthermore, the highest values of $\gamma_{\text{ext}}$ are also usually found in the most elliptical lenses.

In mock data, 84 per cent of external shears are aligned with the mass distribution: their mean offset is 3° with an rms scatter of 5°. 14 per cent of external shears are anti-aligned with the mass distribution, with a mean of 85° and scatter of 6°. Only one lens has a best-fit external shear that is neither aligned nor anti-aligned, but this also has the lowest shear amplitude ($\gamma_{\text{ext}} = 0.0003$), so the angle $\phi_{\text{ext}}$ is ill-defined. The Pearson correlation coefficient between the best-fit axis ratios and external shears is $-0.63$.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Relative orientation of the strong lensing 'external' shear and the major axis of the lens mass, in mock (top panel) and real HST (bottom panel) data. In both cases, most of the shears are suspiciously aligned ($\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{ext}}^{\text{PL+ext}} \le 30^{\circ}$) or anti-aligned ($\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{ext}}^{\text{PL+ext}} \ge 60^{\circ}$) with the lens mass distribution. If the $\gamma_{\text{ext}}$ parameter were measuring true external perturbations, the orientations would be random. Highly elliptical lenses often lead to high values of $\gamma_{\text{ext}}$.</figcaption></figure>

Real *HST* data produce a similar pattern, but with more scatter. Best-fit ‘external’ shears are aligned with the mass distribution in 68 per cent of the lenses, and anti-aligned in 20 per cent. All the remaining 12 per cent have $\gamma_{\text{ext}} < 0.04$, so the angles $\phi_{\text{ext}}^{\text{PL+ext}}$ are noisy. The Pearson correlation coefficient between the best-fit axis ratios and external shears is −0.60. Despite the inferred external shears being an order of magnitude larger in the observations than in the mock data, the mass distributions are similarly elliptical: $\langle q \rangle = 0.77 \pm 0.02$ in the mocks, and $0.69 \pm 0.02$ for *HST* data.

## 5.3 SL shear amplitudes are too large

Along random lines of sight through the universe, a shear amplitude of ~1–3 per cent is expected (Keeton et al. 1997), even accounting for the different scales on which measurements need to be averaged to detect a signal (Valageas, Barber & Munshi 2004; Wong et al. 2011). SLACS and GALLERY strong lenses typically exist in field environments away from galaxy clusters, where shear amplitudes should never exceed ~5 per cent (Treu et al. 2009). However, most of the measured SL shears have amplitude exceeding this, and a large fraction exceed 10 per cent (Fig. 3). Therefore, the SL measurements alone indicate that their ‘external shear’ parameter cannot be tracing external shear, even without WL data or detailed inspection of the SL mass models.

## 5.4 Lenses in clusters

We now discuss the four lenses in clusters, which due to their cluster environment are expected to have external shear magnitudes around 5 per cent to 15 per cent (Treu et al. 2009), for these lenses, the argument presented in Section 5.3 above does not apply, noting their inferred shears are in the expected range (see Table 2).

### 5.4.1 RX J2129-GGL1 (snail)

Of the galaxy–galaxy lenses in clusters, the snail has the closest match between SL and WL measurements of shears. It also has the largest shear amplitude. Such measurements are even independent of analysis details, with Desprez et al. (2018)'s cluster-scale strong lensing analysis confirming a value of shear that agrees with that derived from weak lensing analysis of Cluster Lensing And Supernova survey with Hubble (CLASH) data (fig. 13 of that paper). This agreement provided early optimism in using galaxy–galaxy lensing within clusters as a high signal-to-noise measurement of weak lensing shear (e.g. Hogg et al. 2022); but it appears an outlier in our sample. Indeed, while the strong lensing external shear is anti-aligned with the distribution of mass in the galaxy, this is not surprising because the galaxy is coincidentally orientated with its major axis pointing towards the cluster centre (Table 2).

### 5.4.2 MACS1149-GGL20

The shears measured using the independent probes for MACS1149-GGL20 are also in agreement, although the shear magnitude is much lower than is measured for the snail. In fact, this is one of the few lenses that measures a lower best-fitting value of shear magnitude with strong lensing $\gamma^{\text{ext}} = 0.01$ than it does with weak lensing $\gamma^{\text{WL}} = 0.04$. Desprez et al. (2018) found that measurements of external shear from modelling the GGL alone were underestimated compared to the shears constrained using a full-scale model of the galaxy cluster. However, both of these measurements are larger than either of the shears measured in this work. When modelling the cluster as a double pseudo-isothermal ellipse (dPIE), they infer shear amplitude $\gamma^{\text{ext}} = 0.13_{-0.06}^{+0.08}$ at the location of the galaxy, and axis ratio $q = 0.81$. We instead measure significantly lower shear $\gamma^{\text{ext}} = 0.01_{-0.01}^{+0.02}$ and more elliptical mass distribution $q = 0.51_{-0.01}^{+0.02}$. It is therefore possible that the degeneracy between shear and axis ratio may explain this discrepancy. As with the snail, the galaxy's major axis coincidentally points towards the cluster centre, therefore our observed anti-alignment of external shear with the galaxy is to be expected.

### 5.4.3 Abell370-GGL19

We measure a similar shear magnitude with strong lensing for Abell370-GGL19 as we do with weak lensing, but the strong lensing external shear is suspiciously orientated towards a nearby neighbour galaxy and is aligned with the mass distribution (see Table 2). We therefore repeat the fit including free parameters for a singular isothermal sphere (SIS; $\gamma = 2$ and $q = 1$ in equation 11) fixed at the centre of the neighbour galaxy. The results including the mass of the neighbour galaxy do not change significantly (see the SIS neighbour column of Table 2 compared to the no neighbour column), although the power-law mass distribution does become less elliptical.
<table><caption>**Table 2.** Summary of strong and weak lensing parameters for the four galaxy–galaxy lenses that reside in clusters. All angles are in degrees anticlockwise from West.</caption><thead><tr><th rowspan="2">Lens name</th><th colspan="4">SIS neighbour</th><th colspan="2">No neighbour</th><th rowspan="2">γ<sup>WL</sup><sub>1</sub></th><th rowspan="2">γ<sup>WL</sup><sub>2</sub></th><th rowspan="2">γ<sup>WL</sup></th><th rowspan="2">φ<sup>WL</sup></th><th rowspan="2">φ<sup>BCG</sup></th><th rowspan="2">φ<sup>neighbour</sup></th></tr><tr><th>γ<sup>ext</sup></th><th>φ<sup>ext</sup></th><th>q<sup>PL</sup></th><th>φ<sup>PL</sup></th><th>γ<sup>ext</sup></th><th>φ<sup>ext</sup></th><th>q<sup>PL</sup></th></tr></thead><tbody><tr><td>MACS1149-GGL18</td><td>0.13<sup>+0.04</sup><sub>-0.02</sub></td><td>40<sup>+5</sup><sub>-5</sub></td><td>0.48<sup>+0.04</sup><sub>-0.02</sub></td><td>147<sup>+9</sup><sub>-8</sub></td><td>0.24<sup>+0.02</sup><sub>-0.02</sub></td><td>26<sup>+2</sup><sub>-2</sub></td><td>0.82<sup>+0.02</sup><sub>-0.02</sub></td><td>174<sup>+19</sup><sub>-22</sub></td><td>0.01<sup>+0.01</sup><sub>-0.01</sub></td><td>-0.04<sup>+0.01</sup><sub>-0.01</sub></td><td>0.04</td><td>-40</td><td>85</td><td>95</td></tr><tr><td>Abell370-GGL19</td><td>0.06<sup>+0.01</sup><sub>-0.01</sub></td><td>40<sup>+3</sup><sub>-6</sub></td><td>0.92<sup>+0.01</sup><sub>-0.01</sub></td><td>23<sup>+1</sup><sub>-1</sub></td><td>0.07<sup>+0.03</sup><sub>-0.02</sub></td><td>34<sup>+7</sup><sub>-14</sub></td><td>0.72<sup>+0.03</sup><sub>-0.02</sub></td><td>26<sup>+5</sup><sub>-6</sub></td><td>0.05<sup>+0.02</sup><sub>-0.02</sub></td><td>-0.01<sup>+0.02</sup><sub>-0.02</sub></td><td>0.05</td><td>-8</td><td>-100</td><td>35</td></tr><tr><td>MACS1149-GGL20</td><td>–</td><td>–</td><td>–</td><td>–</td><td>0.01<sup>+0.02</sup><sub>-0.01</sub></td><td>13<sup>+47</sup><sub>-49</sub></td><td>0.51<sup>+0.02</sup><sub>-0.01</sub></td><td>106<sup>+1</sup><sub>-2</sub></td><td>0.01<sup>+0.02</sup><sub>-0.02</sub></td><td>0.04<sup>+0.02</sup><sub>-0.02</sub></td><td>0.04</td><td>40</td><td>-83</td><td>–</td></tr><tr><td>RX J2129-GGL1 (snail)</td><td>–</td><td>–</td><td>–</td><td>–</td><td>0.11<sup>+0.01</sup><sub>-0.01</sub></td><td>35<sup>+2</sup><sub>-2</sub></td><td>0.93<sup>+0.01</sup><sub>-0.01</sub></td><td>-61<sup>+10</sup><sub>-12</sub></td><td>0.08<sup>+0.08</sup><sub>-0.08</sub></td><td>-0.02<sup>+0.02</sup><sub>-0.02</sub></td><td>0.08</td><td>28</td><td>-57</td><td>–</td></tr></tbody></table>

### 5.4.4 MACS1149-GGL18

There is also a neighbour galaxy in close proximity to MACS1149-GGL18. We, therefore repeat the fit for including an SIS as was done for Abell370-GGL19. The shear is significantly overestimated as compared with weak lensing when the neighbour is not included in the fit. Including the neighbour galaxy halves the inferred strong lensing external shear, but this is still inconsistent with the weak lensing measurement. The shear is anti-aligned with the power-law mass distribution which, given that the mass distribution is not aligned with the cluster galaxy in this case, suggests the external shear may be acting internally as discussed in the previous section.

# 6 ANALYSIS

## 6.1 External shear may compensate for boxiness/disciness

Our measurements in Section 5 suggest that $\gamma^{\text{ext}}$ mostly just compensates for the inability of a power-law model to capture the complex distributions of mass (gNFW + MGE for the mocks, and likely more complex for real galaxies). This is consistent with the conclusion of Keeton et al. (1997), who inferred that the $\langle\gamma^{\text{ext}}\rangle \sim 10-15$ per cent required to fit point-source quad lenses, must reflect an inability of the lens model to capture a complex distribution of mass: perhaps misalignment between light and dark matter. Witt & Mao (1997) reached a similar conclusion, and derived an analytical prediction of the shear required by an elliptical potential to fit quad lenses.

If the external shears result from a lack of complexity in the power-law model to describe the underlying distribution of mass, one can ask what type of complexity the data requires. One possible deviation from the symmetry of an elliptical power law is boxiness and disciness (see Section 4.3). We shall now investigate whether spurious external shear could arise to compensate for boxy/discy lens galaxies.

### 6.1.1 External shear creates boxy/discy critical curves

An isothermal elliptical mass distribution has an elliptical critical curve (oriented in the same direction as the mass distribution but at 90° to the elongation of light from the source galaxy; Kochanek, Schneider & Wambsganss 2004). However, changing the power-law slope, $\gamma \neq 2$, or adding an external shear, $\gamma^{\text{ext}} \neq 0$, perturbs the critical curves (left panel of Fig. 4). These perturbations include significant $a_4/a$ moments (right panels of Fig. 4), although they visually appear to be more than pure $m=4$ modes (c.f. Fig. 1).

### 6.1.2 Discy critical curves come from discy mass distributions

The distributions of mass in our mock lenses happen to be almost all discy. We could measure any isodensity contour, but the $\kappa=1$ contour will be near the most sensitive region for lens fitting. These iso-convergence contours have mean $\langle|a_4/a|\rangle = 0.01$ and $\langle|b_4/a|\rangle =$

0.0005. Only three lenses are boxy, but not usefully so, with very low values $\langle a_4/a \rangle = -0.0003$.

Critical curves of the best-fit PL + ext models to our mock data show $a_4/a$ moments highly correlated with those of the density contours (Fig. 5). Again, $\langle|b_4/a|\rangle = 0.0001$ is an order of magnitude lower. Studying the same mocks, C22 also noted that ‘external’ shear allowed the best-fit critical curves to better match the true critical curves. We find two systems with boxy critical curves $a_4/a < -0.01$. Subject to some scatter, however, we conclude that the disciness of isodensity contours and critical curves are highly correlated.

Notably, all mock lenses whose best-fit external shear is aligned with the mass distribution have very discy critical curves (red points in Fig. 6), and the three mock lenses with boxy critical curves have anti-aligned shear. Furthermore, $a_4$ typically increases with the external shear (Pearson correlation coefficient 0.45) and with the axis ratio of the lens mass (Pearson correlation coefficient $-0.73$). This may be tentative evidence that (some of) the dichotomy of aligned and anti-aligned shears may be caused by disciness or boxiness.

### 6.1.3 Does boxiness/disciness cause ‘external’ shear?

Scatter in real data is larger than in the mocks. However, for SLACS and GALLERY lenses, the best-fit critical curves have mean $\langle|a_4/a|\rangle = 0.016$, similar to the mocks. Most (79 per cent of) lenses with best-fit ‘external’ shear that is aligned with the mass distribution have discy critical curves; and most (70 per cent) with anti-aligned shear have boxy critical curves (Fig. 7). Moreover, lenses with the largest amplitude of external shear also have critical curves with the largest deviations from elliptical (Pearson correlation coefficient of 0.48 with $a_4$ and 0.65 with $b_4$).

This provides tentative evidence that ‘external’ shear in typical lensing analyses is really caused by the inability of parametric mass models to capture the complex distribution of mass in a lens. A substantial portion of that complexity may be disciness or boxiness of the mass distribution. This creates disciness or boxiness in the critical curves (Section 6.1.2), which leads to a spurious external shear (Section 6.1.1). We have not been able to quantify the relative contributions to external shear from true shear, disciness/boxiness, or other sources.

### 6.1.4 More things probably cause ‘external’ shear too

There are likely more sources of complexity in real mass distributions, which cause (or are compensated by) external shear. These confounding factors would explain the looser correlations and larger scatter than in the mocks. Just the observation that real lenses have external shears with amplitudes six times greater than mocks implies that their distribution of mass deviates more from a power law.

We speculate that the isodensity contours of a lens might be twisted (misaligned as a function of radius), like their isophotoes. Indeed, the critical curves of SLACS and GALLERY lenses have a handedness,
<figure><img src="image_6.png" /></figure>

**Figure 4.** What causes boxiness or disciness? A singular isothermal elliptical (SIE) mass distribution with a horizontal major axis has critical curves that are also elliptical with a horizontal major axis. The critical curves are perturbed if the slope of the density profile $\gamma \neq 2$ (top left panel) or the external shear $\gamma_{\text{ext}} \neq 0$ (bottom left panel). In particular, an aligned shear ($\gamma_{1}^{\text{ext}} > 0$) stretches the critical curves vertically (and the image horizontally); an anti-aligned shear ($\gamma_{1}^{\text{ext}} < 0$) does the opposite. Multipole measurements $a_4/a$ of the critical curve are shown as a function of slope (top right panel) and external shear (bottom right panel), where $a_4/a > 0$ is 'discy' and $a_4/a < 0$ is 'boxy'.

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**Figure 5.** In mock lenses, discy ($a_4/a > 0$) perturbations of the $\kappa = 1$ isodensity contour correlate with discy perturbations of the critical curves – which can also be measured for real lenses. To better visualize the correlation, values have been transformed by $\log[a_4/a + 0.01]$, with dashed lines indicating $a_4 = 0$. Unfortunately, the mock data do not include any lenses with significantly boxy ($a_4/a < 0$) distributions of mass. These points come from fits with Sersic sources, so the uncertainties are not comparable to those from analyses with pixellated sources.

with $\langle |b_4/a| \rangle = 0.012$ two orders of magnitude larger than the mocks. If the critical curves do tell us something about the distribution of mass, this may indicate twisted isodensity contours. Twisting is also

<figure><img src="image_9.png" /></figure>

**Figure 6.** For mock lenses that were simulated without external shear. Angle between the best-fitting values of external shear and the major axis of the lens mass distribution, $|\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{ext}}^{\text{PL+ext}}|$ in degrees, as a function of the amplitude of the best-fit external shear, $\gamma_{\text{ext}}$. Points are coloured by the magnitude of the inferred critical curves deviation from elliptical symmetry $a_4/a$, values of $a_4/a < 0$ correspond to boxy critical curves and $a_4/a > 0$ to discy ones. Values are orders of magnitude different, with the lowest values going down to 0.001.

suggested by the inconsistently measured position angle when real data are fitted with and without external shear (see Fig. 9); this is not present in the mocks. Van de Vyvere et al. (2022) also found that twists in the underlying mass distribution are typically absorbed by changes in orientation of the mass distribution and shear in a
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Same as for Fig. 6 but for the observed SLACS and GALLERY lenses. The inferred external shears have a similar distribution of aligned and anti-aligned shears as the mock data sample, indicating they too may be acting internally. Note the increase in the scale of shear magnitude $\gamma_{\text{ext}}$ and elliptical deviations $a_4/a$ compared to the mock data sample.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** For mock lenses that were simulated without external shear. Angle between the best-fit external shear and the major axis of the lens mass distribution, $|\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{ext}}^{\text{PL+ext}}|$ in degrees, as a function of the change in orientation of the major axis when fitting with and without an external shear, $|\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{mass}}^{\text{PL}}|$. Values are the difference in power-law slope inferred between the models fitted with and without an external shear ($\gamma_{\text{PL+ext}} - \gamma_{\text{PL}}$). Systems with best-fit shear aligned to the mass systematically decrease in power-law slope when the external shear is removed from the model. Anti-aligned shears exhibit the opposite behaviour.</figcaption></figure>

PL+ext model. The maximum amplitude of external shear that can be attributed to these twists is around 5 per cent to 8 per cent (see figs B3 and B4 of Van de Vyvere et al. 2022), but lower values are more common for most of their mock lenses. Given many of our strong lenses have shear amplitudes above ~10 per cent, twists alone cannot explain all of the SL shears measured.

## 6.2 External shear biases other SL model parameters

Since best-fitting values of ‘external’ shear may not (entirely) represent the physical quantity they are imagined to, we now test which other parameters in the mass model are biased by their inclusion, and which are still robustly measured. As an extreme alternative, we repeat all measurements but fix $\gamma_1^{\text{ext}} = \gamma_2^{\text{ext}} = 0$ when we refit the mass distribution.

### 6.2.1 Mock lenses

The orientation of the mass distribution is robustly measured, with mean difference $\langle|\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{mass}}^{\text{PL}}|\rangle \sim 1^{\circ}$ and only $\sim 1^{\circ}$ of scatter when shear is excluded (Fig. 8). The axis ratio is also hardly changed, with $\langle q^{\text{PL+ext}} - q^{\text{PL}}\rangle = 0.01 \pm 0.05$. These values are so robust because all the MGE components share a common axis, which is therefore well-defined.

However, the Einstein radii and slopes of the power-law mass model are systematically biased, by an amount that depends on the relative orientation of the shear and the mass (Einstein radii and power-law slopes are known to be degenerate parameters: see e.g. fig. 5 of E23). For lenses whose shears were aligned with the mass distribution, removing $\gamma^{\text{ext}}$ increases the mean best-fit power-law slope by $0.20 \pm 0.04$ (blue points in Fig. 8; their mean best-fit Einstein radii decrease by $0.2 \pm 0.03$ per cent). For lenses whose shears were anti-aligned, the power-law slopes decrease by $0.07 \pm 0.03$ (red points in Fig. 8; their mean best-fit Einstein radii increase by $0.08 \pm 0.04$). Across the entire sample, Einstein radii had been correctly measured when including $\gamma^{\text{ext}}$ (within $0.05 \pm 0.17$ per cent; C22), but are systematically underestimated by $0.2 \pm 0.05$ per cent if shear is excluded.

Our measurements might be caused by a bias described by Kochanek (2021) in the radial structure of a lens, when a model has too few azimuthal degrees of freedom. Kochanek provides the example of fitting a power-law model to a lens whose ellipticity increases with radius: the density slope is forced to spuriously shallow values to balance the shear inside the Einstein radius relative to the shear outside it. Furthermore, Van de Vyvere et al. (2022) found that decreasing ellipticity outside the Einstein radius spuriously increases $\gamma^{\text{ext}}$ for a power-law model, while ellipticity gradients inside or at the Einstein radius mainly bias the power-law slope. In our mocks, the distribution of mass has an ellipticity that varies as a function of radius. We attempted to measure correlations of measured shear properties (e.g. position angle) with the ellipticity of the mock mass distributions. However, we were unable to detect any statistically significant correlation. We interpret this as the measured shear depending on other properties of each lens (e.g. where the lensed source’s arcs appear in the image-plane). A more detailed investigation is beyond the scope of this work.

### 6.2.2 Real lenses

For SLACS and GALLERY lenses, the orientation of the mass distribution changes considerably when external shear is removed, with $\langle|\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{mass}}^{\text{PL}}|\rangle \sim 27^{\circ}$ and $\sim 27^{\circ}$ scatter (Fig. 9; note the different scale on the horizontal axis to Fig. 8). The best-fit axis ratio (indicated by colour in Fig. 9) also increases by $0.18 \pm 0.016$ for aligned systems, which become more spherical; decreases by $0.08 \pm 0.10$ for anti-aligned systems, which become more elliptical; and is inconsistent for the rest, with mean change $\langle q^{\text{PL+ext}} - q^{\text{PL}}\rangle = -0.02 \pm 0.04$.

Neither the Einstein radii nor slopes of the power-law mass model are systematically biased when external shear is removed. We find $\langle\Delta R_{\text{Ein}}/R_{\text{Ein}}\rangle = -0.013^{+0.030}_{-0.043}$, where $\Delta R_{\text{Ein}} \equiv R_{\text{Ein}}^{\text{PL+ext}} - R_{\text{Ein}}^{\text{PL}}$. However, the best-fitting values scatter about the same mean, such that $\langle\Delta R_{\text{Ein}}^2/R_{\text{Ein}}\rangle = 0.045^2$, which is 2 orders of magnitude larger than the equivalent $0.0003^2$ value for the mocks (Fig. 10). Furthermore, the scatter increases with increasing external shear magnitude (Fig. 11).

The best-fit centre of the mass distribution moves on average by $-0.031 \pm 0.061$ arcsec when external shear parameters are
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Orientation angle offset of the external shear from the PL+ext mass distribution ($\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{ext}}^{\text{PL+ext}}$) as a function of the difference in orientation angle when the mass distribution is fitted with and without an external shear ($\phi_{\text{mass}}^{\text{PL+ext}} - \phi_{\text{mass}}^{\text{PL}}$) for the observed SLACS and GALLERY samples. Values are the difference in axis ratio of the models fitted with a PL+ext and a PL ($q^{\text{PL+ext}} - q^{\text{PL}}$).</figcaption></figure>

introduced, and on average remains $0.05 \pm 0.06$ arcsec from the centre of light. With or without shear, this offset is an order of magnitude greater than in the mock data.

We suspect this indicates that the distribution of mass in real galaxies is more complex than that in the mocks: for example with multiple components that are rotationally offset from one another and no single, well-defined axis of ellipticity.

# 7 INTERPRETATION

## 7.1 Strong lensing external shears are not measuring shear

We find that external shear, as measured by galaxy–galaxy strong lensing, does not correlate well with the true shear along a line-of-sight (as measured independently by weak lensing; Section 5.1 and Fig. 2). Best-fitting values of external shear are also frequently several times higher than that along typical lines of sight through the Universe, so only a small fraction of it could be due to true shear. Rather, external shear tends to align with either the major axis or minor axis of the lens mass distribution (Section 5.2). It appears to be compensating for the inflexibility of typical mass models (here an elliptical power-law) to represent the complex distributions of mass. A substantial portion of that complexity appears to be discinity and boxiness (Section 6.1), especially in mock data. In real data, isophotal twists, elliptical gradients, and offsets in the centres and alignments of dark and stellar matter all increase uncertainty on model parameters. These have all been seen in the stellar mass of SLACS lenses (Nightingale et al. 2019). We have not been able to quantify how much each source contributes to the scatter or bias in a real measurement of external shear.

Hogg et al. (2022) suggest a different, ‘minimal line-of-sight’ way of parametrizing the shear that is less degenerate with lens model parameters, but this is still subject to biases in the shear parameters when simplifying assumptions are made for the lens model. For complex distributions of lens mass, the false inference of external shear will pose a challenge for efforts to use strong lensing to measure cosmic shear (e.g. Birrer et al. 2017; Fleury et al. 2021).

## 7.2 Implications for strong lensing science goals

Including external shear alters the best-fitting values of other parameters (Section 6.2). In mock data, some parameters do move closer to the known truth: the mean error in power-law slopes is $-0.02 \pm 0.10$ with shear, or $-0.14 \pm 0.20$ without it. A $\sim 0.2$ per cent bias in Einstein radii also disappears with shear. However, the shear is not physically real and, contributing zero convergence, does not correspond to a physically meaningful distribution of mass. For studies that rely on accurate reconstructions of the mass distribution (e.g. galaxy evolution, dark matter physics, and the Hubble constant), this degeneracy with key parameters will eventually limit statistical precision. For example a key result of the SLACS survey was that elliptical galaxies have isothermal ($\gamma = 2$) mass profiles (Gavazzi et al. 2007; Vegetti & Koopmans 2009; Auger et al. 2010). Including external shear, E23 measured $\langle\gamma\rangle = 2.0756_{-0.024}^{+0.023}$ for SLACS galaxies; here we measure $\langle\gamma\rangle = 2.0159_{-0.032}^{+0.027}$ (Fig. 10), highlighting the systematic uncertainty that will remain fixed even if the sample size increases. Furthermore, a PL+ext model leaves false detections of sub-haloes in a mock discy galaxy (He et al. 2022c) or real HST data (Nightingale et al. 2024).

The Hubble constant $H_0$ can be measured from the time-delay between multiple images (Suyu et al. 2017; Birrer et al. 2019; Wong et al. 2019). However, assuming specific functional forms for the mass model can artificially break the mass-sheet transformation. Our tests on mock data (Section 6.2.1) demonstrate a coupling, predicted by Kochanek (2021), between angular and radial structure. Oversimple models bias the slope, and hence bias $H_0$. C22 estimated $\sim 9$ per cent bias in $H_0$ when using a PL+ext model to interpret time-delays generated from gNFW + MGE lenses. The angular degrees of freedom added by ‘external’ shear are insufficient to compensate for even this complexity. Given that our analysis of SLACS and GALLERY lenses indicate even more angular degrees of freedom, such as twists in the mass distribution, biases for real lensing systems will likely be closer to the 20–50 per cent suggested by other studies (Schneider & Sluse 2013; Xu et al. 2016; Gomer & Williams 2018, 2020, 2021; Kochanek 2020). More flexible models, such as adding an internal mass sheet to the PL+ext model which is constrained via stellar dynamics, have been introduced to mitigate these systematics.

Combining $H_0$ measurements from a large population of lenses might average away individual biases up to $10 \text{ kms}^{-1} \text{ Mpc}^{-1}$, on the assumption that boxy and discy mass distributions are equally well represented (Van de Vyvere et al. 2022). However, the discinity of the mass distribution correlates with the discinity of the (observable) PL+ext critical curves (Section 6.1.2), and our observationally selected sample contains an overpopulation of 71 per cent discy galaxies (Section 6.1.3). More flexible mass models (e.g. a PL with an internal mass sheet) have also been introduced (Birrer et al. 2020), which infer unbiased $H_0$ values on mock lens samples (Ding et al. 2021). Moreover, Van de Vyvere et al. (2022) did not investigate the effect of $b_4$ ‘twisting’ perturbations (Section 6.1.4) or mis-centring (Section 6.2.2), both of which we observe in the critical curves of real galaxies, and whose effects may not average away. For example the centre of mass in the H0LiCOW model of lens WFI 2033-4723 is offset from the centre of light by $\sim 10 \times$ astrometric uncertainty (Suyu et al. 2010; Barrera et al. 2021), and a similar offset in iPTF16geu increases asymmetry (Diego et al. 2022). Such offsets are unphysical (Schaller et al. 2015), so it is not clear that complexities in the mass distributions of real lenses can be safely ignored by averaging over a population.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Histograms of the difference between inferred PL mass distribution model parameters with external shear as free parameters, minus those without shear, for the observed (wider distributions) and mock data (narrower distributions) samples. From left to right panels the parameters are: fractional Einstein radius, logarithmic slope of density profile, axis ratio of mass distribution, and radial distance of the centre of the mass distribution, in arcsec.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Fractional difference between Einstein radii inferred for models with and without the free parameters for external shear, as a function of the external shear amplitude.</figcaption></figure>

## 7.3 Future work: more complex mass models are needed

External shear does not appear to be sufficient as (the sole) parameter to encode all the complexity in real lenses. Although Einstein radii are expected to be ‘model independent’ within ~2 per cent uncertainty Bolton et al. (2008b) and Sonnenfeld et al. (2013), the 4.5 per cent rms fractional difference we measure with and without shear, suggests an unknown unknown. A single power-law model leads to mass discrepancies with stellar dynamics (Etherington et al. 2022), and spurious false-positives in searches for dark matter subhaloes (Nightingale et al. 2024). Further work (e.g. Van de Vyvere et al. 2020, 2022; C22) to understand the types of asymmetries that must be accounted for in the lens modelling, and the possibility of constraining such models, will be invaluable.

What parametric forms allow sufficient complexity – in a minimum number of parameters that need to be constrained? Even our gNFW + MGE mocks still do not capture the full complexity of real lenses, but their MGEs were forced to be aligned. Perhaps the model could have both $a_4$ and $b_4$ perturbations, varying as a function of radius. Cluster-scale models frequently use a soft core inside some scale radius, and (Limousin et al. 2022) add B-spline functions. The number of free parameters must be balanced with the information available: pixellated mass models are generally underconstrained (and although they might be able to fit observations without external shear; Valls-Gabaud et al. 2006, some line-of-sight shear is *expected* for galaxy-scale lenses). To instead increase the available information, the total mass could be decomposed into dark matter and stellar components, with the latter informed by the lens light (which is otherwise a nuisance). Whatever the parametric form, the azimuthal degrees of freedom must be defined carefully to avoid the bias described by Kochanek (2021) on the inference of the radial mass distribution. We suggest more studies (e.g. Vyvere et al. 2020; Van de Vyvere et al. 2022, C22), embedded deeply in each specific science case, to quantify the impact of simplifying assumptions.

A silver lining is that the sensitivity of galaxy-galaxy strong lensing data may be a new opportunity. Both our results and those of Vyvere et al. (2022) suggest that it might be possible to measure
the disciness or boxiness of galaxies' dark matter haloes. Our results are even more optimistic: if measurable multipole perturbations in critical curves traces those in the distributions of both mass and light, then one could study the dark morphology of galaxies at high-redshift with relative ease.

## SOFTWARE

This work uses the following software packages:

(i) ASTROPY (Astropy Collaboration 2013; Price-Whelan et al. 2018)

(ii) CORNER.PY (Foreman-Mackey 2016)

(iii) DYNESTY (Speagle 2020)

(iv) MATPLOTLIB (Hunter 2007)

(v) NUMBA (Lam, Pitrou & Seibert 2015)

(vi) NUMPY (van der Walt, Colbert & Varoquaux 2011)

(vii) PYAUTOFIT (Nightingale, Hayes & Griffiths 2021a)

(viii) PYAUTOGALAXY (Nightingale et al. 2023a)

(ix) PYAUTOLENS (Nightingale & Dye 2015; Nightingale et al. 2018, 2021b)

(x) PYTHON (Van Rossum & Drake 2009)

(xi) SCIKIT-IMAGE (Van der Walt et al. 2014)

(xii) SCIKIT-LEARN (Pedregosa et al. 2011)

(xiii) SCIPY (Virtanen et al. 2020)

(xiv) SQLite (Hipp 2020)

## ACKNOWLEDGEMENTS

We thank Yiping Shu, Guillaume Desprez, Johan Richard, Mathilde Jauzac, and Keichi Umetsu for providing data from previously published papers, and Liliya Williams and Adi Zitrin for stimulating discussions.

AE was supported by UK STFC via grants ST/R504725/1 and ST/T506047/1. JN and RM are supported by STFC via grant ST/T002565/1, and the UK Space Agency via grant ST/W002612/1. XYC and RL are supported by the National Nature Science Foundation of China (via grants 11988101, 11773032, 12022306), the China Manned Space Project (science research grants CMS-CSST-2021-B01, CMS-CSST-2021-A01) and by the K.C.Wong Education Foundation. QH, AA, SMC, and CSF are supported by the ERC via grant GA 786910. MJ was supported by UKRI via grant MR/S017216/1. This work used both the Cambridge Service for Data Driven Discovery (CSD3) and the DiRAC Data-Centric system, which are operated by the University of Cambridge and Durham University on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). These were funded by BIS capital grant ST/K00042X/1, STFC capital grants ST/P002307/1, ST/R002452/1, ST/H008519/1, ST/K00087X/1, STFC operations grants ST/K003267/1, ST/K003267/1, and Durham University. DiRAC is part of the UK National e-Infrastructure.

## DATA AVAILABILITY

Text files and images of every model-fit performed in this work are available at https://zenodo.org/record/6104823. Full DYNESTY chains of every fit are available upon request.

## REFERENCES

Ade P. A. et al., 2016, A&A, 594, A13

Amorisco N. C. et al., 2022, MNRAS, 510, 2464

Astropy Collaboration, 2013, A&A, 558, A33

Auger M. W., Treu T., Bolton A. S., Gavazzi R., Koopmans L. V., Marshall P. J., Moustakas L. A., Burles S., 2010, ApJ, 724, 511

Barrera B., Williams L. L. R., Coles J. P., Denzel P., 2021, Open J. Astrophysics, 4, 12

Bertin E., Arnouts S., 1996, A&AS, 117, 393

Birrer S., Welschen C., Amara A., Refregier A., 2017, J. Cosmol. Astropart. Phys., 2017, April 2017

Birrer S. et al., 2019, MNRAS, 484, 4726

Birrer S. et al., 2020, A&A, 643, A165

Bolton A. S., Burles S., Koopmans L. V. E., Treu T., Gavazzi R., Moustakas L. A., Wayth R., Schlegel D. J., 2008a, ApJ, 682, 964

Bolton A. S., Treu T., Koopmans L. V. E., Gavazzi R., Moustakas L. A., Burles S., Schlegel D. J., Wayth R., 2008b, ApJ, 684, 248

Bolton A. S. et al., 2012, ApJ, 757, 82

Cao X. et al., 2022, RAA, 22, 30 C22

Cappellari M., 2002, MNRAS, 333, 400

Cappellari M. et al., 2013, MNRAS, 432, 1709

Chaware L., Cannon R., Kembhavi A. K., Mahabal A., Pandey S. K., 2014, ApJ, 787, 102

Chu Z., Lin W. P., Li G. L., Kang X., 2013, ApJ, 765, 134

Despali G., Lovell M., Vegetti S., Crain R. A., Oppenheimer B. D., 2019, MNRAS, 17, 1

Desprez G., Richard J., Jauzac M., Martinez J., Siana B., Clément B., 2018, MNRAS, 479, 2630

Diego J. M., Bernstein G., Chen W., Goobar A., Johansson J. P., Kelly P. L., Mörtsell E., Nightingale J. W., 2022, A&A, 662, A34

Ding X. et al., 2021, MNRAS, 503, 1096

Etherington A. et al., 2022, MNRAS, 517, 3275

Etherington A. et al., 2023, MNRAS, 521, 6005 E22

Faure C. et al., 2008, ApJS, 176, 19

Fleury P., Larena J., Uzan J.-P., 2021, JCAP, 2021, 024

Foreman-Mackey D., 2016, J. Open Source Softw., 1, 24

Gavazzi R., Treu T., Rhodes J. D., Koopmans L. V. E., Bolton A. S., Burles S., Massey R. J., Moustakas L. A., 2007, ApJ, 667, 176

Gomer M. R., Williams L. L. R., 2018, MNRAS, 475, 1987

Gomer M., Williams L. L. R., 2020, JCAP, 2020, 045

Gomer M. R., Williams L. L. R., 2021, MNRAS, 504, 1340

Graham A. W., Driver S. P., 2005, Publ. Astron. Soc. Aust., 22, 118

Harvey D., 2020, MNRAS, 498, 2871

Harvey D., Tam S.-I., Jauzac M., Massey R., Rhodes J., 2019, preprint (arXiv:1911.06333)

He Q. et al., 2022a, MNRAS, 511, 3046

He Q. et al., 2022b, MNRAS, 512, 5862

He Q. et al., 2022c, MNRAS, 518, 220

Hezaveh Y. D. et al., 2016, ApJ, 823, 37

Hilbert S., White S. D. M., Hartlap J., Schneider P., 2007, MNRAS, 382, 121

Hipp R. D., 2020, SQLite, https://www.sqlite.org/index.html

Hogg N. B., Fleury P., Larena J., Martinelli M., 2022, MNRAS, 520, 5982

Hunter J. D., 2007, Comput. Sci. Eng, 9, 90

Keeton C. R., Kochanek C. S., Seljak U., 1997, ApJ, 482, 604

Kilbinger M., 2015, Rep. Prog. Phys., 78, 086901

Kochanek C. S., 2020, MNRAS, 493, 1725

Kochanek C. S., 2021, MNRAS, 501, 5021

Kochanek C., Schneider P., Wambsganss J., in 2004, Meylan G., Jetzer P., North P., eds, Part 2 of Gravitational Lensing: Strong, Weak and Micro, Proceedings of the 33rd Saas-Fee Advanced Course, Springer, Berlin

Koopmans L. V. et al., 2009, ApJ, 703, L51

Kuhn F. A., Bruderer C., Birrer S., Amara A., Réfrégier A., 2020, JCAP, 04, 24

Lam S. K., Pitrou A., Seibert S., 2015, Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC - LLVM, 1

Leauthaud A. et al., 2007, ApJS, 172, 219

Li R., Frenk C. S., Cole S., Gao L., Bose S., Hellwing W. A., 2016, MNRAS, 460, 363

Li R., Frenk C. S., Cole S., Wang Q., Gao L., 2017, MNRAS, 468, 1426

Li R. et al., 2019, MNRAS, 490, 2124
Limousin M., Beauchesne B., Jullo E., 2022, A&A, 664, A90

Massey R. et al., 2007, ApJS, 172, 239

Massey R., Kitching T., Richard J., 2010, Rep. Prog. Phys., 73, 086901

Narayan R., Bartelmann M., 1996, preprint(astro-ph/9606001)

Natarajan P., Springel V., 2004, ApJ, 617, L13

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Nightingale J. W., Dye S., 2015, MNRAS, 452, 2940

Nightingale J. W., Dye S., Massey R. J., 2018, MNRAS, 478, 4738

Nightingale J. W., Massey R. J., Harvey D. R., Cooper A. P., Etherington A.,
Tam S. I., Hayes R. G., 2019, MNRAS, 489, 2049

Nightingale J. W., Hayes R. G., Griffiths M., 2021a, J. Open Source Softw.,
6, 2550

Nightingale J. et al., 2021b, J Open Source Softw., 6, 2825

Nightingale J. W. et al., 2023a, J. Open Source Softw., 8, 4475

Nightingale J. W. et al., 2023b, MNRAS, 521, 3298

Nightingale J. W. et al., 2024, MNRAS, 527, 10480

O'Riordan C. M., Vegetti S., 2024, MNRAS, 528, 1757

Oldham L. J., Auger M. W., 2018, MNRAS, 476, 133

Pedregosa F. et al., 2011, J Mach. Learn. Res., 12, 2825

Price-Whelan A. M. et al., 2018, AJ, 156, 123

Rhodes J., Refregier A., Groth E. J., 2000, ApJ, 536, 79

Ritondale E., Vegetti S., Despali G., Auger M. W., Koopmans L. V., McKean
J. P., 2019, MNRAS, 485, 2179

Robertson A., Smith G. P., Massey R., Eke V., Jauzac M., Bianconi M.,
Ryczanowski D., 2020, MNRAS, 495, 3727

Ryden B. S., 1992, ApJ, 386, 42

Schaller M., Robertson A., Massey R., Bower R. G., Eke V. R., 2015,
MNRAS, 453, L58

Schneider P., Sluse D., 2013, A&A, 559, A37

Shu Y. et al., 2015, ApJ, 803, 73

Shu Y., Bolton A. S., Moustakas L. A., Stern D., Dey A., Brownstein J. R.,
Burles S., Spinrad H., 2016a, ApJ, 820, 43

Shu Y. et al., 2016b, ApJ, 824, 86

Smail I., Ellis R. S., Fitchett M. J., 1994, MNRAS, 270, 245

Sonnenfeld A., Treu T., Gavazzi R., Marshall P. J., Auger M. W., Suyu S. H.,
Koopmans L. V., Bolton A. S., 2012, ApJ, 752, 163

Sonnenfeld A., Treu T., Gavazzi R., Suyu S. H., Marshall P. J., Auger M. W.,
Nipoti C., 2013, ApJ, 752, 163

Speagle J. S., 2020, MNRAS, 493, 3132

Suyu S. H., 2012, MNRAS, 426, 868

Suyu S. H., Marshall P. J., Auger M. W., Hilbert S., Blandford R. D.,
Koopmans L. V. E., Fassnacht C. D., Treu T., 2010, ApJ, 711, 201

Suyu S. H. et al., 2017, MNRAS, 468, 2590

Tam S. I. et al., 2020, MNRAS, 496, 4032

Treu T., Gavazzi R., Gorecki A., Marshall P. J., Koopmans L. V., Bolton A.
S., Moustakas L. A., Burles S., 2009, ApJ, 690, 670

Valageas P., Barber A. J., Munshi D., 2004, MNRAS, 347, 654

Valls-Gabaud D., Cabanac R., Lidman C., Diego J. M., Saha P., 2006, in
Mamon G. A., Combes F., Deffayet C., Fort B. eds, EAS Publ. Ser. Vol.
20. p. 149

Van de Vyvere L., Sluse D., Mukherjee S., Xu D., Birrer S., 2020, A&A, 644,
A108

Van de Vyvere L., Sluse D., Gomer M. R., Mukherjee S., 2022a, A&A, 663,
A179

Van de Vyvere L., Gomer M. R., Sluse D., Xu D., Birrer S., Galan A.,
Vernardos G., 2022b, A&A, 659, A127

van der Walt S., Colbert S. C., Varoquaux G., 2011, Comput Sci Eng,
13, 22

Van der Walt S., Schönberger J. L., Nunez-Iglesias J., Boulogne F., Warner
J. D., Yager N., Gouillart E., Yu T., 2014, PeerJ, 2, e453

Van Rossum G., Drake F. L., 2009, Python 3 Reference Manual. CreateSpace,
Scotts Valley, CA

Vegetti S., Koopmans L. V., 2009, MNRAS, 400, 1583

Vegetti S., Koopmans L. V., Bolton A., Treu T., Gavazzi R., 2010, MNRAS,
408, 1969

Vegetti S., Koopmans L. V. E., Auger M. W., Treu T., Bolton A. S., 2014,
MNRAS, 442, 2017

Virtanen P. et al., 2020, Nat. Met., 17, 261

Witt H. J., Mao S., 1997, MNRAS, 291, 211

Wong K. C., Keeton C. R., Williams K. A., Momcheva I. G., Zabludoff A. I.,
2011, ApJ, 726, 84

Wong K. C. et al., 2019, MNRAS, 498, 1420

Xu D., Sluse D., Schneider P., Springel V., Vogelsberger M., Nelson D.,
Hernquist L., 2016, MNRAS, 456, 739

Zhang B., Lee K.-G., Krolewski A., Shi J., Horowitz B., Kooistra R., 2023,
ApJ, 954, 49

Zhao H., 1996, MNRAS, 278, 488

APPENDIX A: SYSTEMATIC TESTS OF SHEAR
MEASUREMENT

Real external shear is caused by galaxies or clusters of galaxies
adjacent to the line of sight (Wong et al. 2011). It is possible
that individual galaxies impart a localized external shear that is
unmeasurable by weak lensing, which averages over large spatial
scales. Here, we repeat our analysis of every SLACS and GALLERY
lens, but including adjacent galaxies in the mass model. Crucially,
we use identical lens light subtraction and source pixellization, so
we can directly compare the inferred shears.

We select which adjacent galaxies to include subjectively, based
on their proximity and size. To implement this, we extend this study's
LP<sup>1</sup> pipeline (Etherington et al. 2022) with a GUI where we look at
10 arcsec cut-outs of each lens and click on up to two galaxies. These
are added to the mass model as spherical isothermal spheres (SIS,
an SIE with *q* = 1), fixed to the centre of their brightest pixel and
the same redshift as the main lens. Their Einstein radius *θ*<sub>E</sub><sup>mass</sup> is
given a flat prior from 0.0 to 0.5 arcsec and fitted with all other
parameters. *θ*<sub>E</sub><sup>mass</sup> = 0.5 arcsec corresponds to a mass far greater
than that suggested by each galaxy's luminosity. It is an intentional
choice not to use more informative priors, so that we can investigate
how line-of-sight galaxies change the shear inference with maximal
freedom.

Adjacent galaxies do not significantly change the inferred external
shear (Fig. A1). For almost every lens, the best-fit external shear for
models with and without adjacent galaxies remains consistent. Large
amplitudes of external shear are still preferred for many lenses, and
their position angles still show the same alignment with the mass
model position angles (not shown). The very largest external shears
γ<sub>ext</sub> ≳ 0.3 may be reduced, with a small statistical sample but with
the two highest (~0.40 and ~0.36) becoming ~0.21 and ~0.24
after accounting for adjacent galaxies. Such values remain physically
infeasible. These results confirm that line-of-sight galaxies cannot
explain the conclusions of this paper.
<figure><img src="image_16.png" /></figure>

**Figure A1.** Best-fit external shear amplitude $\gamma^{\text{ext}}$ and position angle $\phi^{\text{ext}}$ for every lens in the sample. Shear values are shown for the lens models presented in the main paper which do not include nearby line-of-sight galaxies and for an additional set of lens model fitted, which include up to two nearby galaxies as SISs. Error bars are at $1\sigma$ confidence intervals. The lens light subtraction and source analysis of every fit are the same, enabling a direct comparison of how the shear values change when nearby galaxies are included in the mass model. The inferred $\gamma^{\text{ext}}$ and $\phi^{\text{ext}}$ values are consistent with one another in the vast majority of lenses, irrespective of whether line-of-sight galaxies are included. In the right figure, models with $\gamma^{\text{ext}} < 0.05$ are shown in black, illustrating that the measured values are only inconsistent for low shear lenses. Lensing due to nearby galaxies therefore is not responsible for the large shears which align with the lens mass model found in the main paper.

<figure><img src="image_17.png" /></figure>

**Figure A2.** Best-fit external shear amplitude $\gamma^{\text{ext}}$ and position angle $\phi^{\text{ext}}$ for every lens in the sample. Shear values are shown for the lens models presented in the main paper which use a Voronoi mesh with luminosity based regularization and for an additional set of model fitted, which use a Delaunay mesh with interpolation regularization (see appendix A of Nightingale et al. 2023b). Error bars are at $1\sigma$ confidence intervals. In the right figure, models with $\gamma^{\text{ext}} < 0.05$ are shown in black, illustrating that the measured values are only inconsistent for low shear lenses. The inferred $\gamma^{\text{ext}}$ and $\phi^{\text{ext}}$ values are consistent with one another in the vast majority of lenses, irrespective of the source model. The source analysis therefore does not impact our main conclusions.

## APPENDIX B: SENSITIVITY TO SOURCE MODEL

We investigate whether the shear measurements presented in the main paper depend on assumptions related to the source model. In the main paper, we followed Etherington et al. (2022) in reconstructing the source using a Voronoi pixel mesh and a regularization scheme that adapted to the unlensed source luminosity (see Nightingale et al. 2018). As an alternative, here, we model the same lenses using

a Delaunay pixel mesh and the interpolation-based regularization scheme described in appendix A of Nightingale et al. (2024).

The best-fit external shear with this alternative source-plane pixelization has amplitude and angle that is consistent with those in the main paper, for nearly every lens (see Fig. A2). Our main conclusions are thus unaffected by even substantial changes in source reconstruction.

This paper has been typeset from a TeX/LaTeX file prepared by the author.