∎

11institutetext:  1\. Department of Astronomy and Astrophysics
22institutetext: University of Chicago
Chicago, IL 60637, USA
22email<EMAIL_ADDRESS>
2\. Kavli Institute for Cosmological Physics 33institutetext: University of
Chicago
Chicago, IL 60637, USA
3\. Institute of Physics, Laboratory of Astrophysics 44institutetext: Ecole
Polytechnique Fédérale de Lausanne (EPFL)
Observatoire de Sauverny, 1290 Versoix, Switzerland
4\. Institute of Cosmology and Gravitation 55institutetext: University of
Portsmouth
Burnaby Rd, Portsmouth PO1 3FX, UK
5\. Instituto de Fisica y Astronomia, Facultad de Ciencias 66institutetext:
Universidad de Valparaiso
Avda. Gran Bretana 1111, Valparaiso, Chile
6\. STAR Institute 77institutetext: Quartier Agora - Allée du six Aout
19c B-4000 Liége, Belgium
7\. School of Physics and Astronomy 88institutetext: University of Minnesota
116 Church Street SE, Minneapolis MN, 55455, USA
8\. Physik-Institut 99institutetext: University of Zurich
Winterthurerstrasse 190, 8057 Zurich, Switzerland
9\. Kavli Institute for Particle Astrophysics and Cosmology and Department of
Physics 1010institutetext: Stanford University
Stanford, CA 94305, USA
10\. SLAC National Accelerator Laboratory
Menlo Park, CA, 94025
11\. Department of Physics and Astronomy 1111institutetext: Stony Brook
University
Stony Brook, NY 11794, USA
12\. Department of Physics 1212institutetext: University of Oxford
Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK
13\. INAF – Osservatorio Astronomico di Capodimonte
Via Moiariello 16, 80131, Naples, Italy
14\. Department of Physics and Astronomy 1313institutetext: University of
California, Los Angeles
430 Portola Plaza, Los Angeles, CA 90095, USA

# Strong Lensing by Galaxies

A. J. Shajib${}^{1,2}\dagger$ $\dagger$ NHFP Einstein Fellow G. Vernardos3 T.
E. Collett4 V. Motta5 D. Sluse6 L. L. R. Williams7 P. Saha8 S. Birrer9,10,11
C. Spiniello12,13 T. Treu14

(Received: date / Accepted: date)

###### Abstract

Strong lensing at the galaxy scale is useful for numerous applications in
Astrophysics and Cosmology. Some of the principal applications are studying
the mass structure, formation, and evolution of elliptical galaxies,
constraining the stellar initial mass function, and measuring cosmological
parameters. Since the first discovery of a galaxy-scale strong lens in the
eighties, this field has come a long way regarding data quality and techniques
to model the data. In this review article, we describe the most common
methodologies to model lensing observables of galaxy-scale strong lenses,
especially the imaging data, as it is the most available and informative
source of lensing observables. We review the main results from the literature
in astrophysical and cosmological applications of galaxy-scale strong lenses.
We also discuss the current limitations of the data and methodologies and
provide a future outlook of the expected development and improvements in both
aspects in the near future.

###### Keywords:

gravitational lensing: strong, galaxies: elliptical

††journal: my journal

###### Contents

  1. 1 Introduction
    1. 1.1 Brief description of lensing phenomenology at the galaxy-scale
    2. 1.2 Unique advantages of lensing as a probe of galaxy structure
  2. 2 Historical background
    1. 2.1 Initial discoveries of strong lensing systems
    2. 2.2 Prominent samples of galaxy-scale lenses
      1. 2.2.1 MG-VLA survey-based samples
      2. 2.2.2 SDSS-based samples
      3. 2.2.3 CFHTLens-based samples
      4. 2.2.4 BOSS-based samples
    3. 2.3 Other samples and ongoing efforts from recent surveys
  3. 3 Observables and analysis methods
    1. 3.1 Lensing observables
      1. 3.1.1 Imaging of the lens system
      2. 3.1.2 Time delay
    2. 3.2 Lens modeling methods
      1. 3.2.1 Likelihood-based inference
      2. 3.2.2 Machine-learning based parameter extraction
    3. 3.3 Lens mass models
      1. 3.3.1 Simply parametrized mass models
      2. 3.3.2 Common degeneracies in simple parametric modeling
      3. 3.3.3 Free-form models
      4. 3.3.4 Model ensembles
    4. 3.4 Bayesian hierarchical framework
    5. 3.5 Incorporating non-strong-lensing observables
      1. 3.5.1 Combining stellar kinematics with strong lensing
      2. 3.5.2 Combining weak lensing with strong lensing
  4. 4 Applications in galaxy properties and evolution
    1. 4.1 Galaxy mass density profile
      1. 4.1.1 Luminous (baryonic) mass profile
      2. 4.1.2 Dark matter profile
    2. 4.2 Galaxy evolution
    3. 4.3 Stellar initial mass function
    4. 4.4 Constraints on the very central densities and super-massive black holes from central images
    5. 4.5 Spiral galaxies
  5. 5 Applications in cosmology
    1. 5.1 Utilizing multiple sources at different redshifts
    2. 5.2 Utilizing stellar kinematics of single source lenses
    3. 5.3 Utilizing galaxy–galaxy lensing statistics
  6. 6 Open problems and future outlook
    1. 6.1 Selection function
      1. 6.1.1 Parametrized selection function
      2. 6.1.2 Posterior predictive test
    2. 6.2 Assumption of self-similarity
    3. 6.3 Degeneracies in strong-lensing observables
    4. 6.4 Degeneracies in non-strong-lensing observables
    5. 6.5 Comparison with galaxy simulations
    6. 6.6 Statistical framework for comparison with simulations
  7. 7 Concluding remarks

##  1 Introduction

This review article discusses applications of galaxy-scale strong lenses to
study properties of the deflector galaxies (i.e., the central lensing
galaxies). These have so far typically been massive elliptical galaxies at
$0.1\lesssim z\lesssim 1$ and strong lensing has been chiefly applied to study
their internal structure and composition. However, this decade, a large sample
of strong lenses with deflectors other than elliptical galaxies is expected to
be discovered with the upcoming deep sky surveys (see Saha et al. (in
preparation)). We can also gain important insights into the formation and
evolution of the deflectors by comparing their structural properties (e.g.,
the logarithmic slope of the density profile, the dark matter fraction, etc.)
across cosmic times. Furthermore, in this review article, we present
cosmological applications of the galaxy-scale lenses that do not require time
delay information – i.e., measuring cosmological parameters such as
$\Omega_{\mathrm{m}}$. Cosmological application involving the time delay
measurements, i.e., measuring primarily the Hubble constant, is presented in
detail in Birrer et al. (2022). Additionally, the application of galaxy-scale
lensing to study sub-galactic structures of dark matter is presented in
Vegetti et al. (in preparation).

In this introductory section, we provide a brief description of the lensing
phenomenology (Sect. 1.1) and discuss the advantages of strong lensing
observables in comparison with other probes of galaxy mass, such as stellar
dynamics (Sect. 1.2). The remainder of this review article is organized as
follows. In Sect. 2, we highlight the significant historical results involving
galaxy-scale lenses and introduce several prominent lens samples. In Sect. 3,
we describe the strong lensing observables and their modeling and analysis
methods. We discuss the application of galaxy-scale strong lensing to study
galaxy properties, formation, and evolution in Sect. 4 and to constrain
cosmological parameters in Sect. 5. Next, in Sect. 6, we discuss open issues –
both in technical aspects and scientific questions – and provide future
outlooks. We conclude the review article in Sect. 7.

###  1.1 Brief description of lensing phenomenology at the galaxy-scale

The background extended source in a galaxy-scale strong-lens system can be
lensed into arcs or a full Einstein ring. Multiple point images will also
appear if there is a point source within the background galaxy, e.g., an
active galactic nucleus (AGN) or a supernova. The former type of system is
called galaxy–galaxy lenses. In contrast, the latter is usually referred to as
“quads” (for the case of 4 detected point images) or “doubles” (for the case
of 2 detected point images). The different manifestations of strong lenses –
the appearance of arcs or a full Einstein ring, or the number of point images
– depend on the position of the source with respect to the lens caustics, as
introduced in Saha et al. (in preparation). Strong lensing of point sources
can provide three types of observables: image positions, image magnification
ratios, and time delays between the images. It follows from the lensing theory
that all three are properties of the Fermat potential, or the arrival-time
surface: images form at the local extrema – minima, saddle points, and maxima
– of this potential, magnification is inversely proportional to the
determinant of its Jacobian matrix, and time delays are the differences in its
height at the locations of the images (for a detailed explanation, see Saha et
al. (in preparation)).

However, not all of these observables are available for every lens. Image
positions are almost always observable for point-source lens systems. Image
flux ratios, while always observable, can be turned to actual magnifications
only when the intrinsic source flux is known. However, interpreting these
magnifications during modeling requires extra care, as they can be affected by
more complex features in the galaxy-scale mass distribution, e.g., baryonic
disks (Hsueh et al., 2016, 2017), microlensing by individual stars and planets
in the lensing galaxy (see Vernardos et al. (in preparation)), some
intermediate-mass-scale structures like dark matter subhalos (see Vegetti et
al. (in preparation)), or dust extinction by the lens galaxy (e.g., Motta et
al., 2002; Mediavilla et al., 2005). Time delays can be obtained only for
variable point sources, like quasars, supernovae, and in the future,
gravitational waves and fast radio bursts (see Birrer et al. (2022)).
Measuring time delays can be a difficult task, especially for quasars, which
require long-term monitoring spanning from a few seasons to years (e.g.,
Eigenbrod et al., 2005; Bonvin et al., 2016; Millon et al., 2020).

The above description of point-source lensing must be modified a little for
extended sources comparable in size to the lens caustics. Instead of point-
like image positions, the light from such a source is spread into an extended
area in the image plane and then further smeared by the point spread function
(PSF). As a result, it is hard to know a priori how the observed flux of the
multiple images of the source on the lens plane is traced back to the same
location on the source plane. This means that, even if the actual source
brightness is known, a lens model is required to compute the magnification
field (i.e., the Jacobian of the lens potential) across the lens plane. Due to
their size, extended sources are not variable on time scales relevant for
lensing. Thus, time delays are not observable without a variable source. Note
that extended arcs from the point-source host galaxy are also usually present
for a lens system including a point-like source. Thus, the lensed arcs’ flux
distribution and point-image positions can be utilized simultaneously to
constrain a lens model.

###  1.2 Unique advantages of lensing as a probe of galaxy structure

Other than lensing, the only commonly used probe of the mass distribution in
galaxies is kinematics, i.e., the velocity dispersion or streaming motions of
stars, gas, or globular clusters. For this to be fully informative of the mass
distribution, however, spatially resolved measurements are necessary
(Cappellari, 2016), which is often limited to nearby galaxies at $z\leq 0.5$
due to the sizeable observational cost, or due to the smaller size on the sky
of objects at higher redshifts. Using individual stars’ motions, e.g., from
Gaia, to map a galaxy’s mass is only limited to our Milky way (e.g., Nitschai
et al., 2020), thus an aperture-integrated velocity dispersion measurement is
the only possible dynamical observable for galaxies other than ours. However,
strong lensing observables are obtained from high-resolution imaging data,
which are much more informative than the dynamical observables for galaxies
beyond the local Universe ($z\gtrsim 0.03$, e.g., Smith et al., 2015). Strong
lensing can provide $\sim$1–2% constraints on the mass enclosed within the
Einstein radius from imaging data alone. To obtain a similarly precise mass
constraint from dynamical observations at high redshift, either integration
times longer by a factor $\mathcal{O}(10)$ are needed on current facilities,
or we must wait for better quality adaptive optics systems planned for future
extremely large telescopes.

Furthermore, the dynamical observables have their own intrinsic degeneracies,
i.e., the mass–anisotropy degeneracy (e.g., Treu and Koopmans, 2002). There is
an important complementarity between lensing and dynamical observables, which
can be used to break their corresponding degeneracies (e.g., Courteau et al.,
2014). Finally, gravitational deflection of light depends only on the
deflecting mass, so lensing responds to both baryonic and dark matter without
any assumptions on their dynamical state, i.e., regardless of whether it is in
equilibrium or not.

##  2 Historical background

This section provides a historical note on the initial discoveries of galaxy-
scale strong lenses (Sect. 2.1). Then, in Sect. 2.2, we briefly introduce
several prominent samples of galaxy-scale lenses that have contributed to the
major science applications described in Sect. 4 and Sect. 5. Additional
references can be found in the review by Treu (2010).

###  2.1 Initial discoveries of strong lensing systems

The first strong lens, the double quasar 0957$+$561A,B, was discovered in 1979
by Walsh et al. (1979). The system consists of two images of the quasar
separated by $5\,\hbox{$.\\!\\!^{\prime\prime}$}{7}$. The authors first
offered a “conventional” interpretation that the two images are different,
individual quasars that happen to be close to each other and share the same
physical characteristics. Since no gravitational lens was known before that,
their less conventional view was that the two are multiple images of the same
source. Subsequent work showed that the lensing hypothesis was correct. For
example, the structure of the radio jets emanating from the two images of the
quasar is consistent with them being mirror imaged, as one would expect for a
minimum and a saddle image (Gorenstein et al., 1988; Garrett et al., 1994).
The discovery of the first quadruply imaged quasar, PG 1115$+$080, was
announced the following year (Weymann et al., 1980). It was initially called a
“triple” because the second arriving minimum and its neighboring saddle point
were too close to be resolved. These two lenses’ discoveries opened up a new
field in astrophysics: observations of multiply-imaged, “strong” gravitational
lenses. Discoveries of other types of lensed sources followed: for example,
the first dust-obscured Seyfert 2 AGN, IRAS F10214$+$4724, was detected in the
infra-red, and later identified as being lensed (Eisenhardt et al., 1996;
Lehar and Broadhurst, 1996).

The increasing number of detections of such lensed systems spurred a lens
modeling effort. The initial studies that used more detailed models beyond a
point-like lens mass appeared in the early 1980s (Young et al., 1980, 1981b,
1981a), and already recognized “that there are several plausible ways to
reproduce the observations”, foreshadowing the importance of lens model
degeneracies. The role of the lens mass granularity due to individual stars in
the lensing galaxy was also recognized very early on (Chang and Refsdal, 1979;
Young, 1981) and later grew into the rich sub-field of extragalactic stellar
microlensing (see Vernardos et al. (in preparation)).

###  2.2 Prominent samples of galaxy-scale lenses

In this subsection, we briefly introduce some of the most prominent samples of
galaxy-scale lenses that had an impact on the science applications presented
in Sect. 4 and Sect. 5. Note that this is not a complete list of all the
discovered galaxy-scale lenses. The specifics of lens searching and discovery
to build samples like these are the topic of McMahon et al. (in preparation),
where we refer the reader for recent developments in the search methods and
discoveries.

####  2.2.1 MG-VLA survey-based samples

