11institutetext:  Department of Astronomy, School of Physics and Astronomy,
Shanghai Jiao Tong University, Shanghai 200240, China
11email<EMAIL_ADDRESS>22institutetext: Leiden Observatory, Leiden
University, P.O. Box 9513, 2300 RA Leiden, The Netherlands33institutetext:
Alma Mater Studiorum - Università di Bologna, Dipartimento di Fisica e
Astronomia ”Augusto Righi”, Via Gobetti 93/2, Bologna, Italy44institutetext:
CNRS and CNES, Laboratoire d’Astrophysique de Marseille, Aix-Marseille
Université 38, Rue Frédéric Joliot-Curie, F-13388 Marseille,
France55institutetext: Department of Astronomy & Astrophysics, University of
Chicago, Chicago, IL 60637, USA66institutetext: Kavli Institute for
Cosmological Physics, University of Chicago, Chicago, IL 60637,
USA77institutetext: NHFP Einstein Fellow88institutetext: Centre for
Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn
3122, Australia

# Strong lensing selection effects

Alessandro Sonnenfeld Strong lensing selection effectsStrong lensing selection
effectsStrong lensing selection effectsStrong lensing selection effects Shun-
Sheng Li Strong lensing selection effectsStrong lensing selection effects
Giulia Despali Strong lensing selection effectsStrong lensing selection
effects Raphael Gavazzi Strong lensing selection effectsStrong lensing
selection effects Anowar J. Shajib Strong lensing selection effectsStrong
lensing selection effectsStrong lensing selection effectsStrong lensing
selection effectsStrong lensing selection effectsStrong lensing selection
effects Edward N. Taylor Strong lensing selection effectsStrong lensing
selection effects

###### Abstract

Context.  Strong lenses are a biased subset of the general population of
galaxies.

Aims.  The goal of this work is to quantify how lens galaxies and lensed
sources differ from their parent distribution, namely the strong lensing bias.

Methods.  We first studied how the strong lensing cross-section varies as a
function of lens and source properties. Then, we simulated strong lensing
surveys with data similar to that expected for Euclid and measured the strong
lensing bias in different scenarios. We focused particularly on two
quantities: the stellar population synthesis mismatch parameter,
$\alpha_{\mathrm{sps}}$, defined as the ratio between the true stellar mass of
a galaxy and the stellar mass obtained from photometry, and the central dark
matter mass at fixed stellar mass and size.

Results.  Strong lens galaxies are biased towards higher stellar masses,
smaller half-mass radii, and higher dark matter masses. The amplitude of the
bias depends on the intrinsic scatter in the mass-related parameters of the
galaxy population and on the completeness in Einstein radius of the lens
sample. For values of the scatter that are consistent with observed scaling
relations and a minimum detectable Einstein radius of $0.5^{\prime\prime}$,
the strong lensing bias in $\alpha_{\mathrm{sps}}$ is $10\%$, while that in
the central dark matter mass is $5\%$. The bias has little dependence on the
properties of the source population: samples of galaxy-galaxy lenses and
galaxy-quasar lenses that probe the same Einstein radius distribution are
biased in a very similar way.

Conclusions.  Given current uncertainties, strong lensing observations can be
used directly to improve our current knowledge of the inner structure of
galaxies, without the need to correct for selection effects. Time-delay
measurements of $H_{0}$ from lensed quasars can take advantage of prior
information obtained from galaxy-galaxy lenses with similar Einstein radii.

###### Key Words.:

Gravitational lensing: strong

##  1 Introduction

Strong gravitational lensing is a powerful tool for studying galaxies and
cosmology. Strong lenses have been used to probe the mass structure of massive
galaxies (Auger et al., 2010; Oguri et al., 2014; Sonnenfeld et al., 2015;
Shajib et al., 2021), to detect substructure (Vegetti et al., 2012; Hezaveh et
al., 2016; Nierenberg et al., 2020), to carry out detailed studies of
magnified star-forming galaxies (Jones et al., 2013), and to measure the
expansion rate of the Universe with time delays (see Treu & Marshall, 2016 for
a review).

Strong lenses, however, are a biased subset of the general population of
galaxies and background sources. In order for a lens-source pair to be
included in a strong lensing survey, the following conditions must be met: (1)
the lens must produce multiple images of the source, (2) the multiple images
must be detected and resolved, and (3) the system must be recognised as a lens
by the survey.

Condition 1 poses constraints on the mass distribution of the lens, as well as
on the geometry of the system. First of all, not all lenses are capable of
producing multiple images, only those whose surface mass density,
$\Sigma(\boldsymbol{\theta})$, is higher than the critical surface mass
density for lensing, $\Sigma_{\mathrm{cr}}$, at at least one position,
$\boldsymbol{\theta}$ (Schneider et al., 1992). This condition excludes
objects with a very diffuse mass distribution from the population of lenses.
Second, objects with a higher concentration of mass allow for a larger portion
of the source plane to be multiply imaged and are therefore over-represented
in lens samples. Third, with fixed lens properties, sources that are located
farther away along the line of sight are more likely to be lensed.

Condition 2 sets limits on the observable properties of the sources and,
again, the lens mass distribution, in a way that depends on the survey
characteristics. On the one hand, brighter sources are more likely to be
detected when lensed. On the other hand, lensing magnification allows the
detection of fainter objects with respect to the field. However, if the
angular separation between the multiple images is too small, the system may
not be classified as a lens. Since the image separation increases with lens
mass, this means that the lens sample is biased against low-mass objects.

Condition 3 can introduce additional biases, depending on the efficiency of
the lens finding process on which the survey is based. Most commonly, lenses
with a smaller image separation are more difficult to identify even if
detected and resolved, because of contamination from the lens galaxy. This
introduces an additional bias towards lenses with a higher or more
concentrated mass.

The cumulative effect of these conditions can be summarised as follows. The
probability distribution, ${\rm P}_{\mathrm{{SL}}}$, of a sample of strong
lenses from a given survey with selection criterion $S$ is given by
(Sonnenfeld, 2022)

${\rm
P}_{\mathrm{{SL}}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|S)\propto{\rm
P}_{\mathrm{g}}(\boldsymbol{\psi}_{\mathrm{g}}){\rm
P}_{\mathrm{s}}(\boldsymbol{\psi}_{\mathrm{s}}){\rm
P}_{\mathrm{sel}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|S).$
(1)

In this equation, $\boldsymbol{\psi}_{\mathrm{g}}$ is the set of parameters
describing the properties of foreground galaxies that are relevant for
lensing, such as their redshift and mass distribution;
$\boldsymbol{\psi}_{\mathrm{s}}$ is the set of parameters describing
background sources; ${\rm P}_{\mathrm{g}}$ and ${\rm P}_{\mathrm{s}}$ describe
the parent distribution of foreground galaxies and background sources in the
absence of lensing; and ${\rm P}_{\mathrm{sel}}$ describes the probability of
selecting a lens-source system with parameters
$\boldsymbol{\psi}_{\mathrm{g}}$ and $\boldsymbol{\psi}_{\mathrm{s}}$ given
the criterion, $S$, used to define a lens. This last factor takes both
physical and survey selection effects into account, that is, whether a lens
with parameters $\boldsymbol{\psi}_{\mathrm{g}}$ produces a strongly lensed
image of a source with parameters $\boldsymbol{\psi}_{\mathrm{s}}$, and
whether such an image can be detected and recognised as a strong lens.

The information needed for the left-hand side of Eq. 1 is directly accessible
from strong lensing observations. If the main goal of a lensing survey is to
characterise the properties of the strong lens population, then it can be
accomplished by directly analysing this term. For many applications of strong
lensing, however, the aim is to constrain the properties of the general galaxy
or source population, ${\rm P}_{\mathrm{g}}$ and ${\rm P}_{\mathrm{s}}$, which
are coupled in a non-trivial way via the lens selection probability, ${\rm
P}_{\mathrm{sel}}$. In order to obtain an unbiased estimate of either ${\rm
P}_{\mathrm{g}}$ or ${\rm P}_{\mathrm{s}}$, it is necessary to invert Eq. 1.
In principle, this can be done with a Bayesian hierarchical formalism
(Sonnenfeld, 2022), but knowledge of the lens selection probability, ${\rm
P}_{\mathrm{sel}}$, is required. This factor can be written as the following
product:

${\rm
P}_{\mathrm{sel}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|S)={\rm
P}_{\mathrm{det}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}){\rm
P}_{\mathrm{find}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|S),$
(2)

where ${\rm P}_{\mathrm{det}}$ is the probability of detecting a strong
lensing event, and ${\rm
P}_{\mathrm{find}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}|S,\mathrm{det})$
is the probability of correctly classifying it as such111Sonnenfeld (2022)
implicitly assumed ${\rm P}_{\mathrm{find}}\equiv 1$; therefore, ${\rm
P}_{\mathrm{det}}$ and ${\rm P}_{\mathrm{sel}}$ can be used interchangeably in
the context of that work.. The detection probability, ${\rm
P}_{\mathrm{det}}$, can be obtained via simulations. The main challenge is
characterising ${\rm P}_{\mathrm{find}}$: in most of the existing strong
lensing surveys, the process of determining whether a system is included in a
strong lens sample is typically a combination of several cuts, usually
involving a non-trivial visual selection step.

For the above reasons, the problem of inverting Eq. 1 is a difficult one to
tackle exactly. A few studies have attempted to explicitly account for strong
lensing selection effects, usually by making ad hoc simplifying assumptions
(Sonnenfeld et al., 2015; Oldham et al., 2017; Sonnenfeld et al., 2019).
Whether it is necessary to invert Eq. 1, however, depends on the severity of
the strong lensing bias that needs to be corrected and on the accuracy
requirements for the key quantities of interest.

In this paper we aim to quantify the strength of the strong lensing bias on a
series of foreground galaxy and background source parameters. In particular,
we aim to determine how strong lenses differ from the parent population of
foreground galaxies and background sources in terms of (a) the radial mass
structure of the lenses (i.e. their stellar and dark matter mass density
profiles), (b) the ellipticity of the lenses, and (c) the size-magnitude
distribution of the lensed sources. This determination depends on (1) how the
lens detection probability, ${\rm P}_{\mathrm{det}}$, varies as a function of
galaxy and source properties, (2) the efficiency of a survey in correctly
classifying detected strong lenses (i.e. ${\rm P}_{\mathrm{find}}$), and (3)
the shape of the galaxy and source parameter distribution, ${\rm
P}_{\mathrm{g}}$ and ${\rm P}_{\mathrm{s}}$. To understand this third point we
can imagine the limiting case in which both ${\rm P}_{\mathrm{g}}$ and ${\rm
P}_{\mathrm{s}}$ are Dirac delta functions (i.e. all lenses and sources are
identical): in this limit, ${\rm P}_{\mathrm{{SL}}}$ simply reduces to the
product ${\rm P}_{\mathrm{g}}{\rm P}_{\mathrm{s}}$ up to a scaling constant;
this corresponds to a case in which there is no lensing bias.

Mandelbaum et al. (2009) carried out a thorough study of point (1): they
quantified how the properties of a lens galaxy determine its probability of
creating a lensing event with a point source. In this work we revisit the
Mandelbaum et al. (2009) study, expanding it to the extended source case. We
simulated individual lenses and examined how the lens detection probability
varies with lens and source properties. In addition, we address points (2) and
(3): we simulated large populations of strong lenses using empirical models,
we simulated the lens detection and finding phase, and we quantified the
lensing bias under various scenarios. We explored how the results change as a
function of the efficiency of a lens survey in discovering small image
separation lenses, and as a function of the scatter in mass parameters at
fixed light.

Finally, we considered galaxy-quasar lenses. The strong lensing bias of lensed
quasars was studied by Collett & Cunnington (2016) in the context of power-law
lens mass models with external convergence. Here we address the question of
how different galaxy-quasar lenses are from galaxy-galaxy lens samples. This
point is relevant for time-delay cosmography studies, in which measurements of
the time delay between the multiple images of a strongly lensed quasar are
used to constrain the expansion rate of the Universe. Galaxy-galaxy lenses can
in principle be used to help break some of the model degeneracies affecting
these measurements (Birrer et al., 2020; Birrer & Treu, 2021), but any
difference between the two lens classes can introduce biases if not corrected
for. With this study we aim to quantify this bias.

The structure of this paper is as follows. In Sect. 2 we introduce the basics
of gravitational lensing. In Sect. 3 we study individual lens systems. In
Sect. 4 we describe our simulations of lens surveys. In Sect. 5 we show the
results of our analysis of the simulated lens survey data. We discuss the
results in Sect. 6 and draw conclusions in Sect. 7.

The Python code used for the simulation and analysis of the lens sample can be
found in a dedicated section of a GitHub
repository222https://github.com/astrosonnen/strong_lensing_tools/tree/main/papers/selection_effects.

##  2 Strong lensing basics

The lensing properties of an object with respect to a source depend solely on
its dimensionless surface mass density distribution,
$\kappa(\boldsymbol{\theta})$ (also referred to as the convergence). This is
the ratio between the surface mass density and the critical surface mass
density for lensing:

$\kappa(\boldsymbol{\theta})=\dfrac{\Sigma(\boldsymbol{\theta})}{\Sigma_{\mathrm{cr}}}.$
(3)

The latter quantity is defined as

$\Sigma_{\mathrm{cr}}=\dfrac{c^{2}D_{\mathrm{s}}}{4\pi
GD_{\mathrm{d}}D_{\mathrm{ds}}},$ (4)

where $c$ is the speed of light, $G$ the gravitational constant, and
$D_{\mathrm{d}}$, $D_{\mathrm{s}}$, and $D_{\mathrm{ds}}$ are the angular
diameter distances between the observer and the lens, the observer and the
source, and the lens and the source, respectively.

Given a source at angular position $\boldsymbol{\beta}$, images of it form at
the positions $\boldsymbol{\theta}$ that are solutions of the lens equation

$\boldsymbol{\beta}=\boldsymbol{\theta}-\boldsymbol{\alpha}(\boldsymbol{\theta}),$
(5)

where $\boldsymbol{\alpha}$ is the deflection angle of the lens. This can be
expressed in terms of the dimensionless surface mass density by means of the
following integral over the whole sky:

$\boldsymbol{\alpha}(\boldsymbol{\theta})=\frac{1}{\pi}\int_{\mathbb{R}^{2}}d^{2}\boldsymbol{\theta}^{\prime}\kappa(\boldsymbol{\theta}^{\prime})\dfrac{\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}}{|\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}|^{2}}.$
(6)

The images of the background source are in general magnified in total flux and
in size, while preserving the original surface brightness.

