11institutetext:  Institute of Physics, Laboratory of Astrophysics, Ecole
Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290
Versoix, Switzerland 22institutetext: AIM, CEA, CNRS, Université Paris-Saclay,
Université de Paris, F-91191 Gif-sur-Yvette, France

# Using wavelets to capture deviations from smoothness
in galaxy-scale strong lenses

A. Galan Using wavelets to capture deviations from smoothness in galaxy-scale
strong lensesUsing wavelets to capture deviations from smoothness in galaxy-
scale strong lenses G. Vernardos Using wavelets to capture deviations from
smoothness in galaxy-scale strong lensesUsing wavelets to capture deviations
from smoothness in galaxy-scale strong lenses A. Peel Using wavelets to
capture deviations from smoothness in galaxy-scale strong lensesUsing wavelets
to capture deviations from smoothness in galaxy-scale strong lenses F. Courbin
Using wavelets to capture deviations from smoothness in galaxy-scale strong
lensesUsing wavelets to capture deviations from smoothness in galaxy-scale
strong lenses J.-L. Starck Using wavelets to capture deviations from
smoothness in galaxy-scale strong lensesUsing wavelets to capture deviations
from smoothness in galaxy-scale strong lenses

Modeling the mass distribution of galaxy-scale strong gravitational lenses is
a task of increasing difficulty. The high-resolution and depth of imaging data
now available render simple analytical forms ineffective at capturing lens
structures spanning a large range in spatial scale, mass scale, and
morphology. In this work, we address the problem with a novel multi-scale
method based on wavelets. We test our method on simulated Hubble Space
Telescope imaging data of strong lenses containing different types of mass
substructures making them deviate from smooth models: (1) a localized small
dark matter subhalo, (2) a Gaussian random field that mimics a non-localized
population of sub-halos along the line of sight, (3) galaxy-scale multipoles
that break elliptical symmetry. We show that wavelets are able to recover all
of these structures accurately. This is made technically possible by using
gradient-informed optimization based on automatic differentiation over
thousands of parameters, also allowing us to sample the posterior
distributions of all model parameters simultaneously. By construction, our
method merges all current modeling paradigms—analytical, pixelated, and deep
learning—into a single modular framework. It is also well-suited for the fast
modeling of large samples of lenses. All methods presented here are publicly
available in our new Herculens package .

###### Key Words.:

Cosmology: dark matter – Galaxies: structure – Gravitation – Gravitational
lensing: strong – Methods: data analysis

##  1 Introduction

The Lambda cold dark matter cosmological model encapsulates our current
understanding of the Universe, accurately explaining a number of observations
at large scales ($>\\!10$ Mpc), such as the cosmic microwave background
temperature and polarization anisotropies (Planck Collaboration et al. 2020),
baryon acoustic oscillations (e.g. Raichoor et al. 2021), and the accelerated
expansion of the Universe (e.g. Riess 2020). The dark matter (DM) component of
this model plays a major role in the hierarchical collapse of matter due to
gravitational instability, which eventually produced the galaxies populated by
stars that we observe today (e.g. Toomre & Toomre 1972; Dubinski 1994;
Springel et al. 2006). Despite its successes at large scales, the effect of DM
on galactic and sub-galactic scales ($<10$ Mpc), where its interplay with
baryons is mediated by non-linear and poorly understood mechanisms (e.g.
Blumenthal et al. 1986; Zubovas & King 2012), poses a challenge in explaining
observations (e.g. Vogelsberger et al. 2020). For instance, predictions from
the cold DM model lead to overly cuspy central density profiles (the so
called“cusp-core problem”, Moore 1994; de Blok 2010), too many low-mass
($\lesssim\\!10^{8}\ {\rm M}_{\odot}$) dwarf galaxies (the “missing satellites
problem”, Moore et al. 1999; Klypin et al. 1999) and intermediate-mass
($\sim\\!10^{10}\ {\rm M}_{\odot}$) galaxies populated by too few stars (the
“too-big-to-fail problem”, Boylan-Kolchin et al. 2011; Papastergis et al.
2015).

One path to reconcile theory with observations is to examine alternative DM
models that assume different properties of the DM particle. For example, warm
DM simulations (e.g. Li et al. 2017) have been able to produce far fewer dwarf
galaxy satellites, in line with observations. Another path is to improve our
ability to detect observational signatures of DM through its gravitational
influence on baryonic matter, like those observed in the local Universe within
the Milky Way (e.g. with stellar streams, Erkal et al. 2016) and its
neighborhood (e.g. within dwarf galaxy satellites, Nadler et al. 2021). At
cosmological distances, gravitational lensing is a direct and more efficient
probe of the mass in galaxies, and has the potential to measure DM effects
down to kpc in size and $10^{7}\ {\rm M}_{\odot}$ in mass (e.g. Hezaveh et al.
2016a; Chatterjee & Koopmans 2018; Gilman et al. 2019). This is a crucial
range of mass were most DM models tend to disagree in their predictions.

Galaxy-galaxy strong lensing occurs when two galaxies at different redshifts
align along the line of sight, with the foreground galaxy (the lens)
deflecting incoming light rays of the background one (the source). From the
observer’s perspective, the source appears magnified, distorted, and split
into multiple images. By modeling the observed lensed features through the
well-established physics of the lens equation, one can constrain the total
mass distribution in galaxies and measure their DM content and formation
history. However, lens modeling is an under-constrained problem, mostly
because both the lens mass and the source light distributions are a priori
unknown, and the lensed source is often blended with the lens light.
Assumptions are therefore needed to reconstruct the lensing mass, the lens
light, and the source light distributions by inverting the lens equation, such
as priors on the shape of mass and light profiles, or on their higher order
statistical properties. Currently used analytical profiles are sufficient to
capture first-order properties of the lens mass distribution (e.g. power-law
profiles), but they lack the degrees of freedom to capture those small scale
features that are critical for determining DM properties (e.g. He et al.
2022).

Increasing the complextiy of a lens model is not trivial due to degeneracies
between the lens potential and the source light: complex structure in the
observed lensed features could be attributed to either local perturbations in
the potential or to an intrinsicly complex source light distribution (e.g.
clumpy star forming galaxies or galaxy mergers). Koopmans (2005) used a Taylor
expansion to address this, deriving a perturbative correction to the lens
equation that combines spatial derivatives of both the source and the
potential. This approach extends the semi-linear inversion technique of Warren
& Dye (2003) by adding a perturbing field to the smooth lens potential, which
does not assume any specific shape and can be solved for in the same way as
the source. The lens potential and the source light are discretized and cast
on two grids of pixels, in the lens and source planes respectively, on which
the two fields are reconstructed (see also Vegetti & Koopmans 2009; Suyu et
al. 2009; Vernardos & Koopmans 2022, for extensions of the technique). For the
source, various priors have been explored, including analytical functions
ranging from the Sérsic profile (Sérsic 1963) to shapelet basis sets (Tagore &
Keeton 2014; Birrer et al. 2015), gaussian processes (Karchev et al. 2022) or
deep generative models (Chianese et al. 2020). More complex sources can be
modeled using grids of pixels, combined with curvature-based (Suyu et al.
2006; Nightingale & Dye 2015), gaussian process (Vernardos & Koopmans 2022),
or wavelet-based (Joseph et al. 2019; Galan et al. 2021) regularizations.
Recently, Vernardos & Koopmans (2022) studied the effect of specific prior
assumptions on both the potential perturbations and the source light
distribution, confirming their degeneracy and that the particular choice of
priors plays an important role in recovering the underlying potential
perturbations.

In a real-world application, we do not know a priori the dominating type of
perturbations, such that we need a method that is both flexible yet robust in
recovering their main properties. These perturbations may be due to isolated
subhalos, or subhalo populations along the line of sight. Both cases allow us
to probe the low-mass end of the dark matter mass function, and possibly the
shape of DM halos, for which different DM model predictions disagree. The
gravitational imaging technique has been used for the detection of subhalos
based on the reconstructed potential perturbations. If the pixelated
reconstruction contains a well-localized mass over-density, it is replaced by
an analytical profile such as a Navarro-Frenk-White halo (NFW, Navarro et al.
1996) whose parameters are further optimized (Vegetti & Koopmans 2009). Using
this method, Vegetti et al. (2012) reported the detection of a $10^{8.28}\
{\rm M}_{\odot}$ dark halo, and Vegetti et al. (2014) used non-detections in
eleven systems of the Sloan Lens ACS Survey (SLACS) sample observed with HST
(Bolton et al. 2006) to constrain the mean projected substructure mass
fraction in the context of cold DM. Using a different method based on
comparing model likelihoods between a model that explicitly includes a subhalo
at a given position with a model that does not, Hezaveh et al. (2016b)
reported the detection of a substructure of mass $10^{8.96}\ {\rm M}_{\odot}$
from ALMA observations. In general, constraints from individual subhalo
detections can mainly be improved with the analysis of larger samples of
lenses (Vegetti et al. 2014), and depend on the shape of the profile used to
describe the subhalo mass distribution (Despali et al. 2022).

These direct detection methods quickly become too computationally expensive
for detecting $\gtrsim 2$ subhalos. Considering populations of subhalos
instead allows one to probe even lower subhalo masses ($\lesssim 10^{7}\ {\rm
M}_{\odot}$) through their collective effects, since low-mass subhalos are
predicted to be more numerous (Hezaveh et al. 2016a). A population of subhalos
is described in a statistical way, but its properties are directly related to
those of DM models, which can thus be constrained (e.g. the thermal relic mass
of the warm DM particle, Birrer et al. 2017). One such statistical description
of the perturbed lensing potential has been introduced in Chatterjee &
Koopmans (2018), using a power-spectrum analysis of the (lensed) source
surface brightness fluctuations, reaching a sensitivity down to a few kpc in
the subhalo mass power-spectrum (see also Bayer et al. 2018). Similarly, Cyr-
Racine et al. (2019) decomposed model residuals into Fourier modes and linked
them to the substructure power-spectrum. More recently, Vernardos & Koopmans
(2022) used the gravitational imaging technique to reconstruct the perturbing
field of a population of subhalos and recovered its power-spectrum properties,
especially the slope, remarkably well. Several studies have also employed deep
learning techniques to infer the presence of subhalo populations (e.g. Brehmer
et al. 2019; Diaz Rivero & Dvorkin 2020; Varma et al. 2020; Coogan et al.
2020; Vernardos et al. 2020; Ostdiek et al. 2022; Adam et al. 2022). While it
is still unclear if these methods are strongly limited by the simplifying
assumptions on their training data (e.g. the absence of lens light, fixed
instrumental properties, etc.), they are a promising path forward to speed up
computations for application to large samples of lenses, possibly in
combination with more classical approaches.

Independently of the presence of subhalos, deviations from smoothness can
occur within the lens galaxy mass distribution itself, manifesting as
additional radial or azimuthal structures. For instance, deviations along the
radial direction can be mass-to-light radial gradients (Oldham & Auger 2018;
Sonnenfeld et al. 2018; Shajib et al. 2021), while angular structures can be
the consequence of ellipticity gradients and twists (Keeton et al. 2000; Van
de Vyvere et al. 2022). These cannot be captured by the typically employed
elliptical power-law mass models. Free-form techniques have been one way to
address these limitations by dismissing the smooth component of the potential
and relying solely on a grid of mass pixels governed by a set of priors
(either physically or mathematically motivated) to prevent the appearance of
un-physical mass distributions (Saha & Williams 1997; Coles et al. 2014).
However, retaining the smooth component already provides a reliable first-
order model, from which to explore higher order deviations. Azimuthal
structures can be described by higher-order multipoles that have been
identified in stellar populations from both real observations and cosmological
simulations (Trotter et al. 2000; Claeskens et al. 2006). These can be
explained by AGN feedback that suppresses the formation of discs in massive
galaxies, as they tend to evolve from disky to elliptical or even boxy shapes
(i.e. quadrupoles, Frigo et al. 2019). Recently, Van de Vyvere et al. (2021)
demonstrated that quadrupole moments of low amplitude, based on results of Hao
et al. (2006), can be detected in current HST observations of strong lenses,
although their detectability depends on numerous factors such as their
alignment with the smooth potential or the degrees of freedom in the source
model. Using very long baseline interferometric observations of a strong lens,
Powell et al. (2022) recently reported the detection of multipole structures
beyond ellipticity in the deflector mass distribution. While not accounting
for those multipoles in the model can bias the measurement of the lens mass
profile up to a few percent (Van de Vyvere et al. 2021; Powell et al. 2022),
their detection most importantly holds valuable information on the formation
history of galaxies, as signatures of their past evolution.

Our goal in this paper is to unify the modeling of generic lens potential
perturbations in a robust framework that includes, but is not limited to, the
three categories presented above: individual or populations of subhalos and
higher order moments in the lens mass distribution. To achieve this, we extend
our previous work in Galan et al. (2021) by including a reconstruction of lens
potential perturbations regularized with wavelets, and demonstrate that our
method can successfully reconstruct perturbed lens potentials of different
origin. Our technique benefits from the multi-scale properties of the wavelet
transform, along with well-motivated sparsity constraints to reconstruct the
various spatial scales over which perturbations to the smooth potential can
occur. We implement our method using a fully differentiable algorithm based on
automatic differentiation, first introduced by Wengert (1964). This
programming framework gives direct access to all the derivatives of the highly
non-linear loss function at a negligible computational cost, which enables the
use of robust gradient-informed algorithms to explore the parameter space and
optimize the model parameters. As a result, convergence to the maximum-a-
posteriori (MAP) solution is fast and first-order error estimates are obtained
through Fisher matrix analysis. An efficient exploration of the parameter
space for estimating posterior distributions is then performed via gradient-
informed Hamiltonian Monte-Carlo (HMC) sampling.

We present our methodology in Sect. 2. The simulated examples used to
demonstrate the capabilities of our method are presented in Sect. 3. We
perform the reconstruction of perturbations to the lens potential by uniformly
modeling these examples and present our results in Sect. 4. We then evaluate
the reconstructions of the perturbations in Sect. 5. We conclude this work and
discuss its future prospects in Sect. 7.

##  2 Methodology

Figure 1: Flow chart of the fully differentiable method used in this work and
implemented in Herculens, which merges analytical profiles with pixelated
components. We indicate the typical number of parameters $N$ for the
analytical potential ${\boldsymbol{\eta}}_{\psi}$ and the pixelated potential
${\boldsymbol{\psi}_{\rm pix}}$. Partial derivatives are computed and
propagated for all parameters. These derivatives are then used to update
parameter values in the direction indicated by the gradient of the loss
function $L$. We note that the blurring operator $\boldsymbol{\mathsf{B}}$
(see Eq. 4) is also used to generate the model ${\boldsymbol{m}}$ but is not
shown in the diagram to avoid clutter. This work focuses on modeling the
lensing potential ${\boldsymbol{\psi}}$, and we refer to Joseph et al. (2019)
and Galan et al. (2021) for the modeling of surface brightness distributions.

In this section we introduce the gravitational lensing formalism and describe
in detail the various aspects of our method. Fig. 1 summarizes the main
components of the method and can be referred to throughout this section. Some
important mathematical notations are also summarized Table 1.