The first systematic search for strongly lensed systems at radio wavelengths
took place in the eighties within the MG-VLA survey (Lawrence et al., 1986).
This survey led to the discovery of a few famous radio-loud lensed quasars
among thousands of radio sources scrutinized at high resolution with the Very
Large Array (VLA). The Jodrell Bank–VLA Astromtric Survey (JVAS; Patnaik et
al., 1992; King et al., 1996), and its successor, the Cosmic Lens All Sky
Survey (CLASS; Myers et al., 1995) was the largest survey carried out for a
long time. This survey targeted the whole northern sky
($0^{\circ}<\textrm{Dec}<75^{\circ}$) for multiple images (separated by
$0\hbox{$.\\!\\!^{\prime\prime}$}{3}<\Delta\theta<6\hbox{$.\\!\\!^{\prime\prime}$}{0}$)
among flat spectrum radio sources brighter than 30 mJy. CLASS yielded the
discovery of 22 new systems, among which ten are doubles, nine are quads, and
one displays six images (Browne et al., 2003).

The CfA-Arizona Space Telescope LEns Survey
(CASTLES111https://lweb.cfa.harvard.edu/castles/) is a follow-up Hubble Space
Telescope (HST) imaging survey of $\sim$100 galaxy-scale lenses222There are a
handful cluster-scale lenses included in this sample as well. known at the
time, some from previous surveys such as CLASS (Muñoz et al., 1998). This
survey collected the first uniform ensemble of high-resolution images of known
galaxy-scale lens systems, including both galaxy–quasar and galaxy–galaxy
lenses. An account of early systematic lens searches, which generally unveiled
small samples of less than half a dozen systems, can be found in Claeskens and
Surdej (2002).

####  2.2.2 SDSS-based samples

The Sloan Lens ACS (SLACS) survey discovered 85 galaxy–galaxy lenses from the
Sloan Digital Sky Survey (SDSS) spectroscopic data and followed them up with
multi-band HST imaging (Bolton et al., 2006; Auger et al., 2009). The sample
was later expanded with 40 new systems having smaller deflector masses – the
SLACS for the Masses (S4TM) sample (Shu et al., 2015). In addition to
galaxy–galaxy lenses, the SDSS Quasar Lens Search (SQLS) discovered a sample
of 28 galaxy–quasar lenses using SDSS multicolor imaging data (Oguri et al.,
2006). More galaxy–quasar lens systems were discovered from joint SDSS and
UKIRT Infrared Deep Sky Survey (UKIDSS) data by the Major UKIDSS-SDSS Cosmic
Lens Survey (MUSCLES; Jackson et al., 2012).

####  2.2.3 CFHTLens-based samples

The Strong Lensing Legacy Survey (SL2S) discovered a sample of $\sim 35$
galaxy–galaxy lenses from the Canada-France-Hawaii Telescope Legacy Survey
data (CFHTLens, Gavazzi et al., 2012). Some newer lens samples have also been
discovered from the CFHTLens data (More et al., 2016a; Paraficz et al., 2016).

####  2.2.4 BOSS-based samples

The Baryon Oscillation Spectroscopic Survey Emission-Line Lens Survey (BELLS)
discovered $\sim$30 galaxy–galaxy systems from the Baryon Oscillation
Spectroscopic Survey (BOSS) and obtained HST imaging for them (Brownstein et
al., 2012). This survey was later expanded into the BELLS for the GALaxy-
Ly$\alpha$ EmitteR sYstems (BELLS GALLERY) survey, where the source galaxies
are specifically Ly$\alpha$ emitters (Shu et al., 2016). A sample of 13
strongly lensed quasars has also been discovered from the BOSS data (More et
al., 2016b).

###  2.3 Other samples and ongoing efforts from recent surveys

Several other lens samples have been recently discovered thanks to the
numerous large-area sky surveys. The STRong-Lensing Insights into the Dark
Energy Survey (STRIDES) collaboration has discovered $\sim$30 quadruply lensed
quasar systems from the Dark Energy Survey (DES) data – often in combination
with data from other large-area sky surveys – and obtained multi-band HST
imaging of them in infra-red (IR), optical, and ultra-violet (UV) bands
(Shajib et al., 2019; Schmidt et al., 2022). Many galaxy–galaxy lens
candidates have also been identified in the DES data (Jacobs et al., 2019,
2019; Rojas et al., 2021; Tran et al., 2022). In addition to the DES, surveys
such as the Canada-France Imaging Survey (CFIS), Gaia, the Hyper Suprime-Cam
(HSC) survey, the Kilo Degree Survey (KiDS), and the Panoramic Survey
Telescope and Rapid Response System (Pan-STARRS) have provided a plethora of
newly discovered galaxy-scale strong lenses (e.g., Petrillo et al., 2017,
2019; Agnello et al., 2018; Krone-Martins et al., 2018; Delchambre et al.,
2018; Lemon et al., 2018, 2020; Cañameras et al., 2021; Savary et al., 2021;
Li et al., 2021; Lemon et al., 2022; Wong et al., 2022). Most of these samples
still contain candidate lenses and require spectroscopic confirmation (by
measuring the redshifts) and high-resolution imaging (to perform lens
modeling) for the science applications described in Sect. 4 and Sect. 5. See
McMahon et al. (in preparation)and the references therein for a detailed
discussion on the ongoing efforts.

Figure 1:  Examples of galaxy-scale lenses with different types of background
sources. These images are false color images created from multi-band HST
imaging. The first two panels show two SLACS systems, where the background
source is a star-forming galaxy without a central AGN. Image credit: NASA/ESA,
A. Bolton, and the SLACS team. The next two panels show quadruply lensed
quasars discovered by the STRIDES collaboration (Shajib et al., 2019). The
fifth panel shows a quadruply imaged supernova (Dhawan et al., 2020). For the
quasar and supernova examples, lensed arcs or a complete Einstein ring from
their extended host galaxy are also noticeable.

##  3 Observables and analysis methods

This section describes the strong lensing observables (Sect. 3.1) and the
analysis techniques to constrain galaxy properties from them. We describe the
lens modeling methods in Sect. 3.2 and briefly discuss their common
degeneracies. We outline the Bayesian hierarchical framework for a sample of
lens systems, which allows inferring population characteristics of galaxies
from the sample in hand in Sect. 3.4. Lastly, additional non-lensing
observables most commonly combined with strong lensing observables are
presented in Sect. 3.5.

###  3.1 Lensing observables

The two types of lensing observables for galaxy-scale lenses are imaging of
the lensing system and the time delay between different point-source images.

####  3.1.1 Imaging of the lens system

The most common and informative lensing observables result from imaging data
with angular resolution much better than the Einstein radius. An extended
source can be lensed into clearly identifiable arcs that can form partial or
complete Einstein rings (see Fig. 1). The conjugate points on the arcs, i.e.,
locations that are traced back to the same location in the source plane, can
be used to simultaneously constrain a lens model (through its deflection
angles and magnification) and reconstruct the surface brightness of the source
that is a priori unknown. Although it is possible to constrain simple models
from any image that displays lensing features with sufficient signal-to-noise
ratio, high-resolution imaging from space- or ground-based telescopes offers
many more observational constraints (i.e., conjugate pixels). This is crucial
for exploring more sophisticated models, which are required for precise
science applications (see Sect. 4 and Sect. 5).

If the background source contains a point-like emitting region – e.g., a
quasar or supernova – the positions of its multiple images (i.e., conjugate
points) can be extracted from the data and used as constraints for a lens
model. Although the simple addition of 2 or 4 conjugate pixels may initially
seem insignificant compared to the hundreds of pixels that correspond to the
extended source, the flux contained in these pixels may be brighter than the
entire host galaxy, e.g., in the case of an AGN. Hence, these conjugate points
have a much stronger effect on the resulting model and must be treated
separately. In addition to the astrometry of the point-like images, their
magnification ratio (or flux ratio) can be used as an observational
constraint. However, image magnifications are susceptible to micro- and milli-
lensing, dust extinction333Lensed arcs from extended sources are also
vulnerable to dust extinction., and the effect of higher order moments in the
mass distribution of the lens usually attributed to the complex, non-linear
physics of baryons, e.g., galactic disks, bars (Hsueh et al., 2016, 2017).
Therefore, many lensing analyses involving systems with such point-like
sources choose not to use magnification ratios as constraints (e.g., Shajib et
al., 2019). However, suppose the microlensing and dust extinction effects can
be incorporated and quantified within a lens model. In that case, any residual
flux ratio anomaly would signal a departure from the smooth macro-model for
the deflector galaxy and thus could indicate the presence of dark subhalos
within the main halo (e.g., Mao and Schneider, 1998; Nierenberg et al., 2017;
Gilman et al., 2020). Such detection of dark subhalos can provide important
insights into the nature of the dark matter properties, as presented in
Vegetti et al. (in preparation).

####  3.1.2 Time delay

If the background source is a variable point source, e.g., a quasar or
supernova, the delay between the arrival times of photons at its multiple
different images can be measured through long-term monitoring that spans from
a few seasons to decades (e.g., Eigenbrod et al., 2005; Bonvin et al., 2016;
Millon et al., 2020). The time delays are variant under the well-known mass-
sheet degeneracy (MSD) in lensing, unlike the imaging observables (Falco et
al., 1985; Schneider and Sluse, 2014). Thus, they can be combined with the
imaging observables to break the MSD when constraining the potential of the
lensing galaxy. However, such a combination of these observables requires a
fiducial cosmological model since the time delays depend on the cosmology,
particularly the Hubble constant. A more detailed discussion on the
measurement of time delays is provided in Birrer et al. (2022).

###  3.2 Lens modeling methods

Lens modeling is the process of constraining properties of the lens galaxy and
the source from the lensing observables. Traditional methods are based on
reconstructing the source light and lens potential to fit the data under some
assumptions, such as regularization, and are described in Sect. 3.2.1. More
recently, machine learning methods have been developed for these purposes,
presented in Sect. 3.2.2, which provide faster models with data-driven instead
of ad-hoc assumptions.

####  3.2.1 Likelihood-based inference

Such methods require a likelihood function to be maximized that leads to an
optimized model, able to reproduce an observed multiply-imaged system down to
the noise level. In general, the model consists of three main components: the
background source’s flux distribution, the mass distribution in the lensing
galaxy (or galaxies), and the flux distribution in the lensing galaxy (or
galaxies). The likelihood function that measures the goodness-of-fit of the
model to the data can be defined as

$\mathcal{L}(\bm{d}\mid\bm{m})\propto\exp\left[-\frac{1}{2}\left(\bm{d}-\bm{m}\right)^{T}\Sigma_{\bm{d}}^{-1}\left(\bm{d}-\bm{m}\right)\right],$
(1)

where $\bm{d}$ is the vector of pixel values in the data,

$\bm{m}\equiv\bm{m}(\bm{\xi}_{\mathrm{mass}},\ \bm{\xi}_{\mathrm{source}},\
\bm{\xi}_{\mathrm{light}})$ (2)

is the model’s estimate of the flux in the data pixels, where $\bm{\xi}_{\rm
mass}$ is the set of model parameters defining the lens mass distribution,
$\bm{\xi}_{\rm source}$ is the set of model parameters defining the source’s
flux distribution, $\bm{\xi}_{\rm light}$ is the set of model parameters
defining the lens galaxy’s flux distribution, and $\Sigma_{\bm{d}}$ is the
covariance matrix of the data. Some studies choose to subtract the lens
galaxy’s flux distribution from the imaging data before optimizing a lens
model without it (e.g., Bolton et al., 2008). Given that the lens flux is
loosely related to the lensing phenomenon – only through mass-follows-light
arguments that are not strictly required – we omit it in the following for
brevity. This likelihood function can be extended in a Bayesian framework to
include prior (regularization) terms on the source (Warren and Dye, 2003; Treu
and Koopmans, 2004; Koopmans, 2005), or used to compute the Bayesian evidence
(Suyu et al., 2006; Vernardos and Koopmans, 2022) and perform model
comparisons.

When the lensed arcs from an extended source are resolved, each pixel is a
constraint for the lens model. Given a deflection field,
$\bm{\alpha}(\bm{\theta})$, which depends on the lens mass distribution
through $\bm{\xi}_{\mathrm{mass}}$, the lens equation

$\bm{\beta}(\bm{\theta})=\bm{\theta}-\bm{\alpha}(\bm{\theta})$ (3)

can be used to map any position on the image plane, $\bm{\theta}$, to the
corresponding one on the source plane, $\bm{\beta}$. The lensed flux at any
location on the image plane can then be easily calculated as

$I(\bm{\theta})=S[\bm{\beta}(\bm{\theta})],$ (4)

where $S$ is the light distribution of the source that depends on
$\bm{\xi}_{\rm source}$, and we use the fact that lensing conserves surface
brightness. The dependence of $I(\bm{\theta})$ on the mass is almost always
non-linear, which means that we cannot directly solve Eq. 3 to obtain the true
parameters $\bm{\xi}_{\rm mass}$, even when $S$ (and therefore
$\bm{\xi}_{\mathrm{mass}}$) is perfectly known.

In practice, $S$ is an unknown that must be solved simultaneously with the
lens potential and requires special attention. The most straightforward choice
for it is an analytic function, e.g., a Sérsic profile (Sersic, 1968), whose
parameters are treated in the same forward-modeling way as for the lens
potential. However, several more advanced techniques have been developed over
the years that allow a free-form reconstruction of the source, each with its
advantages and disadvantages. The semi-linear inversion technique of Warren
and Dye (2003) uses a regular grid of pixels to approximate the source
brightness. Although the degrees of freedom are much higher in this case, the
use of priors explicitly selected to correspond to quadratic terms in the
likelihood allows the direct computation of the best-fit source for given
$\bm{\xi}_{\mathrm{mass}}$ (hence the name ‘inverse’ as opposed to ‘forward’
modeling). This turns the source reconstruction into a linear inversion
problem for a given set of non-linear parameters within $\\{\bm{\xi}_{\rm
mass},\bm{\xi}_{\rm source}\\}$, where $\bm{\xi}_{\rm source}$ refers only to
the ‘hyper-parameters’ of the source that are not solved through the linear
inversion.

This linear inversion problem requires some form of regularization of the
source (i.e., prior) to be solved. Typical choices of source regularization
are gradient and curvature that impose smoothness on the source solution
through its derivatives (a standard approach in this kind of problem, e.g.
Press et al., 1992). Alternative choices are covariance kernel priors
(Vernardos and Koopmans, 2022), which are based on observations of the galaxy
brightness power spectrum and are thus a more physically justifiable choice,
or multi-scale regularization through the use of sparsity constraints on a
wavelet representation of the source (Galan et al., 2021). In fact, choosing a
different basis to represent the source can in itself significantly reduce the
number of degrees of freedom while still having enough flexibility to
represent complex light profiles across different scales (Birrer et al., 2015;
Tagore and Jackson, 2016, see also Fig. 3). This can be similarly achieved by
reconstructing the source on an irregular, adaptive grid that can have
increased resolution in the most magnified areas of the source, i.e., near the
caustics (Vegetti and Koopmans, 2009; Nightingale and Dye, 2015; Vernardos and
Koopmans, 2022). Fig. 2 illustrates an example of lens modeling based on high-
resolution HST imaging of a lensing system with a bright, extended source
component. The source is reconstructed on an adaptive grid using curvature
regularization (see Chapter 4 of Bayer, 2021).

If there is a point source within the extended source galaxy, then the
positions of its point-like images provide additional constraints on the mass
model as conjugate points. One approach to model the point-like source is to
include the location of its multiple images as free parameters and then
require the mass model to trace them back to the same location on the source
plane through the lens equation. Inversely, the point-source location can be
free, and the predicted image positions must match the observed ones. In
practice, the choice between the two approaches is based on computational and
sub-grid effect arguments (see Keeton, 2001, for a more detailed presentation
on the topic). Fig. 3 illustrates an example of lens modeling based on high-
resolution HST imaging of a lensing system with both an extended and a point-
like source component. In this case, the model for the extended source
brightness (the quasar host galaxy) assumes an analytic Sérsic profile and an
additional free-form component described by a basis set of shapelets (Birrer
et al., 2015; Shajib et al., 2019).

Figure 2:  First panel: False-color RGB image combining three HST filters
(Sonnenfeld et al., 2012). Second panel: HST image in the F814W filter of the
system with the galaxy light subtracted. The lensed arcs from the source
galaxies at different redshifts are indicated in solid and dashed contours.
Third panel: Model of the lensed arc from only the brightest source. Fourth
panel: Corresponding source reconstructed on an adaptive grid using every data
pixel and curvature regularization [see Chapter 4 of Bayer (2021) for the
analysis of the data and Vernardos and Koopmans (2022) for the modeling
method]. Figure 3:  Example of modeling a galaxy-scale lens system - here the
lensed quasar DES J2038$-$4008\. First panel: The false-color image of the
system from the DES data, where only the point images are resolved (Agnello et
al., 2018). Second panel: A lens model based only on the positions of the
quasar images (Agnello et al., 2018). The dashed lines show the contours of
the time-delay surface, whose extrema define the locations where images of the
point source appear. Third panel: The false-color image of the system taken by
the Hubble Space Telescope (HST), where the lensed arcs of the quasar host
galaxy can be seen in great detail. Fourth panel: The best-fit model of the
HST data from Shajib et al. (2019). Fifth panel: The corresponding
reconstructed quasar host galaxy (without the central quasar) in the F160W
band. The blue star marks the position of the central quasar.

####  3.2.2 Machine-learning based parameter extraction

Recently, machine learning (ML) based methods have been developed to extract
lensing parameters – e.g., the Einstein radius, power-law index, ellipticity,
and shear parameters – from the imaging data (e.g., Hezaveh et al., 2017;
Morningstar et al., 2019; Adam et al., 2022; Schuldt et al., 2022). Some
studies also explored reconstructing the lens mass and source flux
distributions using machine learning algorithms (e.g., Chianese et al., 2020;
Karchev et al., 2022; Mishra-Sharma and Yang, 2022; Biggio et al., 2022). In
these approaches, a machine learning algorithm – usually a neural network – is
trained using synthetic data since real examples of lens systems are not
adequate in number for such data-intensive training. This leads to the
critical caveat that the mass density profile in the simulated galaxies is
based on empirical priors (e.g., Hezaveh et al., 2017) or cosmological
simulations (e.g., Adam et al., 2022). Thus, the inferred parameters are
prior-dependent – either empirical or physical – in a similar manner that the
forward modeling approach depends on the adopted mass model and the associated
priors.

Although ML-based inferences are yet to demonstrate their potential in science
applications with real data, they are essential for quickly extracting lensing
parameters with minimum human supervision. This is crucial for modeling huge
samples of lenses like those expected from several large surveys upcoming in
the 2020s (e.g., Park et al., 2021; Wagner-Carena et al., 2021). For such
large datasets, the traditional forward modeling techniques would be
unfeasible due to either computational or human time restrictions. Finally,
ML-based inference provides a direct and fast way to constrain quantities of
interest that can be otherwise too cumbersome to infer from the data through a
traditional approach, e.g., detecting individual or populations of dark
subhalos or constraining the subhalo-mass function (Brehmer et al., 2019; Diaz
Rivero and Dvorkin, 2020; Coogan et al., 2020, 2022; Ostdiek et al., 2020,
2022; Vernardos et al., 2020; Wagner-Carena et al., 2021; Zhang et al., 2022;
Anau Montel et al., 2022). Note that this particular science application is
discussed in detail in Vegetti et al. (in preparation).

###  3.3 Lens mass models

Lens mass models can vary in complexity and eventually in the number of free
parameters, reflecting that the mass distribution in galaxies is a non-trivial
problem for which methodologies and algorithms are still evolving. Although it
is possible to extract lensing information with free-form models and directly
connect them to galaxy properties of interest (see Saha et al. (in
preparation)), lens modeling with simple parametrization has been the most
common practice in the literature (e.g., Ritondale et al., 2019; Schmidt et
al., 2022). We briefly introduce some commonly used models with simple
parametrization in Sect. 3.3.1, discuss some degeneracies that impact them in
Sect. 3.3.2, and describe free-form models in Sect. 3.3.3 and model ensembles
in Sect. 3.3.4.

####  3.3.1 Simply parametrized mass models

Saha et al. (in preparation) provides an in-depth discussion about several
popular mass models with simple parametrization, e.g., the singular isothermal
sphere/ellipsoid (SIS/SIE), the softened power-law mass distribution (SPEMD;
Kassiola and Kovner, 1993; Barkana, 1998), and the Navarro–Frenk–White (NFW)
profile (Navarro et al., 1996, 1997). It is also often necessary to include a
constant external shear (XS) lensing potential in addition to the primary lens
mass distribution to accurately model the distortions in the lensed arcs or
the Einstein ring arising from nearby perturbers and large-scale structures
(Keeton et al., 1997). Until the mid-2010s, the SIE+XS model has been the most
popular choice to model large samples of galaxy-scale lenses (e.g., Bolton et
al., 2008; Sonnenfeld et al., 2013), adequate for the required science
applications at the time. Although simple SIE models can sufficiently
constrain the Einstein radius, obtaining other important properties, such as
the radial slope of the mass profile, requires models with additional degrees
of freedom. Improved data quality and analysis techniques have allowed the use
of such models (e.g., Ritondale et al., 2019; Shajib et al., 2021).

Simply parametrized lens models with larger degrees of freedom, e.g., a
superposition of a stellar component with a constant/varying mass-to-light
ratio and an NFW profile for the dark matter, have been adopted by some
studies (e.g., Treu et al., 2010; Sonnenfeld et al., 2015; Shajib et al.,
2021). Although necessary for the addressed science questions, such mass
profiles with more free parameters amplify the impact of degeneracies inherent
to lensing. It is thus often necessary to constrain these additional
parameters by incorporating non-strong-lensing observables such as stellar
kinematics and weak lensing, or by incorporating informative priors (e.g.,
Sonnenfeld, 2018; Shajib et al., 2021).

While the additional degrees of freedom beyond the simple SIE model focus on
the radial profile of the lens in most cases, azimuthal structures such as
discyness, boxyness, ellipticity gradients, or isodensity twists may leave
noticeable imprints in the lensed images. For point-like sources, the image
flux ratios are the most susceptible to perturbations (Möller et al., 2003;
Keeton et al., 2003, 2005), while for extended sources the imprint is
generally more subtle and detectable only from high-resolution imaging data
(Van de Vyvere et al., 2022b, a). Fortunately, most of those structures arise
from baryonic physics and may also be detectable in the luminosity profile of
the lens (Van de Vyvere et al., 2022b).

####  3.3.2 Common degeneracies in simple parametric modeling

The degeneracies in lens modeling – both intrinsic in the data and stemming
from the parametrization scheme – are discussed in detail in Saha et al. (in
preparation). Here, we summarize the common degeneracies that largely impact
simply parametrized lens models for the convenience of the readers.

The MSD, which is intrinsic to imaging observables in lensing, originates from
the mass-sheet transform (MST, Falco et al., 1985; Saha, 2000)

$\kappa\to\kappa^{\prime}=\lambda\kappa+1-\lambda,$ (5)

where $\lambda$ is a constant. This equation implies a source position
transformation (e.g., Schneider and Sluse, 2014), where the unknown source
position is altered as

$\beta\to\beta^{\prime}=\lambda\beta.$ (6)

A more general but approximate degeneracy happens when $\lambda$ is not a
constant anymore but depends on the position, i.e.,
$\lambda\equiv\lambda(\theta)$ (Unruh et al., 2017; Wertz et al., 2018).

A particular choice of the mass profile, e.g., the power-law mass profile,
artificially limits the MSD, as the MST of a power law is mathematically not a
power law anymore. However, it is a common practice in the literature to
assume such a mass profile in galaxy-scale lens modeling, since numerous non-
lensing constraints on the mass profile of elliptical galaxies demonstrate
that this is a good approximation (e.g., Thomas et al., 2007; Tortora et al.,
2014; Bellstedt et al., 2018). The strong lensing imaging observables can only
constrain in a model-independent way $\theta_{\mathrm{E}}$ and the MST-
invariant quantity

$\xi_{\rm
rad}\propto\frac{\theta_{\mathrm{E}}\alpha^{\prime\prime}({\theta_{\mathrm{E}}})}{1-\kappa({\theta_{\mathrm{E}}})},$
(7)

where $\alpha^{\prime\prime}$ is the second derivative of the deflection angle
[see Equation 42 of Birrer (2021) for the full definition of $\xi_{\rm rad}$,
also Kochanek (2021)]. For a power-law convergence profile
$\kappa(\theta)\propto\theta^{-\gamma_{\rm PL}+1}$, this quantity becomes
$\xi_{\rm rad}=\gamma_{\rm PL}-2$. As a result, the choice of a SIE+XS lens
model fixes $\xi_{\rm rad}=0$ and only extracts $\theta_{\mathrm{E}}$ from the
imaging data.

Degeneracies in lens modeling can also arise from particular parametrizations
of the lens model. One such example is the shear–ellipticity degeneracy, as
the total shear can be redistributed between the “internal” shear, arising
from the ellipticity of the central deflector, and the external shear
(Kassiola and Kovner, 1993). Specifically, a SIE+XS model can be modeled with
only an SIE model that has a quadrupole moment equaling 1/3rd of the shear
(An, 2005). This degeneracy is particularly apparent when modeling with point-
like image positions as the only constraints (Witt, 1996).

####  3.3.3 Free-form models

In the most commonly used lens models, the main focus is on capturing the
radial shape of the mass profile, while any azimuthal structure beyond an
elliptical shape with external shear is of secondary importance. Real galaxies
can obviously have more complicated mass profiles, with higher order moments
present, such as discyness, boxyness, or bar- or disc-like structures (e.g.
Trotter et al., 2000; Claeskens et al., 2006; Hsueh et al., 2017; Frigo et
al., 2019), radial dependence of the ellipticity or the orientation of the
isodensity contours (gradients and twists, e.g., Hao et al., 2006; Oldham and
Auger, 2018; Sonnenfeld et al., 2018; Gomer and Williams, 2018; Nightingale et
al., 2019; Shajib et al., 2021), or even features that do not fit into a
simple parametric description, such as merger products that are not yet
completely relaxed and populations of substructures (e.g., satellite subhalos
or perturbers along the line of sight). Detecting such deviations from the
simple parametric profiles depends on numerous factors, such as their
alignment with the smooth potential (Van de Vyvere et al., 2022b), the
complexity of the brightness profile of the source (Vernardos and Koopmans,
2022), the signal-to-noise ratio, etc. However, multipole components in the
mass potential beyond the combined effect of ellipticity and external shear
have been recently detected in very-long-base-line interferometric
observations of a strong lens (Powell et al., 2022). Although not accounting
for such structures can bias the mass model by up to several percent – which
is mostly acceptable except for time-delay cosmography (see Birrer et al.
(2022)) and for detecting dark matter substructure (see Vegetti et al. (in
preparation)) – their detection holds valuable information on the formation
history and evolution of galaxies.

To this extent, free-form techniques have been developed that either entirely
dismiss any parametric mass component and employ a grid of mass pixels to
describe the lens potential (Saha and Williams, 1997), or retain a parametric
model as a first-order smooth component and combine it with a similar pixel
grid that now focuses specifically on capturing higher order deviations
(Koopmans, 2005). In both cases, regularization (prior) assumptions on the
free-form pixel grid are necessary to obtain a solution and to prevent the
appearance of unphysical mass distributions. Existing techniques are based on
forward models or extensions of the semi-linear inversion technique (Vernardos
and Koopmans, 2022).

The specific form of the regularization priors plays an important role in the
quality of the obtained solutions. It has been shown that purely
mathematically motivated priors (e.g., curvature) can lead to biased
potentials as opposed to more physically driven ones, based on the observed
light properties of real galaxies (e.g. Vernardos and Koopmans, 2022).
Recently, Galan et al. (2022) have proposed a wavelet-based regularization
technique that finds solutions that satisfy sparsity constraints. Baggio et
al. (in preparation) completely replaced the pixel grid with a neural network.
These approaches allow the use of purely data-driven regularization that is
the most compatible with the data (without computing and comparing the
Bayesian evidence). It remains to be seen how well these new and promising
techniques can perform on real data and robustly recover deviations from
smooth, parametric models that encode galaxy evolution.

####  3.3.4 Model ensembles

Whereas the goal of most modeling done at present in the literature, including
in free-form models as discussed above, is to find the single best-fitting
mass model with associated statistical uncertainties, an alternative approach
aims to find an ensemble of models, which explores the degenerate “model-
space” consistent with the data. In this approach, there is no single optimal
model and no need for regularization. However, a prior is still necessary, but
it is less informative.

Such model ensembles were first introduced for free-form models made up of
mass tiles or pixels (Saha and Williams, 2004; Coles et al., 2014). In these
and related works, the mass distribution is required to be non-negative and
centrally concentrated in a broad sense, and given these hard constraints, the
prior probability is uniform in an abstract space of pixels. From this prior,
models that correctly reproduce the observed positions of point-like image
features are sampled in a random-walk fashion to form an ensemble. The
ensemble can be further filtered according to how well image data can be
fitted. For an example, see Fig. 4 adopted from Denzel et al. 2021b).