###  2.1 The axisymmetric lens

In the special case of an axisymmetric lens we can simplify the notation by
considering a single coordinate axis with origin at the centre of the lens. We
label the components of the image position, source position and deflection
angle along this axis as $\theta$, $\beta$, and $\alpha$, respectively. The
lens equation for an axisymmetric lens then becomes

$\beta=\theta-\alpha(\theta),$ (7)

and the expression for the deflection angle reduces to

$\alpha(\theta)=\frac{2}{\theta}\int_{0}^{\theta}d\theta^{\prime}\kappa(\theta^{\prime})\theta^{\prime}.$
(8)

This can also be expressed in terms of the total projected mass enclosed
within a circle with angular radius equal to $\theta$:

$\alpha(\theta)=\frac{1}{\pi\theta}\dfrac{M(<\theta)}{\Sigma_{\mathrm{cr}}D_{\mathrm{d}}^{2}}.$
(9)

A very important quantity for describing the strength of a strong lens is the
Einstein radius, $\theta_{\mathrm{Ein}}$. For an axisymmetric lens, this is
the radius corresponding to the solution of Eq. 7 for a source placed at the
same angular position as the lens centre ($\beta=0$). The circle with radius
equal to $\theta_{\mathrm{Ein}}$ is known as the tangential critical curve.
Images that form there have infinite magnification in the tangential
direction. It can be shown that $\theta_{\mathrm{Ein}}$ satisfies the
following condition,

$\bar{\kappa}(<\theta_{\mathrm{Ein}})=1,$ (10)

where $\bar{\kappa}(<\theta)$ is the average surface mass density within a
radius equal to $\theta$:

$\bar{\kappa}(<\theta)\equiv\frac{2}{\theta^{2}}\int_{0}^{\theta}d\theta^{\prime}\kappa(\theta^{\prime})\theta^{\prime}.$
(11)

Axisymmetric lenses of the kind considered in this work typically produce
either one or three images of a point source. Figure 1 helps visualise this
property. Plotted in Figure 1 is the quantity $\theta-\alpha(\theta)$ as a
function of the position in the image plane $\theta$, for a few lens models.
According to Eq. 7, images form at the locations where this quantity equals
the position of the source. Therefore, given a source position $\beta$, the
number of images and their location can be determined by drawing a horizontal
line at the value $\beta$ on the vertical axis, and finding the points where
this line intersects the $\theta-\alpha(\theta)$ curve.

For small values of $\beta$ (for sources close to the lens centre), the number
of images that are produced is three: the source is strongly lensed. For large
values of $\beta$, instead, only one image is formed. The value of $\beta$
where the transition occurs is known as the radial caustic, which is marked by
the horizontal dotted line in Figure 1. As can be seen from Figure 1, a source
at this location is mapped to the stationary point of the
$\theta-\alpha(\theta)$ curve. That location on the image plane is known as
the radial critical curve. Images that form there have a formally infinite
magnification in the radial direction333The slope of the
$\theta-\alpha(\theta)$ curve is the inverse of the radial magnification. This
can be understood by taking the derivative of Eq. 7 with respect to $\theta$..

Figure 1:  Lens equation of axisymmetric lenses. Solid lines represent the
right-hand side of Eq. 7 for five lenses with different density profiles.
Given a source at position $\beta$, its lensed images form at the values of
$\theta$ where the solid line intersects a horizontal line located at $\beta$
on the vertical axis. Dotted lines represent positions of the radial caustics.
Sources located within the radial caustic of a given lens produce three
images. The lens models used in this simulation consist of a stellar component
and a dark matter halo, as described in Sect. 3.1. Their density profile is
plotted in Fig. 4.

###  2.2 The elliptical lens

In this paper we focus mostly on lens galaxies with elliptical isodensity
contours. Given a surface mass density profile $\Sigma(R)$, a lens with an
elliptical mass distribution can be obtained by replacing the radial
coordinate with the circularised radius:

$R\rightarrow\sqrt{qx^{2}+\frac{y^{2}}{q}},$ (12)

where $x$ and $y$ are Cartesian axes centred on the lens centre, with $x$
pointing towards the major axis, and where $q$ is the minor-to-major axis
ratio.

Figure 2 shows the source-plane caustics of elliptical lenses with different
values of the axis ratio. We used the software Glafic (Oguri, 2021) to obtain
the caustic curves. The outermost curves are radial caustics, while the inner
ones are tangential ones. The most striking difference with respect to the
axisymmetric case (blue curves in Figure 2) is the fact that the tangential
caustic is transformed from a point into a diamond-like curve. Sources located
within the diamond produce five images (one of which is usually highly de-
magnified). Sources that lie in the region enclosed between the two caustics
are imaged three times. Sources outside the radial caustic are imaged only
once, as in the axisymmetric case. The fact that the number of images changes
by two at a caustic crossing is a general feature of gravitational lenses with
no singularities (Schneider et al., 1992).

Figure 2:  Caustics of lenses with fixed radial structure and different
ellipticity. Source-plane angular coordinates are in units of the Einstein
radius. The outer curve is the radial caustic, while the inner diamond is the
tangential caustic. Point sources located outside the caustic are not strongly
lensed. Sources that lie in the region enclosed by the two caustics produce
three images, while sources inside the tangential caustic produce five images.
The lens model adopted for this experiment consists of a stellar component and
a dark matter halo, as described in Sect. 3.1. The two components have the
same ellipticity.

###  2.3 Lensing event definition: Point sources

In order to compute the probability of a lensing event we must provide an
exact definition for it. A necessary condition for a lens-source system to
qualify as a strong lens is the presence of multiple images of the source. As
we showed above, this requires the source to lie within the region enclosed by
the radial caustic. In order to recognise a strong lens in a real survey,
however, it is not sufficient for multiple images to exist: they must be
detected. For this reason, given the detection limit for a point source,
$m_{\mathrm{lim}}$, we define as strong lensing event any lens-source system
with at least two images brighter than $m_{\mathrm{lim}}$. This is the same
definition used by (van de Ven et al., 2009), upon which the Mandelbaum et al.
(2009) work is based. Labelling with $m_{2}$ the magnitude of the second-
brightest image, then, in the absence of photometric noise the lens detection
probability ${\rm P}_{\mathrm{det}}$ becomes

${\rm
P}_{\mathrm{det}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}})=\left\\{\begin{array}[]{ll}1&\rm{if}\,m_{2}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}})<m_{\mathrm{lim}}\\\
0&\rm{otherwise}\end{array}\right..$ (13)

###  2.4 Lensing event definition: Extended sources

Defining a strong lensing event in the case of an extended source is less
straightforward. In principle, we should require parts of the source to be
multiply imaged. In practice, it is not always easy to determine whether a
lensed image contains multiple images or not. This is because, when the source
is extended, some of the images can be blended together. In real strong
lensing surveys, it is common to find lens candidates in which the lensed
source consists of only a single arc. In those cases it is difficult to
establish whether the arc is a set of blending images or not, and the decision
of including such systems in a strong lens sample is often arbitrary. Here we
adopt the following working definition: an extended source is strongly lensed
if the number of surface brightness peaks that are detected is larger than its
intrinsic number of peaks in the absence of lensing.

We explain how this definition applies in practice with a few examples. For
simplicity, we focus on the case of a source with a single surface brightness
peak, such as a Sérsic profile. All of the sources considered in this work
belong to this family of surface brightness models. We adopt the following
procedure to determine the number of detected peaks. Given the observed image
of a lensed source, we define its footprint as the ensemble of pixels where
the source is detected with $S/N>2$. The footprint is in general composed of
multiple disconnected regions, corresponding to the different images. In order
to only include images that can be clearly identified, we add the condition
that the integrated $S/N$ of each separate region must be $S/N>10$. This
condition has the effect of removing from the source footprint any isolated
region consisting only of a very small number of pixels. In a real-world
application, it would be very hard to classify such marginal detections as
images. If, after applying this cut, the source footprint is spread over
multiple separated regions, then the system is classified as a lens. If the
source footprint consists instead of a single region, we iteratively increase
the surface brightness threshold used to define the footprint and count the
number of isolated regions with $S/N>10$. The maximum number of detected
regions defines the number of detected peaks.

Figure 3 shows a few examples of how this criterion can be used to classify
lenses. The first column shows the caustic structure and source position. The
second column shows an image of the lensed source. The third column shows the
source footprint defined with the procedure described above. Pink pixels
correspond to the largest footprint that maximises the number of detected
images, while purple pixels belong to the $2\sigma$ detection footprint. In
the first, second and sixth example, these two footprints coincide. In the
third, fourth and fifth example, instead, the $2\sigma$ detection footprint
consists of a single region, but an increase of the surface brightness
threshold leads to the detection of additional images.

Strictly speaking, our lens definition criterion fails for a perfect Einstein
ring. In order to cover such a scenario, we also classify as strongly lensed
any sources that produce a footprint with a hole.

We point out that, although our lens definition relies on peak detections, it
is not necessary for the peak of the source surface brightness to lie within
the caustics in order for the source to be strongly lensed. This can be
observed in the fourth example of Fig. 3: although the source centroid lies
outside of the caustic (as shown in the first column), the outskirts of the
source overlap with the lens, and therefore multiple images are produced.

Figure 3:  Criterion used to classify lensed images of extended sources. Six
examples are shown. First column: Caustics (red curves) and source position
(blue circle). The radius of the circle corresponds to the radius at which the
surface brightness is equal to $2\sigma$ the sky background rms fluctuation
for a single pixel. In other words, the blue circle delimits the area of the
source that can be detected. Second column: Mock image of the lensed source,
with added noise. Third column: Footprint of the source. The purple footprint
is obtained with a $2\sigma$ detection criterion. The coral region is the
largest footprint with the highest number of detected images.

We can identify three different regimes in strong lensing of extended sources,
depending on the relative size of the source and of the lens caustics. In the
limit of a very small source, the image configurations that are produced are
qualitatively similar to those that can be obtained in the point-source case.
When the source and caustic size are comparable, the multiple images tend to
be blended into arcs or rings. Most, if not all, of the galaxy-scale strong
lenses known belong to these two categories. However, there is a third regime
in strong lensing, corresponding to the case in which the source size is much
bigger than the caustics, such as in the fifth example of Fig. 3. In this
regime, the overall size and total flux of the source are roughly preserved,
and the lens produces only a relatively minor perturbation on a localised
region of the image. Strong lenses of this kind can be difficult to detect,
especially if the region of the image subject to strong lensing overlaps with
the light from the foreground galaxy. The ultimate limiting factor to the
detection of strong lenses in the large source size regime, however, is the
ability to spatially resolve the multiple images. This limit is set by the
size of the point spread function (PSF).

All of the examples that we considered were based on sources with a single
surface brightness peak. In general this is not necessarily the case, and the
intrinsic number of surface brightness peaks of a strongly lensed source is
not known a priori. When dealing with a real lens candidate, applying our lens
definition criterion requires showing that the observed number of surface
brightness peaks can be reproduced with a lens model in which the source has a
smaller number of peaks. In samples of lenses defined via visual inspection
this process is typically done implicitly, by identifying multiply imaged
blobs that belong to the same source element.

To our knowledge, we are the first to propose a surface brightness peak-based
definition of a strong lensing event. A popular alternative definition of a
strong lens is one based on magnification: only images that are magnified by
more than a given threshold are considered as strongly lensed (see e.g.
Hilbert et al., 2007). The problem with a definition of this kind is that
magnification cannot be determined unambiguously from observations, unless the
intrinsic properties of the lensed source are known from other means (e.g. if
the source is a standard candle or a standard ruler). Because of the mass-
sheet degeneracy (Falco et al., 1985), it is possible to vary the
magnification of a lensed image while keeping its observed properties fixed:
this could lead to the paradox of two identical-looking lenses that are
classified differently on the basis of the underlying magnification. Although
we could still use a magnification-based definition for the sake of carrying
out our experiments, it would then be difficult to apply our results to real
data. Our definition of a strong lensing event, instead, is robust with
respect to the mass-sheet degeneracy.

##  3 Individual lenses

In this section we study how the probability of a strong lensing event varies
as a function of lens and source properties. In order to do so, it is useful
to introduce the concept of strong lensing cross-section. Given a foreground
galaxy with parameters $\boldsymbol{\psi}_{\mathrm{g}}$, a background source
with parameters $\boldsymbol{\psi}_{\mathrm{s}}$, and a criterion $S$ to
define a strong lensing event, the strong lensing cross-section is defined as
(Sonnenfeld, 2022)

$\sigma_{\mathrm{{SL}}}=\int_{\mathbb{R}^{2}}d\boldsymbol{\beta}{\rm
P}_{\mathrm{det}}(\boldsymbol{\psi}_{\mathrm{g}},\boldsymbol{\psi}_{\mathrm{s}}^{(-\boldsymbol{\beta})},\boldsymbol{\beta}|S),$
(14)

where $\boldsymbol{\beta}$ is the position of the background source,
$\boldsymbol{\psi}_{\mathrm{s}}^{(-\boldsymbol{\beta})}$ is the ensemble of
source parameter except for the position, and the integral is carried out over
the whole sky. The definition above is valid for both a point source and an
extended source: although there is no unique way of defining the position of
an extended source, the integral over the sky ensures that the result is
independent of how the source position is defined. In the limit of low density
of background sources, which is satisfied in all practical cases, the
probability of a lensing event is proportional to $\sigma_{\mathrm{{SL}}}$.
The lensing cross-section defined via Eq. 14 depends solely on the lens
detection probability ${\rm P}_{\mathrm{det}}$ and does not take into account
whether the lens can be correctly classified by a lens finder. This separate
selection step is captured by the term ${\rm P}_{\mathrm{find}}$ in Eq. 2. In
this section we consider exclusively the detectability of a lens, and
therefore focus only on $\sigma_{\mathrm{{SL}}}$. In Sect. 4, when considering
specific lens survey simulations, we introduce ${\rm P}_{\mathrm{find}}$.

We computed $\sigma_{\mathrm{{SL}}}$ in a series of different scenarios of
increasing complexity. The model family adopted to describe the radial density
profile of the lenses was the same in all of our experiments. We describe this
in Sect. 3.1. In Sect. 3.2 we show calculations of the strong lensing cross-
section in the case of axisymmetric lenses and point sources. In Sect. 3.3 we
generalise the lens geometry to elliptical, while in Sect. 3.4 we replace
point sources with extended sources.

###  3.1 Lens density profile