Table 1: List of the main notations used in the paper. Notation | Definition
---|---
Observation
${\boldsymbol{d}}$ | Imaging data
${\rm\bf C}_{\boldsymbol{d}}$ | Data covariance matrix
Model components
${\boldsymbol{m}}$ | Image model
$\boldsymbol{\mathsf{B}}$ | Blurring operator
$\boldsymbol{\mathsf{L}}_{{\boldsymbol{\psi}}}$ | Lensing operator
${\boldsymbol{\psi}}$, ${\boldsymbol{s}}$, ${\boldsymbol{\ell}}$ | Lens potential, source light, lens light
Model parametrization
${\boldsymbol{\eta}}$ | Vector containing all model parameters
${\boldsymbol{\eta}}_{{\boldsymbol{\psi}}}$, ${\boldsymbol{\eta}}_{{\boldsymbol{s}}}$, ${\boldsymbol{\eta}}_{{\boldsymbol{\ell}}}$ |  Parameters for analytical profiles $\tilde{{\boldsymbol{\psi}}}$, $\tilde{{\boldsymbol{s}}}$, $\tilde{{\boldsymbol{\ell}}}$
${\boldsymbol{\psi}_{\rm pix}}$ | Pixel values of the pixelated potential
$\lambda_{\rm\psi,bl}$, $\lambda_{\rm\psi,bl}$ |  Regularization strengths for ${\boldsymbol{\psi}_{\rm pix}}$
${\rm\bf W}_{\rm\psi,bl}$, ${\rm\bf W}_{\rm\psi,bl}$ |  Regularization weights for ${\boldsymbol{\psi}_{\rm pix}}$
Optimization & Inference
$L$, $\nabla L$ | Loss function and its gradient
${\boldsymbol{\eta}}_{\rm MAP}$ | MAP solution
${\rm\bf C}_{\boldsymbol{\eta}}$ | Parameter covariance matrix

###  2.1 Discrete formulation of lensing

From the observation of a strongly lensed source, or data ${\boldsymbol{d}}$,
we aim to recover the lens potential $\psi$ and the (unlensed) source light
${\boldsymbol{s}}$. Additionally, the lens light ${\boldsymbol{\ell}}$, which
is often partially blended with the lensed features, must also be modeled,
either jointly with the lens potential and the source, or carefully subtracted
from the data beforehand.

The mapping between an observed angular position $\@vec{\theta}$ on the sky
and a correpsonding position on the source plane $\@vec{\beta}$ before lensing
is given by the lens equation

$\displaystyle\@vec{\beta}=\@vec{\theta}-\@vec{\alpha}(\@vec{\theta})\ ,$ (1)

where the reduced deflection angle $\@vec{\alpha}$ is given by the first
derivatives of the lens potential

$\displaystyle\@vec{\alpha}(\@vec{\theta})=\nabla\psi(\@vec{\theta})\ .$ (2)

It is possible to infer the (dimensionless) surface mass density of the lens
from the lens potential, called the convergence $\kappa$, which can be
computed from the second derivatives of the lens potential, as

$\displaystyle\kappa(\@vec{\theta})\equiv\frac{\Sigma(\@vec{\theta})}{\Sigma_{\rm
crit}}=\frac{1}{2}\nabla^{2}\psi(\@vec{\theta})\ ,$ (3)

where $\Sigma$ is the surface mass density in physical units and $\Sigma_{\rm
crit}$ the critical density that depends on the cosmology and the redshifts of
the lens and source galaxies. For quoting quantities in physical units (e.g.
masses), we assume a flat $\Lambda$CDM model with Hubble constant $H_{0}=70$
and matter density at present time $\Omega_{{\rm m},0}=0.3$.

The lens equation holds for any light ray from the source to the observer, and
as such describes the lensing effect in a continuous way. However, the data is
pixelated and includes additional nuisance effects due to limiting seeing
conditions and instrumental noise. We thus discretize the problem and write
${\boldsymbol{d}}$ as

$\displaystyle{\boldsymbol{d}}=\boldsymbol{\mathsf{B}}\,\boldsymbol{\mathsf{L}}_{{\boldsymbol{\psi}}}{\boldsymbol{s}}+\boldsymbol{\mathsf{B}}\,{\boldsymbol{\ell}}+{\boldsymbol{n}}\
,$ (4)

where ${\boldsymbol{s}}$ and ${\boldsymbol{\ell}}$ are vectors holding the
true (discretized) surface brightness values of the unlensed source and the
lens, respectively. The lensing operator
$\boldsymbol{\mathsf{L}}_{{\boldsymbol{\psi}}}$ encodes a discrete version of
Eq. 1 that models the lensing effect by mapping source surface brightness
values onto the lens plane, based on the lens potential ${\boldsymbol{\psi}}$
(which can also be discretized) through Eq. 2. The blurring operator
$\boldsymbol{\mathsf{B}}$ models the seeing conditions by convolving the
images of the lens and the lensed source with the point spread function (PSF)
of the instrument. Throughout this paper, we assume that the PSF has been
modeled beforehand (e.g. from stars in the field) and that it remains constant
in Eq. 4. The last term, ${\boldsymbol{n}}$, represents additive noise, and is
usually a combination of instrumental read-out noise (gaussian noise) and shot
noise (Poisson noise that depends on ${\boldsymbol{s}}$ and
${\boldsymbol{\ell}}$). We note that while the model depends linearly on
${\boldsymbol{s}}$ and ${\boldsymbol{\ell}}$, the lensing operator
$\boldsymbol{\mathsf{L}}_{{\boldsymbol{\psi}}}$ depends non-linearly on
${\boldsymbol{\psi}}$ through the lens equation (Eq. 1).

###  2.2 Lens modeling and inversion of the lens equation

The model defined in Eq. 4 cannot easily be inverted to retrieve the lens
potential, the unlensed source light, and the lens light. The problem is ill-
posed and subject to degeneracies that cannot be mitigated based only on the
data: (1) some locations in the source plane do not map to any data pixel
(thus are unconstrained), (2) the lens potential is directly constrained only
at locations where there are lensed features, (3) the lens light is often
blended with these same lensed features. Therefore additional constraints,
based on priors on ${\boldsymbol{\psi}}$, ${\boldsymbol{s}}$ and
${\boldsymbol{\ell}}$, are needed to solve Eq. 4 and obtain physically-
motivated solutions. This is particularly important in situations where the
number of model parameters is comparable to the number of pixels in the data,
which is a common situation when modeling complex sources or perturbations in
the lens potential.

Before specifying the choice of priors, we first simplify the notation by
defining the model ${\boldsymbol{m}}$ and the corresponding full set of
parameters ${\boldsymbol{\eta}}$ that describe ${\boldsymbol{\psi}}$,
${\boldsymbol{s}}$ and ${\boldsymbol{\ell}}$ as

$\displaystyle{\boldsymbol{m}}({\boldsymbol{\eta}})\equiv{\boldsymbol{m}}({\boldsymbol{\psi}},{\boldsymbol{s}},{\boldsymbol{\ell}})=\boldsymbol{\mathsf{B}}\,\boldsymbol{\mathsf{L}}_{{\boldsymbol{\psi}}}{\boldsymbol{s}}+\boldsymbol{\mathsf{B}}\,{\boldsymbol{\ell}}\
.$ (5)

Inverting Eq. 4 boils down to obtaining the set of maximum a posteriori (MAP)
parameters ${\boldsymbol{\eta}}_{\rm MAP}$ that maximizes the posterior
probability distribution
$P\left({\boldsymbol{\eta}}\,|\,{\boldsymbol{d}},{\boldsymbol{m}}\right)$ of
the parameters given the data. From Bayes theorem we have

$\displaystyle
P\left({\boldsymbol{\eta}}\,|\,{\boldsymbol{d}},{\boldsymbol{m}}\right)=\frac{P\left({\boldsymbol{d}}\,|\,{\boldsymbol{\eta}},{\boldsymbol{m}}\right)P\left({\boldsymbol{\eta}}\,|\,{\boldsymbol{m}}\right)}{P\left({\boldsymbol{d}}\,|\,{\boldsymbol{m}}\right)}\
,$ (6)

where the first term of the numerator is the data likelihood, the second term
is the prior, and the denominator is the Bayesian evidence. The evidence does
not change ${\boldsymbol{\eta}}_{\rm MAP}$, but is particularly relevant for
comparing different models (through evidence ratios).

In practice, we obtain ${\boldsymbol{\eta}}_{\rm MAP}$ by minimizing a loss
function $L({\boldsymbol{\eta}})$ defined as

$\displaystyle L({\boldsymbol{\eta}})$ $\displaystyle=-\log
P\left({\boldsymbol{d}}\,|\,{\boldsymbol{\eta}},{\boldsymbol{m}}\right)-\log
P\left({\boldsymbol{\eta}}\,|\,{\boldsymbol{m}}\right)$
$\displaystyle\equiv\mathcal{L}\left({\boldsymbol{d}}\,|\,{\boldsymbol{\eta}},{\boldsymbol{m}}\right)+\mathcal{P}\left({\boldsymbol{\eta}}\right)\
,$ (7)

where the data-fidelity term $\mathcal{L}$ is the negative log-likelihood, and
the regularization term $\mathcal{P}$ is the negative log-prior. The MAP
solution can then be obtained as

$\displaystyle{\boldsymbol{\eta}}_{\rm
MAP}=\operatorname*{arg\,min}_{{\boldsymbol{\eta}}}\ L({\boldsymbol{\eta}})\
.$ (8)

The choice of the data-fidelity term is tightly linked to the statistical
properties of the data noise ${\boldsymbol{n}}$, characterized by the
covariance matrix of the data ${\rm\bf C}_{\boldsymbol{d}}$. The noise is
composed of instrumental readout noise and shot noise due to both the flux
from the observed target and the sky background, which we assume follow
Gaussian distributions. The target shot noise is estimated from the modeled
flux, as estimating it from the data itself can introduce biases (Horne 1986);
therefore we formally have ${\rm\bf C}_{\boldsymbol{d}}\equiv{\rm\bf
C}_{\boldsymbol{d}}({\boldsymbol{\eta}})$ that we implicitly assume throughout
the following equations to avoid clutter. Moreover, we assume uncorrelated
noise as is usually true for charge-coupled device images, hence ${\rm\bf
C}_{\boldsymbol{d}}$ is a diagonal matrix.

Under the assumption of Gaussian noise, the data-fidelity term is the
$\chi^{2}$ of the data given the model

$\displaystyle\mathcal{L}\left({\boldsymbol{d}}\,|\,{\boldsymbol{\eta}},{\boldsymbol{m}}\right)$
$\displaystyle=\frac{1}{2}\,\chi^{2}$
$\displaystyle=\frac{1}{2}\,\Big{[}\,{\boldsymbol{d}}-{\boldsymbol{m}}\left({\boldsymbol{\eta}}\right)\Big{]}^{\top}{\rm\bf
C}_{\boldsymbol{d}}^{-1}\
\Big{[}\,{\boldsymbol{d}}-{\boldsymbol{m}}\left({\boldsymbol{\eta}}\right)\Big{]}\
.$ (9)

The priors and corresponding regularization terms depend on the specific
choices of parametrization of model components ${\boldsymbol{\psi}}$,
${\boldsymbol{s}}$ and ${\boldsymbol{\ell}}$. In general, each of these model
components can be described by a set of analytical profiles, pixelated
profiles, or a combination of both.

###  2.3 Analytical and pixelated components

We model smooth mass and light distributions with a set of analytical
profiles. In this case, regularization terms in Eq. 2.2 are not explicitly
defined, but rather directly encoded in the parametrization of the model
${\boldsymbol{m}}$, which we write as

$\displaystyle{\boldsymbol{\psi}}\equiv\tilde{}{\boldsymbol{\psi}}({\boldsymbol{\eta}}_{\psi}),\
{\boldsymbol{s}}\equiv\tilde{}{\boldsymbol{s}}({\boldsymbol{\eta}}_{s}),\
{\boldsymbol{\ell}}\equiv\tilde{}{\boldsymbol{\ell}}({\boldsymbol{\eta}}_{\ell})\
,$ (10)

where ${\boldsymbol{\eta}}_{\psi}$, ${\boldsymbol{\eta}}_{s}$,
${\boldsymbol{\eta}}_{\ell}$ are the set of parameters for the analytical
profiles $\tilde{}{\boldsymbol{\psi}}$, $\tilde{}{\boldsymbol{s}}$ and
$\tilde{}{\boldsymbol{\ell}}$ that describe the lens potential, the source and
the lens light, respectively.

To describe more complex mass and light distributions that cannot be captured
by analytical functions, we rely on pixelated components, where pixel values
represent the lens potential, the lens light, or source light (in source
plane) at each pixel position. We adopt the following notation for those
components

$\displaystyle{\boldsymbol{\psi}}\equiv{\boldsymbol{\psi}_{\rm pix}},\
{\boldsymbol{s}}\equiv{\boldsymbol{s}_{\rm pix}},\
{\boldsymbol{\ell}}\equiv{\boldsymbol{\ell}_{\rm pix}}\ .$ (11)

Such pixelated components typically imply a much larger number of parameters
(i.e. each single pixel value). The model inversion is then highly under-
constrained, and the choice of regularization plays a central role in the
success of the method in recovering the underlying lens potential and light
distributions (see e.g. Vernardos & Koopmans 2022).

###  2.4 Multi-scale regularization in source plane

In Galan et al. (2021), we used a multi-scale regularization technique to
reconstruct the source light distribution in the source plane. This
regularization is based on wavelet transforms, which decompose an image into a
set of wavelet coefficients organized by spatial scale. Each of these scales
is a filtered version of the signal (i.e. same number of pixels) that contains
emphasized features at a given spatial scale, similar to a frequency
decomposition using the Fourier transform (see fig. 1 of Galan et al. 2021).

The regularization term in Eq. 2.2 used to reconstruct ${\boldsymbol{s}_{\rm
pix}}$ is

$\displaystyle\mathcal{P}_{{\boldsymbol{s}}}\left({\boldsymbol{s}_{\rm
pix}}\right)$ $\displaystyle=i_{\geq 0}({\boldsymbol{s}_{\rm
pix}})+\lambda_{s}\,\left\lVert\,{\rm\bf
W}_{s}\circ{\boldsymbol{\Phi}}^{\top}\,{\boldsymbol{s}_{\rm
pix}}\,\right\rVert_{1}\ .$ (12)

The first term in the above equation is a positivity constraint111The
indicator function $i_{\geq 0}(\,\cdot\,)$ is formally equal to 0 if its
argument contains only non-negative values, and $+\infty$ otherwise. that
enforces pixel values to be non-negative. The second term combines the
$\ell_{1}$-norm $\left\lVert\,\,\cdot\,\,\right\rVert_{1}$ with the wavelet
transform operator ${\boldsymbol{\Phi}}^{\top}$. The effect of the
$\ell_{1}$-norm is to impose a sparsity constraint on wavelet coefficients,
which is effectively equivalent to a soft-thresholding222Soft-thresholding is
the proximal operator of the $\ell_{1}$-norm. of the coefficients (Starck et
al. 2015). The threshold level depends directly on the hyper-parameter
$\lambda_{s}$, which is further adapted to each wavelet scale through the
weight matrix ${\rm\bf W}_{s}$ to efficiently regularize features that span
different spatial scales in the source plane (the operation $\circ$ represents
the element-wise product). We note that a similar regularization term
$\mathcal{P}_{{\boldsymbol{\ell}}}\left({\boldsymbol{\ell}_{\rm pix}}\right)$
can also be written for reconstructing the lens light (Joseph et al. 2019).

The success of this regularization strategy thus relies on our ability to
correctly estimate the regularization weights ${\rm\bf W}_{s}$. We follow a
data-driven approach that relies on the noise, in the goal to control the
statistical significance of the reconstructed source light distribution via
$\lambda_{s}$ (Paykari et al. 2014). We achieve this by estimating the
standard deviation of the noise in the source plane for each wavelet scale.
The details of this procedure are given in Joseph et al. (2019). The only
remaining hyper-parameter is the overall strength of the regularization
$\lambda_{s}$, which is scalar and usually set between 3 and 5 to ensure high
enough statistical significance of the reconstructed source light distribution
(e.g. Starck et al. 2007).

###  2.5 Multi-scale regularization in the lens potential

Figure 2: Example noise maps used to properly scale the multi-scale
regularization of our pixelated lens potential model ${\boldsymbol{\psi}_{\rm
pix}}$. Each panel shows the standard deviation for each potential pixel, per
wavelet scale of the starlet transform. The resulting maps are similar for the
Battle-Lemarié wavelet transform. Note how the standard deviation strongly
decreases as a function of the wavelet scale. More details are given Sect.
2.5.  source