Free-form modeling is not, however, the only way to generate lens model
ensembles; parametric models, with the number of parameters somewhat exceeding
the number of data constraints, have also been used recently (Williams and
Zegeye, 2020; Barrera et al., 2021).

Figure 4:  An example of a model ensemble for the lensing system SW05
J1434+5228, showing the cumulative mass profiles from the lens model (purple)
and the MassiveFIRE simulation (A1, A4, C1, D7; Feldmann et al., 2016, 2017).
The solid lines show the total enclosed mass, and the dashed lines show the
enclosed stellar mass. The shaded regions encompass 99.7% confidence region.
The models are best constrained at the image positions marked by vertical
dashed lines.

###  3.4 Bayesian hierarchical framework

The Bayesian hierarchical framework can be used to constrain the population
properties from a sample of strong lenses (e.g., Sonnenfeld et al., 2015).
This framework also allows one to incorporate a selection function of the
lensing galaxies and generalize the sample properties to the population of all
galaxies that are of the same type as the lensing ones (e.g., Sonnenfeld et
al., 2019). Within the hierarchical analysis, there are two levels of
parameters: hyper-parameters that dictate the distribution of the parent
population of the lens galaxies and parameters pertaining to individual lens
galaxies sampled from the parent population. The hierarchical framework
connects the “top”-level population hyper-parameters to the observed data of
individual galaxies through the “bottom”-level parameters. According to the
Bayes’ theorem, the posterior probability distribution of the hyper-parameters
$v$ is given by

$p(v\mid\mathcal{D})\propto p(\mathcal{D}\mid v)p(v),$ (8)

where $\mathcal{D}$ is the dataset, $p(\mathcal{D}\mid v)$ is the likelihood,
and $p(v)$ is the prior. For model parameters $w_{i}$ pertaining to individual
galaxies, the above equation can be expressed as