In this work we focus on massive early-type galaxies as lenses, as these make
up the vast majority of known strong lenses. We describe their mass
distribution with a model consisting of two concentric components, one
describing the baryons and one for the dark matter. We assume that the
baryonic component consists entirely of stars, thereby neglecting gas, which
is known to contribute very little to the mass of early-type galaxies in the
inner regions that are probed by strong lensing. We then assume that the stars
follow a Sérsic profile, with projected surface mass density given by

$\Sigma(R)=\Sigma_{0}\exp{\left\\{-b\left(\frac{R}{R_{\mathrm{e}}}\right)^{1/n}\right\\}}.$
(15)

In the above equation,

$\Sigma_{0}=\frac{M_{*}b^{2n}}{2\pi nR_{\mathrm{e}}^{2}\Gamma(2n)},$ (16)

$M_{*}$ is the total mass, $R_{\mathrm{e}}$ is the radius enclosing half of
the total mass, $n$ is the Sérsic index, $\Gamma$ is the incomplete Gamma
function, and $b$ is given by (Ciotti & Bertin, 1999)

$b(n)\approx 2n-\frac{1}{3}+\frac{4}{405n}+\frac{46}{25515n^{2}}+O(n^{-3}).$
(17)

Throughout this paper we indicate with $\theta_{\mathrm{e}}$ the angular size
of the half-light radius.

We fix the Sérsic index of the lenses to $n=4$, corresponding to a de
Vaucouleurs model. Although early-type galaxies are often described with a
free Sérsic index, a de Vaucouleurs profile is able to reproduce their surface
brightness distribution to a few percent in the radial range
$1\rm{kpc}<R<30\rm{kpc}$ (see e.g. Sonnenfeld, 2020), which is the region that
is most relevant for strong lensing. Finally, we assume that the light
distribution of the stellar component follows its mass distribution exactly.
That is, we do not allow for the presence of gradients in the stellar mass-to-
light ratio. In Sect. 6.5 we discuss qualitatively what the implications of
relaxing this assumption would be.

We describe the dark matter halo with a generalised Navarro, Frenk, and White
(gNFW) profile. We first define it by its three-dimensional distribution,
which for a spherically symmetric profile is

$\rho(r)=\dfrac{\rho_{0}}{(r/r_{s})^{\gamma_{\mathrm{DM}}}\left(1+r/r_{s}\right)^{3-\gamma_{\mathrm{DM}}}}.$
(18)

The parameter $\gamma_{\mathrm{DM}}$ is the inner density slope, $\rho_{0}$ is
a normalisation parameter, while $r_{s}$ is the scale radius. The logarithmic
slope of the density profile transitions from $-\gamma_{\mathrm{DM}}$ to $-3$
at a radius $r\approx r_{s}$. The projected surface mass density of a gNFW
profile can be expressed in terms of the following integral (Wyithe et al.,
2001):

$\Sigma(R)=2r_{s}\rho_{0}\left(\frac{R}{r_{s}}\right)^{1-\gamma_{\mathrm{DM}}}\int_{0}^{\pi/2}dx\sin{x}(\sin{x}+R/r_{s})^{\gamma_{\mathrm{DM}}-3}.$
(19)

Figure 4 shows the dimensionless surface mass density profile of Sérsic + gNFW
models with various values of the inner dark matter density slope
$\gamma_{\mathrm{DM}}$ and of the fraction of projected dark matter mass
within the half-light radius, $f_{\mathrm{DM}}$. All of the profiles have a
dark matter scale radius equal to ten times $R_{\mathrm{e}}$, and are
normalised in such a way that the Einstein radius is equal to the half-mass
radius of the stellar component (these assumptions are dropped later). Two
main features emerge from Fig. 4. First, the baryons generally dominate the
total density in the inner regions $(\theta<\theta_{\mathrm{Ein}})$, while the
dark matter is the main component at large radii. Second, models with
different dark matter fractions and inner dark matter slopes can conspire to
produce very similar total density profiles. This is the case, for example, of
the ($f_{\mathrm{DM}}=0.5,\gamma_{\mathrm{DM}}=1.5$) and the
($f_{\mathrm{DM}}=0.2,\gamma_{\mathrm{DM}}=2.0$) models (green and blue lines
in Fig. 4).

Figure 4:  Dimensionless surface mass density profile of Sérsic + gNFW
composite models. The Sérsic index of the baryonic component is fixed to
$n=4,$ and the scale radius of the dark matter component is fixed to ten times
the half-light radius. All of the profiles are normalised in such a way that
the Einstein radius is equal to the half-light radius. Solid lines show the
total density profile, dotted lines the dark matter density profile, and
dashed lines the baryonic density profile. The dashed blue, orange, and green
lines are identical, as they correspond to profiles with the same fraction of
baryonic mass within the half-light radius.

###  3.2 Axisymmetric lenses, point sources

Axisymmetric lenses with a density profile of the kind introduced in Sect. 3.1
can produce either one or three images of a point source. This can be seen in
Fig. 1, which shows the lens equation for various values of the dark matter
fraction and the inner dark matter slope. All of the lenses shown in this
figure have the same Einstein radius, which is equal in size to the half-mass
radius of the stellar component. The radius of the radial caustic, marked by
the dotted lines in Fig. 1, is a strong function of the lens properties: it is
largest in lenses with a smaller dark matter fraction or steeper dark matter
slopes. As a result, the source plane area that is subject to strong lensing
is an even stronger function of these properties, since it scales with the
square of the radial caustic radius. In order to compute the lensing cross-
section, however, we must take the magnification into account because it
determines whether multiple images can be detected or not.

Figure 5 shows the magnification of the secondary image as a function of
source position. The secondary image is the one located in the region between
the radial and tangential critical curves, opposite to the source with respect
to the lens centre. In most practical cases this is the second brightest
image. As Fig. 5 shows, the magnification is very large for sources close to
the lens centre (small values of $\beta$), decreases with increasing source
position and then increases in close proximity to the radial caustic. While
for the model with $f_{\mathrm{DM}}=0.8$ (purple curve) the magnification is
above unity everywhere, other lens models can produce highly de-magnified
secondary images for large values of $\beta$. Depending on the intrinsic
brightness of the lensed source, these images may or may not be detected.

Figure 5:  Magnification of the secondary image as a function of the source
position. The lenses are axisymmetric composite models. Their density profile
is plotted in Fig. 4. Vertical dotted lines mark the position of the radial
caustic.

Using the definition of Eq. 14, we computed the lensing cross-section of a set
of axisymmetric lenses, with respect to point sources with different
brightnesses. In particular, we considered model lenses with fixed Einstein
radius, with angular half-mass radius $\theta_{\mathrm{e}}$ fixed to
$\theta_{\mathrm{Ein}}$, and varying values of the dark matter fraction and
dark matter inner slope. The results are shown in Fig. 6. Each line
corresponds to a source with a given intrinsic (i.e. unlensed) magnitude,
$m_{s}$. The difference between this magnitude and the limiting magnitude of
the survey is indicated as follows:

$\Delta m\equiv m_{s}-m_{\mathrm{lim}}.$ (20)

We can see a clear trend between $\Delta m$ and the lensing cross-section, in
both panels of Fig. 6: $\sigma_{\mathrm{{SL}}}$ is larger for brighter
sources. The trend saturates below a certain $\Delta m$, for sufficiently
small values of $\gamma_{\mathrm{DM}}$ or for sufficiently large values of
$f_{\mathrm{DM}}$. In these cases the lensing cross-section coincides with the
full area enclosed within the caustic, and further increasing the source
brightness does not result in an increased value of $\sigma_{\mathrm{{SL}}}$.

Figure 6:  Strong lensing cross-section of an axisymmetric lens and a point
source. The lens galaxy is a composite model, introduced in Sect. 3.1, with an
angular half-light radius equal to the Einstein radius. The left panel shows
the effect of varying the slope $\gamma_{\mathrm{DM}}$ for fixed
$f_{\mathrm{DM}}$; conversely, the right panel shows variations as a function
of $f_{\mathrm{DM}}$ with $\gamma_{\mathrm{DM}}$ fixed. The system is defined
as a strong lens if at least two images are detected. Different lines
correspond to the difference, $\Delta m$, between the source magnitude and the
survey detection limit for a point source.

At fixed source brightness, trends with the dark matter inner slope or dark
matter fraction are generally weak. The lines of Fig. 6, however, have been
computed by keeping the Einstein radius fixed while varying
$\gamma_{\mathrm{DM}}$ or $f_{\mathrm{DM}}$. This is achieved by adjusting
other properties of the lens, such as the total mass of the baryonic or the
dark matter component (see Fig. 4). In practice, when varying one ingredient
of the lens density profile, the Einstein radius varies in response. To get a
complete view of how the lensing cross-section depends on different lens
properties, we also computed how $\sigma_{\mathrm{{SL}}}$ responds in absolute
terms by varying one lens parameter at a time. Figure 7 shows
$\sigma_{\mathrm{{SL}}}$ as a function of stellar mass, half-light radius,
inner dark matter slope, and projected dark matter mass enclosed within an
aperture of $5$ kpc, $M_{\mathrm{DM},5}$.

Figure 7:  Absolute value of the strong lensing cross-section as a function of
various lens properties. The reference lens is a galaxy at $z=0.3$ with
$\log{M_{*}}=11.5$, $R_{\mathrm{e}}=7\,{\rm kpc}$, $\gamma_{\mathrm{DM}}=1.5$,
$\log{M_{\mathrm{DM},5}}=11.0$, $r_{s}=100$ kpc, and a source at $z=1.5$ In
each panel, only one property of the lens is varied, as indicated on the label
of the horizontal axis. Each curve corresponds to a different value of $\Delta
m$, in accordance with Fig. 6. The dashed line in each panel shows the
quantity $\pi\theta_{\mathrm{Ein}}^{2}$.

The lensing cross-section increases with increasing stellar mass and dark
matter mass, decreases with increasing $R_{\mathrm{e}}$ for bright sources,
while is only a weak function of $\gamma_{\mathrm{DM}}$. The lack of a clear
trend between $\sigma_{\mathrm{{SL}}}$ and $\gamma_{\mathrm{DM}}$ appears to
be in contradiction with the result of Mandelbaum et al. (2009), who found a
strong positive correlation between $\sigma_{\mathrm{{SL}}}$ and
$\gamma_{\mathrm{DM}}$. The origin of this discrepancy lies in the different
ways in which $\gamma_{\mathrm{DM}}$ is varied in the two experiments. While
we varied $\gamma_{\mathrm{DM}}$ at fixed $M_{\mathrm{DM},5}$, Mandelbaum et
al. (2009) kept fixed the virial mass of the dark matter halo. At fixed virial
mass, increasing the inner dark matter slope results in a correspondingly
larger dark matter mass in the inner regions, which naturally results in a
larger lensing cross-section.

Figures 6 and 7 show that the trends between lens properties and the strong
lensing cross-section can be different for sources with different
brightnesses. The net effect in a strong lensing survey is an average over the
source population, weighted by the source luminosity function. This implies
that surveys that target different families of sources, with different
luminosity functions (for instance, galaxies or quasars), can have different
strong lensing biases. We explore this possibility in Sect. 5.

###  3.3 Elliptical lenses, point sources

We measured $\sigma_{\mathrm{{SL}}}$ for lenses with a fixed radial density
profile and different ellipticities, with respect to point sources of
different brightnesses. In particular, we set $f_{\mathrm{DM}}=0.5$,
$\gamma_{\mathrm{DM}}=1.5$, $r_{s}=10R_{\mathrm{e}}$,
$\theta_{\mathrm{Ein}}=\theta_{\mathrm{e}}$, and set the ellipticities of both
the baryonic and dark matter components to be the same, with the same
orientation of the major axis. This is the lens model used to produce the
caustics plot of Fig. 2. Figure 2 suggests that the size of the source plane
area subject to strong lensing, the area enclosed within the outermost
caustic, does not vary strongly with the ellipticity of the lens. Therefore,
we expect the strong lensing cross-section to be a weak function of
ellipticity.

We carried out the computation of $\sigma_{\mathrm{{SL}}}$ by means of
simulation: we generated a large number of point sources randomly distributed
over a given area that includes the caustic, then used Glafic to solve the
lens equation, find the number of images and their magnification. We then
measured the fraction of sources that are strongly lensed according to the
criterion of Eq. 13 and multiplied this value by the area over which sources
are located. The resulting $\sigma_{\mathrm{{SL}}}$ is shown in Fig. 8 as a
function of the minor-to-major axis ratio, $q$.

Figure 8:  Point source strong lensing cross-section as a function of lens
axis ratio. Solid lines show the cross-section based on the lens event
definition of Eq. 13 and dashed lines the cross-section based on the detection
of four images (quad cross-section). Lines of different colours correspond to
sources of different intrinsic magnitudes. The dashed blue line, corresponding
to the brightest source magnitude, overlaps completely with the dashed orange
line. The dashed purple line is zero: very faint sources cannot produce any
detectable quads. The parameters of the lens density profile are
$f_{\mathrm{DM}}=0.5$, $\gamma_{\mathrm{DM}}=1.5$, $r_{s}=10R_{\mathrm{e}}$,
and $\theta_{\mathrm{e}}=\theta_{\mathrm{Ein}}$. The baryonic and dark matter
components have the same ellipticity and direction of the major axis.

For bright sources the strong lensing cross-section is approximately constant
with axis ratio. This is because, as pointed out earlier, in the bright source
regime the cross-section is determined by the area enclosed within the radial
caustic, which does not vary much with lens ellipticity. For faint sources we
observe a larger variation with $q$, with a factor of two difference between
the largest and smallest value of $\sigma_{\mathrm{{SL}}}$ at fixed source
brightness.

Most of the sources that result in detectable lenses produce two detectable
images. These are sources that are located in the region enclosed between the
radial and the tangential caustic. If the source is located within the
tangential caustic, however, four detectable images are usually created.
Lenses with four visible images, usually referred to as quad lenses, are
sometimes given a high priority in certain lensing studies, because they offer
more constraints compared to double lenses. For instance, quads make up the
majority of the lenses used so far in time-delay studies (Millon et al.,
2020). For this reason we also computed an alternative lensing cross-section,
in which the definition of a lensing event requires the detection of four
images, instead of two. This is plotted in Fig. 8 with dashed lines. The
cross-section for quads is a strong function of lens ellipticity, for bright
sources. This is a consequence of the fact that the area enclosed within the
tangential caustic, which is where a source needs to be in order to produce
four images, increases with increased lens ellipticity, as Fig. 2 shows. For
sources that are intrinsically fainter than the detection limit, however, the
quad cross-section is extremely small, regardless of ellipticity.

###  3.4 Elliptical lenses, extended sources