We use a pixelated component ${\boldsymbol{\psi}_{\rm pix}}$ to model more
complex lens potential distributions not captured solely with smooth
analytical profiles. The regularization of ${\boldsymbol{\psi}_{\rm pix}}$
presented here is very general, and could be used to model either the full
lens potential, or only perturbations to the smooth potential. In this work,
we focus on the latter case and use ${\boldsymbol{\psi}_{\rm pix}}$ to
reconstruct various types of perturbations to the underlying smooth potential.

Similarly to the pixelated source, we need to define a regularization term for
${\boldsymbol{\psi}_{\rm pix}}$. This regularization needs to be flexible
enough to allow for a variety of different features to be reconstructed,
ranging at least from localized subhalos to populations of subhalos and
multipolar moments. Thus, similarly to ${\boldsymbol{s}_{\rm pix}}$, we apply
on ${\boldsymbol{\psi}_{\rm pix}}$ a multi-scale regularization based on
wavelets and sparsity constraints. However we do not impose a positivity
constraint on the values of ${\boldsymbol{\psi}_{\rm pix}}$, as negative
potential pixels can be physically interpreted as local decrease of the smooth
lens potential $\psi_{\rm m}$. The full regularization term for
${\boldsymbol{\psi}_{\rm pix}}$ is

$\displaystyle\mathcal{P}_{{\boldsymbol{\psi}}}\left({\boldsymbol{\psi}_{\rm
pix}}\right)$ $\displaystyle=\lambda_{\psi,\rm st}\,\left\lVert\,{\rm\bf
W}_{\psi,\rm st}\circ{\boldsymbol{\Phi}}_{\rm
st}^{\top}\,{\boldsymbol{\psi}_{\rm pix}}\,\right\rVert_{1}$
$\displaystyle+\lambda_{\psi,\rm bl}\,\left\lVert\,{\rm\bf W}_{\psi,\rm
bl}\circ{\boldsymbol{\Phi}}_{\rm bl}^{\top}\,{\boldsymbol{\psi}_{\rm
pix}}\,\right\rVert_{1}\ ,$ (13)

where we use two different wavelet transforms: the starlet transform
(${\boldsymbol{\Phi}}_{\rm st}^{\top}$) and the Battle-Lemarié wavelet
transform of order 3 (${\boldsymbol{\Phi}}_{\rm bl}^{\top}$). As in Eq. 12,
the scalars $\lambda_{\psi,\rm s}$ and $\lambda_{\psi,\rm bl}$ are the
regularization strengths, and the elements of the matrices ${\rm\bf
W}_{\psi,\rm s}$ and ${\rm\bf W}_{\psi,\rm bl}$ are the regularization
weights. The starlet transform is the transform used for modeling the source
light distribution in Galan et al. (2021), and is well-suited for the
reconstruction of multiple locally isotropic features at different spatial
scales. The Battle-Lemarié transform is introduced to reduce the appearance of
spurious isolated potential pixels, inspired by a similar use for mass-mapping
studies from weak lensing observations (Lanusse et al. 2016).

Similarly to the regularization of the source, we aim at finding the
regularization weights ${\rm\bf W}_{\psi,\rm st}$ and ${\rm\bf W}_{\psi,\rm
bl}$ based on the noise in the data. The situation is however more complicated
than for the pixelated source reconstruction discussed in Sect. 2.5, because
the relation between the lens potential and the lensed source light
distribution is highly non-linear due to the lens equation. Nevertheless, in
the limit of small perturbations $\delta{\boldsymbol{\psi}}$ to the smooth
potential ${\boldsymbol{\psi}}$, it is possible to linearize the lens equation
and relate changes in the lens potential to changes in the lensed light
distribution of the source. We use this first-order treatment for estimating
the covariance—more specifically the standard deviation—in the lens potential,
and scaling the regularization strengths accordingly.

A first-order Taylor expansion of the lens equation around $\@vec{\beta}$
leads to the following equation (Blandford et al. 2001; Koopmans 2005)

$\displaystyle\delta{\boldsymbol{d}}(\@vec{\theta})\approx-\nabla_{\beta}\,{\boldsymbol{s}}(\@vec{\beta})\cdot\nabla_{\theta}\,\delta{\boldsymbol{\psi}}(\@vec{\theta})\
,$ (14)

where residuals
$\delta{\boldsymbol{d}}\equiv{\boldsymbol{d}}-\tilde{{\boldsymbol{m}}}({\boldsymbol{\psi}},{\boldsymbol{s}})$
are based on a preliminary model $\tilde{{\boldsymbol{m}}}$ that does not
include any perturbations $\delta{\boldsymbol{\psi}}$ in the lens potential.
Gradients are computed with respect to the coordinates indicated by the
subscripts $\beta$ and $\theta$ (defined in Eq. 1). The above equation
provides a linear relation that connects individual pixels in potential space
to individual pixels in data space, which we can use to propagate noise levels
for tuning the regularization strengths. We give all the remaining details of
the computation in Appendix A. The resulting weight matrices ${\rm\bf
W}_{\psi,\rm st}$ and ${\rm\bf W}_{\psi,\rm bl}$ contain the standard
deviation of the noise for each wavelet scale in the lens potential. We show
an example of weights with respect to the starlet transform (i.e. ${\rm\bf
W}_{\psi,\rm st}$) on Fig. 2.

As for the source, the remaining hyper-parameters are the two regularization
strengths $\lambda_{\psi,\rm st}$ and $\lambda_{\psi,\rm bl}$, which directly
control the statistical significance of the starlet and Battle-Lemarié wavelet
coefficients, respectively. These scalars are set in practice set between 3
and 5, depending on how strongly certain features need to be regularized.

###  2.6 Optimization with differentiable programming

So far we have only defined how the different model components
(${\boldsymbol{\psi}}$, ${\boldsymbol{s}}$, ${\boldsymbol{\ell}}$) can be
parametrized and cast into an optimization problem. However, combining
components which are described either with analytical profiles or on pixelated
grids is a challenging task, as the standard methods to minimize the loss
function can be fundamentally different in each case. With analytical model
components, the MAP solution ${\boldsymbol{\eta}}_{\rm MAP}$ can be approached
via stochastic algorithms such as particle swarm optimization (PSO, Kennedy &
Eberhart 2001). With wavelet-regularized pixelated components (Eqs. 12 and
2.5), convergence to ${\boldsymbol{\eta}}_{\rm MAP}$ usually requires the use
of carefully-chosen iterative algorithms relying on the formalism of proximal
operators to apply constraints such as sparsity and positivity to the solution
(so-called proximal splitting algorithms, see e.g. Starck et al. 2015).

In Galan et al. (2021), we used a hybrid scheme that first optimizes the
smooth lens potential described analytically using a PSO, followed by a source
reconstruction step (at fixed lens potential parameters) using an iterative
proximal algorithm. However, this strategy does not scale well with model
complexity. Each additional pixelated model component would require similar
hybrid schemes, which rapidly become inefficient at converging to the MAP
solution. In addition, estimating the joint posterior distribution of model
parameters using traditional sampling techniques (e.g. Markov chain Monte
Carlo, MCMC) becomes computationally too expensive. In this work, we overcome
this issue by implementing a fully differentiable loss function, meaning that
we can obtain its full gradient and higher order derivatives analytically.
This allows us to simultaneously optimize all parameters, both analytical
parameters and individual pixel values, using robust gradient descent
algorithms which remain efficient even in a large parameter space. This
naturally replaces proximal algorithms that are usually necessary to solve for
pixelated components, and leads to a self-consistent combination of both
analytical and pixelated model components.

Gradient descent optimization guarantees convergence to a minimum of the loss
function, which is typically not the case for stochastic optimization
algorithms that do not use the gradient (or higher order derivatives) of the
loss function. However there is still no definitive guarantee to converge to
the global minimum (i.e. ${\boldsymbol{\eta}}_{\rm MAP}$), which can depend on
the parameter initial values. It is possible to address this limitation by
using a multi-start gradient descent optimization, which runs the same
minimization multiple times for different parameter initializations (e.g. Gu
et al. 2022). We leave this optimization improvement to future work.

We use the automatic differentiation Python library JAX (Bradbury et al. 2018)
to construct a fully differentiable model and loss function. Wavelet
transforms are implemented using convolutions, making them straightforwardly
differentiable (as in convolutional neural networks). We also take advantage
of efficient compilation features of JAX to speed up computations during
evaluation of the loss function and its derivatives. All of our modeling
methods and algorithms are implemented in the Python software package
Herculens, which we make publicly
available333https://github.com/austinpeel/herculens. The code structure and
part of the modeling routines of Herculens are based on the open-source
modeling software package Lenstronomy (Birrer & Amara 2018; Birrer et al.
2021).

###  2.7 Estimating the parameter covariance matrix

Estimating parameter uncertainties and covariances is crucial to reliably
interpret the model that corresponds to the MAP solution. Sampling techniques
such as MCMC, which draw samples from the full posterior distribution function
(Eq. 6), are often used for this purpose. However, for large parameter spaces,
such stochastic techniques become inefficient, especially when parameters
depend non-linearly on each other, as is the case in lensing.

In this regime, Hamiltonian Monte Carlo sampling (HMC, introduced by Duane et
al. 1987; Neal 2011, for a review) is particularly efficient as each new
sample is drawn based on the gradient of the loss function, resulting in high
acceptance rates. While individual HMC steps might be more expensive to
compute, the number of samples required for a reliable estimation of the
parameters’ posterior distributions is largely reduced compared to stochastic
techniques like MCMC, which tend to deliver noisier distributions. Moreover,
we note that gradient-informed nested sampling for calculating the Bayesian
evidence also exists (see e.g. Albert 2020a), which is an important tool for
model comparison (although not explored in this work).

In addition to HMC sampling, we also explore the possibility of using the
Fisher information matrix (FIM) and its relation to the second-order
derivatives of the loss function to estimate the parameters’ covariance
matrix. This is usually referred to as a Fisher analysis, which has been used
in various studies including the forecast of constraints on cosmological
parameters (e.g. Philcox et al. 2021), the mapping of large-scale structures
(e.g. Abramo 2012), the analysis of gravitational waves (e.g. Belgacem et al.
2019), and the substructure power-spectrum from lensing (Cyr-Racine et al.
2019).

Second-order partial derivatives of the loss function $L$ define the Hessian
matrix, which reflects the local shape of the loss function at any point of
the parameter space. Each entry of the Hessian matrix ${\rm\bf H}_{L}$ is
defined as

$\displaystyle{\rm\bf
H}_{L}({\boldsymbol{\eta}})=\left\\{\frac{\partial^{2}L({\boldsymbol{\eta}})}{\partial\eta_{i}\,\partial\eta_{j}}\right\\}\
,$ (15)

where the indices $i$ and $j$ indicate different model parameters from the
entire parameter set ${\boldsymbol{\eta}}$. The expected Fisher information
matrix, denoted as $\mathcal{I}$, is defined as

$\displaystyle\mathcal{I}({\boldsymbol{\eta}})=\mathbb{E}\Big{[}{\rm\bf
H}_{L}({\boldsymbol{\eta}})\Big{]}\ ,$ (16)

where $\mathbb{E}$ is the expectation operator over the distribution of all
the realizations of the data for a fixed set of parameters
${\boldsymbol{\eta}}$ (i.e.
$P\left({\boldsymbol{d}}\,|\,{\boldsymbol{\eta}}\right)$). In practice,
however, we cannot compute the expected value
$\mathcal{I}({\boldsymbol{\eta}})$, because we only have access to a single
observation ${\boldsymbol{d}}$, with a specific realization of the noise.
Instead, we can compute the observed FIM, ${\rm\bf I}({\boldsymbol{\eta}})$,
evaluated at the MAP solution

$\displaystyle{\rm\bf I}({\boldsymbol{\eta}}_{\rm MAP})={\rm\bf
H}_{L}\left({\boldsymbol{\eta}}_{\rm MAP}\right)\ .$ (17)

The FIM can then be used to approximate the covariance matrix of the MAP
solution through matrix inversion

$\displaystyle{\rm\bf C}_{\boldsymbol{\eta}}\approx{\rm\bf
I}({\boldsymbol{\eta}}_{\rm MAP})^{-1}\ .$ (18)

This covariance directly gives a lower bound on the uncertainty for each
parameter444This is more formally called the Cramér-Rao bound, which states
that the variance of an unbiased parameter estimator is at least as high as
the inverse of the Fisher information (e.g. Cramér 1999).. We note that this
(first-order) approximation of ${\rm\bf C}_{\boldsymbol{\eta}}$ becomes exact
if the loss function locally behaves quadratically. We expect this to be the
case for simple, smooth models described with analytical profiles, as shown in
Vernardos & Koopmans (2022).

We find that while for fully analytical models we can rely on the above first-
order approximation of parameter covariance matrix, a sampling-based
exploration of the highly non-linear parameter space using HMC is warranted.

##  3 Experimental setup

Figure 3: Overview of the simulated HST observations. The left-most column
shows the simulated data without perturbations to the smooth lens potential.
The dashed lines in the bottom-left panel are isopotential contours. The
remaining columns show the mock data used in this work (top row), the three
different input potential perturbations (middle row), and the difference
between the unperturbed and perturbed lenses (bottom row). The range (min and
max) of the perturbations varies according to the different nature of the
perturbers. The solid lines enclose the region where the signal-to-noise of
the lensed source is higher than 5. The isopotential contours are almost
indistinguishable between the unperturbed and perturbed cases, therefore we
omit them.  source

We test our method on mock Hubble Space Telescope (HST) observations that
include realizations of the three perturbation categories introduced in Sect.
1: a localized DM subhalo, a population of DM subhalos, and higher-order
multipoles. We simulate these perturbing fields based on physical
considerations and results of previous works (detailed in the sections below),
which result in each case in very different levels of perturbations when
compared to the underlying smooth lens potential.

We focus on the strong lensing of a source with smooth surface brightness,
resembling early-type/early-type lenses (EELs, Oldham et al. 2017). We
simulate both the source and the lens light as single elliptical Sérsic
profiles, as early-type galaxies are often well-fitted by these profiles (e.g.
Shu et al. 2017). The smooth lens potential is described by a singular
isothermal ellipsoid (SIE) for the main deflector, embedded in an external
shear to simulate the net influence of neighboring galaxies. All analytical
profiles used in this work are defined in Appendix B.

To create the mock data, we use an algorithm that is entirely independent of
our modeling code, in order to prevent the occurrence of artificially
advantageous minima during optimization. This also has the advantage to
closely mimic a real-world situation, as the data never exactly corresponds to
any model generated by the modeling code itself. We use the software package
MOLET (Vernardos 2021) to simulate typical observations of strongly lensed
galaxies as observed with HST and the Wide Field Camera 3 (WFC3) instrument,
in the near infrared (F160W filter). The pixel size is
$0\aas@@fstack{\prime\prime}08$, and the field of view is
$8^{\prime\prime}\times 8^{\prime\prime}$. In this work we do not explore
effects due to incorrect PSF modeling, hence for simplicity we use a gaussian
PSF with $0\aas@@fstack{\prime\prime}3$ FWHM. This results in a simulated data
set that is sufficiently realistic to evaluate our method. The first column of
Fig. 3 shows the simulated image without perturbations to the smooth
potential. The remaining three columns show the different perturbation cases
to which we apply our method. Full details of all instrumental settings and
input model parameters are listed in Appendix C.

###  3.1 Localized subhalo (LS)