$p(v\mid\mathcal{D})\propto
p(v)\prod_{i}\int\mathrm{d}w_{i}p(\mathcal{D}_{i}\mid w_{i})p(w_{i}\mid v),$
(9)

where $\mathcal{D}_{i}$ is the data from to the $i$-th individual lens galaxy.
In this expression, the likelihood function $p(\mathcal{D}_{i}\mid w_{i})$ for
each galaxy can be more easily computed than the top-level likelihood function
$p(\mathcal{D}\mid v)$. This approach can infer any property at the population
level with the associated mean and scatter values. See, for example,
Sonnenfeld and Cautun (2021) for a detailed presentation of the hierarchical
framework with specific examples of hyper-parameters $v$.

###  3.5 Incorporating non-strong-lensing observables

Incorporating non-strong-lensing observables can be used to break the
degeneracies in strong lensing analysis. For example, stellar kinematics data
can be used to break the MSD (e.g., Romanowsky and Kochanek, 1999; Treu and
Koopmans, 2002, 2004; Shajib et al., 2018; Birrer et al., 2020; Shajib et al.,
2021), whereas spectroscopic stellar population analysis (e.g., Spiniello et
al., 2011), weak lensing (e.g., Gavazzi et al., 2007; Sonnenfeld et al., 2018;
Shajib et al., 2021), and microlensing (e.g., Schechter et al., 2014; Oguri et
al., 2014) information can help mitigate the degeneracy between the stellar
and dark matter distributions.

####  3.5.1 Combining stellar kinematics with strong lensing

Imaging observables probe the mass distribution of the lens projected on the
plane of the sky, whereas stellar kinematics probe its full three-dimensional
mass distribution. Thus, a combination of the two helps break the MSD to
constrain the mass distribution in galaxies robustly. Although elliptical
galaxies are triaxial, assuming spherical symmetry has been a standard
practice for the case of a single aperture-integrated stellar velocity
dispersion measurements (see Sonnenfeld et al., 2012, for a discussion on the
impact of this assumption). Then, the stellar velocity dispersion is obtained
by solving the spherical Jeans equation

$\frac{{\rm d}\left(l(r)\ \sigma_{\rm r}^{2}\right)}{{\rm
d}r}+\frac{2\beta_{\rm ani}(r)\ l(r)\ \sigma_{\rm r}^{2}}{r}=-l(r)\ \frac{{\rm
d}\Phi}{{\rm d}r}.$ (10)

Here, $l(r)$ is the three-dimensional luminosity density of the stars,
$\sigma_{\rm r}$ is the intrinsic radial velocity dispersion, and $\beta_{\rm
ani}(r)$ is the anisotropy parameter relating $\sigma_{\rm r}$ with the
tangential velocity dispersion $\sigma_{\rm t}$ given by

$\beta_{\rm ani}(r)\equiv 1-\frac{\sigma_{\rm t}^{2}}{\sigma_{\rm r}^{2}}.$
(11)

By solving the Jeans equation, the line-of-sight velocity dispersion, which is
the kinematic observable, is obtained as

$\sigma_{\rm
los}^{2}(R)=\frac{2G}{I(R)}\int_{R}^{\infty}\mathcal{K}_{\beta}\left(\frac{r}{R}\right)\frac{l(r)\
M(r)}{r}\ {\rm d}r$ (12)

(Mamon and Łokas, 2005). Here, $M(r)$ is the enclosed mass within radius $r$.
The function $\mathcal{K}_{\beta}(u)$ depends on the parameterization of
$\beta(r)$ [see Mamon and Łokas (2005) for specific forms of
$\mathcal{K}_{\beta}(u)$ corresponding to different $\beta(r)$]. Thus, the
observed velocity dispersion in Eq. 12 can be written as a function of the
lens model parameters as

$\sigma_{\rm los}^{2}=\frac{D_{\mathrm{s}}}{D_{\mathrm{ds}}}\ c^{2}\
J(\bm{\xi}_{\rm mass},\ \bm{\xi}_{\rm light},\ \beta_{\rm ani},\ \lambda),$
(13)

where $\bm{\xi}_{\rm mass}$ are the deflector’s mass model parameters,
$\bm{\xi}_{\rm light}$ are the deflector’s light model parameters. In this
form, the function $J$ is independent of cosmology. It only depends on the
lens model parameters $\\{\bm{\xi}_{\rm mass},\ \bm{\xi}_{\rm light}\\}$, the
anisotropy profile $\beta_{\rm ani}(r)$, and the MST parameter $\lambda$. All
the cosmological dependence of $\sigma_{\rm los}$ is contained in the distance
ratio $D_{\mathrm{s}}/D_{\mathrm{ds}}$ (Birrer et al., 2016). To reproduce the
observed velocity dispersion integrated within an aperture, the computed
luminosity-weighted velocity dispersion needs to be blurred with the seeing
$\mathcal{S}$ by performing the convolution

$\sigma^{2}_{\rm ap}=\frac{\int_{\rm ap}\left[I(R)\sigma^{2}_{\rm
los}(R)*\mathcal{S}\right]\ R\ \mathrm{d}R\mathrm{d}\theta}{\int_{\rm
ap}\left[I(R)*\mathcal{S}\right]\ R\ \mathrm{d}R\mathrm{d}\theta},$ (14)

where the $*$ symbol denotes the convolution operation. To obtain analytic
solutions with specific choices for mass, light, and anisotropy profiles, see
Koopmans (2004) in the case of power-law mass and light profiles with constant
anisotropy, and Suyu et al. (2010) in the case with power-law mass profile,
the Hernquist light profile (Hernquist, 1990), and isotropic stellar orbits.

The constraints from the velocity dispersion measurement can be folded in the
lens model posterior with a multiplicative likelihood term

$\mathcal{L}_{\rm kin}\propto\exp\left[-\frac{(\sigma_{\rm ap}^{\rm
obs}-\sigma_{\rm ap}^{\rm model})^{2}}{2\sigma^{2}_{\sigma_{\rm ap}^{\rm
obs}}}\right],$ (15)

where $\sigma_{\sigma_{\rm ap}^{\rm obs}}$ is the uncertainty in the observed
velocity dispersion $\sigma_{\rm ap}^{\rm obs}$. Whereas the imaging
observables from strong lensing cannot constrain the MST parameter $\lambda$,
the stellar velocity dispersion constrains $\lambda$ through the likelihood
$\mathcal{L}_{\rm kin}$ (see Birrer et al. (2022)for a detailed discussion
within a hierarchical Bayesian framework).

Although integrated velocity dispersion from single slit spectra is the most
commonly used kinematic observable in strong lensing studies, a few studies
have also incorporated spatially resolved velocity dispersion, mainly from IFU
spectra (e.g., Barnabè and Koopmans, 2007; Barnabè et al., 2011; Czoske et
al., 2012; Spiniello et al., 2015).

####  3.5.2 Combining weak lensing with strong lensing

Weak lensing measures the excess shear quantity $\Delta\Sigma$ from tidal
distortions of galaxies far away ($\gtrsim 10$ arcsec) from the central
lensing galaxy. Thus, weak lensing provides information on the mass
distribution of the lensing galaxy’s outer region, i.e., where the dark matter
halo dominates. As a result, the weak lensing information constrains the
normalization of a dark matter component in the lens mass model. Strong
lensing observables provide information on the total mass distribution (e.g.,
the enclosed mass). Still, there remains a degeneracy between the luminous and
dark matter fractions within the total enclosed mass. If a specific profile is
assumed for the dark matter distribution, weak lensing data can break this
degeneracy between luminous and dark component normalizations (Sonnenfeld,
2018; Shajib et al., 2021).

The weak lensing information needs to be incorporated within the hierarchical
framework. Since the weak lensing signal from one single lensing galaxy does
not usually have enough constraining power on the dark matter normalization,
it is often required to stack weak lensing signals from a large sample of
elliptical galaxies. If it is justified to assume that this large sample of
elliptical galaxies and the lens galaxy sample under consideration are
subsamples of the same parent population, then the weak lensing observables
$\mathcal{D}_{\rm weak}$ and the strong lensing observables $\mathcal{D}_{\rm
strong}$ can be jointly considered in the posterior probability function of
the population hyper-parameters $v$ as

$p(v\mid\mathcal{D}_{\rm strong},\mathcal{D}_{\rm weak})\propto
p(\mathcal{D}_{\rm strong}\mid v)p(\mathcal{D}_{\rm weak}\mid v)p(v).$ (16)

Here, the strong lensing likelihood $p(\mathcal{D}_{\rm strong}\mid v)$ can
then be expanded using the single system likelihoods similar to the product
term in Eq. 9. Since from Bayes’ theorem, we have

$p(\mathcal{D}_{\rm weak}\mid v)p(v)\propto p(v\mid\mathcal{D}_{\rm weak}),$
(17)

it is also valid for numerical convenience to first obtain the posterior of
$v$ from weak lensing observables only and then fold this posterior as the
prior of $v$ in the hierarchical analysis of strong lensing data only.

##  4 Applications in galaxy properties and evolution

This section describes what we can learn about galaxy structure and evolution
using the lensing galaxy properties. Since strong lensing galaxies are
typically massive ellipticals, most of the strong-lensing studies in the field
relate to this type of galaxies (Sect. 4.1–Sect. 4.4). However, we briefly
discuss strong lensing by spiral galaxies at the end of this section in Sect.
4.5.

###  4.1 Galaxy mass density profile

All galaxies are believed to form and grow inside their dark matter halos.
Thus, a massive galaxy’s total mass density profile comprises two components:
the baryonic matter distribution, which includes stars and gas, and the dark
matter distribution. As seen from the comparison of the observed galaxy
luminosity function and the distribution of simulated dark-matter halos (i.e.,
abundance matching, see Moster et al., 2010), the stellar-mass fraction
decreases with mass, for the mass range of typical lensing galaxies. This is
attributed to AGN feedback and is evident from simple comparisons of the
stellar and total mass in lensing galaxies (Auger et al., 2009; Küng et al.,
2018).

It is observed from several lensing and non-lensing observables that the total
density profile of massive galaxies is well approximated by a power law
$\rho(r)\propto r^{-\gamma_{\mathrm{pl}}}$ that is close to the isothermal
profile, i.e., $\gamma_{\mathrm{pl}}=2$, e.g., from strong lensing only or in
combination with stellar dynamics (Kochanek, 1995; Treu and Koopmans, 2004;
Dye and Warren, 2005; Koopmans et al., 2006, 2009; Barnabè et al., 2009; Auger
et al., 2010; Dutton and Treu, 2014; Ritondale et al., 2019; Powell et al.,
2022), from the combination of strong and weak lensing (Gavazzi et al., 2007),
from stellar dynamics only (Bertin and Stiavelli, 1993; Gerhard et al., 2001;
Thomas et al., 2007; Tortora et al., 2014; Bellstedt et al., 2018; Derkenne et
al., 2021), and from X-ray luminosity (Humphrey and Buote, 2010). This
phenomenon – that the total mass profile in ellipticals approximately follows
the power law, whereas neither the baryonic nor the dark matter component
individually follows the power law – is referred to as the “bulge–halo
conspiracy” (Treu and Koopmans, 2004), similar to the “disk–halo conspiracy”
in spiral galaxies (van Albada and Sancisi, 1986). Fig. 5 illustrates the
power-law mass profile for seven SLACS lens galaxies constrained from lensing
imaging data and the individual dark and luminous components that were
constrained by combining kinematic and weak lensing information with strong
lensing (Shajib et al., 2021). In most cases, the total density profile is
very close to a power-law form, with deviations only being prominent far from
the half-light, or effective, radius. Galaxy formation simulations suggest
that the close-to-isothermal nature of the total density profile originates
from rearranging the mass distribution through collisionless accretion in gas-
poor mergers (Johansson et al., 2009; Remus et al., 2013).

Analysis of the SLACS lenses finds the mean logarithmic slope
$\langle\gamma_{\mathrm{pl}}\rangle=2.08\pm 0.03$ with an intrinsic scatter of
$0.16\pm 0.02$ from a sample of 85 galaxy–galaxy lenses (Koopmans et al.,
2009; Auger et al., 2010). The median redshift of the SLACS lenses is $\langle
z_{\rm SLACS}\rangle\simeq 0.19$. The SIE lens model was adopted in this
analysis to constrain $\theta_{\mathrm{E}}$ from the imaging data. Then, the
power-law index $\gamma_{\mathrm{pl}}$ was obtained from the stellar velocity
dispersion measured by the SDSS. Shajib et al. (2021) re-analyze 23 systems
from the SLACS sample with a power-law model instead of the SIE model and
obtained the mean logarithmic slope
$\langle\gamma_{\mathrm{pl}}\rangle=2.08\pm 0.03$ with an intrinsic scatter of
$0.13\pm 0.02$. Ritondale et al. (2019) analyze 17 galaxy–galaxy lens systems
from the BELLS sample with a power-law mass model to find the average
logarithmic slope $\langle\gamma_{\mathrm{pl}}\rangle=2.00\pm 0.01$. The mean
redshift of the BELLS sample is approximately $\langle z_{\rm
BELLS}\rangle\simeq 0.5$.

The surface mass in elliptical galaxies constrained from strong lensing can
explain the origin of the so-called tilt of the fundamental plane, i.e., the
tight correlation between the effective radius, the effective surface
brightness, and the central velocity dispersion (Ciotti et al., 1996). If the
surface mass is used instead of the surface brightness, then the resulting
mass plane (in place of the fundamental plane) is not tilted (Bolton et al.,
2008). This result implies that the tilt of the fundamental plane stems from
the increase in the dark matter fraction with increasing velocity dispersion
or dynamical mass (Auger et al., 2009).

In summary, the total mass profile in galaxies seems to be well described by a
power law. The implications of this finding for the individual baryonic and
dark matter components are discussed separately below.

Figure 5:  3D mass distribution in dark matter and baryonic matter in lensing
elliptical galaxies constrained by Shajib et al. (2021). The power law mass
distribution (orange) refers to an imaging-data-only-based constraint that
uses a power-law lens model. Shajib et al. (2021) combined strong lensing
imaging with the stellar kinematics and weak lensing information to decouple
baryonic (blue) and dark matter (grey) contributions to the total mass density
profile (purple). In most cases, the total density profile is close to the
power-law profile, with occasional deviations appearing near the center or far
outside the effective radius $R_{\rm eff}$.

####  4.1.1 Luminous (baryonic) mass profile

In recent years spatially resolved spectroscopic surveys of nearby ellipticals
using integral field units (IFUs) and high-resolution HST imaging have
considerably advanced our understanding of the structure and evolution of
these galaxies (Cappellari, 2016). Such detailed studies of intermediate
redshift galaxies ($z\sim 0.2-0.7$) using direct observations are not
possible. However, the mass of galaxies acting as strong lenses at these
redshifts can be mapped out in detail, thereby providing critical information
about the progenitors of present-day galaxies and galaxy evolution.

In the last decade, imaging and kinematic data have led to a revision of
galaxy classification. The modern analyses classify ellipticals into fast and
slow rotators, and those with and without central cores, resulting in a total
of 4 classes (Cappellari, 2016; Krajnović et al., 2020). Massive ellipticals
($M\gtrsim 10^{11}\mbox{M}_{\odot}$) in the local Universe tend to have small
cores with size 0.02–0.50 kpc (Krajnović et al., 2020), i.e., much smaller
than the Einstein radius of typical lenses. In contrast, less massive ones
appear to have cuspy central light profiles. Cores are believed to be the
result of in-spiraling super-massive black hole (SMBH) binaries, which
transfer their angular momentum outward to stars and leave a flatter density
core at the center. The population of massive elliptical galaxies is also more
likely to be composed of slow rotators, being also morphologically rounder
than the lower mass fast rotator counter-part (e.g., Weijmans et al., 2014;
van de Sande et al., 2017).