In the case of an extended source, the complexity of the problem is increased
due to the addition of a series of features: the source surface brightness
distribution, with its radial profile, shape and orientation, and the PSF.
Moreover, as we discussed in Sect. 2.4, there are different regimes in strong
lensing of extended sources, depending on the relative size of the source and
the caustics of the lens. For this reason, we split our analysis into two
parts. First, we explore the small source size regime, where the source size
is comparable to or smaller than the lens caustics. Then, we consider cases in
which the source size is bigger than the lens caustics, which we refer to as
the large source regime.

####  3.4.1 Small source sizes

For the sake of reducing the dimensionality of the analysis, we focused on
circularly symmetric sources. We also fixed the surface brightness profile to
an exponential disk (i.e. a Sérsic model with $n=1$). We then took a lens
model with the same parameters used in Sect. 3.3 and a minor-to-major axis
ratio of $q=0.7$. We simulated a large number of images of extended sources
with Glafic and measured the fraction of them that results in a strong lens,
according to the definition of Sect. 2.4. For the small source size
experiment, we used pixels with a size equal to $1/20\theta_{\mathrm{Ein}}$
and convolved the images with a Moffat PSF with a full width at half maximum
(FWHM) of two pixels and a $\beta$ parameter of $5.0$. We also assumed that
the background noise is an uncorrelated Gaussian field. We carried out
experiments with sources with different values of the total unlensed flux,
$f$, and half-light radius, $\theta_{\mathrm{e,s}}$. The results are shown in
Fig. 9. The total flux of the source, indicated in the legend, is measured in
units of the sky background rms fluctuation within an area equal to the square
of the Einstein radius.

Figure 9:  Strong lensing cross-section of an extended source: small source
size regime ($\theta_{\mathrm{e,s}}<\theta_{\mathrm{Ein}}$). The cross-section
is plotted as a function of the ratio between the half-light radius of the
source and the Einstein radius of the lens. The lens model is fixed to be the
same as in Fig. 8, with axis ratio $q=0.7$. The source is a circular
exponential profile. Each solid line corresponds to a different value of the
total unlensed flux of the source. The flux, $f$, is expressed in terms of the
background noise rms fluctuation measured over an area equal to
$\theta_{\mathrm{Ein}}^{2}$. The vertical dashed lines correspond to the
maximum size for which a galaxy with a given flux can be detected in the
absence of lensing. The horizontal dotted line indicates the cross-section for
a very bright point source (i.e. the area enclosed within the caustics).

As in the point source case, the lensing cross-section increases with
increasing total flux, at fixed source size. At fixed flux,
$\sigma_{\mathrm{{SL}}}$ stays approximately constant with increasing source
size until a given value, then drops rapidly for larger sizes. From a
qualitative point of view, this behaviour can be observed also in the absence
of lensing: increasing the size of a galaxy while keeping its flux fixed
lowers its average surface brightness. If the surface brightness drops below
the sky rms fluctuation level, then it becomes very difficult to detect it. In
order to determine whether there are lensing-specific features in the
$\sigma_{\mathrm{{SL}}}-\theta_{\mathrm{e,s}}$ relation of Fig. 9, we
computed, for each source flux, the maximum half-light radius for which it can
be detected in the absence of lensing. We used the same criterion as that of
Sect. 2.4 to define a detection: we defined the source footprint as the
ensemble of pixels that are $2\sigma$ above the background and required the
total signal-to-noise ratio within the footprint to be larger than ten. The
resulting limiting sizes are shown as vertical lines in Fig. 9. For each given
total flux, the non-lensing size limit is similar to the value of
$\theta_{\mathrm{e,s}}$ at which the lensing cross-section drops.

This result suggests that, to the first approximation, lensing does not
introduce a strong selection in source size. While perhaps surprising, this
follows from the fact that gravitational lensing preserves surface brightness:
a source that can be detected in the absence of lensing will produce images
with the same surface brightness when lensed, which can be detected as well.
In order to classify a source as strongly lensed, however, we require that
multiple images are observed. Only sources that lie within a well-defined
region give rise to a strong lensing configuration. If part of the source
extends outside of this region, then only a fraction of its flux contributes
to creating a set of strongly lensed images. This lowers the signal-to-noise
ratio of the multiple images compared to the point source case. The result is
that $\sigma_{\mathrm{{SL}}}$ starts to decrease with increasing source size
at values of $\theta_{\mathrm{e,s}}$ that are smaller than the no-lensing
detection limit, as observed in Fig. 9.

At the brightest flux explored in the experiment (green line in Fig. 9), the
lensing cross-section at small source sizes is larger than the area enclosed
within the radial caustic (black dotted line in Fig. 9). This is because, when
the source is very bright, it can give rise to multiple images even while its
centroid lies outside of the radial caustic, as long as its surface brightness
distribution extends into it. This also explains why $\sigma_{\mathrm{{SL}}}$
increases with increasing source size, before dropping to zero: the more
extended the source, the farther away from the lens centre it can be while
still producing multiple images.

####  3.4.2 Large source sizes

For the large source size case we fixed the surface brightness profile of the
source and varied the Einstein radius of the lens. In particular, we set the
source half-light radius to $20$ pixels and adjusted its total flux in such a
way that the $2\sigma$ detection footprint in the absence of lensing extends
out to the half-light radius. Then, starting from the lens model used in the
previous section, we progressively increased the critical surface mass density
to reduce the Einstein radius down to values comparable to the pixel size.
Figure 10 shows the resulting strong lensing cross-section as a function of
the ratio between lens Einstein radius and source half-light radius. For
values of the Einstein radius close to the size of the PSF, the lensing cross-
section (blue solid line) is very small: this is because the perturbations
caused by lensing are not well resolved. For larger values,
$\sigma_{\mathrm{{SL}}}$ stays approximately constant, around values that are
comparable to the source size (horizontal dotted line) and much larger than
the area enclosed within the caustics (red dashed line). We conclude that, in
the large source size regime, the main factor that determines the lensing
cross-section is the area of the background source, provided that the Einstein
radius of the lens is larger than the PSF.

Figure 10:  Strong lensing cross-section of an extended source: large source
size regime ($\theta_{\mathrm{e,s}}\gtrsim\theta_{\mathrm{Ein}}$). The cross-
section is plotted as a function of the ratio between the Einstein radius of
the lens and the half-light radius of the source. The source model is fixed to
a circular Sérsic profile with $n=1$, detected out to the half-light radius.
The source half-light radius measures $20$ pixels, and the PSF has an FWHM of
two pixels. The lens model is varied: starting from the model of Fig. 9, the
critical surface mass density is increased to produce progressively smaller
values of $\theta_{\mathrm{Ein}}$. The solid line shows the lensing cross-
section, in units of the footprint area of the source in the absence of
lensing, as a function of the ratio between the lens Einstein radius and the
source half-light radius. The dashed line delimits the area enclosed within
the lens caustics, and the dotted line shows the value of the cross-section
corresponding to the area of the source. The dash-dotted line is the FWHM of
the PSF.

##  4 Lens population simulations

In the previous section we showed how the lensing cross-section, which is
closely related to the lens detection probability ${\rm P}_{\mathrm{det}}$,
varies as a function of lens and source properties. From here on we focus on
the effect that those trends have on populations of lenses. We addressed this
question by simulating populations of foreground galaxies and background
sources, selecting strong lenses among them, and comparing the properties of
the strong lens sample with the parent galaxy population. We simulated a lens-
based search (as opposed to a source-based one), in which strongly lensed
images are searched among a stellar mass-selected sample of galaxies. Our
simulations are based on empirical models, in which existing observations of
the baryonic component of galaxies are complemented with a set of assumptions
on the mass distribution of the lenses. In Sect. 4.1 we explain how we built
our foreground galaxy sample, while in Sect. 4.2 we describe the simulation of
the background source population. In Sect. 4.3 we describe how our mock
observations of lenses are generated. In Sect. 4.4 we apply a further
selection step, based on the angular size of the Einstein radius.

###  4.1 Foreground galaxies

Our foreground galaxy population consists of a volume-limited sample of early-
type galaxies, complete above a minimum observed stellar mass of
$10^{11}M_{\odot}$. We chose to focus on early-type galaxies because many
strong lensing surveys have preferentially targeted this class of objects in
their lens-finding phase (e.g. Bolton et al., 2006; Gavazzi et al., 2012;
Sonnenfeld et al., 2018). This, in turn, had a dual motivation: first, early-
type galaxies are among the most massive objects in the Universe, and
therefore they are more likely to be lenses; second, their smooth surface
brightness distribution and red colour makes it easier to detect arcs from
strongly lensed star forming galaxies around them.

We describe lenses with elliptical versions of the two-component model
introduced in Sect. 3.1. In the following sections we explain how their
parameters are generated.

####  4.1.1 Stellar mass and redshift distribution

We generated lens galaxies over a finite redshift range, $0.1<z<0.7$. We chose
these lower and upper limits because the value of the critical surface mass
density $\Sigma_{\mathrm{cs}}$ becomes very large outside of this range, and
the fraction of galaxies that act as strong lenses drops substantially as a
result.

We drew stellar masses from the stellar mass function of quiescent galaxies
measured by Muzzin et al. (2013). In particular, we used the following
comoving number density distribution:

$\Phi(M_{*}^{(\mathrm{obs})})=\Phi^{*}\left(\frac{M_{*}^{(\mathrm{obs})}}{M_{*}^{*}}\right)^{\alpha}\exp{\left\\{-\frac{M_{*}^{(\mathrm{obs})}}{M_{*}^{*}}\right\\}}.$
(21)

We set $\Phi^{*}=1.009\times 10^{-3}\,{\rm Mpc}^{-3}$, $\alpha=-0.92$ and
$\log{M_{*}^{*}}=11.21$: these are the best-fit values measured by Muzzin et
al. (2013).

We refer to the stellar masses drawn from this distribution as the observed
stellar masses, $M_{*}^{(\mathrm{obs})}$, as opposed to the true stellar
masses. The observed stellar mass is meant to represent an estimate of $M_{*}$
based on stellar population synthesis modelling, which is the method used by
Muzzin et al. (2013) to measure the galaxy stellar mass function. Stellar
population synthesis measurements, however, are subject to systematic
uncertainties, since they have not been calibrated on galaxies whose stellar
mass is known by other means. We quantify the discrepancy between the observed
and true stellar mass by means of the stellar population synthesis mismatch
parameter $\alpha_{\mathrm{sps}}$, defined as follows:

$\alpha_{\mathrm{sps}}\equiv\frac{M_{*}}{M_{*}^{(\mathrm{obs})}}.$ (22)

The most important source of systematic uncertainty in the measurement of the
stellar mass is the assumption of the stellar initial mass function (IMF).
Muzzin et al. (2013) assumed a Kroupa IMF (Kroupa, 2001) to derive their
measurements. This means that, in the absence of other systematic effects, a
galaxy with a Kroupa IMF has a value of $\alpha_{\mathrm{sps}}=1$.

For each galaxy in the sample, we randomly drew a value of
$\log{\alpha_{\mathrm{sps}}}$ from the following distribution:

${\rm
P}(\log{\alpha_{\mathrm{sps}}})\sim\mathcal{N}(0.1,\sigma_{\alpha}^{2}),$ (23)

where the notation $\mathcal{N}(\mu,\sigma^{2})$ indicates a Gaussian with
mean $\mu$ and variance $\sigma^{2}$. We set the mean of
$\log{\alpha_{\mathrm{sps}}}$ to $0.1$, as this is an intermediate value among
estimates of $\alpha_{\mathrm{sps}}$ from the literature (Conroy & van Dokkum,
2012; Cappellari et al., 2013; Smith et al., 2015; Sonnenfeld et al., 2015,
2019). We adopted a few different values for the scatter $\sigma_{\alpha}$. We
explain in Sect. 4.1.4 how these were chosen.

####  4.1.2 Stellar mass density distribution

We described the surface mass density distribution of each galaxy as an
elliptical de Vaucouleurs profile (a Sérsic profile with $n=4$). Given the
observed stellar mass of a galaxy, we assigned a half-mass radius by drawing
it from the following distribution in $\log{R_{\mathrm{e}}}$:

${\rm
P}(\log{R_{\mathrm{e}}})\sim\mathcal{N}(1.20+0.63(\log{M_{*}^{(\mathrm{obs})}}-11.4),0.14^{2}).$
(24)

Then, we assigned an axis ratio $q$ by drawing it from the following beta
distribution:

${\rm P}(q)\propto q^{\alpha-1}(1-q)^{\beta-1},$ (25)

with $\alpha=6.28$ and $\beta=2.05$. The choice for these distributions was
motivated by observations of a sample of early-type galaxies. In Appendix A we
explain how this sample was defined and how the coefficients of Eq. 24 and Eq.
25 were determined.

####  4.1.3 Dark matter distribution

We modelled the dark matter distribution of each lens galaxy with an
elliptical gNFW halo (see Sect. 3.1 for its definition). The parameters of the
radial density profile were assigned as follows. Given the stellar mass of a
lens, we first determined the value of its halo virial mass, $M_{\mathrm{h}}$.
We defined this as the mass enclosed within a spherical shell with average
density equal to $200$ times the critical density of the Universe. Using
current weak lensing constraints on the halo mass of elliptical galaxies as a
reference (Sonnenfeld et al., 2022), we drew halo masses from the following
distribution:

${\rm
P}(\log{M_{\mathrm{h}}})\sim\mathcal{N}(13.0+1.0(\log{M_{*}}-11.5),\sigma_{h}^{2}),$
(26)

which is a Gaussian in $\log{M_{\mathrm{h}}}$ with a mean that scales with
stellar mass and scatter $\sigma_{h}$. The intrinsic scatter in halo mass is
not well constrained observationally; therefore, we ran simulations with
different values of $\sigma_{h}$, as for the stellar population synthesis
mismatch parameter.

To determine the density profile of the halo we relied on a theoretically
motivated model that takes into account the effect of baryons on the dark
matter. We assumed the dark matter distribution to be initially described by
an NFW profile with a concentration444The concentration is the ratio between
the halo virial radius and the scale radius, $r_{s}$. of five, then used the
prescription of Cautun et al. (2020) to model the response of the halo to the
infall of baryons. In the Cautun et al. (2020) model the halo response is
approximated with an analytical function that depends on the present stellar
mass distribution, and typically results in a more concentrated and steeper
density profile compared to the original NFW model. Finally, we fitted a gNFW
profile to the surface mass density of the contracted halo. By doing this, we
were able to fully describe the dark matter density profile with three
parameters: $M_{\mathrm{h}}$, $\gamma_{\mathrm{DM}}$ and $r_{s}$. The dark
matter density profile defined in this way is determined uniquely by the halo
mass, the stellar mass and the half-light radius (the more concentrated the
stellar distribution, the stronger the halo response and the steeper the dark
matter density profile). In principle we could have allowed for additional
degrees of freedom, for instance by relaxing the assumption of a fixed initial
halo concentration. In practice, as we explain in Sect. 6.5, our main results
are not affected by this choice.