We simulate a single localized DM subhalo, with a spherical isothermal profile
(SIS) and Einstein radius $\theta_{\rm E,halo}=0\aas@@fstack{\prime\prime}07$
as input perturbing potential (for reference, the Einstein radius of the main
deflector is $1\aas@@fstack{\prime\prime}6$). For simplicity, we assume that
the subhalo is in the same redshift plane as the main lens, although it can be
located at any redshift (through a simple scaling of its mass). We assume
typical lens and source redshifts of EEL lenses $z_{\rm d}=0.3$ and $z_{\rm
s}=0.7$ respectively (Oldham & Auger 2018), from which we get a mass of
$10^{9}\ {\rm M}_{\odot}$ within the SIS Einstein radius, which is comparable
to previous dark halo detections based on HST observations (Vegetti et al.
2010). The resulting simulated observation is shown in the second column of
Fig. 3. The mean perturbation level, computed within the region containing the
lensed arcs, is $6.4\%$ relative to the smooth potential.

###  3.2 Population of subhalos (PS)

We follow recent work and simulate the net effect of a population of DM
subhalos along the line-of-sight with a gaussian random field (GRF)
(Chatterjee & Koopmans 2018; Bayer et al. 2018; Vernardos et al. 2020). GRF
perturbations are random fluctuations that have a Fourier power-spectrum
following a power-law. We parametrize this power-law relation as

$\displaystyle{\rm PS}\big{(}\psi_{\rm
GRF}\big{)}(k)=C_{\beta}\,\sigma^{2}_{\rm GRF}\,k^{-\beta_{\rm GRF}}\ ,$ (19)

where $k$ is the wave-number, $\beta_{\rm GRF}$ is the power-law slope, and
$C_{\beta}$ is a normalizing factor that depends on $\beta_{\rm GRF}$ and the
size of the field of view, which ensures that $\sigma^{2}_{\rm GRF}$ is the
variance of the GRF (for the exact formula, see Chatterjee 2019). The power-
spectrum is then converted to a specific random realization of direct-space
perturbations using the inverse Fourier transform. The value of $\beta_{\rm
GRF}$ determines the distribution of power at each length scale: a large value
leads to extended and smooth variations, whereas a small value creates a large
number of localized and grainy structures.

Typical ranges for $\sigma^{2}_{\rm GRF}$ and $\beta_{\rm GRF}$ have been
explored in the literature and are justified in Chatterjee (2019) and
Vernardos et al. (2020). In these works, the authors explored ranges
$\sigma^{2}_{\rm GRF}\in[10^{-5},10^{-2}]$ and $\beta_{\rm GRF}\in[3,8]$.
Additionally, Bayer et al. (2018) excluded GRF variance larger than
$\sigma^{2}_{\rm GRF}=10^{-2.5}$, based on HST observations of the strong lens
system SDSS J0252$+$0039\. Hence, we set $\sigma^{2}_{\rm GRF}=10^{-3}$ and
$\beta_{\rm GRF}=4$, such that it leads to a GRF that is not unphysically
large, and contains both small and large-scale features. The specific GRF
realization and the corresponding simulated observation are shown in the third
column of Fig. 3. The corresponding relative mean perturbation level is
$0.83\%$.

###  3.3 Higher-order multipoles (HM)

We introduce higher-order deviations to the smooth potential as a multipole of
order 4 (octupole). We use the same definition of multipole as in Van de
Vyvere et al. (2021)

$\displaystyle\psi_{\rm
multipole}(\@vec{\theta})=\frac{r}{1-m^{2}}\,a_{m}\cos\\!\left(m\phi-m\phi_{m}\right)\
,$ (20)

where $r$ and $\phi$ are polar coordinates transformed from $\@vec{\theta}$.
In our case, we fix the multipole order to $m=4$ (octupole), and the
orientation $\phi_{4}$ to be aligned with the main axis of the SIE component
of the smooth potential (in this case the full potential becomes more disky).
The octupole strength is set to $a_{4}=0.06$, which corresponds to the high
end of the distribution found by Hao et al. (2006), from isophote measurements
over a large sample of elliptical and lenticular galaxies in SDSS data. The
resulting simulated observation is shown in the last column of Fig. 3. The
corresponding relative mean perturbation level is $0.16\%$.

##  4 Modeling the full lens potential

###  4.1 Baseline model, parameter optimization and sampling

We model the simulated data set presented above, with the goal of retrieving
the full lens potential, including the perturbations. The lens potential is
modeled as
${\boldsymbol{\psi}}=\tilde{}{\boldsymbol{\psi}}({\boldsymbol{\eta}}_{\psi})+{\boldsymbol{\psi}_{\rm
pix}}$. The smooth component $\tilde{}{\boldsymbol{\psi}}$ is parameterized as
a SIE and external shear, and the a priori unknown perturbations are captured
in the pixelated component ${\boldsymbol{\psi}_{\rm pix}}$, regularized with
sparsity and wavelets as discussed in Sect. 2.5. The deflection angles at each
position in the image plane are computed based on bicubic interpolation of
${\boldsymbol{\psi}_{\rm pix}}$. We use bicubic instead of bilinear
interpolation in order to compute the the surface mass density (Eq. 3)
corresponding to the pixelated model, as it requires second-order spatial
derivatives of the potential (these derivatives are always zero with bilinear
interpolation). We model its surface brightness using a Sérsic profile, that
is ${\boldsymbol{s}}=\tilde{}{\boldsymbol{s}}({\boldsymbol{\eta}}_{s})$.

The lens light profile is modeled only once with a Sérsic profile for one of
the system, assuming the other model components are known, then it is then
fixed during the rest of the modeling. Fixing the lens light is not identical
to subtracting it from the observation, as it still contributes to the noise
model and reduces the contrast of the lensed features. We assume the lens
light centroid is a good tracer of the mass centroid, which is a realistic
scenario for fairly isolated lens galaxies (see e.g. Shajib et al. 2019).
Therefore, we join the center of the SIE profile to that of the lens light
profile. We note that it is however outside the scope of this work to assess
the impact of inaccurate lens light modeling.

We model instrumental and seeing effects by assuming perfect knowledge of the
PSF, background noise level (read-out noise and shot noise from sky
brightness), however we estimate the shot noise from the modeled lens and
source light distributions to estimate the covariance matrix ${\rm\bf
C}_{\boldsymbol{d}}$.

When optimizing only analytical profile parameters
(${\boldsymbol{\eta}}_{\psi}$, ${\boldsymbol{\eta}}_{s}$)—which we do before
including pixelated components in the model—, we find that the quasi-Newton
optimization method
BFGS555https://docs.scipy.org/doc/scipy/tutorial/optimize.html (Nocedal &
Wright 2006) is sufficient to reach convergence to the MAP solution. However,
the optimization of both analytical and pixelated model components is more
challenging, as the dynamic range of analytical and pixel parameters can vary
significantly during optimization according to their impact on the loss
function and the different regularization terms. Therefore, in this case,
parameter updates are performed using the adaptive gradient descent algorithm
AdaBelief (Zhuang et al. 2020), which is extremely efficient for optimizing a
large number parameters (typically as large as for convolutional neural
networks). The initial learning rate is set in order to obtain a smooth
decrease of the loss function until convergence, coupled with an exponential
decay of the learning rate. We use the optimization library Optax (Hessel et
al. 2020) that implements the AdaBelief algorithm.

For fully analytical models, we find that using the fast estimation from the
FIM leads to a parameter covariance matrix almost indistinguishable from the
one obtained via sampling methods. Therefore we only rely on the FIM for these
simple models. However, for the more complex models that include a pixelated
component, we find that HMC sampling of the parameter space is warranted to
obtain reliable estimates of the posterior distributions. We use the python
package BlackJAX666https://github.com/blackjax-devs/blackjax that provides an
implementation of HMC well-integrated with JAX, that we run using the No
U-Turn Sampler algorithm (NUTS, Hoffman & Gelman 2011) to dynamically adapt
the step size, and their “Stan’s adaptation window” feature to improve
sampling efficiency.

Our fiducial model is defined with a ${\boldsymbol{\psi}_{\rm pix}}$ pixel
size set to three times the data pixel size, leading to a resolution of
$0\aas@@fstack{\prime\prime}24$. The choice of pixel size is based on
preliminary models comparing the best-fit reconstructions and model residuals
for different ${\boldsymbol{\psi}_{\rm pix}}$ resolutions. While the residuals
do not vary significantly for a pixel scale between 4 and 1.5, best-fit
reconstructions obtained with pixel scales larger than 3 display artifacts at
the scale of individual pixels. Although most of these artifacts can be
efficiently reduced using our multi-scale regularization strategy by
increasing further the regularization strength for small wavelet scales, we
find that using a larger number of model parameters for only marginal
improvements in terms of residuals is not necessary for this work. We refer to
Appendix D for further discussion on the choice of pixel scale for
${\boldsymbol{\psi}_{\rm pix}}$. The total number of parameters for our
fiducial model is thus $1101$ ($1089$ for ${\boldsymbol{\psi}_{\rm pix}}$, $5$
for ${\boldsymbol{\eta}}_{\psi}$ and $7$ for ${\boldsymbol{\eta}}_{s}$),
jointly optimized and constrained by $10^{4}$ data pixels.

The regularization strengths for ${\boldsymbol{\psi}_{\rm pix}}$ are set to
$\lambda_{\psi,\rm st}=3$ and $\lambda_{\psi,\rm bl}=4$, that is in the range
of values discussed in Sect. 2.5. We use a $1\sigma$-higher strength for the
Battle-Lemarié regularization in order to penalize more the appearance of
spurious pixels in the solution.

###  4.2 Modeling the perturbations only

Figure 4: Best-fit pixelated potential reconstruction of our simulated data
set, assuming all other model components are perfectly known. Each column
represents the different types of perturbations considered in this work. From
top to bottom: simulated data, image model, normalized residuals, input
perturbations, ${\boldsymbol{\psi}_{\rm pix}}$ model. The outlined annular
region on the potential panels corresponds to the solid lines in the panels of
Fig. 3, and is where the signal-to-noise of the lensed source is higher than
5. Perturbations outside this area are cropped only to ease the visual
comparison with the input perturbations (there is no cropping or masking
during modeling). Additionally, for lens LS only, the minimum value of the
perturbations is subtracted from the model, again to ease the comparison with
the input SIS profile (this does not affect the lensing observables). An
animated visualization of the gradient descent optimization is available on
this link.  source

We first reconstruct the perturbations in the idealized case where all smooth
components are perfectly known and fixed. This unrealistic scenario allows us
to assess the best level of perturbations that can be recovered, given the
quality of the data set. As the source and the smooth potential are fixed, the
regularization weights ${\rm\bf W}_{\psi,\rm st}$ and ${\rm\bf W}_{\psi,\rm
bl}$ can be pre-computed and kept constant throughout the gradient descent
optimization.

The resulting models, corresponding to the MAP parameters
${\boldsymbol{\eta}}_{\rm MAP}$, are shown in Fig. 4. We also note that the
pixelated model ${\boldsymbol{\psi}_{\rm pix}}$ is expected to differ from the
input by a uniform offset, as a uniform value in the potential corresponds to
zero deflection of light rays (Eq. 2). This constant offset is thus never
constrained by the data alone (and depends on the initialization). Hence, for
better visualizing the reconstruction for lens LS (for which we know that the
input potential is positive), we shift the ${\boldsymbol{\psi}_{\rm pix}}$
model such that the minimum pixel value displayed on the figure is zero. We do
not add such an offset to our reconstructions for lenses PS and HM shown on
Fig. 4, as the underlying perturbations have roughly zero-mean.

Overall the characteristic features of each type of perturbations are well
recovered. For lens LS, both the subhalo position and the shape of the
underlying SIS profile are captured by the model. For lens PS, the clear over-
density region on the bottom right part of the arc is recovered, although the
correspondence with the input perturbations is less than for lens LS. For lens
HM, the reconstruction displays imprints of azimuthal periodicity between
over- and under-density regions, despite the low level (0.2%) of perturbation.
Lastly, we notice that the amplitude of the modeled perturbations is
systematically lower, by a factor of $\sim$2, than the input perturbations. We
note that a similar result can be seen in some of the pixelated
reconstructions of Vernardos & Koopmans (2022, see their fig. 8).

###  4.3 Modeling the full potential and the source

Figure 5: Best-fit models on the simulated data set, after full modeling of
the lens potential and the source. The first four rows are as in Fig. 4. The
last two rows show, in order: the standard deviation for each
${\boldsymbol{\psi}_{\rm pix}}$ pixel obtained from HMC sampling of the
posterior, and the SNR maps that correspond to the absolute value of the
modeled potential divided by the standard deviation.  source

Our method is then applied to the more realistic situation in which the full
lens potential and the source light are unknown as well. We do so by uniformly
modeling the simulated data set using a series of steps involving gradient
descent optimization. After each step, the model degrees of freedom are
gradually increased, in order to prevent the optimizer to be trapped in local
minima. These optimization steps are as follows:

  1. 1.

The pixelated component ${\boldsymbol{\psi}_{\rm pix}}$ of the lens potential
is initialized to zero (in practice to $10^{-8}$, to prevent gradients to be
evaluated at zero, see also Sect. E) and kept fixed. The initial model is thus
fully smooth, with corresponding parameters
${\boldsymbol{\eta}}_{\boldsymbol{\psi}}$ and
${\boldsymbol{\eta}}_{\boldsymbol{s}}$. We use a multi-start gradient descent
(as advocated by Gu et al. 2022) with 30 runs, for which we verified that it
leads to convergence to the global minimum.

  2. 2.

The pixelated potential component is released and regularization strengths are
deliberately set to large values, such that only the most significant
potential pixels enter the solution. The regularization weights are computed
based on the smooth model from the previous steps, and fixed throughout the
gradient descent. This prevents the pixelated component from fitting all model
residuals from the previous step, which can strongly bias the recovered
perturbations. We find that setting $\lambda_{{\boldsymbol{\psi}},\rm st}=10$
and $\lambda_{{\boldsymbol{\psi}},\rm st}=20$ leads to an intermediate
solution that contains the main features of the perturbations, reducing
slightly model residuals but preventing the model from getting trapped in a
spurious minimum. All model parameters
(${\boldsymbol{\eta}}_{\boldsymbol{\psi}}$,
${\boldsymbol{\eta}}_{\boldsymbol{s}}$ and ${\boldsymbol{\psi}_{\rm pix}}$)
are simultaneously optimized.

  3. 3.

The regularization strengths are set to their fiducial values
($\lambda_{{\boldsymbol{\psi}},\rm st}=3$, $\lambda_{{\boldsymbol{\psi}},\rm
st}=4$), and regularization weights are re-computed based on the previous
models. All model parameters are simultaneously further optimized.

  4. 4.

HMC sampling is performed to estimate the error and posterior distributions
for all parameters. We find that, thanks to the gradient information, only 400
samples (after a warm-up phase of 100 samples) are necessary to obtain well-
sampled posteriors. For posterior distributions that are have multiple modes,
which is not the case in our models, the number samples should be increased.

In Fig. 5 we summarize, for each lens, the models corresponding to the best-
fit parameters ${\boldsymbol{\eta}}_{\rm MAP}$. The panels are identical to
those of the ideal case (Fig. 4), although we added two additional panels to
help the interpretation of the results. The first of these panels is the
standard deviation for each ${\boldsymbol{\psi}_{\rm pix}}$ pixel, based on
the HMC samples, which can be interpreted as an error map for the best-fit
model. These error maps are almost identical for each modeled system, but are
necessary to compute the second additional panel, which shows a measure of the
SNR of the reconstruction. We define the SNR for ${\boldsymbol{\psi}_{\rm
pix}}$ as the absolute value of the best-fit model divided pixel-wise by the
standard deviation. These maps reveal the regions where the reconstruction is
the most statistically significant.

For lens LS, the position of the subhalo is clearly and accurately recovered.
However, the overall shape of the underlying smooth subhalo profile is not
recovered equally well as in the ideal case. We also observe some features on
the other side of the Einstein ring, although those are less significant than
at the subhalo position. The SNR map clearly shows that the feature at the
position of the subhalo has a higher likelihood compared to other features
present in the solution.