Although it is well known that the baryons and the dark matter do not follow
the same radial profile, it is not known a priori whether they have the same
angular structure or not. Strong lensing allows us to compare the azimuthal
distributions of the light and the total mass distributions, thus detecting
any possible difference between the angular structures of the dark matter and
the baryons. Several studies report a strong correlation in the ellipticity
between the matter and light distributions (e.g., Koopmans, 2006; Gavazzi et
al., 2012; Sluse et al., 2012; Kostrzewa-Rutkowska et al., 2014). In contrast,
some other studies only report weak or no correlation (e.g., Keeton et al.,
1998; Ferreras et al., 2008; Rusu et al., 2016; Shajib et al., 2019, 2021).
Some differences can be attributed to data quality, modeling procedure, or
selection effects (Shajib et al., 2021). The major axes are usually well
aligned ($\lesssim 10^{\circ}$) between the mass and light distributions
(e.g., Keeton et al., 1998; Kochanek, 2002; Koopmans, 2006; Treu et al., 2009;
Sluse et al., 2012; Bruderer et al., 2016; Shajib et al., 2019, 2021). The
cases where the major axes of the mass and light do not align within $\sim
10^{\circ}$ also generally have large external shear. These findings suggest
that stellar orbits highly misaligned with the potential can only be sustained
in non-isolated galaxies (indicated by the large external shear), which is
consistent with what is found in simulations (Heiligman and Schwarzschild,
1979; Martinet and de Zeeuw, 1988; Adams et al., 2007; Debattista et al.,
2015). Although the stellar IMF pertains to the luminous structure, we present
the strong lensing results on the stellar IMF separately in Sect. 4.3.

####  4.1.2 Dark matter profile

Numerical simulations show that the dark matter halos and the baryonic mass
within them initially follow the Navarro–Frenk–White profile (Navarro et al.,
1996, 1997) before star formation begins. However, the baryonic gas has to
cool down and fall inward for star formation to begin. The contraction in the
baryonic matter deepens the gravitational potential. Thus, the dark matter
distribution also contracts adiabatically in response (Blumenthal et al.,
1986). In the process of adiabatic contraction, the initial radius $r_{\rm i}$
of a dark matter particle and its final radius $r_{\rm f}$ is related as

$r_{\rm i}M_{\rm i}(r_{\rm i})=r_{\rm f}M_{\rm f}(r_{\rm f}),$ (18)

where $M(r)$ is the enclosed 3D mass within radius $r$ (Blumenthal et al.,
1986). However, numerical simulations find that dark matter does not fully
respond to the baryonic infall according to the theoretical model of
Blumenthal et al. (1986) (e.g., Gnedin et al., 2004; Abadi et al., 2010).
Dutton et al. (2007) prescribe a formalism defining a halo response parameter
$\nu$ to adjust the degree of contraction (i.e., the response to the baryonic
infall) as

$r_{\rm i}\equiv\Gamma^{-\nu}(r_{\rm f})\ r_{\rm f},$ (19)

where $\Gamma(r_{\rm f})\equiv r_{\rm f}/r_{\rm i}$ is the contraction factor.
In this formalism, $\nu=0$ corresponds to no contraction, and $\nu=1$
corresponds to fully responsive contraction according to the model of
Blumenthal et al. (1986). The simulations of Gnedin et al. (2004) point to
$\nu\sim 0.8$, and those of Abadi et al. (2010) to $\nu\sim 0.4$ (Dutton and
Treu, 2014).

Several studies have combined additional information, e.g., stellar kinematics
and weak lensing, with the strong lensing data to decouple the dark matter
distribution and the baryonic distribution from the total density profile that
strong lensing can constrain. Dutton and Treu (2014) adopt multiple dark
matter contraction models with fixed $\nu$ to values between -0.5 and 1 and
find that $\nu=0$, i.e., no contraction, provides the best match to the data
from the SLACS sample among all the adopted models. Shajib et al. (2021) allow
for a fully variable $\nu$ parameter within -0.5 and 1 in their dark matter
model and find the average contraction in their sample to be
$\langle\nu\rangle=-0.03^{+0.04}_{-0.05}$, which is consistent with no
contraction and rules out the contraction results from simulations with high
statistical significance (Fig. 6).

Figure 6:  The halo response parameter $\nu$ constrained by Shajib et al.
(2021) is consistent with no adiabatic contraction in the dark matter halo
with an initial NFW density profile. The blue-shaded distribution corresponds
to a model with a spatially constant mass-to-light ratio, and the red-shaded
distribution corresponds to a baryonic distribution with a mass-to-light ratio
gradient.

###  4.2 Galaxy evolution

Observations of elliptical galaxy structure at different redshifts – using
strong lensing constraints from HST imaging, or IFU kinematics, or both – have
provided crucial insights into formation pathways of elliptical galaxies.
However, observations of distant galaxies at different redshifts are, by
nature, snapshots of different galaxies at different times. To form a
statement on the assembly history of these galaxies, it is necessary first to
identify evolutionary trends (or the lack thereof) in galaxy properties.
Numerical simulations have been an essential tool to interpret the presence
(or absence) of such evolutionary trends in terms of astrophysical processes.
Cosmological and hydrodynamical simulations adopt different prescriptions for
gravitational and baryonic interactions. By comparing the properties of the
simulated galaxies with those of the observed ones, we can identify the
appropriate interactions or processes at play to assemble the elliptical
galaxies under consideration. From cosmological $N$-body simulations, it has
been evident that the elliptical galaxies embedded in their dark matter halos
are the end results of hierarchical mergers of smaller dark matter halos
(e.g., White and Rees, 1978; Navarro et al., 1996). As a result, the structure
of the baryonic and dark matter density profiles observed in these galaxies
carry information on their assembly history.

Here, we state some key observations from strong lensing studies that shed
light on the formation and evolution of massive ellipticals. The observed
“bulge–halo conspiracy” described in Sect. 4.1 requires a particular mixture
of adiabatic contraction and feedback-driven expansion along the galaxies’
formation history. The near isothermality of the total mass distribution can
be explained as the end result of collisionless accretion of matter in gas-
poor mergers, where after achieving isothermality, further mergers preserve
this feature (e.g., Nipoti et al., 2009; Johansson et al., 2009; Remus et al.,
2013; Wang et al., 2019). As a result, gas-poor mergers appear to be the
primary formation mechanism for massive ellipticals (van der Wel et al.,
2009). However, dissipational gas-rich mergers cannot be completely ruled out
due to the observed existence of contracted central dark matter profiles
(Sonnenfeld et al., 2014). Furthermore, simulations predict that the dark
matter distribution would adiabatically contract due to the baryonic infall
through cooling at the onset of star formation ($z\sim 2$) (Blumenthal et al.,
1986; Dubinski, 1994; Kazantzidis et al., 2004; Gnedin et al., 2004; Abadi et
al., 2010). However, strong lensing studies find the dark matter halos
consistent with the “un-contracted” NFW profile at $z\sim 0.2$ (Dutton and
Treu, 2014; Shajib et al., 2021). To reconcile the prediction from simulations
for the dark halo shape at $z\sim 2$ with the observed one at $z\sim 0.2$, the
presence of baryonic feedback mechanisms is necessary, which can expand the
dark halos by transferring baryons outward. Although the relative importance
between the stellar, supernova, and AGN feedbacks in the evolutionary track of
ellipticals have not been conclusive yet, comparing structural properties of
simulated and observed ellipticals suggests that AGN feedback has been non-
negligible in addition to stellar and supernova feedbacks (Remus et al., 2017;
Shajib et al., 2021). This scenario is consistent with the findings from the
IllustrisTNG hydrodynamic simulation, where the density profile gets steeper
than isothermality until $z\sim 2$ during the initial baryonic infall.
Primarily AGN feedback makes the profile close to isothermal again until
$z\sim 1$, after which the elliptical galaxies passively evolve mainly through
gas-poor mergers that preserve their near-isothermal nature with little to no
change in the logarithmic slope (Wang et al., 2019, 2020).

We can summarize the above discussion on the formation history of ellipticals
as follows (see also Koopmans et al., 2006). The initial mass density profile
of the dark matter and baryons in halos follow the NFW profile. The baryons
fall inward through cooling and trigger star formation, and as a result of the
baryonic infall, the dark matter halos contract up to $z\sim 2$. Then, several
feedback mechanisms – i.e., stellar, supernova, and AGN – counteract the
initial adiabatic contraction by transferring baryons outward and thus
expanding the dark halo shape back to the NFW profile. The dark matter halo
and the embedded galaxy grow in size through a few major mergers and multiple
minor mergers, most of which are gas-poor and thus dissipationless (Nipoti et
al., 2009; Newman et al., 2012). These mergers lead to the galaxy being
pressure-supported and tri-axial in shape, thus forming ellipticals out of
disky or spiral galaxies. Dissipationless mergers lead to isothermality in the
total density profile, and such mergers conserve the isothermality after it
has been achieved.

Nevertheless, many open issues and potential tensions exist in investigating
galaxy evolution by comparing simulations with the observed ones using strong
lensing. We discuss these open issues in Sect. 6.5.

###  4.3 Stellar initial mass function

When combined with ancillary data, strong lensing provides a method to infer
the stellar mass-to-light ratio ($M_{\star}/L$) (Treu et al., 2010), which can
then be related to the low-mass end of the stellar initial mass function (IMF)
slope in the lensing galaxies (for a review, see Smith, 2020). This is mainly
because low-mass stars ($M_{\star}<0.5M_{\odot}$) contribute only by a few
percent to the integrated light in the optical, but they give a much larger
contribution to the mass (Conroy and van Dokkum, 2012).

As described in Sect. 4.1, the total mass distribution in the lensing galaxy
can be decomposed into the stellar and dark contribution by combining the
lensing and dynamical observables. Thus, the total stellar mass
$M_{\star}^{\rm LD}$ can be obtained. An independent method to get the stellar
mass $M_{\star}^{\rm SPS}$ of the galaxy is the stellar population synthesis
(SPS) method applied on the lens galaxy’s photometric or spectroscopic data
(Spiniello et al., 2011, 2012). The $M_{\star}/L$ computed via SPS analysis
depends on the choice of the IMF slope. In particular, a bottom-heavier IMF
implies a larger $M_{\star}/L$, because dwarf stars contributed more to the
mass than to the light. Hence, the mismatch parameter: $\alpha_{\rm IMF}\equiv
M_{\star}^{\rm LD}/M_{\star}^{\rm SPS}$. can be used to infer the lightness or
heaviness of the IMF (Treu et al., 2010). For example, if a light IMF – such
as the Chabrier IMF (Chabrier, 2003) – is adopted in the SPS method, then
$\alpha_{\rm IMF}\sim 1$ would point to the Chabrier IMF being consistent with
the lensing and dynamical observables. However, $\alpha_{\rm IMF}\sim 2$ would
indicate that the IMF in the lens galaxies is bottom-heavier (i.e.,
characterized by a larger number of dwarf stars), with a slope more similar to
that of a Salpeter IMF (Salpeter, 1955), or even steeper. Note that
$M_{\star}^{\rm LD}$ may depend on lens modeling assumptions, such as the
choice of dark matter density profile, the stellar mass-to-light ratio being
spatially constant or varying, and the assumed anisotropy profile of the
stellar orbits in the dynamical modeling (e.g., see the discussions in
Sonnenfeld et al., 2018, 2019).

One of the greatest advantages of lensing, when combined with dynamics and
stellar population analysis, is that it allows not only the inference of the
IMF slope but also to put constraints on the low-mass cutoff
($M_{\text{low}}$). Barnabè et al. (2013) has shown this by studying two
strong lenses from the X-Shooter Lens Survey (Spiniello et al., 2011) for
which both HST imaging (for precise lens modeling) and X-Shooter spectra (for
stellar population analysis) were available. Chromatic, microlensing-induced
flux anomalies in a galaxy–quasar strong lens can also be used to constrain
the stellar IMF (Schechter et al., 2014). This technique is described in
Vernardos et al. (in preparation), which discusses the theory and applications
of microlensing.

Whereas the IMF within the Milky Way is light – i.e., consistent with the
Chabrier IMF regardless of the stellar population age and environment
(Chabrier, 2003; Bastian et al., 2010) – the majority of strong lensing
studies on elliptical galaxies report consistency with a heavier IMF slope
(e.g. Spiniello et al., 2011; Sonnenfeld et al., 2012; Oldham and Auger,
2018). For instance, the SLACS analysis – by combining lensing and dynamics –
finds the IMF in elliptical galaxies at mean redshift $\langle z\rangle\sim
0.2$ to be consistent with the Salpeter IMF, i.e., $\alpha_{\rm IMF}\sim 2$
(see Fig. 7; Treu et al., 2010, also, Grillo et al. 2009). This result is
reproduced with more flexible models for the same SLACS systems or a subset of
them (Auger et al., 2010; Shajib et al., 2021). This is also in agreement with
the more general findings, based on dynamics or SPS analysis only, that the
IMF is bottom-heavier for more massive galaxies in general (Cappellari et al.,
2012; La Barbera et al., 2013; Spiniello et al., 2014).The general picture is
that the low-mass end of the IMF slope might not be universal across all
galaxies, as generally assumed in the last 30 years. However, a consensus on
the physical mechanisms that could be responsible for its variation has not
yet been reached. According to theoretical work (e.g., Hopkins 2013; Chabrier
et al. 2014), high density, temperature, and turbulence of the gas are key
parameters that drive the fragmentation of molecular clouds (higher density
and temperature make the fragmentation easier, forming more dwarf stars, i.e.,
a bottom-heavier IMF).

In contrast with the studies mentioned above on massive, local elliptical
galaxies, Sonnenfeld et al. (2019) find $\log_{10}\alpha_{\rm IMF}=-0.04\pm
0.11$ with respect to the Chabrier IMF in the SPS method for a sample of
strong lensing galaxies at $z\sim 0.6$, which is in tension with Shetty and
Cappellari (2014) that reports consistency with a Salpeter IMF for a galaxy at
$z\sim 1$ from a dynamical analysis. However, allowing spatial gradients in
the $M/L$ can alleviate this tension (Sonnenfeld et al., 2018, 2019). Indeed,
a radial gradient in the $M_{\star}/L$, or equivalently a radially varying
IMF, has been reported by studies based on the SPS method applied on local
massive ellipticals (Martín-Navarro et al., 2015; van Dokkum et al., 2017;
Sarzi et al., 2018; Barbosa et al., 2021a, b). These authors find that the
central region within $\sim 2$ kpc has a heavy IMF (even super-Salpeter), and
the IMF in the outer regions gradually becomes light (i.e., Milky-Way-like).
The current belief is that the IMF is bottom-heavy for stars formed very early
in cosmic time via a quick and violent SF burst. These stars usually form a
“red nugget” ($z\sim 2$, Damjanov et al. 2011): an ultra-compact red-and-dead
massive core. Then, with a second and more time-extended phase, red nuggets
merge, interact with other structures in the Universe, and accrete gas. This
causes a growth in size up to a factor of $\sim 5$, and (only slightly) in
mass, transforming them into giant local, massive ellipticals. Depending on
the merger history of each single galaxy, the red nugget can remain (almost)
untouched in the innermost region, dominating the light there. In this case,
one would measure a bottom-heavy IMF in this region. This seems to be the case
for NGC 3311, the central galaxy in the Hydra cluster (Barbosa et al., 2021a),
M87 (Sarzi et al., 2018), and many other very massive low-z early-type
galaxies. However, the red nugget can also be destroyed or contaminated by
accreted or lately formed stars. In that case, one would measure a Milky Way
like IMF, which is the characteristic for stars formed later on and through
more time-extended star formation channels. This “two-phase formation
scenario” (Naab et al., 2009) is also supported by the discovery of “relic
galaxies" (Trujillo et al., 2014; Spiniello et al., 2021), the local
counterparts of red nuggets that somehow wholly missed the size-growth and
evolved passively and undisturbed across cosmic time. The IMF slope for these
peculiar and rare objects has been measured to be steep everywhere up to at
least one effective radius (Ferré-Mateu et al., 2017). Finally, simulations
have lately shown that not all elliptical galaxies formed via this two-phase
formation scenario (e.g., Pulsoni et al., 2021), although this becomes more
and more common with increasing stellar mass.