Given the radial profile of the dark matter halo, we obtained an elliptical
version of it by applying a transformation of the kind of Eq. 12 to its
projected surface mass density. We assumed that the axis ratio and orientation
of the halo is the same as that of the stellar component.

####  4.1.4 Intrinsic scatter

The distribution in stellar mass, halo mass and dark matter inner slope of our
sample of simulated foreground galaxies depends on parameters describing the
intrinsic scatter in these properties, namely $\sigma_{\mathrm{sps}}$ and
$\sigma_{h}$. Direct observational constraints on these quantities are poor.
However, we can derive upper limits on them on the basis of observed scaling
relations.

Early-type galaxies lie on the stellar mass fundamental plane, a scaling
relation between stellar mass, half-light radius and central velocity
dispersion (Hyde & Bernardi, 2009; de Graaff et al., 2021):

$\sigma_{\mathrm{e}}\propto
M_{*}^{(\mathrm{obs})\beta_{\sigma}}R_{\mathrm{e}}^{\xi_{\sigma}}.$ (27)

The existence of a fundamental plane relation is a consequence of the virial
theorem: the velocity dispersion of a galaxy in dynamical equilibrium is
directly related to the three-dimensional mass distribution in its inner
regions, which is typically dominated by the stellar component. At fixed
observed stellar mass distribution, however, the central velocity dispersion
can vary depending on the stellar population synthesis mismatch parameter, on
the dark matter mass and on the dark matter density profile. This gives rise
to a spread in the values of the velocity dispersion given
$M_{*}^{(\mathrm{obs})}$ and $R_{\mathrm{e}}$. Therefore, we can use the
observed scatter in velocity dispersion to put an upper limit on the intrinsic
scatter in the stellar population synthesis mismatch parameter, halo mass and
dark matter slope.

We used Jeans modelling for this purpose. We generated samples of early-type
galaxies with the recipes described above and with different values of
$\sigma_{\mathrm{sps}}$ and $\sigma_{\mathrm{h}}$. Using the spherical Jeans
equation under the assumption of isotropic orbits, we predicted the central
velocity dispersion of each galaxy in the sample. Then, we fitted a
fundamental plane relation to this mock sample and measured the predicted
scatter in velocity dispersion at fixed $M_{*}^{(\mathrm{obs})}$ and
$R_{\mathrm{e}}$. We then varied $\sigma_{\mathrm{sps}}$ and $\sigma_{h}$ to
match the observed and the predicted scatter. We did not attempt to match the
other parameters of the fundamental plane (i.e. the constant of
proportionality of Eq. 27 and the power-law indices $\beta_{\sigma}$ and
$\xi_{\sigma}$), because these are sensitive to the orbital anisotropy of the
galaxies, which we are asserting to be zero. The details of this procedure are
given in Appendix B.

We settled on three different sets of intrinsic scatter parameters, as
indicated in Table 1. We label them the fiducial, the low-scatter and the
high-scatter scenarios. The high-scatter scenario is ruled out both by the
fundamental plane and by weak lensing constraints (Sonnenfeld et al., 2022).
Moreover, our dynamical model is quite simplistic, as it neglects the effects
of orbital anisotropy and departures from spherical symmetry, which are
additional sources of scatter. Nevertheless, we use it in our experiment in
order to obtain a more conservative upper limit on the amplitude of the strong
lensing bias.

Table 1: Intrinsic scatter scenarios.

Model name | $\sigma_{\mathrm{sps}}$ | $\sigma_{\mathrm{h}}$ | Predicted FP
---|---|---|---
|  |  | scatter
Fiducial | 0.08 | 0.20 | 0.034
Low scatter | 0.05 | 0.10 | 0.021
High scatter | 0.10 | 0.30 | 0.043
555 Adopted values of the intrinsic scatter parameters $\sigma_{\mathrm{sps}}$
and $\sigma_{\mathrm{h}}$ in different simulations. For each set of values,
the fourth column indicates the scatter around the fundamental plane predicted
via the spherical Jeans analysis of Appendix B. The observed fundamental plane
scatter is $0.035$.

###  4.2 Background sources

####  4.2.1 Galaxies

Our background galaxy population is taken from the surfs-based KiDS-Legacy-
Like Simulation (SKiLLS) input catalogue, a hybrid simulation catalogue
integrating cosmological simulation with high-quality imaging observations (Li
et al., 2023). The cosmological simulation is obtained from the Synthetic
UniveRses For Surveys (surfs) simulations, a set of $N$-body simulations from
Elahi et al. (2018). The galaxy properties, including the star formation
history and the metallicity history, are from an open-source semi-analytic
model named Shark666https://github.com/ICRAR/shark (Lagos et al., 2018). The
original photometry is drawn from a stellar population synthesis technique
using stellar synthesis libraries with physically motivated dust attenuation
and re-emission models (Robotham et al., 2020). Li et al. (2023) further
applied an empirical correction to the original synthetic photometry to better
agree with the COSMOS2015 observations (Laigle et al., 2016). The galaxy
morphology is described by a Sérsic profile with three parameters: the half-
light radius in angular units $\theta_{\mathrm{e,s}}$, the Sérsic index
$n_{\mathrm{s}}$, and the axis ratio $q_{\mathrm{s}}$. These structural
parameters are learned from the imaging data obtained with the Advanced Camera
for Surveys (ACS) instrument on the Hubble Space Telescope (Griffith et al.,
2012). We refer to Li et al. (2023) for details on the learning algorithm and
validation.

The complete SKiLLS catalogue contains ${\sim}108\leavevmode\nobreak\ {\rm
deg}^{2}$ of galaxies with redshift up to $2.5$ and $r$-band apparent
magnitude down to $27$. We applied a lower limit to the source redshift, by
selecting only sources with $z_{s}>0.8$. This ensures that all of the sources
lie behind all of the lenses. A similar cut could be applied to a real survey
using photometric redshifts, to reduce the incidence of false positives (e.g.
arc-like features physically associated with the lens galaxy) in the lens
finding phase. The resulting number density of sources is $70$ arcmin-2. We
approximated their spatial distribution as uniform in the sky, that is, we
neglected clustering of the sources.

####  4.2.2 Quasars

We described the population of background quasars with the following double
power-law luminosity function in the rest-frame UV absolute magnitude, $M$:

$\Phi(M,z_{\mathrm{qso}})=\dfrac{\Phi(M^{*})}{10^{0.4(\alpha_{Q}+1)(M-M^{*})}+10^{0.4(\beta_{Q}+1)(M-M^{*})}}.$
(28)

Following Manti et al. (2017), we set $\alpha_{Q}=-1.35$ (faint-end slope),
$\beta_{Q}=-3.23$ (bright-end slope), and adopted a redshift-evolving
normalisation

$\log{\Phi^{*}}=-6.0991+0.0209z_{\mathrm{qso}}+0.0171z_{\mathrm{qso}}^{2}$
(29)

and characteristic magnitude

$M^{*}=-22.5216-1.6510z_{\mathrm{qso}}+0.2869z_{\mathrm{qso}}^{2}.$ (30)

Given the redshift and rest-frame UV luminosity of a quasar, we then computed
the apparent magnitude in the observed $i-$band, $m_{\mathrm{qso}}$, using a
quasar spectral template777https://archive.stsci.edu/hlsps/reference-
atlases/cdbs/grid/comp_qso/ built from optical and near-infrared spectra
obtained by Vanden Berk et al. (2001); Glikman et al. (2006).

For the sake of consistency with the population of background extended
sources, we limited the redshift distribution of quasars to the range
$0.8<z_{\mathrm{qso}}<2.5$. We then truncated the distribution in
$m_{\mathrm{qso}}$ at two magnitudes fainter than the detection limit (which
is specified in Sect. 4.3). Finally, we randomly placed quasars in the source
plane with a projected number density of $70$ arcmin-2. This is a much larger
number density than observed in the real universe, but we are allowed to do so
because we are not interested in predicting the absolute number of strong
lenses, so this is a legitimate choice. The advantage of boosting the number
density of quasars is that it allows us to produce a large number of lenses
without the need for generating too big a population of foreground galaxies.

###  4.3 Observations

For each of the three intrinsic scatter scenarios, we drew a population of
foreground galaxies covering $1000$ square degrees. The expectation value of
the number of galaxies in the corresponding volume, given the redshift and
stellar mass cuts described in Sect. 4.1.1, is around $300\,000$. We
approximated the foreground galaxies as isolated: when determining whether a
galaxy acts as a strong lens, we only modelled the contribution to the lensing
signal from the galaxy itself, and neglected that of the environment. We
discuss the possible implications of this approximation in Sect. 6. For each
lens, we determined its caustics relative to the highest source redshift,
$z_{\mathrm{s}}=2.5$. We then placed sources randomly behind the lens. If at
least one source fell within a circular region enclosing the caustics, we
proceeded to compute its lensed images.

For the simulation with extended sources, we produced images with properties
similar to those expected for the Euclid Wide survey (Euclid Collaboration et
al., 2022). We used a pixel size of $0.1^{\prime\prime}$ and we applied a
Moffat PSF with an FWHM of $0.2^{\prime\prime}$ and a $\beta$ parameter of
$5.0$. Finally, we assumed a background noise level such that an extended
source with half-light radius $0.5^{\prime\prime}$ and an apparent magnitude
in the absence of lensing of $m_{\mathrm{s}}=25$ is detected with $S/N=10$. We
then applied the peak detection-based lens selection criterion introduced in
Sect. 2.4 to find the strong lenses.

The fiducial scatter simulation with extended background sources produced a
sample of $2113$ lenses, corresponding to a number density of $2.1$ deg-2.
This is about a factor of five smaller than the number density predicted by
Collett (2015) for the Euclid survey. This is a result of differences in the
description of the source population, in the criteria used to define a strong
lensing event, and in the redshift and stellar mass cuts that we applied to
define the foreground galaxy population.

For the lensed quasars we did not simulate pixel-level data, but simply
computed the observed magnitudes of the multiple images. Following the
definition of Sect. 2.3, we included in the sample of strong lenses only
systems with at least two images brighter than a limiting magnitude
$m_{\mathrm{lim}}$. We set $m_{\mathrm{lim}}=23.3$, which corresponds to the
$10\sigma$ detection limit of the Legacy Survey of Space and Time (LSST) in a
single visit (Oguri & Marshall, 2010). This is motivated by the fact that the
quasar lenses are meant to simulate a sample assembled for the purpose of
carrying out time-delay measurements, which in turn require combining single-
visit detections over many epochs. The scenario that we are simulating, then,
is that of a lens search in a Euclid-like survey, followed-up with LSST time-
domain observations.

The resulting number of quasar lenses in our simulation with the fiducial
scatter is $1621$. This number is meaningless, given that the simulation was
created with an unrealistically large number density of quasars. More
interesting is the relative number of quad lenses with respect to the total,
which is about $9\%$. This is a slightly smaller value than the fraction of
quads predicted by Oguri & Marshall (2010). The reason for this discrepancy
lies in the differences between the lens mass models in the two simulations.

###  4.4 Lens finding probability

Figure 11 shows the distribution in Einstein radius of the simulated lens
samples. In all cases, the distribution peaks at $\theta_{\mathrm{Ein}}\approx
0.7^{\prime\prime}$. Although all of the lenses in these samples are detected,
it does not necessarily follow that they would all be included in a strong
lensing study. There can be a few reasons for excluding certain lenses from a
sample. One is the low accuracy of lens finders: current automated lens
finding algorithms tend to produce lens candidates samples with low purity
(see e.g. Sonnenfeld et al., 2018; Petrillo et al., 2019; Savary et al.,
2022). Such samples are then visually inspected, and only those candidates
that can be clearly distinguished from false-positives are kept. This visual
inspection step tends to disfavour lenses with a small image separation,
because of the contamination from the light of the lens galaxy.

Figure 11:  Einstein radius distribution of the simulated lens samples. These
are: galaxy-galaxy lenses in the fiducial, high scatter, and low scatter
scenarios and galaxy-quasar lenses with fiducial scatter, considering all
lenses or only lenses that produce four images.

Another possible reason for refining a sample of lens candidates is the
availability of redshift measurements. Redshifts of both the lens and the
source are needed in order to convert a lens model into a measurement of mass.
When working with large samples of lenses, obtaining spectroscopic
measurements is not a viable option, and photometric redshifts are a
necessity. Measuring the photometric redshift of a strongly lensed source,
however, is challenging, especially when the Einstein radius is small and the
source light is blended with the light from the lens (Langeroodi et al.,
2023).

Both of these scenarios can result in samples that are incomplete below a
certain value of the Einstein radius. We simulate this situation via the
following Einstein radius-dependent lens finding probability:

${\rm
P}_{\mathrm{find}}(\theta_{\mathrm{Ein}}|S,\mathrm{det})=\left\\{\begin{array}[]{ll}1&\rm{if}\,\theta_{\mathrm{Ein}}>\theta_{\mathrm{Ein,min}}\\\
0&\rm{otherwise}\end{array}\right..$ (31)

In words, all lenses with Einstein radii larger than
$\theta_{\mathrm{Ein,min}}$ are included in the sample, while all those with
smaller Einstein radius are excluded. We refer to $\theta_{\mathrm{Ein,min}}$
as the completeness limit: our simulated lens samples are complete down to
$\theta_{\mathrm{Ein}}=\theta_{\mathrm{Ein,min}}$. We explored scenarios with
different values of $\theta_{\mathrm{Ein,min}}$. As the next section shows,
the larger the minimum Einstein radius, the higher the strong lensing bias.

##  5 Results

In this section we present the results of the lens population simulations.
Section 5.1 shows the results of the galaxy-galaxy lens experiment, while
Sect. 5.2 focuses on the population of galaxy-quasar lenses. Given the number
of parameters that are needed to describe our model, providing a complete
characterisation of the strong lensing bias is a problem with relatively high
dimensionality, and is beyond the scope of this paper. For the sake of
conciseness, we focus instead on the quantities that we consider most
important. Nevertheless, the output of our simulations is available
online888https://github.com/astrosonnen/strong_lensing_tools/tree/main/papers/selection_effects.
We encourage readers who are interested in studying aspects of the strong
lensing bias that are not covered in this section to download our data and
analyse them directly.

###  5.1 Galaxy-galaxy lenses

In this section we show the results of the experiments with populations of
galaxy-galaxy lenses. We first present the results in a qualitative way, then
proceed to quantify the amplitude of the strong lensing bias in various
quantities of interest.

Figure 12 shows the distribution in the parameters of the foreground galaxies
of the lens systems, compared to those of the parent population, for the
fiducial scatter scenario with two different values of the minimum Einstein
radius: $0.5^{\prime\prime}$ and $1.0^{\prime\prime}$. A completeness limit of
$0.5^{\prime\prime}$ is close to what can currently be achieved via visual
inspection of high-resolution space-based images (Garvin et al., 2022), while
the $1.0^{\prime\prime}$ limit can be seen as a more conservative case.

Figure 12:  Comparison between the properties of lens samples and the parent
population: distribution in foreground galaxy parameters. Filled contours
represent the distribution of the parent sample, solid green lines the
distribution of the lenses with Einstein radii larger than
$0.5^{\prime\prime}$, and solid red lines the distribution of the lenses with
Einstein radii larger than $1.0^{\prime\prime}$.

The most striking difference between the samples is in the stellar mass:
strong lensing selects preferentially galaxies with larger values of
$M_{*}^{(\mathrm{obs})}$. Lenses tend to also have a larger halo mass, a
smaller half-light radius at fixed stellar mass, and a larger stellar
population synthesis mismatch parameter. The distribution in ellipticity and
inner dark matter slope of the lenses instead look very similar to that of the
parent population. Additionally, we can see that the amplitude of the strong
lensing bias appears to be always larger in the lens sample with the more
restrictive selection on Einstein radius. We quantify the amplitude of these
biases later in this section.

Figure 13 shows the distribution in parameters describing the background
source population, for the same simulations of Fig. 12. Additionally, Fig. 13
shows the subset of the parent population that consists of detectable sources.
These are background galaxies that, in the absence of lensing, can be detected
according to the same criterion used for the lensed sources (i.e. the $S/N$
over their $2\sigma$ footprint is larger than ten). The detection limit of the
survey is at $m_{\mathrm{s}}\approx 25$ (the actual limit varies depending on
the surface brightness distribution parameters). Because our simulated
background source population extends to much fainter magnitudes, the
distribution of detectable sources differs substantially from that of the
parent population. We can then consider two different strong lensing bias
definitions: one that quantifies the difference in lensed source properties
with respect to the parent population, and one that describes the difference
with respect to the detectable source population. We are mostly interested in
the second definition.

Figure 13:  Comparison between the properties of lens samples and the parent
population: distribution in background source parameters. Filled contours show
the distribution of the parent sample, solid black lines the distribution of
the detectable sources, solid green lines the distribution of the lenses with
Einstein radii larger than $0.5^{\prime\prime}$, and solid red lines the
distribution of the lenses with Einstein radii larger than
$1.0^{\prime\prime}$.

Strong lensing tends to preferentially select sources at higher redshift,
especially in the more restrictive case with
$\theta_{\mathrm{Ein}}>1.0^{\prime\prime}$. This is because, at fixed lens
properties, increasing the source redshift lowers the critical surface mass
density and, consequently, increases the size of the caustics and the Einstein
radius. The distribution in the magnitude of the lensed sources is also very
different from that of the parent distribution, as it drops rapidly for values
larger than $m_{\mathrm{s}}\approx 25$. Interestingly, however, there does not
seem to be a large difference with respect to the distribution of detectable
sources in the absence of lensing. This result can appear to be somewhat
counter-intuitive: lensing magnification should allow the detection of sources
that are intrinsically fainter than the detection limit. To some extent, this
is the case: the distribution of lensed sources shows a slight excess of
fainter galaxies compared to the unlensed case. However, as we quantify later
in this section, the difference is far from large. The reason for this
behaviour lies in the fact that, in the detection of both lensed and unlensed
sources, the most important quantity is surface brightness, which is preserved
by gravitational lensing when the source is larger than the PSF size. As we
showed in Sect. 3.4, the lensing cross-section drops to zero once the surface
brightness of the source reaches a value that would make it undetectable in
the absence of lensing. For this reason, also the distribution in half-light
radius is very similar between the lensed sources and the detectable source
population.

Qualitatively, the results shown in Figs. 12 and 13 match our expectations
from Sect. 3. In the rest of this section we quantify the lensing bias. We
present the results in three different parts. First, we focus on the
properties of the lens galaxies that can be observed directly. These are
quantities that can be derived from photometry and spectroscopy with minimal
assumptions: the lens redshift, the observed stellar mass, the half-light
radius and the axis ratio999Strictly speaking, the axis ratio in mass is not
necessarily observable, but in the context of our simulations it is, since
light and mass have the same ellipticity.. Any bias in these quantities can be
determined relatively easily in a real strong lens survey. In the second part,
we consider the parameters related to the mass distribution: the stellar
population synthesis mismatch parameter and the dark matter distribution
parameters. Determining the bias on these parameters is much more difficult,
but these are quantities of great interest from a galaxy science point of
view. In the third part, we focus on source distribution parameters.

####  5.1.1 Bias in observable lens parameters

Figure 14 shows the median redshift, median $M_{*}^{(\mathrm{obs})}$, median
size for a given $M_{*}^{(\mathrm{obs})}$, and median axis ratio of various
lens population simulations, as a function of the minimum Einstein radius. We
defined the median size at fixed $M_{*}^{(\mathrm{obs})}$ by fitting the
following mass-size relation to the $R_{\mathrm{e}}-M_{*}^{(\mathrm{obs})}$
distribution:

$\log{R_{\mathrm{e}}}\sim\mu_{R,0}+\beta_{R}(\log{M_{*}^{(\mathrm{obs})}}-11.4).$
(32)

The quantity shown in the third panel of Fig. 14 is the parameter $\mu_{R,0}$,
which is the average $\log{R_{\mathrm{e}}}$ at an observed stellar mass of
$\log{M_{*}^{(\mathrm{obs})}}=11.4$. The horizontal dashed line in each panel
shows the value of the parent population: the larger the distance between the
curve of a simulation and this line, the higher the strong lensing bias.

Figure 14:  Bias on lens observable properties as a function of the minimum
Einstein radius. First panel: Median redshift. Second panel: Median observed
stellar mass. Third panel: Median half-light radius at an observed mass of
$\log{M_{*}^{(\mathrm{obs})}}=11.4$. Fourth panel: Median axis ratio. In each
panel, the dashed black line indicates the value of the parent population.
Error bars indicate the standard deviation of the mean of the lens sample.

As was already visible in Fig. 12, there is a clear bias towards lower
redshift, higher stellar mass and smaller sizes, with the bias becoming
stronger for more restrictive cuts in Einstein radius. Biases in stellar mass
and size are stronger for the simulation with low scatter in
$\alpha_{\mathrm{sps}}$ and $M_{\mathrm{h}}$. This can be explained as
follows. For a given value of the observed stellar mass, strong lensing
selection favours galaxies with a larger $\alpha_{\mathrm{sps}}$ or
$M_{\mathrm{h}}$, or with a smaller half-light radius. Since the galaxy
stellar mass function of the parent population is steep, lenses tend to have a
relatively small $M_{*}^{(\mathrm{obs})}$ and large values of
$\alpha_{\mathrm{sps}}$ or $M_{\mathrm{h}}$ for their observed stellar mass
(this is shown more clearly in the next section). If the intrinsic scatter in
$\alpha_{\mathrm{sps}}$ and $M_{\mathrm{h}}$ is low, however, the number of
galaxies with a small $M_{*}^{(\mathrm{obs})}$ and a large stellar or halo
mass is greatly reduced. Only galaxies with a large $M_{*}^{(\mathrm{obs})}$
or a small size can therefore act as strong lenses. This is an interesting
result, because it suggests that it is in principle possible to use the strong
lensing bias on $M_{*}^{(\mathrm{obs})}$ or $R_{\mathrm{e}}$, which is
observable, as a way to constrain the amplitude of the intrinsic scatter in
the mass parameters, which is poorly known.

####  5.1.2 Bias in lens mass parameters

Figure 15 shows the median of the distribution in various mass-related
quantities, as a function of the minimum Einstein radius. The top panel is the
median stellar population synthesis mismatch parameter. In all simulations,
strong lenses are biased towards larger values than the parent population. The
bias is larger the higher the intrinsic scatter in $\alpha_{\mathrm{sps}}$ and
$M_{\mathrm{h}}$, and increases with increasing $\theta_{\mathrm{Ein,min}}$.
For the fiducial model, the bias on $\alpha_{\mathrm{sps}}$ can be as small as
$0.03$ dex (7%), if no cut on Einstein radius is applied. However, it can rise
up to $0.09$ dex (23%) when considering only lenses with
$\theta_{\mathrm{Ein}}>2^{\prime\prime}$. For comparison, we also show two
reference values of $\alpha_{\mathrm{sps}}$, corresponding to a Chabrier IMF
(Chabrier, 2003) and a Salpeter IMF (Salpeter, 1955). These values roughly
bracket the current systematic uncertainty on $\alpha_{\mathrm{sps}}$.

The second panel of Fig. 15 shows the median halo mass of galaxies with an
observed stellar mass of $\log{M_{*}^{(\mathrm{obs})}}=11.4$. We measured this
quantity by fitting the following relation to the
$M_{\mathrm{h}}-M_{*}^{(\mathrm{obs})}$ distribution,

$\log{M_{\mathrm{h}}}\sim\mu_{\mathrm{h},0}+\beta_{\mathrm{h}}(\log{M_{*}^{(\mathrm{obs})}}-11.4),$
(33)

and taking the resulting value of the parameter $\mu_{\mathrm{h},0}$.
Similarly to the $\alpha_{\mathrm{sps}}$ case, the strong lenses are biased
towards larger values of the halo mass, with the bias being larger for higher-
scatter simulations and more restrictive cuts on $\theta_{\mathrm{Ein}}$. In
the fiducial scatter scenario with a completeness limit of
$\theta_{\mathrm{Ein,min}}=1.0^{\prime\prime}$, the halo masses of lenses are
on average $0.16$ dex larger than those of their parent population.

The third panel of Fig. 15 shows the median projected dark matter mass
enclosed within an aperture of $5$ kpc, $M_{\mathrm{DM},5}$, at fixed observed
stellar mass and half-light radius. We obtained this quantity by fitting the
following relation to the
$M_{\mathrm{DM},5}-M_{*}^{(\mathrm{obs})}-R_{\mathrm{e}}$ distribution:

$\log{M_{\mathrm{DM},5}}\sim\mu_{\mathrm{DM},0}+\beta_{\mathrm{DM}}(\log{M_{*}^{(\mathrm{obs})}}-11.4)+\xi_{\mathrm{DM}}(\log{R_{\mathrm{e}}}-1.2).$
(34)

Figure 15 shows the value of $\mu_{\mathrm{DM},0}$. The bias on this quantity
is qualitatively similar to that on the total halo mass, but much smaller in
amplitude. There are two reasons for this. First, the dependence of
$M_{\mathrm{DM},5}$ on $M_{\mathrm{h}}$ is shallower than linear. This is
because, as the virial mass increases, the virial radius increases as well:
the extra mass is spread over a larger volume, and therefore the mass within
the inner region does not increase proportionally. Second, while
$\mu_{\mathrm{h},0}$ is the halo mass at fixed stellar mass,
$\mu_{\mathrm{DM},0}$ is measured at fixed half-light radius as well. The
density profile of the dark matter halos in our simulation have a dependence
on galaxy size: the response of dark matter to baryons is stronger for more
concentrated stellar distributions (see Sect. 4.1.3). By capturing this
dependence in the model of Eq. 35, the residual scatter in $M_{\mathrm{DM},5}$
around the mean is reduced, and so is the strong lensing bias. This, however,
is a minor effect: we repeated the analysis while setting the dependence on
size to zero and found minimal differences on the derived values of
$\mu_{\mathrm{DM},0}$. In the fiducial scenario, the bias on
$M_{\mathrm{DM},5}$ is as small as $0.02$ dex for
$\theta_{\mathrm{Ein,min}}<0.5^{\prime\prime}$, and is negligible in the low-
scatter scenario.

In the fourth panel of Fig. 15 we show the average inner dark matter slope at
fixed observed stellar mass and half-light radius, which we measured by
fitting the following model

$\gamma_{\mathrm{DM}}\sim\mu_{\gamma,0}+\beta_{\gamma}(\log{M_{*}^{(\mathrm{obs})}}-11.4)+\xi_{\gamma}(\log{R_{\mathrm{e}}}-1.2),$
(35)

and taking the resulting value of $\mu_{\gamma,0}$. Different simulations show
different trends in the average $\gamma_{\mathrm{DM}}$. These are the results
of underlying correlations between the dark matter density profile and the
stellar and halo mass. In our model, galaxies of a given size with a larger
stellar mass have a steeper dark matter slope, because the response of the
dark matter to the infall of baryons is stronger. Vice versa, galaxies with a
larger halo mass have a shallower density slope. At fixed
$M_{*}^{(\mathrm{obs})}$ and $R_{\mathrm{e}}$, strong lenses have both a
larger stellar mass (because their $\alpha_{\mathrm{sps}}$ is larger) and a
larger halo mass. In the simulation with low scatter, the bias in stellar mass
is more important; therefore, $\gamma_{\mathrm{DM}}$ has a positive bias. In
the simulation with large scatter, the bias in halo mass dominates, and
therefore $\gamma_{\mathrm{DM}}$ is negatively biased.

Figure 15:  Bias on lens mass properties as a function of the minimum Einstein
radius. First panel: Median $\log{\alpha_{\mathrm{sps}}}$. Dotted lines
indicate values of $\alpha_{\mathrm{sps}}$ corresponding to a Chabrier and a
Salpeter IMF ($\alpha_{\mathrm{sps}}=1$ corresponds to a Kroupa IMF). Second
panel: Median $M_{\mathrm{h}}$ at $\log{M_{*}^{(\mathrm{obs})}}=11.4$. Third
panel: Median $M_{\mathrm{DM},5}$ at $\log{M_{*}^{(\mathrm{obs})}}=11.4$ and
$\log{R_{\mathrm{e}}}=1.2$. Fourth panel: Median $\gamma_{\mathrm{DM}}$ at
$\log{M_{*}^{(\mathrm{obs})}}=11.4$ and $\log{R_{\mathrm{e}}}=1.2$. In each
panel, the dashed line indicates the value of the parent population. Error
bars indicate the standard deviation of the mean of the lens sample.

####  5.1.3 Bias in source parameters

Figure 16 shows the lensing bias in source-related parameters, with respect to
the population of sources that are detectable without lensing. As previously
seen in Fig. 13, lensed sources are biased towards higher redshift (top panel
of Fig. 16), with the trend being larger for more restrictive cuts on the
Einstein radius. The second panel shows the bias in the source magnitude. This
bias is negative for small values of $\theta_{\mathrm{Ein,min}}$, that is,
lensed sources tend to be brighter than their field counterparts. Naively one
might have expected the opposite trend, as strong lensing magnification allows
the detection of sources that are intrinsically fainter than the detection
limit. However, as the analysis of Sect. 3 clearly shows, the lensing cross-
section is always larger for brighter sources, and that explains the sign of
the bias. At the same time, the median does not capture the whole picture of
the bias in $m_{\mathrm{s}}$. For instance, in the fiducial simulation with no
cut on the Einstein radius, the $90$th percentile of the $m_{\mathrm{s}}$
distribution is $26.05$, while that of the population of detectable sources is
$0.14$ mag brighter: indeed, strong lensing allows the detection of fainter
sources.

The third and fourth panels of Fig. 13 show the bias in half-light radius and
Sérsic index, respectively, for sources within a magnitude bin centred on
$m_{\mathrm{s}}=25$ and with a $0.4$ mag width. We detected no clear sign of
bias. Finally, the fifth panel shows the bias in source axis ratio. Also in
this, case no obvious sign of bias was detected, except in the largest values
of $\theta_{\mathrm{Ein,min}}$.

Figure 16:  Bias on source properties as a function of the minimum Einstein
radius. First panel: Median source redshift. Second panel: Median source
magnitude. Third panel: Median half-light radius in the magnitude bin
$24.8<m_{\mathrm{s}}<25.2$. Fourth panel: Median Sérsic index in the magnitude
bin $24.8<m_{\mathrm{s}}<25.2$. Fifth panel: Median axis ratio. In each panel,
the dashed line indicates the value of the population of detectable sources.
Error bars indicate the standard deviation of the mean of either the full
sample (first, second, and fifth panel) or the bin (third and fourth panel).

###  5.2 Galaxy-quasar lenses

####  5.2.1 Bias in lens parameters

The main goal of the experiment with lensed quasars is to check whether there
are any differences in the strong lensing bias with respect to the extended
source case, at fixed properties of the foreground galaxy population. For this
reason, we only ran simulations with lensed quasars in the fiducial scatter
scenario and compared the results with those from the extended source
simulation. We are also interested in understanding how the subset of quad
lenses differs from the entire population of quasar lenses; therefore, we also
analysed that subsample on its own.

The first four rows of Fig. 17 show the strong lensing bias as a function of
minimum Einstein radius in the following quantities: stellar population
synthesis mismatch parameter, halo mass at fixed stellar mass, enclosed dark
matter mass at fixed stellar mass and half-light radius, and axis ratio. The
biases of the population of lensed quasars (red curves) are very similar to
that of lensed galaxies (blue curves), especially for values of
$\theta_{\mathrm{Ein},min}<1^{\prime\prime}$, When considering only quad
lenses, however, there are some differences, the most remarkable of which is
the bias in the axis ratio: quad lenses tend to be on average galaxies with a
higher ellipticity. This is a well-known effect (see e.g. Keeton et al., 1997)
that can be explained in the context of the results of Sect. 3.2: the higher
the ellipticity, the larger the area enclosed within the inner caustic, which
is where a source needs to lie in order to produce four or more images. Quad
lenses tend to also have a slightly larger halo mass at fixed stellar mass,
for $\theta_{\mathrm{Ein},min}<1^{\prime\prime}$.

Figure 17:  Bias in the population of quasar lenses (all and quads only),
compared to the extended source simulation. First panel: Median
$\log{\alpha_{\mathrm{sps}}}$. Second panel: Median $M_{\mathrm{h}}$ at
$\log{M_{*}}=11.5$. Third panel: Median $M_{\mathrm{DM},5}$ at
$\log{M_{*}}=11.5$ and $\log{R_{\mathrm{e}}}=1.2$. Fourth panel: Median lens
axis ratio. Fifth panel: Median mass-sheet transformation parameter, defined
in Appendix C. In each panel, the dashed line indicates the value of the
parent population. Error bars indicate the standard deviation of the mean of
the lens sample.

####  5.2.2 Impact on $H_{0}$

A measurement of the time delay between the images of a lensed quasar can be
converted into an estimate of the expansion rate of the Universe, $H_{0}$. In
general, though, the value of $H_{0}$ is degenerate with the mass distribution
of the lens (and of matter along the line of sight), which can be difficult to
constrain on an individual lens basis. One possible strategy to break this
degeneracy is to use information gathered from a larger set of lenses as a
prior for the time-delay lenses. This was the strategy adopted by Birrer et
al. (2020) in the most recent measurement of $H_{0}$ by the TDCOSMO
collaboration: they combined a larger sample of extended source lenses with a
smaller set of time-delay lenses. In this section we want to quantify the
possible bias on $H_{0}$ that can result from this approach.

The lens mass model used by Birrer et al. (2020) consists in a power-law
density profile modified by a mass-sheet transformation. The corresponding
projected density profile is

$\kappa(\theta)\propto\lambda_{\mathrm{mst}}\theta^{1-\gamma}+(1-\lambda_{\mathrm{mst}}),$
(36)

where $\lambda_{\mathrm{mst}}$ is called the mass-sheet transformation
parameter. Strong lensing data can constrain the power-law index $\gamma$ very
well, but not at all $\lambda_{\mathrm{mst}}$. Given the time delay, the
inferred value of the Hubble constant $H_{0}$ scales linearly with
$1/\lambda_{\mathrm{mst}}$. Birrer et al. (2020) measured the population
distribution of $\lambda_{\mathrm{mst}}$ on the sample of extended source
lenses, by combining strong lensing and stellar kinematics data, and then used
it as a prior in the time-delay analysis. Although the mass model used in our
work is different from that of Eq. 36, we can estimate the bias associated
with this approach by computing the value of $\lambda_{\mathrm{mst}}$ for each
lens and then comparing the strong lensing bias on $\lambda_{\mathrm{mst}}$ of
extended source and quasar lenses. Appendix C explains how we obtained
$\lambda_{\mathrm{mst}}$ for each lens. The strong lensing bias on the median
value of $\lambda_{\mathrm{mst}}$ is shown in the fifth panel of Fig.
17101010Since $\lambda_{\mathrm{mst}}$ is only defined for strong lenses, the
median value of $\lambda_{\mathrm{mst}}$ of the general population of galaxies
is not defined..

At fixed minimum Einstein radius, there is a small difference between the
median $\lambda_{\mathrm{mst}}$ in extended source and quasar lenses,
typically of order 1%. The associated bias on $H_{0}$ is, therefore, equally
small. The value of $\lambda_{\mathrm{mst}}$, however, increases with
increasing minimum $\theta_{\mathrm{Ein}}$. This implies that if the two
samples probe different ranges in Einstein radius, the bias on $H_{0}$ can be
larger than that. The bias can also be larger if the properties of the general
population of galaxies (as opposed to the strong lenses) are used to inform
the structure of time-delay lenses (see Collett & Cunnington, 2016). This is
simply because quasar lenses are on average more similar to galaxy-galaxy
lenses than they are to the underlying galaxy population (the blue and red
curves in Fig. 17 are closer to each other than they are to the dashed
horizontal lines).

##  6 Discussion

###  6.1 Key results

In order for a galaxy-source pair to be included as a strong lens in a survey,
three conditions must be met. First of all, at least part of the source must
be multiply imaged. Second, the multiple images must be detectable. Third, the
lens must be recognised as such. Each one of these conditions introduces a
bias with respect to the parent population of foreground galaxies and
background sources. Together, they define the lens selection probability term
${\rm P}_{\mathrm{sel}}$ of Eq. 1. The first two points are intrinsic to a
strong lensing survey and constitute an unavoidable source of bias. The best
case scenario occurs when the efficiency of including a detected strong lens
in a survey is always one, and the sample is $100\%$ complete. In this case,
the third condition does not introduce any further selection and the strong
lensing bias is minimised. We explore this scenario in Sect. 5 when setting
the minimum Einstein radius to zero. If only lenses with Einstein radii larger
than a given threshold are selected, however, the bias generally increases.

The strong lensing bias affects all quantities that are related to the mass
distribution of the lens, as well as the redshifts of lens and source and the
source surface brightness parameters. Some of these quantities, such as the
observed stellar mass and half-light radius of the lenses, can be directly
measured, and it is straightforward to quantify their lensing bias. Other
quantities, however, such as the stellar population synthesis mismatch
parameter or the dark matter content, are difficult to obtain via traditional,
non-lensing, observations. Our simulations are particularly useful to quantify
the bias on these properties.

In the fiducial scatter scenario, the one consistent with the observed scatter
around the fundamental plane, the bias on $\alpha_{\mathrm{sps}}$ varies from
$0.03$ dex, when no restrictions on the lens Einstein radius are applied, to
$0.09$ dex, corresponding to the extreme case in which only lenses with
$\theta_{\mathrm{Ein}}>2.0^{\prime\prime}$ are selected (see Fig. 15). A
reasonable value for the minimum Einstein radius in a space-based survey like
Euclid or the Chinese Space Station Telescope is
$\theta_{\mathrm{Ein,min}}=0.5^{\prime\prime}$. In this case, the bias on
$\alpha_{\mathrm{sps}}$ is slightly smaller than $0.04$ dex. The current
systematic uncertainty on $\alpha_{\mathrm{sps}}$ is $0.2-0.3$ dex: this is
roughly the difference in measurements of the stellar mass of a galaxy
obtained with a Chabrier or a Salpeter IMF. Compared with this uncertainty,
the amplitude of the strong lensing bias on $\alpha_{\mathrm{sps}}$ is small:
strong lensing observations can be used directly to discriminate between these
two alternative choices of IMF, without the need to correct for selection
effects. Such a goal could be reached with a sample size of a thousand lenses
and a statistical study of the kind proposed by Sonnenfeld & Cautun (2021).

Strong lenses are also biased towards larger halo masses. Nevertheless, when
focusing on the dark matter content in the inner regions, the amplitude of the
strong lensing bias is relatively small, especially when controlling for the
stellar distribution of the lens galaxies. For instance, at fixed observed
stellar mass and half-light radius, the bias on the projected dark matter mass
enclosed within $5$ kpc, $M_{\mathrm{DM},5}$, is only a few percent in the
fiducial scatter scenario with $\theta_{\mathrm{Ein,min}}<1.0^{\prime\prime}$.
This means that strong lenses can indeed be used to understand the inner dark
matter distribution of galaxies, as long as the dependence of the dark matter
distribution on the properties of the baryonic component is accurately
modelled (for example, by following the approach of Sonnenfeld & Cautun,
2021).

We also looked at the bias on source-related parameters. Strong lensing causes
the luminosity function of background sources to be broadened, compared to the
distribution of detectable sources. On the one hand, it preferentially selects
brighter sources, because the lensing cross-section increases with increasing
source brightness. On the other hand, it allows for the detection of sources
that are intrinsically fainter than the detection limit in the absence of
lensing. Interestingly, we did not find any significant bias on the source
size, at fixed magnitude. This result follows from the fact that our lens
detection criterion relies on a surface brightness threshold, and surface
brightness is preserved by lensing. However, it appears to be in contradiction
with the work of Oldham et al. (2017), who argued that their strong lens
sample selected preferentially compact sources. The origin of this discrepancy
probably lies in the differences between the criteria used to define a strong
lens in the two studies. The Oldham et al. (2017) sample was selected
primarily via spectroscopy, by looking for signatures of two galaxies at
different redshifts in the Sloan Digital Sky Survey (SDSS; York et al., 2000)
data. SDSS spectra were taken in fibres with a $1.5^{\prime\prime}$ radius.
Lensed galaxies that are comparable in size to this scale, or larger, are less
likely to be detected, because part of their flux extends outside of the
fibre. Compact galaxies with a large overall magnification, instead, are more
likely to be detected (see the discussion in Sect. 5.3 of Oldham et al.,
2017). We refer the reader to Arneson et al. (2012) for a thorough study of
selection effects associated with spectroscopy.

There is also an indirect way in which strong lensing could preferentially
probe more compact sources than possible in the field. Sources that are
smaller than the size of the PSF can be easily confused as stars in our own
Galaxy. Since these stars cannot be strongly lensed into multiple images, the
detection of a strong lens automatically confirms the extragalactic nature of
a lensed source. This type of selection effect could explain, for example, the
detection of an apparent outlier in the magnitude-size relation of Lyman-break
galaxy by means of strong lensing, by Jaelani et al. (2020). It is, however,
improper to refer to this effect as strong lensing bias, since the star-galaxy
separation is a bias that primarily affects field observations. Our analysis
showed that the tendency to select compact sources is not a general feature of
strong lens samples.

###  6.2 Combining lensed quasars with lensed galaxies

Our experiments were also useful for understanding the possible biases that
might incur when combining information from samples of galaxy-galaxy lenses
with samples of galaxy-quasar lenses. On the one hand the biases in the mass-
related quantities ($\alpha_{\mathrm{sps}}$ and $M_{\mathrm{DM},5}$) are very
similar in simulations with the same foreground galaxy population and
different background sources (see Fig. 17). On the other hand, these biases
are a function of the completeness limit. In order to use a sample of strong
lenses as a prior for another sample, then, it is important to make sure that
1) the parent population of foreground galaxies among which lenses are
searched for is the same in both surveys, and 2) the two surveys can probe the
same distribution in Einstein radius.