For lens PS, we see a relatively good agreement with the input field, where
both over- and under-dense regions are recovered. However, again, the
amplitude of the perturbation is well below the input value, by a factor
between 1.5 and 3. Contrary to the reconstruction of the ideal case, the over-
dense region on the lower part of the ring is now the most prominent. This
feature is well aligned with a lower-intensity region of the arc, which can
partly explain why the model favors a correction to the smooth potential at
this location. The SNR map confirms the significance of this feature. All the
reconstructed features have a SNR of approximately 3 and above. This system is
arguably the most complex to model, as the perturbations affect the lensed
features at many different locations with different strengths and
orientations, which is likely to translate to a larger set of possible
solutions.

For lens HM, the reconstruction reveals a clear azimuthal periodicity, with no
regions significantly different from the input perturbations. Interestingly,
the reconstruction is very similar to the ideal case, despite the input level
of perturbations (relative to the smooth potential) being lower than in the
other systems. Similar to lens PS, we notice a region where the model over-
predicts a correction to the smooth potential. This region is better revealed
in the SNR map and corresponds to the same low-flux region of the arc. This is
not surprising as lower SNR in the data leads to looser constraints on the
model parameters. For this lens we also notice that the model does not predict
any perturbation at the position of left-most image of the source. At this
location, which is closer to the lens center compared to other parts of the
arc, the input perturbation level is lower while the data noise is higher
(shot noise), which leads to almost no detection of perturbations.

Figure 6: MAP values and estimated uncertainties for the smooth potential
parameters ${\boldsymbol{\eta}}_{\psi}$ (from left to right: Einstein radius,
complex ellipticities of the SIE, external shear components). The dashed lines
indicate the input values. The ‘smooth’ models correspond to models without
including pixelated perturbations; ‘pixelated’ models correspond to our
baseline models shown on Fig. 5; ‘high-res’ (lens PS) is similar to the
fiducial model but with a smaller ${\boldsymbol{\psi}_{\rm pix}}$ pixel size;
‘refined’ corresponds to models where the pixelated perturbations have been
replaced with an analytical profile (see text for more details). For several
parameters, error bars are smaller than the marker size.  source

###  4.4 Smooth potential parameters

In the previous section, we have seen that the main features of the
perturbations are overall well recovered, but the reconstructions are not
perfect despite model residuals almost at the noise level. Therefore, we
expect that some of the inaccuracies in the pixelated potential model are
absorbed by the smooth analytical component of the lens potential, or vice
versa.

We check this hypothesis by comparing different models of the smooth
potential, with and without including potential perturbations. As discussed
earlier in this section, parameter uncertainties are either computed using the
fast approximation of the FIM (fully smooth models), or based on HMC sampling
of the parameter space (models including pixelated perturbations). The MAP
values and uncertainties of smooth potential parameters are shown on Fig. 6,
and compared against the input values. We discuss here only lens potential
parameters, but the conclusions are similar for the source parameters as well,
that we present in Appendix F.

We observe that fully smooth models display strong biases in almost all
parameters as expected. Interestingly, models including a pixelated model in
the potential, while having larger error bars, still lead to statistically
significant biases. One particularly informative parameter is the inferred
Einstein radius, as the value of $\theta_{\rm E}$ is slightly more accurate
after including the potential perturbations, but a substantial bias still
remains. We attribute those biases to a manifestation of the degeneracies that
exist between the smooth component of the lens potential and the pixelated
perturbations, where some adjustments of one component can compensate for the
other, still leading to comparable model residuals (see Vegetti et al. 2014,
who also observed re-absorption of residuals in the macro potential model).

###  4.5 Analytically refined models

To assess if the biases discussed above can be mitigated, we reduce the
parameter space by replacing the pixelated component with analytical profiles.
In a real-world scenario, where the underlying type of perturbations is
unknown, this would correspond to imposing stronger priors on the model,
motivated by the characteristic features observed in the pixelated model. This
strategy is similar to previous studies based on the gravitational imaging
technique (e.g. Vegetti & Koopmans 2009; Vegetti et al. 2012), where an
analytical profile is optimized at the position of a tentative detection of a
subhalo, in order to better characterize its properties (position, radial
profile, mass).

Among the different systems modeled in this work, the most characteristic
features that can be noticed in the pixelated models are for the LS and HM
cases, with either a very localized decrease of the potential, or azimuthal
periodicity centered on the lens galaxy, respectively. Therefore, for these
two lenses, we start from the MAP solution obtained with our fiducial model
and include a SIS profile (lens PS) or an multipole component (lens HM)
instead of the pixelated component. We refer to these additional fully
analytical models as ‘refined’ in the remaining of the text.

The smooth potential parameters inferred from these refined models are
compared to the smooth models on Fig. 6. We see that the significant biases
observed with models with too few or too many degrees of freedom have been
correctly mitigated. This result is not surprising as the refined model is now
parametrized identically to the simulated data. Nevertheless, this allows us
to confirm that even after optimizing a more complex and possibly inaccurate
model, it does not prevent us from accessing the optimal global solution after
correctly identifying the underlying type of perturbations. In Sect. 5 we use
these refined models to retrieve the properties of the underlying
perturbations.

###  4.6 Higher resolution pixelated model

Contrary to the localized subhalo and multipolar structures, the underlying
perturbations for lens PS are described by a GRF, which does not have a
specific analytical profile as it is a random realization. Instead, we use
this system to test if using a higher resolution grid for the pixelated
component ${\boldsymbol{\psi}_{\rm pix}}$ (i.e. more parameters) allows us to
reduce the biases on smooth model parameters. The resulting model, named
“high-res”, is compared to the fully smooth and fiducial models on Fig. 6.
Interestingly, we observe as a global trend that smooth parameters tend to be
closer to the input values. While this model is insufficient to recover
unbiased smooth potential parameters, we will see in Sect. 5 that the
recovered power-spectrum of the perturbations is closer to the input one.

Nevertheless, we note that recovering the smooth potential parameters in the
presence of perturbations such as a subhalo population modeled via a GRF is
achievable by imposing informed priors on the pixelated potential model. This
has been shown in the recent work of Vernardos & Koopmans (2022), which
extended the gravitational imaging technique using a covariance-based
regularization of the pixelated potential model. The covariance matrix
governing the regularization term can be specifically adapted to GRF-like
perturbations, leading to an effective regularization of the solution if the
assumption matches the underlying perturbations. We plan to implement the
strategy presented in Vernardos & Koopmans (2022) in the Herculens package,
and leave for future works its comparison with the method presented here.

##  5 Constraints on the underlying perturbations

In the previous section we discussed in a qualitative manner the reconstructed
perturbations. Here we seek to quantify the properties that can be recovered
from these models, and discuss the robustness of those measurements as well as
their applicability to real data sets. All the inferred quantities discussed
in the next subsections are summarized in Table 2.

Table 2: Recovered properties of potential perturbations, from the MAP ${\boldsymbol{\psi}_{\rm pix}}$ models shown on Fig. 5. The quoted values for the subhalo mass is given in solar units and computed within the region with SNR(${\boldsymbol{\psi}_{\rm pix}}$) ¿ 3 (see Fig. 7 and text for more details). The input GRF parameters are given with uncertainties as they were fitted on the power-spectrum of the input perturbing field. Lens | Parameters | Input values | Model | Measured values
---|---|---|---|---
LS | $\big{(}\theta_{x,\rm sub},\theta_{y,\rm sub}\big{)}$ | $\big{(}1.90,-0.40\big{)}$ | pixelated | $\big{(}1.94\pm 0.12,-0.49\pm 0.12\big{)}$
refined | $\big{(}1.899\pm 0.008,-0.408\pm 0.005\big{)}$
$\log_{10}M_{\rm sub}$ | $8.91$ | pixelated | $8.71\pm 0.05$
refined | $8.90\pm 0.01$
PS | $\big{(}\log_{10}\sigma^{2}_{\rm GRF},\beta_{\rm GRF}\big{)}$ | $\big{(}-3.11\pm 0.15,3.23\pm 0.05\big{)}$ | pixelated, ideal | $\big{(}-3.35\pm 0.20,3.40\pm 0.15\big{)}$
pixelated, fiducial | $\big{(}-3.89\pm 0.25,4.15\pm 0.19\big{)}$
pixelated, high-res | $\big{(}-3.49\pm 0.15,3.26\pm 0.09\big{)}$
HM | $\big{(}a_{m},\phi_{m}\big{)}$ | $\big{(}0.06,82.50\big{)}$ | pixelated | $\big{(}0.024\pm 0.002,83.78\pm 0.88\big{)}$
refined | $\big{(}0.061\pm 0.002,83.15\pm 0.39\big{)}$
$m$ | $4$ | refined | $4.01\pm 0.02$

###  5.1 Subhalo mass and position

We consider the two models of lens LS to quantify the properties of the
underlying DM subhalo: our fiducial model including a pixelated component in
the lens potential, and the refined model assuming the detected subhalo mass
distribution follows a SIS profile.

For the pixelated model we assign the position of the detected subhalo to the
minimum of the pixelated potential ${\boldsymbol{\psi}_{\rm pix}}$. We show
the location of the pixel in the top left panel of Fig. 7. For related
uncertainties, we compute the minimum of each HMC sample of
${\boldsymbol{\psi}_{\rm pix}}$, but we note that the minimum remains in the
same potential pixel, leading to error bars smaller than its size. We thus
turn to a more conservative estimate of the uncertainty and simply set it to
half the pixel size (i.e. $0\aas@@fstack{\prime\prime}12$). For the refined
model, we take the optimized position of the SIS as the position of the
subhalo, with uncertainties estimated from the FIM. The resulting positions
and error bars are listed in the top row of Table 2.

Figure 7: Characterization of the position and mass of the localized subhalo
based on our model of lens LS. Top left: ${\boldsymbol{\psi}_{\rm pix}}$
model, from which we assign the detected subhalo position to the minimum of
the clear negative feature, in good agreement with the input subhalo position.
Top right: pixelated convergence model ${\boldsymbol{\kappa}_{\rm pix}}$
obtained from the ${\boldsymbol{\psi}_{\rm pix}}$ model. We indicate the
region within which we compute the subhalo mass, defined as the pixels with
SNR(${\boldsymbol{\psi}_{\rm pix}}$) ¿ 3 (see also the bottom panel of Fig.
5). Bottom: posterior distributions of the subhalo mass, as estimated from the
pixelated model, and the refined model that replaces ${\boldsymbol{\psi}_{\rm
pix}}$ with a SIS profile for the subhalo. The dotted lines indicate the
median value and the $1\sigma$ confidence region of the distributuon. The
black dashed line indicates the subhalo mass as computed from the input
convergence.

The mass of the subhalo is more difficult to estimate from our pixelated
model. Nevertheless, achieving this would be powerful, as it does not require
the choice of a specific shape for the subhalo mass distribution. We start by
computing the surface mass density (i.e. the convergence) corresponding to
each potential pixel using Eq. 3. This results in the pixelated convergence
model ${\boldsymbol{\kappa}_{\rm pix}}$ shown in the top right panel of Fig.
7. Next, we need to define a region in which to integrate
${\boldsymbol{\kappa}_{\rm pix}}$ pixels before converting the surface mass
density to proper solar mass units using $\Sigma_{\rm crit}$. As discussed in
Sect. 3.1, input parameters of the subhalo correspond to a subhalo mass of
about $10^{9}$ within the Einstein radius of $0\aas@@fstack{\prime\prime}07$.
As this scale is much smaller than a pixel of our ${\boldsymbol{\psi}_{\rm
pix}}$ model, we cannot rely on summing the convergence pixels.

We address this issue by considering a larger region within which the subhalo
mass can be inferred for both the analytical and pixelated models. We select a
high-significance region of the reconstruction that contains all pixels with
SNR(${\boldsymbol{\psi}_{\rm pix}}$) ¿ 3 (see bottom left panel of Fig. 5).
This region contains 15 convergence pixels, that we sum and convert to proper
units to estimate the subhalo mass. We repeat the same procedure for each
sample of the joint posterior distribution and find the distribution shown in
the bottom panel of Fig. 7. To compare the inferred value with that of the
input subhalo, we disctretize the input SIS profile by evaluating it on the
same grid of pixels, and compute the mass as for the pixelated model. We note
that the resulting “input” mass is lower than the one computed analytically
within Einstein radius, because of the discretization of the SIS profile that
diverges in the center.

We find that the inferred mean value of the subhalo is lower than the measured
mass on the input perturbations. In addition, we find that the amplitude of
the disagreement is significantly affected by the choice of the region in
which we integrate the convergence. For instance, considering only
${\boldsymbol{\kappa}_{\rm pix}}$ pixels with SNR(${\boldsymbol{\psi}_{\rm
pix}}$) ¿ 5 instead of 3 (6 pixels instead of 15) leads in fact to a very nice
agreement with the input value. The main reason is that this smaller region
essentially excludes the few pixels with negative convergence, seen in brown
color on Fig. 7 (negative convergence pixels do have a physical meaning, as
they indicate a local decrease of the lens mass relative to the smooth
component). On one hand, this leads to a higher mass inferred from the model;
on the other hand, the input value used as a reference is smaller, because it
is computed within a smaller region. These two effects combined lead to an
overall better agreement between the model and the input.

Measuring the mass directly from the pixelated model is therefore challenging
and possibly depends on the set of assumptions. Currently, a more robust
approach is to infer the subhalo mass from our refined model, which is based
on an analytical profile for the subhalo. After applying the same procedure,
we obtain the resulting posterior distribution shown in purple in Fig. 7,
which is in perfect agreement with the input. Again, this is not surprising,
as the SIS profile reflects well the underlying shape of the subhalo.
Nevertheless, these results showcase the requirement of stronger priors on the
shape of the subhalo in order to infer unbiased properties. Additionally, for
proper inference on real data, several works have demonstrated the need to
carefully compare different choices of subhalo profiles (e.g. Çağan Şengül et
al. 2021; Despali et al. 2022).

###  5.2 Statistics of the subhalo population

We analyze our pixelated reconstruction for lens PS following a Fourier power-
spectrum analysis, motivated by the assumption of GRF perturbations (Eq. 19).
We show on Fig. 8 the azimuthally averaged power-spectra from the three
different pixelated models explored in this work: the ‘ideal’ (Fig. 4), the
‘fiducial’ (Fig. 5), and the ‘high-res’ models (i.e. finer
${\boldsymbol{\psi}_{\rm pix}}$ pixels, twice the data pixels). We compute the
power-spectra inside the region of Fig. 5, in order to consider only features
in the region of interest. We then compare the obtained power-spectra with
those from the input perturbations by fitting linear relations in log-space
and list the resulting best-fit values for $\sigma_{\rm GRF}^{2}$ and
$\beta_{\rm GRF}$ in Table. 2. We note that the first bin is excluded from the
linear fit because it corresponds to a wavenumber that translates to the
roughly the size of the region used for computing the power-spectra, hence it
is no informative. The quoted uncertainties are estimated from the least-
square fit.

Figure 8: Azimuthally averaged power spectra of input and modeled
(${\boldsymbol{\psi}_{\rm pix}}$) perturbations for lens PS (computed within
the same mask as in Fig. 5). Linear relations used to estimate the variance
$\sigma^{2}_{\rm GRF}$ and slope $\beta_{\rm GRF}$ are overlayed (dashed
lines), assuming the underlying perturbing field is a GRF (uncertainties are
estimated from the least-squares fit). The crosses indicate the bin ($2\times
10^{-1}$) that is ignored for the linear fit, as it corresponds to the size of
the mask.  source