The finding of the Salpeter IMF from a combination of lensing and dynamics in
the near-by Universe is consistent with this scenario, as strong lensing
information is sensitive to the galaxy’s inner region (typically $\lesssim 6$
kpc). However, at higher redshift, lensing probes larger and larger region,
which explain the results presented in Sonnenfeld et al. (2019). A radially
decreasing $M_{\star}/L$ also explains or alleviates the reported tension
between lensing-based studies themselves (Sonnenfeld et al., 2018, 2019;
Shajib et al., 2021).

Furthermore, studies of few local massive lenses for which the Einstein radius
is much smaller than the effective radius, and hence where the stellar mass
dominates the lensing inference, indicate that bottom-heavy IMFs are excluded
by lensing (Ferreras et al., 2010; Smith et al., 2015; Leier et al., 2016).
The only way that has been proposed so far to alleviate this disagreement is
to allow for a variable cut-off in the low mass end of the IMF (Newman et al.,
2017).

In conclusion, the currently preferred scenario sees the majority of massive
galaxies having a bottom-heavy IMF in their innermost region, where the
pristine stellar population formed at $z>2$ through a star formation burst
dominate, while stars in outskirt are distributed by a Milky-Way like IMF.
However, depending on the single galaxy’s detailed merger tree and cosmic
evolution, the IMF can differ from system to system.

In the future, IFU-based stellar kinematics and population analysis of strong
lensing galaxies in combination with lensing constraints can shed light on the
presence/absence of IMF variations and spatial gradients in ellipticals at
intermediate redshifts. However, to properly track any evolution in the IMF
properties of elliptical galaxies across redshift, it would be essential to
mitigate systematic impacts through a uniform choice of models and to account
for selection differences between samples.

Figure 7:  Summary of measurements on the stellar IMF from lensing and
dynamical observations (Smith et al., 2015).

###  4.4 Constraints on the very central densities and super-massive black
holes from central images

Gravitational lensing theory predicts that the number of multiple images must
always be odd. In systems with 3 (or 5) images, 2 (or 4) are formed roughly at
the Einstein radius from the lens center, which is
$<\\!2\mbox{${}^{\prime\prime}$}$ for galaxies. The odd 3rd (or 5th) image is
formed very near the center of the lens. It is always demagnified, usually
significantly, and superimposed on the light of the lensing galaxy, making it
hard to detect. Most existing searches for central images of strongly lensed
quasars (given that they are much brighter sources) rely on optical, or radio
wavelengths. The radio wavelengths are most favorable as the lensing galaxy is
generally transparent in that range, but these investigations are limited by
the fact that quasars are generally radio-quiet. Out of $\sim\\!200$ doubles
discovered to date, only 2 have observed central images where the lens is a
single galaxy444There are at least 5 cases of central images where the lens
has two or more main lensing galaxies.: PMN J1632$-$0033, demagnified to
0.004, or by 6 mag compared to the brightest image (Winn et al., 2004), and
PKS 1830$-$211, demagnified to 0.007, or by 5.4 mag (Muller et al., 2020). No
reliable detection exists of the $\sim\\!50$ known quads with a single lensing
galaxy. Upper limits have also been measured. For example, Quinn et al. (2016)
place an upper limit of $\sim 10^{-4}$ on the magnification of the central
image with respect to the brightest visible image in the double B1030$+$074\.
Upper limits on the flux of the central image have also been placed at X-ray
wavelengths for several quads and doubles. The stacking of X-ray monitoring
data allows effective exposures of several hundreds of kilo-seconds without
any contamination from the lens. However, the data are still too shallow to
strongly constrain the lens galaxy’s inner density profile. One of the deepest
upper limits has been achieved for HE 0435-1123 (e.g., Chen et al., 2012;
Guerras et al., 2017).

Detection of central images, or the lack thereof, has been used to place
constraints on the central mass density of the lensing galaxies and the mass
of the central SMBH (Mao et al., 2001; Rusin and Ma, 2001; Quinn et al.,
2016). Prospects for detecting supermassive binary black holes due to their
lensing effects are discussed in Li et al. (2012); Hezaveh et al. (2015).

###  4.5 Spiral galaxies

Due to their lower mass and presence of a disc, spiral galaxies have a
substantially lower lensing cross-section than elliptical ones (Keeton and
Kochanek, 1998). Only a handful of spiral galaxies lensing quasars have been
discovered to date. The most well-known and the first example of such a system
is the Einstein cross, which is lensed by the bulge of a nearby spiral galaxy
(Huchra et al., 1985). Targeted searches for galaxies lensed by a spiral have
been carried out, allowing the discovery of several dozens of systems (Treu et
al., 2011, and reference therein). Complications in studying those systems
arise from the dust and the disc mass component. Correcting for the reddening
by dust is needed to model extended lensed images, which rely on conserving
the surface brightness between the lens and the source plane. On the other
hand, the disc component yields strong discontinuities in the gravitational
potential and needs to be explicitly modeled using, e.g., an exponential mass
density profile. Whereas disc-like features can yield flux anomalies in lensed
quasars (Hsueh et al., 2016, 2017; Gilman et al., 2017), they can be
disentangled from the dust using multi-band data (Möller et al., 2003). Once
these ingredients are accounted for, the combination of kinematic and lensing
information can be used to break the degeneracy – that exists in dynamical
studies alone – between the disc and halo components (Maller et al., 2000;
Dutton et al., 2011; Suyu et al., 2012).

##  5 Applications in cosmology

In this section, we present applications of galaxy-scale lens systems to
measure cosmological parameters without using time-delay information.
Measurement of the Hubble constant ($H_{\mathrm{0}}$) and other cosmological
parameters based on the time delays are discussed in detail in Birrer et al.
(2022).

Here, we briefly establish some necessary definitions for use in this section.
More detailed explanations of the cosmological connection with strong lensing
formalism are given in Saha et al. (in preparation). The angular diameter
distance between two objects at redshifts $z_{1}$ and $z_{2}$ (with
$z_{1}<z_{2}$) for a flat Universe555See Saha et al. (in preparation), Section
2.1.1 for the general definition in the non-flat case. is given by

$D_{\rm
ang}(z_{1},z_{2})=\frac{c}{H_{0}(1+z_{2})}\int_{z_{1}}^{z_{2}}\frac{\mathrm{d}z^{\prime}}{E(z^{\prime})},$
(20)

where $E(z)\equiv H(z)/H_{0}$ is the dimensionless Friedman equation given by

$E(z)=\sqrt{{\Omega_{\mathrm{m}}}(1+z)^{3}+{\Omega_{\mathrm{r}}}(1+z)^{4}+{\Omega}_{\rm
de}(1+z)^{3(1+w_{\rm de})}}.$ (21)

Here, ${\Omega_{\mathrm{m}}}$, ${\Omega_{\mathrm{r}}}$, and ${\Omega}_{\rm
de}$ represent the matter, radiation, and dark energy density parameters,
respectively, at $z=0$. The parameter $w_{\rm de}$ is the equation-of-state
parameter of the dark energy given by $w_{\rm de}\equiv p_{\rm de}/\rho_{\rm
de}c^{2}$, where $p_{\rm de}$ and $\rho_{\rm de}$ denote the pressure and
density of the dark energy, respectively. In the $\Lambda$CDM model, $w_{\rm
de}=-1$ is assumed, i.e., the dark energy density stays constant through
cosmic time. The $w$CDM model is one natural extension of the $\Lambda$CDM
model, where $w_{\rm de}\neq-1$ is allowed; however, $w_{\rm de}<-1/3$ should
still be satisfied to reproduce an accelerated Universe.

In the following subsections, we discuss methods to estimate cosmological
parameters (primarily, $\Omega_{\mathrm{m}}$, $\Omega_{\rm de}$, and $w_{\rm
de}$) using multiple-source-plane lenses (Sect. 5.1), using the stellar
kinematics of strong lenses (Sect. 5.2), and using galaxy–galaxy lensing
statistics (Sect. 5.3).

Figure 8:  Example extraction of cosmological parameters from the double
source-plane lens system SDSS J0946$+$1006 (shown in Fig. 2). The 68, 95, and
99.7% confidence intervals on the $\Omega_{\rm m}$–$w$ plane from this double-
source plane system are shown in red and are from Collett et al. (2012), while
the black contours are the constraints from Planck 2013.

###  5.1 Utilizing multiple sources at different redshifts

Strong lens systems with multiple sources (i.e., compound lenses) at different
redshifts can be used as cosmographic probes. Currently, only a tiny sample of
galaxy-scale compound lenses are known (e.g., Lewis et al., 2002; Gavazzi et
al., 2008; Collett and Smith, 2020), but at the HST snapshot depth they are
expected to occur in about 1% of galaxy-scale strong lenses (Gavazzi et al.,
2008). Indeed, if one were to stare deeply at any single plane lens, other
sources would almost inevitably be discovered, as is spectacularly
demonstrated by the discovery of a third source (a $z\approx 6$ Lyman-$\alpha$
emitter) behind the Jackpot lens in deep Multi Unit Spectroscopic Explorer
(MUSE) data from the Very Large Telescope (Collett and Smith, 2020).
Forthcoming surveys are expected to discover $\mathcal{O}(1000)$ compound
lenses.

Since the Einstein radius is a function of the lens mass and the cosmological
distances, the ratio of Einstein radii in a compound lens with sources at two
or more redshifts is independent of the mass (Gavazzi et al., 2008; Collett et
al., 2012). In practice, the method also requires a complete understanding of
the lens density profile and additional lensing by the source galaxies and
other perturbing masses along the line of sight.

For a two-source-plane system with one primary lens, the lens equation can be
written as

$\displaystyle\mathbf{y}=\mathbf{x}-\beta_{12}\bm{\alpha}_{\rm
d}(\mathbf{x}),$ (22) $\displaystyle\mathbf{z}=\mathbf{x}-\bm{\alpha}_{\rm
d}(\mathbf{x})-\bm{\alpha}_{{\rm
s}1}\left(\mathbf{x}-\beta_{12}\bm{\alpha}_{\rm d}(\mathbf{x})\right),$

where $\bf x$ are positions on the image plane, $\bf y$ and $\bf z$ are the
unlensed positions of the first and second source, respectively,
$\bm{\alpha}_{\rm d}$ is the deflection caused by the primary lens,
$\bm{\alpha}_{{\rm s}1}$ is the deflection caused by the closer of the two
sources, and $\beta_{ij}$ is the cosmological scaling factor

$\beta_{ij}=\frac{D_{ij}D_{\mathrm{s}}}{D_{j}D_{i\rm s}}.$ (23)

For realistic redshifts and cosmologies, $\beta_{ij}$ is sensitive to the
matter fraction and the equation-of-state parameter $w_{\rm de}$ but has no
dependence on the Hubble constant $H_{\mathrm{0}}$.

Amongst galaxy-scale compound lenses, only the Jackpot lens (Gavazzi et al.,
2008), shown in Fig. 2, has been used to precisely constrain cosmology since
this system has a favorable redshift configuration – other compound lenses
have the multiple sources that are at similar redshifts, thus having
$\beta\approx 1$ regardless of the cosmology. Collett and Auger (2014) model
the HST imaging performing a pixellated reconstruction of both sources, and
make a 1.1% measurement on $\beta^{-1}$. Converting this into constraints on
the dark energy, this single compound lens with a CMB prior from Planck
Collaboration et al. (2014) constrains $w_{\rm de}$ to 0.2 precision, as shown
in Fig. 8. With hundreds of compound lenses expected in the Vera Rubin
Observatory Legacy Survey of Space and Time (LSST) and the Euclid surveys,
constraints on both $w_{\rm de}$ and its redshift derivative are expected to
be competitive with established cosmological probes.

###  5.2 Utilizing stellar kinematics of single source lenses

The enclosed projected mass inside the Einstein radius is independent of the
mass profile (Schneider et al., 1992). Combining this with the finding that
the density profile of elliptical galaxies is very close to isothermal (see
Sect. 4.1 and references therein), it is sufficient to model lens galaxies as
singular isothermal ellipsoids (SIE) when only a constraint on the Einstein
radius is needed. For a singular isothermal sphere (SIS), we have

$\theta_{\mathrm{E}}=4\pi\left(\frac{\sigma_{\rm
SIS}}{c}\right)^{2}\frac{D_{\mathrm{ds}}}{D_{\mathrm{s}}},$ (24)

where $\rm\sigma_{\rm los}$ is the velocity dispersion of the deflector. Using
this equation, one can constrain the distance ratio $d^{\rm obs}\equiv
D_{\mathrm{ds}}/D_{\mathrm{s}}$ using the observed Einstein radius
$\theta_{\mathrm{E}}$ and assuming that the observed velocity dispersion
$\sigma_{\rm los}=\sigma_{\rm SIS}$. Thus, a sample of single-source lenses
with measured $\sigma_{\rm los}$ can be used to estimate cosmological
parameters (Grillo et al., 2008) by maximizing the likelihood function

$\log\mathcal{L}(\mathbf{\Theta})=-\frac{1}{2}\sum_{i=1}^{N_{\rm
SL}}\frac{\left[d^{\rm
th}_{i}\left(z_{\mathrm{d}},z_{\mathrm{s}};\bf{\Theta}\right)-d^{\rm
obs}_{i}(\theta_{\mathrm{E}},\sigma_{\rm los})\right]^{2}}{(\delta
d^{\rm{obs}}_{i})^{2}},$ (25)

where $\mathbf{\Theta}$ is the set of free parameters in the assumed
cosmological model, $N_{\rm SL}$ is the number of lens systems, and $\delta
d^{\rm{obs}}$ is the uncertainty of each $d^{\rm obs}$ that depends on the
$\sigma_{\rm los}$ and $\theta_{\mathrm{E}}$ uncertainties.