The above argument applies to samples of quasars selected regardless of the
number of multiple images. When dealing with quad lenses, however, the
situation is more complicated. First of all, because our experiment revealed
differences in the bias on the dark matter distribution between quads and the
entire sample of lensed quasars. Second, because quad lenses tend to have a
preferentially higher ellipticity. The ellipticity is important in the context
of stellar dynamical analyses, which are often used in combination with strong
lensing to constrain lens mass parameters (see e.g. Yıldırım et al., 2020). In
order to correctly interpret stellar dynamics data, assumptions on the three-
dimensional structure of a lens must be made. This is directly related to the
projected ellipticity: galaxies with an axis ratio close to one tend to be
preferentially elongated in the line-of-sight direction, and vice versa. When
using stellar dynamics-based mass measurements to inform the properties of a
sample of quad lenses, it is therefore important to take into account possible
biases due to the different three-dimensional structure of the two samples.

###  6.3 The importance of the source population properties

Strictly speaking, the results shown in Sect. 5 apply only to samples of
lenses and background sources with the same properties as our simulations.
Real lens samples can have more complex source distributions. For example,
many existing lens-finding algorithms are tuned to identify sources with
different colours from the lenses, introducing a selection that can modify the
source distribution in redshift-magnitude space. However, the comparison of
Sect. 5.2 shows that, at fixed foreground galaxy population, the strong
lensing bias for the sample of extended sources is indistinguishable from that
of the lensed quasar population, at least for values of
$\theta_{\mathrm{Ein,min}}<1.0^{\prime\prime}$. This result suggests that the
dependence of the strong lensing bias on the details of the background source
population is very weak, as long as no additional selection on the image
configuration is applied (e.g. by selecting only quad lenses). This is not a
coincidence, but is a consequence of the fact that the dependence of the
strong lensing cross-section on the lens parameters is a weak function of
source magnitude111111The cross-section itself is a strong function of source
magnitude, but the trends between $\sigma_{\mathrm{{SL}}}$ and the lens
parameters are not., especially for magnitudes close to the detection limit
(see Sect. 3.2). We can then conclude that the details of the properties of
the background source population play a secondary role in determining the
strong lensing bias. A selection aimed at specific kinds of sources, such as
blue star-forming galaxies or quasars, results primarily in different numbers
of lenses being discovered, but does not change how biased the lenses are with
respect to the general galaxy population.