As expected, the model in the ideal case (i.e. with fixed smooth potential and
source light) agrees very well with the input power-spectra, which translates
in a good agreement for GRF parameters as well. Regarding our fiducial model,
the amplitude is overall lower than the input, consistent with we what we
discussed in Sect. 4. At intermediate wavenumbers, the recovered power-
spectrum is close to the input one, however it is strongly attenuated at large
wavenumbers. This leads to an overall steeper slope, and translates to a
$\sim$ 5-sigma difference in $\beta_{\rm GRF}$ with respect to the input. This
attenuation of small spatial scales is fully mitigated by modeling the
perturbation on a higher-resolution grid, that allows us to better model small
scale features. Indeed, the power-spectrum of the high-res model exhibits an
excellent agreement with the input for all wavenumbers $k>4\cdot 10^{-1}$
arcsec-1. The inferred slope $\beta_{\rm GRF}$ is within 1-sigma with respect
to the input value (see Table 2).

Overall, our pixelated method correctly retrieves the locations of the main
perturbations that mimic a subhalo population, and allows us to obtain a
first-order estimate of its statistics. Our results suggest that using a
higher-resolution grid for the pixelated component allows to better recover
the full power-spectrum. However, the precise characterization of the power-
spectrum of the perturbations under the assumption of GRF is a challenging
task that requires additional priors in the model (e.g. Bayer et al. 2018;
Vernardos & Koopmans 2022).

###  5.3 Properties of multipolar structures

Based on our models of lens HM, we can recover the underlying octupole using
two different methods: fitting an octupole directly on the
${\boldsymbol{\psi}_{\rm pix}}$ model, or using the refined model that
includes explicitly an octupole profile in addition to the other components of
the lens potential.

Figure 9: Measurement of the multipole orientation from the
${\boldsymbol{\psi}_{\rm pix}}$ model of lens HM. An analytical octupole
profile (plus a constant offset) is fitted directly on the reconstructed
perturbations, and compared to the input one. Each pair of lines shows two of
the main octupole axes, to ease the visualization. The gray dotted lines are
corresponds to $\pm 3\sigma$ fits, where $\sigma$ is estimated from the FIM.

We perform the octupole fitting (i.e. we fix $m=4$) on the
${\boldsymbol{\psi}_{\rm pix}}$ model via gradient descent with three free
parameters, namely the amplitude, orientation and an additional constant
offset (remember that this offset in the potential is not constrained by the
data). After converging to the MAP solution, parameter uncertainties are
estimated from the FIM. The recovered octupole orientation, reported in Table
2, agrees extremely well with the input. However, the recovered amplitude is
lower, as expected from the pixelated reconstructions and discussed in Sect.
4.

With the refined model it is possible to measure the constraints on the
multipole order $m$, in addition to the orientation and amplitude. The MAP
values of the multipole component obtained with this refined model are
reported in Table 2, which all agree very well with the input. The parameter
$m$ is expected be more challenging to optimize, as it has non-linear effects
on the profile shape. We tested the robustness of the optimization to
different initial values $m$, and found that initializing $m$ to a value of 3
or higher leads to the correct MAP value, but setting it closer to 2 drives
the model towards a quadrupole ($m=2$), which is degenerate with the shear and
ellipticity of the smooth potential.

##  6 Computation time

Herculens uses JAX to exactly differentiate the loss function and
significantly decrease runtime (Sect. 2.6). The entire analysis of this work,
including parameters optimization and sampling, was performed on a personal
computer. We give average timings of the main modeling steps for a personal
computer 777Note that most of the timings quoted here also include an overhead
time of about 2 to 4 seconds, due to JAX “just-in-time” compilation feature.:

  * •

Optimizing the smooth analytical models (12 parameters), takes one minute for
a multi-start gradient descent with 30 starts (i.e., $\sim 2$ seconds for a
single gradient descent).

  * •

Computing the regularization weights (Eq. 2.5) takes about 20 seconds. This
step can be accelerated at least by a factor of two, by pre-computing some of
the operators involved in Eq. 14; we leave this for future improvements.

  * •

Optimizing the idealized models of Fig. 4 in which only the pixelated
potential component is optimized (1089 parameters) takes 30 seconds. This is
for $10^{3}$ iterations, which is sufficient to reach convergence.

  * •

Compared to these idealized models, the run time is virtually identical for
optimizing the full models of Fig. 5, due to the marginal increase in the
number parameters (1101 parameters).

  * •

Computing the FIM and its inverse typically takes 20 seconds.

  * •

Performing HMC sampling for a total of 500 samples takes about 1.1 hours for a
single chain.

These numbers can be extended to the modeling of a typical HST observation of
a strong lens. For instance, modeling the lens SDSS J$1630+4520$ (from the
SLACS sample, Bolton et al. 2006), with a smooth lens potential and a
pixelated source regularized with wavelets (as in Galan et al. 2021) showed
that convergence to the MAP takes about 1.5 minutes for a single gradient
descent (still on a personal computer). This includes pre-optimization steps
with a smooth model for source whose complexity is progressively improved with
a pixelated source. Fitting the lens light with analytical profiles would only
lead to a marginal increase of the total run time ($\sim 10$ seconds). Next,
modeling deviations to the smooth lens potential assuming the initial model
fits reasonably well the data would require about 30 seconds for a single
gradient descent, similar to the models presented in this work. Finally,
sampling the full parameter space using HMC would take from one to two hours
for $\sim\mathcal{O}(10^{3})$ samples, which we expect to be sufficient to
ensure well-sampled posterior distributions.

While the timings quoted above demonstrate that our code can be conveniently
run on a single CPU, they do not reveal the full potential of the method.
Herculens is fully based on JAX so it supports large scale parallelization and
GPU capabilities, which can lead to dramatic improvements in terms of
computation time (see e.g. Gu et al. 2022). This will be explored in future
analyses of real data sets.

##  7 Summary and conclusion

In this work we developed and applied a novel lens modeling method that is
able to recover perturbations to a mostly smooth lens potential with minimal
assumptions about their nature. This is made possible by modeling the
perturbations on a grid of pixels regularized using a well-established multi-
scale technique based on sparsity and wavelets. This grid of pixels is
superimposed on other analytical profiles for the joint modeling of the full
potential and the source light. We show that merging the two main state-of-
the-art modeling paradigms (analytical and pixelated) is possible within the
framework of differentiable programming. This enables us to seamlessly
optimize lens models with over one thousand individual parameters and obtain
their uncertainties, either via Fisher information analysis or gradient-
informed HMC sampling.

We summarize the key results of this work as follows:

  * •

We extend our previous work in Joseph et al. (2019) and Galan et al. (2021),
by introducing a pixelated mass component in addition to the smooth lens
potential, and recover three different types of perturbations: a localized DM
subhalo, a population of such subhalos, and high-order moments in the lens
potential.

  * •

Sparse regularization is usually performed via iterative algorithms to
converge to the solution, which makes it challenging to incorporate within
standard lens modeling codes. In this work we demonstrate that the solution
can also be obtained via gradient descent, capitalizing on the access to
derivatives of the loss function with automatic differentiation.

  * •

Differentiable programming enables the simultaneous optimization of analytical
and pixelated components. One benefit is that the perturbative approach of
Koopmans (2005) for pixelated perturbations to the lens potential is no longer
warranted, although we can still use it to estimate regularization weights.
This is because the full inverse problem can now be solved from the explicit
superposition of smooth analytical and pixelated components, jointly optimized
via gradient descent.

  * •

For each type of perturbation explored in this work, the main features are
correctly recovered by the pixelated potential component. In particular, the
signature of a localized subhalo and octupolar structures are accurately
captured. The subhalo population, simulated here as a GRF, clearly represents
a more challenging situation, although the main over- and under-density
regions can still be recovered.

  * •

We test for a model-independent recovery of the DM subhalo mass, directly from
the reconstructed pixelated potential model. We find that the resulting
measurement of the mass is sensitive to the region where the surface density
is integrated, leading to either an under-estimation of the mass, or a value
in agreement with the true subhalo mass. Nevertheless, switching to an
analytical profile for the subhalo, as is standard practice, allows to
robustly infer its mass.

  * •

The statistical properties of the subhalo population can be recovered via the
power spectrum of the pixelated model. The underlying GRF parameters, which
effectively act as parameters of the subhalo mass function, are recovered with
our higher resolution model. Nevertheless, in a real-world scenario, we
advocate for a comparison of different model variants, typically with
different grid resolutions, and a cross-checking with other methods relying on
more informative priors (e.g. as in Vernardos & Koopmans 2022) to improve the
robustness of the inference.

  * •

High-order moments in the lens potential (here as an octupole) are remarkably
well recovered, despite being small in amplitude compared to lower-order
moments such as quadrupoles (e.g. external shear). While the amplitude is
biased low in the pixelated reconstruction, the octupole orientation is
accurately measured, either from the pixelated model, or using a refined model
including a multipole component in the lens potential explicitly.

Our method is readily applicable to real HST observations of EELs, such as the
systems presented in Oldham & Auger (2018), as the source surface brightness
is smooth. This assumption of smoothness, although well motivated by real
observations, is arguably the strongest assumption of this work. Indeed, there
are many situations in which the source galaxy is more complex, featuring for
instance spiral arms and localized clumpy star forming regions. This requires
the joint modeling of deviations from smoothness both in the source and the
lens potential, which is challenging due to degeneracies between those two
components, as some of the lensing features may be equally well explained by
underlying features in the source or in the potential.

Nevertheless, the modeling methods presented in this work and implemented in
Herculens enable the joint modeling of more complex sources on a grid of
pixels in the source plane, which we will apply in future analyses. Moreover,
the flexible framework provided by differentiable programming allows one to
implement a large fraction (if not all) of the modeling methods explored in
the literature so far. The careful comparison between their different
assumptions is absolutely critical to mitigate degeneracies and robustly infer
key physical quantities, including constraints on the subhalo mass function
and the mass of DM particles. Recently, there have been a few works along
these lines, including the re-analysis of a subhalo detection based on
different assumptions such as the shape of the subhalo mass profile and its
redshift (Çağan Şengül et al. 2021), or the thorough comparison of different
modeling codes on the same data set (Shajib et al. 2022).

In this work we limit ourselves to three categories of perturbations in the
lens potential, however there exist others. One particular departure from the
simple elliptical profiles is the presence of a disk component in the lens
galaxy. In particular, such a disk structure can have observational effects
that are similar to flux ratio anomalies in multiply imaged lensed quasars.
These anomalies are usually considered as a tracer for DM subhalos, hence
unique probes of the subhalo mass function. In a series of papers based on
cosmological simulations and real observations, Hsueh et al. (2018, and
references therein) showed that not taking into account these baryonic effects
on the lens potential can bias the inferred constraints on DM properties. The
detection and modeling of a disk component in the lens potential is therefore
a situation where our pixelated reconstruction method could be successful.

Herculens is built on top of differentiable programming libraries that readily
enable massive parallelization and GPU acceleration. This important aspect has
been extensively discussed in Gu et al. (2022), where authors demonstrated the
efficiency of gradient-informed methods for modeling large samples of strong
lenses. Their results remain relevant and are applicable in the context of our
work as well. This is a promising path towards the modeling of the tens of
thousands strong lenses that will be discovered in the era of Vera Rubin and
Euclid survey telescopes (Collett 2015). Moreover, highly flexible yet
computationally efficient models will be essential to handle the even higher
resolution data sets soon delivered by the James Webb Space Telescope and
future extremely large telescopes.

###### Acknowledgements.

The authors thank Luca Biggio, Martin Millon, Eric Paic and Simona Vegetti for
useful feedback and discussion on the present work. The authors also thank
Simon Birrer for useful discussion, and for making the modeling software
package Lenstronomy open-source, which made the development of Herculens much
easier. This programme is supported by the Swiss National Science Foundation
(SNSF) and by the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (COSMICLENS: grant agreement No
787886). GV has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodovska-Curie grant
agreement No 897124. This research has also made use of SciPy (Virtanen et al.
2020), NumPy (Oliphant 2006; Van Der Walt et al. 2011), Matplotlib (Hunter
2007), Astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018),
Jupyter (Kluyver et al. 2016) and GetDist (Lewis 2019).