Similar to the multiple-source lens systems, this method is also independent
of $H_{\mathrm{0}}$. However, the assumptions in this method should be
carefully scrutinized: (i) the measured $\theta_{\mathrm{E}}$ is independent
of the choice of the lens mass profile (e.g., Cao et al., 2015), and (ii) the
measured line-of-sight velocity dispersion is equal to that for an SIS
profile. In the early years, Ofek et al. (2003) estimated that systematics
might be up to 20% in $\theta_{\mathrm{E}}$, however more recent modeling
methods provide robust $\theta_{\mathrm{E}}$ measurement within a few per
cents regardless of the mass profile choice (Birrer, 2021). Treu et al. (2006)
argue that $\sigma_{\rm los}\simeq\sigma_{\rm SIS}$ for those lenses that are
elliptical galaxies with velocity dispersions in the range 200–300 km s-1.
These results are further confirmed by analyzing other samples (Bolton et al.,
2006; Auger et al., 2010).

Some authors (Leaf and Melia, 2018; Amante et al., 2020) find that certain
sub-samples provide nonphysical values (i.e., $d^{\rm obs}>1$), providing
cosmological parameters that are inconsistent with other observations (such as
SNIa, CMB). Amante et al. (2020) suggest that those nonphysical values are
related to different kinds of uncertainties in the lens system, such as
complex lens substructure, uncertain redshifts, etc. When these systems are
excluded, the estimated cosmological parameters are consistent with those
expected from other observations (Magaña et al., 2014; Betoule et al., 2014;
Aghanim et al., 2018; Motta et al., 2021).

###  5.3 Utilizing galaxy–galaxy lensing statistics

The statistics of strong lensing were initially expected to be powerful for
probing the cosmological parameters (Fukugita et al., 1990). At the most basic
level, matter clusters can collapse to the densities required to form strong
lenses, whereas the cosmological constant does not cluster and cannot form
lenses. Therefore, it was expected that lensing rates should be suppressed for
larger values of $\Omega_{\Lambda}$. Further information is contained in the
Einstein radius distribution and the lens and source redshift distributions.
In practice, this topic has fallen out of fashion due to the cosmological
sensitivity being overwhelmed by astrophysical uncertainties of the unlensed
source population, the lens discovery selection function, and the lensing
properties of typical galaxies (e.g., Mitchell et al., 2005; Chae, 2010).

##  6 Open problems and future outlook

In this section, we discuss the current open problems that are expected to be
tackled in this decade: the selection function (Sect. 6.1), the self-
similarity assumption (Sect. 6.2), degeneracies in strong-lensing (Sect. 6.3)
and non-strong-lensing observables (Sect. 6.4), comparison with galaxy
simulations (Sect. 6.5) and statistical methods to perform such comparisons
(Sect. 6.6). We also provide future outlooks on these issues whenever
appropriate.

###  6.1 Selection function

Although strong lensing is a powerful probe to study galaxy properties, the
lensing phenomenon is a rare occurrence requiring a serendipitous alignment of
two line-of-sight objects separated by a large cosmological distance. Thus,
samples of strong lensing galaxies inherently occupy a tiny fraction of the
population of all galaxies. When strong lensing studies aim to infer
properties of the general population of galaxies based on such small samples
of galaxies, the lens sample’s selection function must be considered. Is
strong lensing selection inherently biased toward more massive galaxies, given
that larger magnification leads to easier discovery? Does the particular
search algorithm (i.e., based on spectroscopic lines or based on image
morphology) performed over a survey area limit the selection to a particular
sub-population of galaxies? These questions need to be taken into account
before generalizing properties inferred in samples of strong lenses to the
entire population of galaxies. Methodologically, selection functions have been
incorporated by some studies in either of the following two ways: (i) a
parametrized selection function in Bayesian hierarchical analysis (e.g.,
Sonnenfeld et al., 2019), and (ii) posterior predictive tests based on mock
samples applying the same selection procedure as the actual lens sample (e.g.,
Sonnenfeld et al., 2018).

####  6.1.1 Parametrized selection function

The selection function affects the prior probability $p(w_{\rm i}\mid\nu)$ of
the observed lens model parameters $w_{\rm i}$ for the $i$-th lens in a
sample, given population hyper-parameters $\nu$, through the selection
criteria applied on a set of lens and source galaxy properties $u_{i}$ as

$p(w_{\rm i}\mid\nu)=p_{\rm intr}(w_{\rm i}\mid\nu)\ p_{\rm det}(w_{\rm
i},u_{i}\mid\nu),$ (26)

where $p_{\rm intr}(w_{\rm i}\mid\nu)$ is the intrinsic distribution of the
lens model parameters $w_{\rm i}$ for a randomly selected galaxy, and $p_{\rm
det}$ is the probability of detection. Here, the specific properties in the
set $u_{i}$ – that $p_{\rm det}$ is sensitive to – depend on the particular
search algorithm. For example, most search algorithms discover brighter
systems more easily. Hence, the stellar mass $M_{\star}$ of the lens and
source galaxies, which is a proxy of the system’s total flux (magnitude), can
be included in $u_{i}$. The parameter sets $w_{i}$ and $u_{i}$ can have common
members depending on the particular parametrization of the lens model. The
reader is referred to Sonnenfeld et al. (2015, 2019) for specific examples of
parametrized selection functions.

####  6.1.2 Posterior predictive test

Suppose the selection function is too complicated to express in a functional
form. In that case, a posterior predictive test can be performed to evaluate
the appropriateness of the adopted model in retrieving accurate
hyperparameters. In this test, the same selection criteria are applied to the
parent sample of elliptical galaxies to create multiple realizations of mock
strong lensing samples. Then, the posteriors $p(v\mid\mathcal{D})$ are
compared between the real sample and the mock samples. If the posteriors
agree, then the adopted model’s appropriateness in recovering the correct
hyper-level parameters is validated. The reader is referred to (Sonnenfeld et
al., 2018) for a specific example of a posterior predictive test on a lens
sample.

Large lens samples curated from multiple search algorithms (e.g., the sample
used in Shajib et al., 2019) may not allow for a simple parametrization of the
selection function. In such cases, the posterior predictive test can be
applied. However, that comes with an increased computational burden in the
analysis. In the future, lenses discovered from a single survey data – e.g.,
from the Vera Rubin, Euclid, and Roman observatories – through a uniform
searching algorithm can provide a large sample with sufficiently high
statistical constraining power. However, the necessity to properly consider a
selection function cannot be circumvented.

###  6.2 Assumption of self-similarity

When hierarchically inferring population properties of lenses by combining
multiple strong lensing systems (Sect. 3.4), it is crucial to describe and
quantify potential differences between subsets of the considered sample for
accurate population-level inference. Differences in the population may arise
from the sample selection due to search criteria and techniques (Sect. 6.1) or
intrinsic differences in the sources. Such differences in secondary selection
might impact and bias population level constraints when assuming that two
different samples can be described with an identical underlying population.
For example, Birrer et al. (2020) present a hierarchical analysis under the
assumptions that the same population level parameters describe the SLACS and
the Time-Delay COSMOgraphy (TDCOSMO) lenses to constrain the mass density
profiles of the time-delay lenses better. There are different ways to mitigate
such assumptions. The first is to select a purified sample of lenses as self-
similar as possible when performing hierarchical analyses. However, although
the currently small number of known lenses will increase through future
surveys and enable this approach, it is still unclear in which parameter space
these lenses can be considered self-similar and whether there are important
latent variables to consider. The second way is to model and describe
differences between populations on a first-principle level, folding in
differential selection effects and other aspects into the analysis (as stated
in Sect. 6.1).

###  6.3 Degeneracies in strong-lensing observables

The lensing data’s ability to constrain the lens’s radial structure,
particularly in disentangling dark and luminous mass components, strongly
depends on the data quality. As discussed in Saha et al. (in preparation), the
modeling of point-source astrometry primarily encodes information on the
lens’s quadrupole moment. Still, it provides limited constraints on the
monopole (the total mass). This, however, does not mean that any monopole
model combined with a quadrupole component will accurately reproduce a set of
lensed image positions. For some choice of the monopole, extreme (unphysical)
values of the quadrupole may be required, naturally excluding some
mathematically correct solutions. This explains why, e.g., a single component
constant mass-to-light ratio model generally yields large shear amplitudes to
reproduce the observed astrometry of lensed systems to high accuracy (e.g.,
Sluse et al., 2012).

On the other hand, while extended lensed images may provide detailed azimuthal
information and allow one to constrain the ratio of radial magnifications at
different galactocentric distances, their effective constraint on the density
profile remains sensitive to the MSD (e.g., Sonnenfeld, 2018). While the MSD
provides mathematically large leverage to modify the results, its impact may
remain in practice generally small, with a typical change on the total density
profile that may not exceed a few tens percent. The prior on the mass profile
(i.e., choice of families of mass distribution or free-form with some
regularization) may further limit the impact of degeneracies. One may,
however, need to be careful when choosing a mass distribution as a too rigid
model compared to the true mass density may yield too biased or underestimated
parameter uncertainties (e.g., Sonnenfeld, 2018; Kochanek, 2021). Finally, it
is important to note that point images with measured time delays limit the
impact of degeneracies in lensing-only observables for an assumed cosmology
(Saha and Williams, 2001; Kochanek, 2002).

The constraint on the power-law slope of galaxies described in Sect. 4.1
results from the combination of lensing and kinematics data. The lensing data
either consist of an estimate of the Einstein radius (e.g. Bolton et al.,
2006) or is based on models of extended images. By definition, the power-law
index requires the assumption of a power-law density model for both lensing
and kinematics data. For that reason, it is insensitive to mass-sheet
degeneracy. The comparison of power-law indices of lensing galaxies as a
function of redshift requires, however, to be handled with care. Indeed, as
the lens and source redshifts vary, lensed images are formed at different
galactocentric radii, such that minor deviations from the power-law model may
mimic an (absence of) evolution of the galaxy total density profile with
redshift. To account for this effect, (Dutton and Treu, 2014) suggest using a
mass-weighted slope.

###  6.4 Degeneracies in non-strong-lensing observables

The modeling of stellar kinematics data requires an assumption on the
anisotropy profile of stellar orbits. Integrated velocity dispersions obtained
from single-slit spectra cannot constrain the anisotropy profile, which leads
to the so-called mass–anisotropy degeneracy (Treu and Koopmans, 2002) —
typically adopted anisotropy profiles are either isotropic or the
Osipkov–Merritt profile (Osipkov, 1979; Merritt, 1985a, b). Whereas the
isotropy assumption does not entail any free parameter, the Osipkov–Merritt
profile depends on a scale radius. Due to the mass–anisotropy degeneracy, the
posterior of the anisotropy scale radius is dominated by the adopted prior
(e.g., Shajib et al., 2018). Particular choices of the anisotropy profile and
the associated prior may lead to systematic differences between studies
involving strong lensing and kinematics data. For example, Sonnenfeld et al.
(2018) find the existence of a mass-to-light ratio gradient in the SLACS
lenses assuming isotropic orbits, whereas Shajib et al. (2021) find
consistency with a constant mass-to-light ratio assuming the Osipkov–Merritt
anisotropy profile for a subsample of SLACS.

Spatially resolved velocity dispersion measurements can better constrain the
anisotropy profile by breaking the mass–anisotropy degeneracy. However,
spatially resolved kinematics data from integral field unit (IFU) spectroscopy
are more expensive to obtain than an integrated measurement from single-slit
spectroscopy. Thus usage of such data in lensing studies has been limited
(see, e.g., van de Ven et al., 2010; Barnabè et al., 2011, for example of IFU
data being combined with strong lensing).

###  6.5 Comparison with galaxy simulations

The past decade has seen significant progress in understanding galaxies’
structure and formation. Within the $\Lambda$CDM paradigm, there is general
agreement regarding the gravitational part of galaxy formation. The
‘gastrophysics’ is less well understood and requires subgrid models to
simulate, but still, the galaxies formed in simulations like Illustris
(Vogelsberger et al., 2014), FIRE (Hopkins et al., 2014), and EAGLE (Crain et
al., 2015) are much more credible than previous generations of simulated
galaxies. There are also equilibrium galaxy models, including stars, gas, and
dark matter, of which the AGAMA (Vasiliev, 2019) simulations are arguably the
most sophisticated.

Despite the great advances in the fidelity of the simulations, there have been
discrepancies between the simulated predictions and the observed properties of
galaxies. In particular, simulations have not been successful yet in
reproducing the observed distributions of the logarithmic slope
$\gamma_{\mathrm{pl}}$ and the dark matter fraction $f_{\rm dm}$
simultaneously. A no or weak feedback prescription was required in some of the
simulations to reproduce the $\gamma_{\mathrm{pl}}$ distribution, which,
however, led to underestimating $f_{\rm dm}$ compared to the observations due
to overestimating the star formation efficiency (e.g., Naab et al., 2007;
Duffy et al., 2010; Johansson et al., 2009). Similarly, matching the $f_{\rm
dm}$ distribution required strong feedback prescriptions, but these produce
too shallow $\gamma_{\mathrm{pl}}$ compared to the observed distribution. More
recently, the IllustrisTNG simulation reproduced a $f_{\rm
dm}$–$\gamma_{\mathrm{pl}}$ distribution that is consistent with the strong
lensing observations if the stellar IMF corresponds to the Salpeter IMF (Wang
et al., 2020; Shajib et al., 2021). In contrast, Mukherjee et al. (2022), find
that the EAGLE and SLACS $f_{\rm dm}$–$\gamma_{\mathrm{pl}}$ distributions
agree while using a Chabrier IMF, supporting the important role played by
feedback and sub-grid physics in reproducing this relation.

Furthermore, there has been an apparent tension in the redshift evolution of
the logarithmic slope $\gamma_{\mathrm{pl}}$ between observations and
simulations. Strong lensing observations report a steepening of
$\gamma_{\mathrm{pl}}$ with decreasing redshift at $z<1$ (Ruff et al., 2011;
Sonnenfeld et al., 2013). Such a steepening would require dissipative
processes through wet mergers along the evolutionary track of elliptical
galaxies. Simulations instead find no evidence for redshift evolution of
$\gamma_{\mathrm{pl}}$ below $z=1$, or a slightly shallowing trend in
$\gamma_{\mathrm{pl}}$ with decreasing redshift (see Fig. 9, Xu et al., 2017;
Remus et al., 2017; Wang et al., 2020). This tension vanishes if the same
strong lensing analysis is applied to the simulated galaxies from Illustris
and Magneticum to extract $\gamma_{\mathrm{pl}}$ by combining lensing and
kinematic information (Xu et al., 2017; Remus et al., 2017). Therefore, the
modeling systematics in the joint analysis of lensing and dynamical
observables cannot be ruled out as the source of the above tension.

The specific importance of various baryonic processes in the evolution of
ellipticals is still an open question. It would require further fine-tuning in
the simulations to consistently reproduce all of their observable properties.
However, modeling systematics and the selection function must be appropriately
considered for accurate comparison between simulations and observations. The
first results of the statistical framework outlined in the following
subsections for inference of the galaxy evolution scenario, when jointly
considering simulations and observations, indicate that AGN feedback is
essential (Denzel et al., 2021a) but application to a large sample of lenses
is required.

Figure 9:  Comparison of measured logarithmic slopes of the total density
profile at different redshifts between strong lensing measurements (grey
points), dynamical measurements (red points), and simulation results (lines
and shaded region). The figure is taken from Wang et al. (2020). Although
strong lensing observations find steeping of $\gamma_{\mathrm{pl}}$ with time
below $z=1$, simulations find no trend or slightly decreasing trend in
$\gamma_{\mathrm{pl}}$ with decreasing redshift. However, modeling systematics
cannot be ruled out as the source of this tension.

###  6.6 Statistical framework for comparison with simulations