###  6.4 Mitigation strategies

As we discussed throughout this paper, the amplitude of the strong lensing
bias depends on the intrinsic scatter in the mass parameters of the foreground
galaxy population. One way to minimise the bias, therefore, is to identify
scaling relations between observable quantities and mass-related properties
that can account for part of the scatter. Describing the inner dark matter
distribution as a function of stellar mass and half-light radius, as we did in
Sect. 5.1, is the first step in this direction: the bias in
$M_{\mathrm{DM},5}$ at fixed $M_{*}$ and $R_{\mathrm{e}}$ is smaller than the
overall shift in the median $M_{\mathrm{DM},5}$ of the lens population (albeit
by a marginal amount, as we pointed out earlier). This description could be
extended by including the central velocity dispersion as an additional control
parameter. The velocity dispersion is directly related to the mass
distribution of the lens galaxy, which means that, for example, the
distribution in $\alpha_{\mathrm{sps}}$ at fixed velocity dispersion should be
narrower than its global distribution marginalised over the whole population.
The central velocity dispersion, however, is very sensitive to the orbital
anisotropy, to the three-dimensional structure, and to gradients in stellar
mass-to-light ratio, which are not well known. Therefore, it is difficult to
quantitatively estimate the benefit of including it in the description of the
lensing bias.

If one wishes to directly account for the strong lensing bias, the formally
correct procedure is to explicitly model all of the selection steps in a
Bayesian hierarchical formalism, as explained by Sonnenfeld (2022). Although
this can be computationally challenging, machine learning can offer an
efficient alternative (Legin et al., 2023). In order for either of these
approaches to work, however, it is essential that the lens selection procedure
can be simulated. This, in turn, requires having an objective definition of a
strong lensing event. A peak-based definition such as that introduced in Sect.
2.4 could be used for this purpose.

Nevertheless, it can still be difficult to fully forward model strong lensing
selection if visual inspection by humans is involved in the definition of the
lens sample. In that case, a possible alternative is to approximate the lens
finding probability ${\rm P}_{\mathrm{find}}$. For example, we can make the
assumption that ${\rm P}_{\mathrm{find}}$ depends purely on the Einstein
radius. This is a reasonable assumption, as long as the lens finding procedure
does not selectively pick lenses with different image configurations depending
on their Einstein radius. The dependence of ${\rm P}_{\mathrm{find}}$ on
$\theta_{\mathrm{Ein}}$ could then be described empirically and inferred
during the analysis. This is essentially the approach adopted by Sonnenfeld et
al. (2019) in the analysis of strong and weak lensing data from the Hyper
Suprime-Cam Survey.

###  6.5 Limitations of our analysis

The simulations on which our analysis is based are as complex as required by
the goal of the analysis itself, which is to estimate the amplitude of the
strong lensing bias in a few key quantities. We did, however, make some
simplifying assumptions. One such assumption consisted in neglecting the
contribution from line-of-sight structure and the environment to the lensing
signal, which typically introduce an external shear and convergence. The
effect of including external shear in a lens model is similar to that of
changing the ellipticity of the lens. External convergence mimics the effect
of adding or removing a constant sheet of invisible mass, producing an effect
similar to varying the distribution of dark matter. Typical values of the
external shear and external convergence are $|\gamma|<0.1$ and
$|\kappa_{\mathrm{ext}}|<0.1$ (Millon et al., 2020): these values are small
compared to the typical ellipticities and dark matter fractions of our
simulated lenses. Therefore, while including them might modify slightly the
amplitude of the strong lensing bias in the axis ratio and dark matter
distribution parameters, the conclusions of our analysis would not be
affected.

We also assumed that the dark matter density profile of the galaxies is
completely determined by the halo mass and by the stellar mass distribution
(see Sect. 4.1.3). In reality we expect there to be a range of density
profiles, for example as a result of halos having a non-zero scatter in their
initial (i.e. before baryonic infall) concentration parameter. We could in
principle carry out an experiment with such an additional source of variation
in the dark matter profile. Adding a scatter in concentration while keeping
the halo mass scatter parameter $\sigma_{h}$ fixed results in a larger spread
in the inner dark matter distribution, which would increase the amplitude of
the lensing bias. However, the value of $\sigma_{h}$ that we chose for our
fiducial experiment was tuned on the basis of the predicted scatter in
velocity dispersion around the fundamental plane. That scatter would also be
increased by adding flexibility to the dark matter density profile. Then, in
order to be consistent with the analysis carried out in this paper, the value
of $\sigma_{h}$ would need to be lowered accordingly. This, in turn, would
reduce the amplitude of the strong lensing bias and introduce a bias on the
halo concentration. Since lensing and dynamics are sensitive to the mass
distribution at comparable scales (the half-light radius and Einstein radius
are similar for most of the lenses), we expect these two effects to cancel out
to first approximation. The lensing bias on quantities directly related to the
inner mass distribution, such as $\alpha_{\mathrm{sps}}$ and
$M_{\mathrm{DM},5}$, would be left roughly unchanged. The main effect of
adding scatter to the concentration is to reduce the correlation between the
total halo mass and the dark matter mass enclosed within the Einstein radius.
This would then reduce the amplitude of the strong lensing bias on the halo
mass.

Another simplifying assumption that we made was that of adopting a constant
mass-to-light ratio for the stellar component. Massive galaxies are known to
have colour gradients, the main effect of which is to cause the stellar half-
mass radius to be smaller than the half-light radius (Szomoru et al., 2013;
Suess et al., 2019). Since, as shown in Sect. 5, strong lensing bias tends to
select galaxies with a more compact stellar distribution, it would also
preferentially select galaxies with a steeper mass-to-light ratio gradient.
However, testing for the amplitude of the strong lensing bias on the mass-to-
light ratio gradient is beyond the scope of this paper.

###  6.6 Lens versus source selection

The simulations described in Sect. 4, on which our experiment is based, are
examples of so-called lens-selected surveys: the lens samples are built by
first defining a population of possible lens galaxies and then searching for
strong lenses among them. Photometric lens searches are typically carried out
in this way (see e.g. Gavazzi et al., 2012; Sonnenfeld et al., 2018; Petrillo
et al., 2019). An alternative approach is to build a source-selected sample,
by first searching for signatures of strong lensing within a population of
possible background sources and then obtaining information on the lens
galaxies (see e.g. Hezaveh et al., 2013). In principle, these two approaches
lead to different strong lensing biases. For example, from Fig. 12 we can see
that the stellar mass distribution of the lenses extends below the minimum
mass cut that we applied to the population of foreground galaxies: a source-
based search that is able detect those lenses would result in a different
sample and a different associated strong lensing bias. However, when
considering only lenses with $\theta_{\mathrm{Ein}}>1^{\prime\prime}$, lens-
selection and source-selection yield the same sample (see the red histogram in
Fig. 12): in the limit of large Einstein radius, the two operations
(identifying suitable lens galaxies and detecting a strongly lensed source)
commute.

In general, there is no such thing as a purely lens- or source-selected
sample. A lens-based search is always also source-based, because a detection
of the lensed source is required to identify a lens. A source-based search is
also lens-based, because the lens properties determine important observational
features of the lensed source, such as the number of images, the image
separation and the total flux.

###  6.7 Alternative lens definitions

Our extended source analysis was based on a peak detection criterion to define
strong lensing events, as explained in Sect. 2.4. However, many previous
studies aimed at simulating populations of lenses adopted a magnification-
based definition, in which a lens-source system is classified as a strong lens
if the magnification of the source is larger than a minimum value (see e.g.
Collett, 2015; Robertson et al., 2020). In order to check whether our results
are robust to the choice of definition of a strong lensing event, we re-
analysed the fiducial sample with this alternative definition. In particular,
we defined as strong lenses all systems with a detected background source and
$\mu_{\mathrm{tot}}>3$, where $\mu_{\mathrm{tot}}$ is the total magnification
(obtained by integrating the flux of all of the lensed images and dividing it
by the intrinsic source flux). The total number of strong lenses increased by
$34\%$ with this new definition, mostly due to the inclusion of fainter
sources with no detected multiple images. Figures 18 and 19 show the strong
lensing bias on a series of lens- and source-related properties. Besides the
change in the source redshift and magnitude (first two panels in Fig. 19), we
see very little differences in the bias with respect to the peak-based
definition of the fiducial analysis. We then conclude that the exact
definition of a strong lensing event does not have a big impact on the strong
lensing bias of lens-related properties. We stress, however, that
magnification is not a directly observable quantity, and therefore we advise
against applying magnification-based lens definitions to define lens samples,
for the reasons discussed in Sect. 2.4.

Figure 18:  Lens strong lensing bias in the case of a magnification-based lens
definition (brown curves), compared to that of the fiducial analysis (blue
curves). First panel: Median $\log{\alpha_{\mathrm{sps}}}$. Second panel:
Median $M_{\mathrm{DM},5}$ at $\log{M_{*}}=11.5$ and
$\log{R_{\mathrm{e}}}=1.2$. Third panel: Median $\gamma_{\mathrm{DM}}$ at
$\log{M_{*}^{(\mathrm{obs})}}=11.4$ and $\log{R_{\mathrm{e}}}=1.2$. Fourth
panel: Median axis ratio. Fifth panel: Median mass-sheet transformation
parameter, defined in Appendix C. In each panel, the dashed line indicates the
value of the parent population.  Figure 19:  Source strong lensing bias in the
case of a magnification-based lens definition (brown curves), compared to that
of the fiducial analysis (blue curves). First panel: Median source redshift.
Second panel: Median source magnitude. Third panel: Median half-light radius
in the magnitude bin $24.8<m_{\mathrm{s}}<25.2$. The dashed line indicates the
median value of the population of sources detectable in the absence of
lensing.

##  7 Conclusions

Strong lensing is a very active line of research, but a comprehensive
understanding of the selection effects associated with it has so far been
lacking. This work takes a major step towards filling that knowledge gap.
After a thorough investigation, we have learned several lessons regarding the
strong lensing bias in photometrically selected lens samples. The following
are the most important ones.

  1. 1.

The strong lensing cross-section increases primarily with increasing lens
mass, with decreasing half-mass radius, and with increasing source brightness.
At fixed stellar distribution and fixed dark matter mass enclosed within a
given aperture, varying the inner dark matter slope has little impact on the
strong lensing cross-section (as long as the aperture within which the dark
matter mass is normalised is comparable to the Einstein radius).

  2. 2.

The strong lensing cross-section has little dependence on the size of the
source if the former is smaller than the Einstein radius and is detectable in
the absence of lensing. Sources with a surface brightness that is too low to
be detected are still undetected when strongly lensed.

  3. 3.

Lens galaxies tend to be more massive and more compact than their non-lens
counterparts. Their redshift distribution also differs from that of the
general galaxy population. At fixed observed stellar mass (i.e. inferred by
means of stellar population synthesis), lens galaxies have a higher intrinsic
stellar mass (i.e. a larger stellar population synthesis mismatch parameter,
$\alpha_{\mathrm{sps}}$) and a higher dark matter halo mass. At fixed stellar
mass and size, lens galaxies are still biased towards a larger dark matter
content.

  4. 4.

The amplitude of the strong lensing bias depends on how broad the distribution
of the parameters describing the lens population is. Important quantities in
this regard are the intrinsic scatter in the stellar population synthesis
mismatch parameter, $\sigma_{\mathrm{sps}}$, and in the dark matter halo mass,
$\sigma_{\mathrm{h}}$. Increasing the values of $\sigma_{\mathrm{sps}}$ and
$\sigma_{\mathrm{h}}$ results in a stronger bias on $\alpha_{\mathrm{sps}}$
and on the dark matter mass, and a weaker bias on the observed stellar mass
and half-light radius. This implies that, in principle, we could constrain the
intrinsic scatter parameters, which are currently poorly known, by measuring
the amplitude of the strong lensing bias on $M_{*}^{(\mathrm{obs})}$ and
$R_{\mathrm{e}}$, which is easily observable.

  5. 5.

The strong lensing bias varies depending on the completeness of the lens
sample as a function of the Einstein radius. Surveys that can discover lenses
with a smaller Einstein radius have a smaller associated strong lensing bias
in all quantities.

  6. 6.

Under reasonable assumptions regarding the intrinsic scatter parameters, for a
Euclid-like survey that is complete down to
$\theta_{\mathrm{Ein}}=0.5^{\prime\prime}$ the bias on $\alpha_{\mathrm{sps}}$
is smaller than $0.04$ dex (10%). This bias is much smaller than the current
systematic uncertainty on the stellar population synthesis-based stellar
masses. Therefore, strong lensing measurements could be used directly to
calibrate stellar mass measurements of massive galaxies to $10\%$ accuracy,
without the need to correct for selection effects. Under the same assumptions,
the strong lensing bias on the average halo mass at fixed stellar mass is
$0.07$ dex, while that on the inner dark matter distribution at fixed stellar
mass and size is $0.02$ dex.

  7. 7.

Strong lensing selection broadens the magnitude distribution of background
sources, compared to the population of objects that are detectable without
lensing. At the same time, we did not find any evidence for a bias in the size
distribution of background sources at fixed magnitude.

  8. 8.