## References

  * Abramo (2012) Abramo, L. R. 2012, MNRAS, 420, 2042
  * Adam et al. (2022) Adam, A., Perreault-Levasseur, L., & Hezaveh, Y. 2022, arXiv e-prints, arXiv:2207.01073
  * Albert (2020a) Albert, J. G. 2020a, arXiv e-prints, arXiv:2012.15286
  * Albert (2020b) Albert, J. G. 2020b, arXiv e-prints, arXiv:2012.15286
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Bayer et al. (2018) Bayer, D., Chatterjee, S., Koopmans, L. V. E., et al. 2018, arXiv e-prints, arXiv:1803.05952
  * Belgacem et al. (2019) Belgacem, E., Calcagni, G., Crisostomi, M., et al. 2019, J. Cosmology Astropart. Phys., 2019, 024
  * Birrer & Amara (2018) Birrer, S. & Amara, A. 2018, Physics of the Dark Universe, 22, 189
  * Birrer et al. (2015) Birrer, S., Amara, A., & Refregier, A. 2015, ApJ, 813, 102
  * Birrer et al. (2017) Birrer, S., Amara, A., & Refregier, A. 2017, J. Cosmology Astropart. Phys., 2017, 037
  * Birrer et al. (2021) Birrer, S., Shajib, A., Gilman, D., et al. 2021, The Journal of Open Source Software, 6, 3283
  * Blandford et al. (2001) Blandford, R., Surpi, G., & Kundić, T. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 237, Gravitational Lensing: Recent Progress and Future Go, ed. T. G. Brainerd & C. S. Kochanek, 65
  * Blumenthal et al. (1986) Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, ApJ, 301, 27
  * Bolton et al. (2006) Bolton, A. S., Burles, S., Koopmans, L. V. E., Treu, T., & Moustakas, L. A. 2006, ApJ, 638, 703
  * Boylan-Kolchin et al. (2011) Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs
  * Brehmer et al. (2019) Brehmer, J., Mishra-Sharma, S., Hermans, J., Louppe, G., & Cranmer, K. 2019, ApJ, 886, 49
  * Çağan Şengül et al. (2021) Çağan Şengül, A., Dvorkin, C., Ostdiek, B., & Tsang, A. 2021, arXiv e-prints, arXiv:2112.00749
  * Chatterjee (2019) Chatterjee, S. 2019, PhD thesis, University of Groningen
  * Chatterjee & Koopmans (2018) Chatterjee, S. & Koopmans, L. V. E. 2018, MNRAS, 474, 1762
  * Chianese et al. (2020) Chianese, M., Coogan, A., Hofma, P., Otten, S., & Weniger, C. 2020, MNRAS[arXiv:1910.06157]
  * Claeskens et al. (2006) Claeskens, J. F., Sluse, D., Riaud, P., & Surdej, J. 2006, A&A, 451, 865
  * Coles et al. (2014) Coles, J. P., Read, J. I., & Saha, P. 2014, MNRAS, 445, 2181
  * Collett (2015) Collett, T. E. 2015, ApJ, 811, 20
  * Coogan et al. (2020) Coogan, A., Karchev, K., & Weniger, C. 2020, arXiv e-prints, arXiv:2010.07032
  * Cramér (1999) Cramér, H. 1999, Mathematical Methods of Statistics, Princeton Landmarks in Mathematics and Physics (Princeton University Press)
  * Cyr-Racine et al. (2019) Cyr-Racine, F.-Y., Keeton, C. R., & Moustakas, L. A. 2019, Phys. Rev. D, 100, 023013
  * de Blok (2010) de Blok, W. J. G. 2010, Advances in Astronomy, 2010, 789293
  * Despali et al. (2022) Despali, G., Vegetti, S., White, S. D. M., et al. 2022, MNRAS, 510, 2480
  * Diaz Rivero & Dvorkin (2020) Diaz Rivero, A. & Dvorkin, C. 2020, Phys. Rev. D, 101, 023515
  * Duane et al. (1987) Duane, S., Kennedy, A., Pendleton, B. J., & Roweth, D. 1987, Physics Letters B, 195, 216
  * Dubinski (1994) Dubinski, J. 1994, ApJ, 431, 617
  * Erkal et al. (2016) Erkal, D., Belokurov, V., Bovy, J., & Sanders, J. L. 2016, MNRAS, 463, 102
  * Frigo et al. (2019) Frigo, M., Naab, T., Hirschmann, M., et al. 2019, MNRAS, 489, 2702
  * Galan et al. (2021) Galan, A., Peel, A., Joseph, R., Courbin, F., & Starck, J. L. 2021, A&A, 647, A176
  * Gilman et al. (2019) Gilman, D., Birrer, S., Treu, T., Nierenberg, A., & Benson, A. 2019, MNRAS, 487, 5721
  * Gu et al. (2022) Gu, A., Huang, X., Sheu, W., et al. 2022, arXiv e-prints, arXiv:2202.07663
  * Hao et al. (2006) Hao, C. N., Mao, S., Deng, Z. G., Xia, X. Y., & Wu, H. 2006, MNRAS, 370, 1339
  * He et al. (2022) He, Q., Nightingale, J., Massey, R., et al. 2022, arXiv e-prints, arXiv:2202.10191
  * Hessel et al. (2020) Hessel, M., Budden, D., Viola, F., et al. 2020, Optax: composable gradient transformation and optimisation, in JAX!
  * Hezaveh et al. (2016a) Hezaveh, Y., Dalal, N., Holder, G., et al. 2016a, J. Cosmology Astropart. Phys., 2016, 048
  * Hezaveh et al. (2016b) Hezaveh, Y. D., Dalal, N., Marrone, D. P., et al. 2016b, ApJ, 823, 37
  * Hoffman & Gelman (2011) Hoffman, M. D. & Gelman, A. 2011, arXiv e-prints, arXiv:1111.4246
  * Horne (1986) Horne, K. 1986, Publications of the Astronomical Society of the Pacific, 98, 609
  * Hsueh et al. (2018) Hsueh, J.-W., Despali, G., Vegetti, S., et al. 2018, MNRAS, 475, 2438
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Joseph et al. (2019) Joseph, R., Courbin, F., Starck, J. L., & Birrer, S. 2019, A&A, 623, A14
  * Karchev et al. (2022) Karchev, K., Coogan, A., & Weniger, C. 2022, MNRAS[arXiv:2105.09465]
  * Keeton (2001) Keeton, C. R. 2001, arXiv e-prints, astro
  * Keeton et al. (2000) Keeton, C. R., Falco, E. E., Impey, C. D., et al. 2000, ApJ, 542, 74
  * Kennedy & Eberhart (2001) Kennedy, J. & Eberhart, R. C. 2001, Swarm Intelligence (San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.)
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt, IOS Press, 87 – 90
  * Klypin et al. (1999) Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82
  * Koopmans (2005) Koopmans, L. V. E. 2005, MNRAS, 363, 1136
  * Lanusse et al. (2016) Lanusse, F., Starck, J. L., Leonard, A., & Pires, S. 2016, A&A, 591, A2
  * Lee et al. (2020) Lee, W., Yu, H., Rival, X., & Yang, H. 2020, arXiv e-prints, arXiv:2006.06903
  * Lewis (2019) Lewis, A. 2019, arXiv e-prints, arXiv:1910.13970
  * Li et al. (2017) Li, R., Frenk, C. S., Cole, S., Wang, Q., & Gao, L. 2017, MNRAS, 468, 1426
  * Moore (1994) Moore, B. 1994, Nature, 370, 629
  * Moore et al. (1999) Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJ, 524, L19
  * Nadler et al. (2021) Nadler, E. O., Birrer, S., Gilman, D., et al. 2021, ApJ, 917, 7
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
  * Neal (2011) Neal, R. 2011, in Handbook of Markov Chain Monte Carlo, 113–162
  * Nightingale & Dye (2015) Nightingale, J. W. & Dye, S. 2015, MNRAS, 452, 2940
  * Nocedal & Wright (2006) Nocedal, J. & Wright, S. J. 2006, Numerical Optimization, 2nd edn., Springer Series in Operations Research and Financial Engineering (New York: Springer)
  * Oldham et al. (2017) Oldham, L., Auger, M. W., Fassnacht, C. D., et al. 2017, MNRAS, 465, 3185
  * Oldham & Auger (2018) Oldham, L. J. & Auger, M. W. 2018, MNRAS, 476, 133
  * Oliphant (2006) Oliphant, T. E. 2006, A guide to NumPy, Vol. 1 (Trelgol Publishing USA)
  * Ostdiek et al. (2022) Ostdiek, B., Diaz Rivero, A., & Dvorkin, C. 2022, A&A, 657, L14
  * Papastergis et al. (2015) Papastergis, E., Giovanelli, R., Haynes, M. P., & Shankar, F. 2015, A&A, 574, A113
  * Paykari et al. (2014) Paykari, P., Lanusse, F., Starck, J. L., Sureau, F., & Bobin, J. 2014, A&A, 566, A77
  * Philcox et al. (2021) Philcox, O. H. E., Sherwin, B. D., Farren, G. S., & Baxter, E. J. 2021, Phys. Rev. D, 103, 023538
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Powell et al. (2022) Powell, D. M., Vegetti, S., McKean, J. P., et al. 2022, arXiv e-prints, arXiv:2207.03375
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123
  * Raichoor et al. (2021) Raichoor, A., de Mattia, A., Ross, A. J., et al. 2021, MNRAS, 500, 3254
  * Riess (2020) Riess, A. G. 2020, Nature Reviews Physics, 2, 10
  * Saha & Williams (1997) Saha, P. & Williams, L. L. R. 1997, MNRAS, 292, 148
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Shajib et al. (2019) Shajib, A. J., Birrer, S., Treu, T., et al. 2019, MNRAS, 483, 5649
  * Shajib et al. (2021) Shajib, A. J., Treu, T., Birrer, S., & Sonnenfeld, A. 2021, MNRAS, 503, 2380
  * Shajib et al. (2022) Shajib, A. J., Wong, K. C., Birrer, S., et al. 2022, arXiv e-prints, arXiv:2202.11101
  * Shu et al. (2017) Shu, Y., Brownstein, J. R., Bolton, A. S., et al. 2017, ApJ, 851, 48
  * Sonnenfeld et al. (2018) Sonnenfeld, A., Leauthaud, A., Auger, M. W., et al. 2018, MNRAS, 481, 164
  * Springel et al. (2006) Springel, V., Frenk, C. S., & White, S. D. M. 2006, Nature, 440, 1137
  * Starck et al. (2015) Starck, J., Murtagh, F., & Fadili, J. 2015, Sparse Image and Signal Processing: Wavelets and Related Geometric Multiscale Analysis (Cambridge University Press)
  * Starck et al. (2007) Starck, J.-L., Fadili, J., & Murtagh, F. 2007, IEEE Transactions on Image Processing, 16, 297
  * Suyu et al. (2009) Suyu, S. H., Marshall, P. J., Blandford, R. D., et al. 2009, ApJ, 691, 277
  * Suyu et al. (2006) Suyu, S. H., Marshall, P. J., Hobson, M. P., & Blandford, R. D. 2006, MNRAS, 371, 983
  * Tagore & Keeton (2014) Tagore, A. S. & Keeton, C. R. 2014, MNRAS, 445, 694
  * Toomre & Toomre (1972) Toomre, A. & Toomre, J. 1972, ApJ, 178, 623
  * Trotter et al. (2000) Trotter, C. S., Winn, J. N., & Hewitt, J. N. 2000, ApJ, 535, 671
  * Van de Vyvere et al. (2021) Van de Vyvere, L., Gomer, M. R., Sluse, D., et al. 2021, arXiv e-prints, arXiv:2112.03932
  * Van de Vyvere et al. (2022) Van de Vyvere, L., Sluse, D., Gomer, M. R., & Mukherjee, S. 2022, arXiv e-prints, arXiv:2206.00022
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Varma et al. (2020) Varma, S., Fairbairn, M., & Figueroa, J. 2020, arXiv e-prints, arXiv:2005.05353
  * Vegetti & Koopmans (2009) Vegetti, S. & Koopmans, L. V. E. 2009, MNRAS, 392, 945
  * Vegetti et al. (2014) Vegetti, S., Koopmans, L. V. E., Auger, M. W., Treu, T., & Bolton, A. S. 2014, MNRAS, 442, 2017
  * Vegetti et al. (2010) Vegetti, S., Koopmans, L. V. E., Bolton, A., Treu, T., & Gavazzi, R. 2010, MNRAS, 408, 1969
  * Vegetti et al. (2012) Vegetti, S., Lagattuta, D. J., McKean, J. P., et al. 2012, Nature, 481, 341
  * Vernardos (2021) Vernardos, G. 2021, arXiv e-prints, arXiv:2106.04344
  * Vernardos & Koopmans (2022) Vernardos, G. & Koopmans, L. V. E. 2022, arXiv e-prints, arXiv:2202.09378
  * Vernardos et al. (2020) Vernardos, G., Tsagkatakis, G., & Pantazis, Y. 2020, MNRAS, 499, 5641
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2020) Vogelsberger, M., Marinacci, F., Torrey, P., & Puchwein, E. 2020, Nature Reviews Physics, 2, 42
  * Warren & Dye (2003) Warren, S. J. & Dye, S. 2003, ApJ, 590, 673
  * Wengert (1964) Wengert, R. E. 1964, Commun. ACM, 7, 463
  * Zhuang et al. (2020) Zhuang, J., Tang, T., Ding, Y., et al. 2020, arXiv e-prints, arXiv:2010.07468
  * Zubovas & King (2012) Zubovas, K. & King, A. 2012, ApJ, 745, L34

##  Appendix A Multi-scale regularization weights

In this section we give details about the estimation of the weights in Eq. 2.5
used to scale the regularization strength of the pixelated potential
${\boldsymbol{\psi}_{\rm pix}}$. Let us assume that we have in hand a model
fitted to the data without ${\boldsymbol{\psi}_{\rm pix}}$, such that we can
assume that a sufficiently good model of the source ${\boldsymbol{s}}$ and the
smooth lens potential are known.

The minimization problem corresponding to the full lens potential,
parametrized as
${\boldsymbol{\psi}}=\tilde{{\boldsymbol{\psi}}}({\boldsymbol{\eta}}_{\boldsymbol{\psi}})+{\boldsymbol{\psi}_{\rm
pix}}$ can be written as

$\displaystyle\small\operatorname*{arg\,min}_{{\boldsymbol{\eta}}_{{\boldsymbol{\psi}}},{\boldsymbol{\psi}_{\rm
pix}}}\ $
$\displaystyle\frac{1}{2}\,\Big{[}\,{\boldsymbol{d}}-\boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{L}}_{{\boldsymbol{\psi}}}{\boldsymbol{s}}\,\Big{]}^{\top}{\rm\bf
C}_{\boldsymbol{d}}^{-1}\
\Big{[}\,{\boldsymbol{d}}-\boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{L}}_{{\boldsymbol{\psi}}}{\boldsymbol{s}}\,\Big{]}$
$\displaystyle+\lambda_{{\boldsymbol{\psi}}}\,\left\lVert\,{\rm\bf
W}_{\psi}\circ{\boldsymbol{\Phi}}^{\top}\,{\boldsymbol{\psi}_{\rm
pix}}\,\right\rVert_{1}\ ,$ (21)

where $\boldsymbol{\mathsf{L}}_{{\boldsymbol{\psi}}}$ depends on the the
smooth potential parameters ${\boldsymbol{\eta}}_{{\boldsymbol{\psi}}}$ and
the pixelated component ${\boldsymbol{\psi}_{\rm pix}}$. The regularization
term is identical to Eq. 2.5.

In contrast to the source reconstruction problem, namely where the variable of
interest is ${\boldsymbol{s}}$, the model we optimize here does not depend
linearly on ${\boldsymbol{\eta}}_{{\boldsymbol{\psi}}}$ and
${\boldsymbol{\psi}_{\rm pix}}$. Consequently, the gradient of the data-
fidelity term in Eq. A does not have a simple closed-form expression, which
prevents us from propagating the noise from the data to the lens potential as
easily as to the source plane.

We address this issue by considering the perturbative approach of Koopmans
(2005) as a means to relate the data space to the potential space. We argue
that we can use the weights obtained with this model for the minimization
problem above.

Based on a Taylor expansion of the lens equation with respect to small
deviations $\delta{\boldsymbol{\psi}}$ to an underlying smooth potential
$\tilde{}{\boldsymbol{\psi}}({\boldsymbol{\eta}}_{{\boldsymbol{\psi}}})$,
Koopmans (2005) established the following linear relation

$\displaystyle\delta{\boldsymbol{d}}$
$\displaystyle=\underbrace{-\boldsymbol{\mathsf{B}}\,\boldsymbol{\mathsf{D}}_{s}\\!\left({\boldsymbol{s}}\right)\boldsymbol{\mathsf{D}}_{\delta\psi}}_{\mbox{$\boldsymbol{\mathsf{D}}$}}\,\delta{\boldsymbol{\psi}}\
,$ (22)

which corresponds to Eq. 14 rewritten with linear operators. The term
$\delta{\boldsymbol{d}}\equiv{\boldsymbol{d}}-\tilde{}{\boldsymbol{m}}$
represents residuals between the data and a model
$\tilde{}{\boldsymbol{m}}=\boldsymbol{\mathsf{B}}\boldsymbol{\mathsf{L}}_{\tilde{}{\boldsymbol{\psi}}}{\boldsymbol{s}}$
without perturbations to the smooth potential (i.e.
$\delta{\boldsymbol{\psi}}=0$). The operator $\boldsymbol{\mathsf{D}}_{s}$
contains spatial derivatives of the source light model with respect to source
plane coordinates mapped from the data grid via
$\boldsymbol{\mathsf{L}}_{\tilde{}{\boldsymbol{\psi}}}$. The operator
$\boldsymbol{\mathsf{D}}_{\delta\psi}$ combines bilinear interpolation and
finite difference coefficients to compute the spatial derivatives of
$\delta{\boldsymbol{\psi}}$ on the data grid. We refer the reader to Koopmans
(2005) for more details about the exact structure of these operators, that we
implement as matrices.

From the linear relation of Eq. 22, we can formulate yet a new minimization
problem for $\delta{\boldsymbol{\psi}}$ as

$\displaystyle\small\operatorname*{arg\,min}_{\delta{\boldsymbol{\psi}}}\ $
$\displaystyle\frac{1}{2}\,\Big{[}\,\delta{\boldsymbol{d}}-\boldsymbol{\mathsf{D}}\delta{\boldsymbol{\psi}}\Big{]}^{\top}{\rm\bf
C}_{\delta{\boldsymbol{d}}}^{-1}\
\Big{[}\,\delta{\boldsymbol{d}}-\boldsymbol{\mathsf{D}}\delta{\boldsymbol{\psi}}\Big{]}$
$\displaystyle+\lambda_{\delta{\boldsymbol{\psi}}}\,\left\lVert\,{\rm\bf
W}_{\delta\psi}\circ{\boldsymbol{\Phi}}^{\top}\,\delta{\boldsymbol{\psi}}\,\right\rVert_{1}\
,$ (23)

where ${\rm\bf C}_{\delta{\boldsymbol{d}}}$ is the covariance matrix
associated to $\delta{\boldsymbol{d}}$. From the definition of model
residuals, we have