An exciting prospect to circumvent the potential systematics associated with
lens modeling would be comparing the simulated galaxies directly with lensing
observations. In general terms, the idea is straightforward. Let $F$ stand for
some formation scenario of galaxies. This could be a particular feedback
prescription in hydrodynamic simulation or some family of galaxy models in
phase space. Scenario $F$ produces simulated galaxies $g$ with a distribution
$P(g\,|\,F)$. Let $G$ denote the data on an observed galaxy. The quantity

$P(G\,|\,F)={\textstyle\sum_{g}}\,P(G\,|\,g)\,P(g\,|\,F)$ (27)

is the likelihood of $G$ in scenario $F$. The posterior probability of the
formation scenario given the observed galaxy is

$P(F\,|\,G)=\frac{P(G\,|\,F)\,P(F)}{P(G)}.$ (28)

The product of this quantity

$P(F\,|\,{\\{G\\}})={\textstyle\prod_{G}}\,P(F\,|\,G)$ (29)

over an entire catalog of observed lensing galaxies ${\\{G\\}}$ is the
ultimate quantity of interest.

The observations of a lensing galaxy will not be only the lensing data (say
$G_{\rm lens}$). There will typically be data on the stellar light (say
$G_{\rm light}$) and on the kinematics of the stars (say $G_{\rm kin}$), and
possibly X-ray emission and other kinds of data as well. Let us consider
$G_{\rm lens}$ first. To assess whether some simulated galaxy $g$ can fit some
lensing data, we need to allow for different galaxy orientations and different
source brightness distributions behind it. Orienting the simulated galaxy by
some $\iota$ (which denotes the Euler angles), we compute the deflection-angle
field $\bm{\alpha}(\bm{\theta}\,)$. Then we consider some unlensed brightness
distribution (say $u$) that will get lensed. The likelihood contribution has
dependencies as

$P(G_{\rm lens}\,|\,g,\iota,u)=P(G_{\rm
lens}\,|\,\bm{\alpha},u)\,P(\bm{\alpha}\,|\,g,\iota)\,.$ (30)

As the source distribution is not of particular interest in the present
context, one would like to marginalize it out

$P(G_{\rm lens}\,|\,g,\iota)=\int\\!P(G_{\rm lens}\,|\,g,\iota,u)\,du$ (31)

with a prior on $u$ if desired. The orientation is a nuisance parameter and
can be likewise marginalized out

$P(G_{\rm lens}\,|\,g)=\int\\!P(G_{\rm lens}\,|\,g,\iota)\,d\iota\,.$ (32)

The marginalization Eq. 32 over orientations is somewhat computationally
intensive but straightforward. The marginalization Eq. 31 over sources, on the
other hand, is complicated. It is necessary to approximate the integral by
fitting for sources that give non-negligible contributions to the integral. A
too-small sample of simulated galaxies is likely to yield no plausible $g$ at
all, or equivalently, that $P(G_{\rm lens}\,|\,g)$ might be vanishingly small
for all $g$, no matter what the orientations and sources. However, Denzel et
al. (2021a) find that plausible matches for observed lensing galaxies can be
found in (at least) the EAGLE simulations, and it is possible to compute
$P(F\,|\,{\\{G\\}})$ for some different feedback prescriptions. The next stage
would be to extend to stellar population and kinematics, using

$\int\\!P(G_{\rm kin}\,|\,g,\iota)\,P(G_{\rm light}\,|\,g,\iota)\,P(G_{\rm
lens}\,|\,g,\iota)\,d\iota$ (33)

to compute $P(G_{\rm kin},G_{\rm light},G_{\rm lens}\,|\,g)$ and hence
$P(F\,|\,{\\{G\\}})$ using multiple data sets.

##  7 Concluding remarks

In this review article, we have provided a review of the applications of
galaxy-scale strong lensing in astrophysics and cosmology. Inevitably, some
special sub-topics within the field of galaxy-scale strong lensing have
evolved into proper research fields, having acquired a methodology and
literature extensive enough to warrant a dedicated review. These are:
detecting dark matter substructures and linking their properties to the dark
matter particle, and measuring $H_{\mathrm{0}}$ and other cosmological
parameters through time delays, examined in Vegetti et al. (in preparation)
and Birrer et al. (2022) respectively.

We started with a brief historical overview in Sect. 2. Then, we have
discussed both strong-lensing observables and complementary non-strong-lensing
observables and methodologies to model and extract meaningful results from
such data in Sect. 3. The most available and informative data for galaxy-scale
strong lenses come from imaging. We have reviewed the most common modeling
methods found in the literature to model such data and constrain galaxy
properties from lensed arcs, and in some cases with the inclusion of multiple
images of a point-like source. Next, we have reviewed the main scientific
results from the literature on the astrophysics of galaxies in Sect. 4 and to
cosmological models in Sect. 5. We then discussed the currently open questions
and provided future outlook in Sect. 6.

The open questions presented in Sect. 6 provide exciting opportunities for the
near future. Several large-area sky surveys – namely the Vera Rubin, Euclid,
and Roman observatories – will discover thousands of new galaxy-scale lensing
systems (Oguri and Marshall, 2010; Collett, 2015). Various automated pipelines
are being developed to tackle the computational aspect of modeling such large
lens samples (Chan et al., 2015; Nightingale et al., 2018; Shajib et al.,
2021; Schmidt et al., 2022). These treasure troves of data will provide the
necessary statistical power to shed light on the open questions on galaxy
evolution and cosmology.

###### Acknowledgements.

We thank the International Space Science Institute in Bern (ISSI) for their
hospitality and the conveners for organizing the stimulating workshop on
$``$Strong Gravitational Lensing$"$. Support for this work was provided by
NASA through the NASA Hubble Fellowship grant HST-HF2-51492 awarded to AJS by
the Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., for NASA, under contract
NAS5-26555. This work was also supported by the U.S. Department of Energy
(DOE) Office of Science Distinguished Scientist Fellow Program. GV has
received funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
897124. This research was made possible by the generosity of Eric and Wendy
Schmidt by recommendation of the Schmidt Futures program. VM acknowledges
partial support from Centro de Astrofísica de Valparaíso and project
"Fortalecimiento del Sistema de Investigación e Innovación de la Universidad
de Valparaíso" (UVA20993).

## References

  * Abadi et al. (2010) M.G. Abadi, J.F. Navarro, M. Fardal, A. Babul, M. Steinmetz, Galaxy-induced transformation of dark matter haloes. Mon. Notices of the Royal Astron. Soc. 407(1), 435–446 (2010). doi:10.1111/j.1365-2966.2010.16912.x
  * Adam et al. (2022) A. Adam, L. Perreault-Levasseur, Y. Hezaveh, Pixelated Reconstruction of Gravitational Lenses using Recurrent Inference Machines. arXiv e-prints, 2207–01073 (2022)
  * Adams et al. (2007) F.C. Adams, A.M. Bloch, S.C. Butler, J.M. Druce, J.A. Ketchum, Orbital instabilities in a triaxial cusp potential. Astrophys. J. 670(2), 1027–1047 (2007). doi:10.1086/522581
  * Aghanim et al. (2018) N. Aghanim, et al., Planck 2018 results. VI. Cosmological parameters (2018)
  * Agnello et al. (2018) A. Agnello, H. Lin, N. Kuropatkin, E. Buckley-Geer, T. Anguita, P.L. Schechter, T. Morishita, V. Motta, K. Rojas, T. Treu, A. Amara, M.W. Auger, F. Courbin, C.D. Fassnacht, J. Frieman, A. More, P.J. Marshall, R.G. McMahon, G. Meylan, S.H. Suyu, K. Glazebrook, N. Morgan, B. Nord, T.M.C. Abbott, F.B. Abdalla, J. Annis, K. Bechtol, A. Benoit-Lévy, E. Bertin, R.A. Bernstein, D. Brooks, D.L. Burke, A.C. Rosell, J. Carretero, C.E. Cunha, C.B. D’Andrea, L.N. da Costa, S. Desai, A. Drlica-Wagner, T.F. Eifler, B. Flaugher, J. García-Bellido, E. Gaztanaga, D.W. Gerdes, D. Gruen, R.A. Gruendl, J. Gschwend, G. Gutierrez, K. Honscheid, D.J. James, K. Kuehn, O. Lahav, M. Lima, M.A.G. Maia, M. March, F. Menanteau, R. Miquel, R.L.C. Ogando, A.A. Plazas, E. Sanchez, V. Scarpine, R. Schindler, M. Schubnell, I. Sevilla-Noarbe, M. Smith, M. Soares-Santos, F. Sobreira, E. Suchyta, M.E.C. Swanson, G. Tarle, D. Tucker, R. Wechsler, DES meets Gaia: discovery of strongly lensed quasars from a multiplet search. Mon. Notices of the Royal Astron. Soc. 479(4), 4345–4354 (2018). doi:10.1093/mnras/sty1419
  * Agnello et al. (2018) A. Agnello, H. Lin, N. Kuropatkin, E. Buckley-Geer, T. Anguita, P.L. Schechter, T. Morishita, V. Motta, K. Rojas, T. Treu, A. Amara, M.W. Auger, F. Courbin, C.D. Fassnacht, J. Frieman, A. More, P.J. Marshall, R.G. McMahon, G. Meylan, S.H. Suyu, K. Glazebrook, N. Morgan, B. Nord, T.M.C. Abbott, F.B. Abdalla, J. Annis, K. Bechtol, A. Benoit-LÃ©vy, E. Bertin, R.A. Bernstein, D. Brooks, D.L. Burke, A. Carnero Rosell, J. Carretero, C.E. Cunha, C.B. D’Andrea, L.N. da Costa, S. Desai, A. Drlica-Wagner, T.F. Eifler, B. Flaughe, J. GarcÃa-Bellido, E. Gaztanaga, D.W. Gerdes, D. Gruen, R.A. Gruendl, J. Gschwend, G. Gutierrez, K. Honscheid, D.J. James, K. Kuehn, O. Lahav, M. Lima, M.A.G. Maia, M. March, F. Menanteau, R. Miquel, R.L.C. Ogando, A.A. Plazas, E. Sanchez, V. Scarpine, R. Schindler, M. Schubnel, I. Sevilla-Noarbe, M. Smith, M. Soares-Santos, F. Sobreira, E. Suchyta, M.E.C. Swanson, G. Tarle, D. Tucker, R. Wechsler, Des meets gaia: discovery of strongly lensed quasars from a multiplet search. Mon. Notices of the Royal Astron. Soc. 479(4), 4345–4354 (2018). doi:10.1093/mnras/sty1419. http://adsabs.harvard.edu/abs/2018MNRAS.479.4345A
  * Amante et al. (2020) M.H. Amante, J. Magaña, V. Motta, M.A. García-Aspeitia, T. Verdugo, Testing dark energy models with a new sample of strong-lensing systems. Mon. Notices of the Royal Astron. Soc. 498(4), 6013–6033 (2020). doi:10.1093/mnras/staa2760
  * An (2005) J.H. An, Gravitational lens under perturbations: symmetry of perturbing potentials with invariant caustics. Mon. Notices of the Royal Astron. Soc. 356(4), 1409–1428 (2005). doi:10.1111/j.1365-2966.2004.08581.x
  * Anau Montel et al. (2022) N. Anau Montel, A. Coogan, C. Correa, K. Karchev, C. Weniger, Estimating the warm dark matter mass from strong lensing images with truncated marginal neural ratio estimation. arXiv e-prints, 2205–09126 (2022)
  * Auger et al. (2009) M.W. Auger, T. Treu, A.S. Bolton, R. Gavazzi, L.V.E. Koopmans, P.J. Marshall, K. Bundy, L.A. Moustakas, The Sloan Lens ACS Survey. IX. Colors, Lensing, and Stellar Masses of Early-Type Galaxies. Astrophys. J. 705(2), 1099–1115 (2009). doi:10.1088/0004-637X/705/2/1099
  * Auger et al. (2010) M.W. Auger, T. Treu, R. Gavazzi, A.S. Bolton, L.V.E. Koopmans, P.J. Marshall, Dark Matter Contraction and the Stellar Content of Massive Early-type Galaxies: Disfavoring “Light” Initial Mass Functions. Astrophys. J. Letters 721(2), 163–167 (2010). doi:10.1088/2041-8205/721/2/L163
  * Auger et al. (2010) M.W. Auger, T. Treu, A.S. Bolton, R. Gavazzi, L.V.E. Koopmans, P.J. Marshall, L.A. Moustakas, S. Burles, The sloan lens acs survey. x. stellar, dynamical, and total mass correlations of massive early-type galaxies. Astrophys. J. 724, 511–525 (2010). doi:10.1088/0004-637X/724/1/511
  * Barbosa et al. (2021a) C.E. Barbosa, C. Spiniello, M. Arnaboldi, L. Coccato, M. Hilker, T. Richtler, A preserved high-z compact progenitor in the heart of NGC 3311 revealed with MUSE 2D stellar population analysis. Astron. Astrophys. 649, 93 (2021a). doi:10.1051/0004-6361/202039809
  * Barbosa et al. (2021b) C.E. Barbosa, C. Spiniello, M. Arnaboldi, L. Coccato, M. Hilker, T. Richtler, What does (not) drive the variation of the low-mass end of the stellar initial mass function of early-type galaxies. Astron. Astrophys. 645, 1 (2021b). doi:10.1051/0004-6361/202039810
  * Barkana (1998) R. Barkana, Fast Calculation of a Family of Elliptical Mass Gravitational Lens Models. Astrophys. J. 502(2), 531–537 (1998). doi:10.1086/305950
  * Barnabè and Koopmans (2007) M. Barnabè, L.V.E. Koopmans, A Unifying Framework for Self-consistent Gravitational Lensing and Stellar Dynamics Analyses of Early-Type Galaxies. Astrophys. J. 666(2), 726–746 (2007). doi:10.1086/520495
  * Barnabè et al. (2009) M. Barnabè, O. Czoske, L.V.E. Koopmans, T. Treu, A.S. Bolton, R. Gavazzi, Two-dimensional kinematics of SLACS lenses - II. Combined lensing and dynamics analysis of early-type galaxies at z = 0.08-0.33. Mon. Notices of the Royal Astron. Soc. 399(1), 21–36 (2009). doi:10.1111/j.1365-2966.2009.14941.x
  * Barnabè et al. (2011) M. Barnabè, O. Czoske, L.V.E. Koopmans, T. Treu, A.S. Bolton, Two-dimensional kinematics of SLACS lenses - III. Mass structure and dynamics of early-type lens galaxies beyond z $\simeq$ 0.1. Mon. Notices of the Royal Astron. Soc. 415(3), 2215–2232 (2011). doi:10.1111/j.1365-2966.2011.18842.x
  * Barnabè et al. (2013) M. Barnabè, C. Spiniello, L.V.E. Koopmans, S.C. Trager, O. Czoske, T. Treu, A low-mass cut-off near the hydrogen burning limit for Salpeter-like initial mass functions in early-type galaxies. Mon. Notices of the Royal Astron. Soc. 436(1), 253–258 (2013). doi:10.1093/mnras/stt1727
  * Barrera et al. (2021) B. Barrera, L.L.R. Williams, J.P. Coles, P. Denzel, Bridging the Gap Between Simply Parametrized and Free-Form Pixelated Models of Galaxy Lenses: The Case of WFI 2033-4723 Quad. arXiv e-prints, 2108–04348 (2021)