$\displaystyle{\rm\bf C}_{\delta{\boldsymbol{d}}}$ $\displaystyle\equiv{\rm
cov}\big{(}\delta{\boldsymbol{d}},\delta{\boldsymbol{d}}\big{)}$
$\displaystyle={\rm
cov}\big{(}{\boldsymbol{d}}-\tilde{{\boldsymbol{m}}},{\boldsymbol{d}}-\tilde{{\boldsymbol{m}}}\big{)}$
$\displaystyle={\rm\bf C}_{\boldsymbol{d}}+{\rm\bf
C}_{\tilde{{\boldsymbol{m}}}}-2\,{\rm
cov}\big{(}{\boldsymbol{d}},\tilde{{\boldsymbol{m}}}\big{)}$
$\displaystyle\approx{\rm\bf C}_{\boldsymbol{d}}\ .$ (24)

The cross-covariance matrix ${\rm
cov}\big{(}{\boldsymbol{d}},\tilde{{\boldsymbol{m}}}\big{)}$ is zero, since
the data and model $\tilde{{\boldsymbol{m}}}$ are uncorrelated (nature does
not correlate with our model). Additionally, the covariance term ${\rm\bf
C}_{\tilde{{\boldsymbol{m}}}}$ is much smaller than ${\rm\bf
C}_{\boldsymbol{d}}$ in the case of smooth analytical profiles for the
potential and source (we checked this numerically). This is because the fully
smooth model $\tilde{{\boldsymbol{m}}}$ is strongly limited by the shape of
the analytical profiles described by a small number of degrees of freedom,
constrained by a much larger number of data pixels. Although not explicitly
specified for conciseness, the elements of ${\rm\bf C}_{\boldsymbol{d}}$
depend on the model $\tilde{}{\boldsymbol{m}}$ to estimate Poisson noise due
to the lensed features.

We now need to establish which operators are necessary to propagate the noise
from $\delta{\boldsymbol{d}}$ space to wavelet coefficients of
$\delta{\boldsymbol{\psi}}$, in order to compute the coefficients of the
matrix ${\rm\bf W}_{\delta{\boldsymbol{\psi}}}$. We do so by considering a
gradient descent update based on the data-fidelity term of Eq. A (e.g. Lanusse
et al. 2016), which reads

$\displaystyle\delta{\boldsymbol{\psi}}^{(n+1)}=\delta{\boldsymbol{\psi}}^{(n)}+\overbrace{\boldsymbol{\mathsf{D}}^{\top}{\rm\bf
C}_{\delta{\boldsymbol{d}}}^{-1}\left(\delta{\boldsymbol{d}}-\boldsymbol{\mathsf{D}}\delta{\boldsymbol{\psi}}^{(n)}\right)}^{\text{negative
data-fidelity\ gradient}}\ .$

The above gradient translates model residuals (the term in parenthesis) into
the coefficients of $\delta{\boldsymbol{\psi}}$. Therefore, we see that the
propagation of a given noise realization $\delta{\boldsymbol{d}}_{N}$ to the
potential pixels is $\boldsymbol{\mathsf{D}}^{\top}{\rm\bf
C}_{\delta{\boldsymbol{d}}}^{-1}\delta{\boldsymbol{d}}_{N}$. We simply apply
the wavelet transform operator to obtain wavelet coefficients instead, which
leads to ${\boldsymbol{\Phi}}^{\top}\boldsymbol{\mathsf{D}}^{\top}{\rm\bf
C}_{\delta{\boldsymbol{d}}}^{-1}\delta{\boldsymbol{d}}_{N}$. Using Monte-Carlo
simulations of $\delta{\boldsymbol{d}}_{N}$, we can compute the standard
deviation $\sigma_{\delta{\boldsymbol{\psi}}}^{(i,j,k)}$ for each potential
pixel $(i,j)$ in each wavelet scale $k$, and populate the matrix ${\rm\bf
W}_{\delta{\boldsymbol{\psi}}}$ accordingly.

As a final step, we return to the original problem described by Eq. A. In this
problem, the pixelated potential parameters are ${\boldsymbol{\psi}_{\rm
pix}}$, instead of $\delta{\boldsymbol{\psi}}$. First, let us assume that both
${\boldsymbol{\psi}_{\rm pix}}$ and $\delta{\boldsymbol{\psi}}$ are defined on
the same grid of potential pixels. Second, ${\boldsymbol{\psi}_{\rm pix}}$ is
used in combination with a smooth lens potential component
$\tilde{{\boldsymbol{\psi}}}({\boldsymbol{\eta}}_{\psi})$, such that it
captures deviations to the smooth potential in a similar way as
$\delta{\boldsymbol{\psi}}$. Third, from Eq. A, we see that the noise in the
data (characterized by ${\rm\bf C}_{\boldsymbol{d}}$) is a good approximation
of the noise in the model residuals $\delta{\boldsymbol{d}}$ (characterized by
${\rm\bf C}_{\delta{\boldsymbol{d}}}$). Therefore, the standard deviation of
the noise in potential space, at the location of potential pixels and in each
wavelet scale, can be approximated by the values
$\sigma_{\delta{\boldsymbol{\psi}}}^{(i,j,k)}$ computed above, namely ${\rm\bf
W}_{\psi}\approx{\rm\bf W}_{\delta\psi}$.

In summary, we have used the linear expression of Eq. 22 as a way to set the
weights of our regularization term (Eq. A) based on the noise levels from the
data. In this work, we use two different wavelet transforms, namely the
starlet and Battle-Lemarié wavelet transforms. Thus we compute the
corresponding weights following the above procedure to obtain ${\rm\bf
W}_{\psi,\rm st}$ and ${\rm\bf W}_{\psi,\rm bl}$ respectively (only the
wavelet operator ${\boldsymbol{\Phi}}^{\top}$ is different). We show in Fig. 2
the corresponding weights for each scale of the starlet transform.

##  Appendix B Analytical profiles

The smooth models used throughout this work are based on analytical profiles,
which we define in detail here. Elliptical profiles are described with
coordinates $(\theta_{1},\theta_{2})$, obtained by rotating the original
coordinates $\@vec{\theta}\equiv(\theta_{x},\theta_{y})$ along the major axis
of the ellipse with position angle $\phi$.

To describe the smooth component of the lens potential, we use the singular
isothermal ellipsoid (SIE) profile. This profile is originally defined in
surface mass density (convergence) but also has analytical formulae for the
potential and deflection angles. The SIE potential is given by (e.g. Keeton
2001)

$\displaystyle\psi_{\rm SIE}(\theta_{1},\theta_{2})$
$\displaystyle=\frac{\theta_{\rm E}\sqrt{q_{\rm m}}\leavevmode\nobreak\
\theta_{1}}{\sqrt{1-q_{\rm m}^{2}}}\,{\rm arctan}\left(\frac{\sqrt{1-q_{\rm
m}^{2}}\leavevmode\nobreak\ \theta_{1}}{\sqrt{q_{\rm
m}^{2}\theta_{1}^{2}+\theta_{2}^{2}}}\right)$ $\displaystyle+\frac{\theta_{\rm
E}\sqrt{q_{\rm m}}\leavevmode\nobreak\ \theta_{2}}{\sqrt{1-q_{\rm
m}^{2}}}\,{\rm arctanh}\left(\frac{\sqrt{1-q_{\rm m}^{2}}\leavevmode\nobreak\
\theta_{2}}{\sqrt{q_{\rm m}^{2}\theta_{1}^{2}+\theta_{2}^{2}}}\right)\ ,$ (25)

where $\theta_{\rm E}$ is the Einstein radius, $q_{\rm m}$ is the axis ratio
of the elliptical profile, and the corresponding position angle is $\phi_{\rm
m}$. In practice, we do not optimize the axis ratio and position angle, but
rather ellipticity components $\\{e_{1,\psi},e_{2,\psi}\\}$ to prevent
sampling issues with the periodicity of $\phi_{\rm m}$ (particularly for small
values of $q_{\rm m}$). These components are defined as

$\displaystyle\begin{cases}e_{1,\psi}&=\frac{1-q_{\rm m}}{1+q_{\rm
m}}\cos\left(2\phi_{\rm m}\right)\\\ e_{2,\psi}&=\frac{1-q_{\rm m}}{1+q_{\rm
m}}\sin\left(2\phi_{\rm m}\right)\ .\end{cases}$ (26)

The singular isothermal sphere (SIS) is a particular case of an SIE with no
ellipticity ($q_{\rm m}=1$).

The influence of other galaxies in the vicinity of the main deflector—whose
lensing effects remain linear with respect to the lensed source—is modeled as
a uniform external shear. It has amplitude $\gamma_{\rm ext}$ and orientation
$\phi_{\rm ext}$. The corresponding lensing potential in polar coordinates
$(r,\phi)$ is

$\displaystyle\psi_{\rm ext}(r,\phi)=\frac{1}{2}\,\gamma_{\rm
ext}\,r^{2}\cos\left(2\phi-2\phi_{\rm ext}\right)\ .$ (27)

The corresponding shear ellipticity parameters that we used for optimisation
are defined as

$\displaystyle\begin{cases}\gamma_{1,\rm ext}&=\gamma_{\rm
ext}\cos\left(2\phi_{\rm ext}\right)\\\ \gamma_{2,\rm ext}&=\gamma_{\rm
ext}\sin\left(2\phi_{\rm ext}\right)\ .\end{cases}$ (28)

The elliptical Sérsic profile, suitable for modeling smooth galaxy light
distributions, is defined as (Sérsic 1963)

$\displaystyle I_{\rm Sersic}\left(\theta_{1},\theta_{2}\right)=I_{\rm
eff}\,\mathrm{exp}\\!\left[-b_{n}\left(\frac{\sqrt{q_{\rm
l}\theta_{1}^{2}+\theta_{2}^{2}/q_{\rm l}}}{\theta_{\rm
eff}}\right)^{\frac{1}{n_{\rm s}}}+b_{n}\right],$ (29)

where $I_{\rm eff}$ is the amplitude of the profile at the effective radius
$\theta_{\rm eff}$, and $n_{\rm s}$ is the Sérsic index which defines the
slope of the profile. The axis ratio is $q_{\rm l}$, and the corresponding
position angle is $\phi_{\rm l}$. Ellipticities are defined as for the
potential (see Eq. 28). The term $b_{n}$ is not a free parameter and is
computed such that $\theta_{\rm eff}$ is always equal to the half-light radius
of the profile.

##  Appendix C Simulated HST observations

In Table 3 we summarize all instrumental settings and specific model
assumptions that were used to generate HST observations of EELs systems. We
used the simulation software package MOLET (Vernardos 2021) to perform high-
resolution ray-tracing and to add instrumental effects (instrumental noise and
PSF convolutions).

Table 3: Instrument settings and model assumptions used for generating HST/WFC3/F160W mock observations of EELs. The coordinates are oriented such as north is up and east is right. The orientation angle is East-of-North. Observation |
---|---
Pixel size [arcsec] | 0.08
Exposure time [s] | $4\times 2400=9600$
Zero-point [mag] | 25.9463
Sky brightness [mag/arcsec2] | 22
Read noise [e-] | 4
PSF (gaussian FWHM) [arcsec] | $0.3$
Noise | readout + shot noise
Deflector, $z_{\rm d}=0.3$ |
Singular isothermal ellipsoid (SIE) |
Einstein radius, $\theta_{\rm E}$ [arcsec] | 1.6
Ellipticity, $q_{\rm m}$ | 0.73
Orientation, $\phi_{\rm m}$ | 82.5
Position, ($x$, $y$) [arcsec] | $(0,0)$
External shear |
Strength, $\gamma_{\rm ext}$ | 0.032
Orientation, $\phi_{\rm ext}$ | 144.2
Sérsic (light) |
Ellipticity, $q_{\rm l}$ | $0.73\ (=q_{\rm m})$
Orientation, $\phi_{\rm l}$ | $82.5\ (=\phi_{\rm m})$
Half-light radius, $\theta_{\rm eff}$ | 2
Sérsic index, $n_{\rm Sersic}$ | 2
Position, ($x$, $y$) [arcsec] | $(0,0)$
Magnitude | 19
Perturbations to the potential |
LS: Localized subhalo (SIS) |
Einstein radius, $\theta_{\rm E,halo}$ [arcsec] | 0.07
Position, ($x$, $y$) [arcsec] | $(1.9,-0.4)$
PS: Population of subhalos (GRF) |
Variance, $\sigma^{2}_{\rm GRF}$ | $10^{-3}$
Slope, $\beta_{\rm GRF}$ | $4$
HM: High-order moments (octupole) |
Strength, $a_{4}$ | $0.06$
Orientation, $\phi_{\rm 4}$ | $82.5=\phi_{\rm m}$
Source, $z_{\rm s}=0.7$ |
Sérsic |
Ellipticity, $q_{\rm l}$ | 0.82
Orientation, $\phi_{\rm l}$ | 170.8
Half-light radius, $\theta_{\rm eff}$ [arcsec] | 0.8
Sérsic index, $n_{\rm Sersic}$ | 2
Position ($x$, $y$) [arcsec] | $(0.4,0.15)$
Magnitude | 21

##  Appendix D Choice of pixel size for ${\boldsymbol{\psi}_{\rm pix}}$

The pixelated model ${\boldsymbol{\psi}_{\rm pix}}$ used to capture deviations
from the smooth lens potential requires a choice of pixel size. In Herculens,
this can be set to any multiplicative factor (not necessarily an integer) of
the data pixel size. In this work we use a pixel scale factor of 3 as our
fiducial model, meaning the grid on which the ${\boldsymbol{\psi}_{\rm pix}}$
component is defined has a pixel size of $3\times
0\aas@@fstack{\prime\prime}08=0\aas@@fstack{\prime\prime}24$. This resolution
is sufficient to accurately model the simulated data set as well as to
characterize the reconstructed features in the lens potential.

The multi-scale aspect of the regularization strategy detailed in Sect. 2.5
enables the independent treatment of different spatial scales in the lens
potential. This avoids the need to choose a specific pixel size, as long as it
is small compared to the lensed features (otherwise it would lead to a poor
fit to the data). Let us take the example of a pixel size much smaller than
that of the data; if there are features below a given spatial scale that are
not supported by the data (typically below the smallest detectable deflection
angle), the multi-scale wavelet decomposition combined with sparsity
constraints will suppress all wavelet coefficients below this scale.
Therefore, these features should not be visible in the reconstruction. This
means that in comparing two models that fit the data equally well, we should
select the one with the smaller pixel size. In practice, however, and from the
perspective of parameter optimization, a smaller pixel size translates into a
larger number of parameters, for which the regularization becomes even more
crucial to balance the lack of constraints provided by the data alone.
Moreover, the principle of Occam’s razor, at the basis of Bayesian approaches,
advocates for fewer parameters if it does not bring significant improvements
to the fit.

Figure 10: Best-fit models for different choices of pixel size of the
${\boldsymbol{\psi}_{\rm pix}}$ component, corresponding to 4, 3, 2, and 1.5
times the data pixel size. The modeling assumptions are otherwise identical to
the results presented in Fig. 5 (the second row is exactly the same). From
left to right, the columns show: model residuals with reduced $\chi^{2}$
value, best-fit ${\boldsymbol{\psi}_{\rm pix}}$ model, followed by each scale
of the starlet transform of ${\boldsymbol{\psi}_{\rm pix}}$. Each of these
scales are images of the same size as ${\boldsymbol{\psi}_{\rm pix}}$, above
which we indicate the scale index $j$ and the characteristic scale
$\tilde{\beta}_{j}$ (in arcsec) of the features. We see that a pixel scale of
3 is sufficient to model the full dynamical range of features supported by the
data. Finer pixels are not required, as the smaller starlet scales reveal
mainly low-significance artifacts over the level of individual pixels.  source

We investigate the impact of the pixelated model resolution further by
modeling lens LS (see Fig. 3) using four different pixel sizes for
${\boldsymbol{\psi}_{\rm pix}}$: 4, 3, 2 and 1.5 times the data pixel size
(leading to 625, 1089, 2500 and 4489 parameters, respectively). For this