# Find the haystacks, then look for needles: The rate of strongly lensed
transients in galaxy-galaxy strong gravitational lenses

Ana Sainz de Murieta1, Thomas E. Collett1, Mark R. Magee2, Justin D. R.
Pierel3, Wolfgang J. R. Enzi1, Martine Lokken4, Alex Gagliano5,6 and Dan
Ryczanowski1,
1 Institute of Cosmology and Gravitation, University of Portsmouth, Burnaby
Rd, Portsmouth PO1 3FX, UK
2 Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4
7AL, UK
3 Space Telescope Science Institute, Baltimore, MD 21218, USA
4 Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain
5 The NSF AI Institute for Artificial Intelligence and Fundamental
Interactions
6 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The time delay between appearances of multiple images of a gravitationally
lensed supernova (glSN) is sensitive to the Hubble constant, $H_{0}$. As well
as time delays, a lensed host galaxy is needed to enable precise inference of
$H_{0}$. In this work we investigate the connection between discoverable
lensed transients and their host galaxies. We find that LSST will discover 88
glSNe per year, of which 54% will also have a strongly lensed host. The rates
can change by approximately 30 percent uncertainty depending primarily on the
choice of unlensed SN population and uncertainties in the redshift evolution
of the deflector population, but the fraction of glSNe with a lensed host is
consistently around a half. LSST will discover 20 glSNe per year in systems
that could plausibly have been identified by Euclid as galaxy-galaxy lenses
before the discovery of the glSN. Such systems have preferentially longer time
delays and therefore are well suited for cosmography. We define a golden
sample of glSNe Ia with time delays over 10 days, image separations greater
than 0.8 arcseconds, and a multiply imaged host. For this golden sample, we
find $91\%$ occur in systems that should already be discoverable as galaxy-
galaxy lenses in Euclid. For cosmology with glSNe, monitoring Euclid lenses is
a plausible alternative to searching the entire LSST alert stream.

###### keywords:

gravitational lensing: strong – transients:supernovae

††pubyear: 2024††pagerange: Find the haystacks, then look for needles: The
rate of strongly lensed transients in galaxy-galaxy strong gravitational
lenses–Find the haystacks, then look for needles: The rate of strongly lensed
transients in galaxy-galaxy strong gravitational lenses

##  1 Introduction

Strong gravitational lensing is the phenomenon by which the intense
gravitational potential from a foreground galaxy or cluster deflects the light
coming from a background source and causes multiple images to appear. It
serves as a powerful probe in astrophysics and cosmology. This rare effect can
provide valuable insights into high-redshift sources, the distribution of
matter in large galaxies and clusters, and the expansion rate of the Universe
(Treu, 2010). Using strong gravitational lensing to measure the expansion rate
of the Universe was first suggested by Refsdal (1964), who proposed using the
time difference between the appearance of lensed images of a supernova,
combined with an accurate model of the lens potential, to measure the time
delay distance and hence the expansion rate, or the Hubble constant ($H_{0}$).
Local measurements of $H_{0}$ combining Cepheids and Type Ia supernovae (SNe
Ia) yield a value of $73.30\pm 1.04$ km s-1 Mpc-1 (Riess et al., 2022), and
are in disagreement with "early-time" measurements from the Cosmic Microwave
Background, for which the latest analysis yields a value of $H_{0}=67.4\pm
0.5$ km s-1 Mpc-1 (Planck Collaboration et al., 2020). Additional independent
measurements, such as those from strong lensing time-delays, must be obtained
if we are to test whether this tension is due to systematics or physics beyond
our current cosmological model (Di Valentino et al., 2021).

Gravitationally lensed quasars have dominated the field of time-delay
cosmography ever since their discovery (Walsh et al., 1979). The discovery of
numerous multiply-imaged quasars along with continuous light curve monitoring
allowed for the first time delay measurements (e.g. Kundić et al. 1997;
Schechter et al. 1997). Subsequent advancements in both light curve analysis
and lens modelling throughout the 21st century have yielded a value of
$H_{0}=73.3^{+1.7}_{-1.8}$ km s-1 Mpc-1 (Wong et al. 2020, Millon et al. 2020)
from 6 lensed quasars, consistent with the local Cepheid and SNe Ia
measurement.

Despite their maturity as a cosmological probe, gravitationally lensed quasars
pose certain challenges. Firstly, their light curves are stochastic and
require years of monitoring to measure precise time delays. Secondly, they are
consistently bright, outshining their host and making it hard to use the
morphology of the host as constraints on the lens model. These challenges can
be overcome by using gravitationally lensed supernovae (glSNe) as they vary on
shorter timescales and they fade away within a few months. Unfortunately,
glSNe are extremely rare. Only a handful of multiply-imaged glSNe have been
found to date, most of which were lensed by galaxy clusters (Kelly et al.,
2015; Rodney et al., 2021; Kelly et al., 2022; Chen et al., 2022; Frye et al.,
2023; Pierel et al., 2024). The first one of these, SN Refsdal, allowed for
the first time-delay measurement from glSNe, which yielded a value of
$H_{0}=64.8^{+4.4}_{-4.3}$km s-1 Mpc-1 (Kelly et al., 2023). A large fraction
of the error budget was linked to underlying uncertainties in modelling the
complex mass structure of the galaxy cluster lens, but statistical errors
remained the largest source of uncertainty (e.g. Grillo et al. 2024). The
simpler profiles of isolated galaxy lenses leads to a smaller level of
modelling uncertainty, making them ideal for precision $H_{0}$
measurementsMillon et al. (2020). But they have shorter time delays, suffer
worse microlensing and potentially larger differential dust extinction. These
combine to give larger fractional time-delay measurement uncertainties (e.g.
Birrer et al. 2022; Treu et al. 2022).

To date, two glSNe lensed by single galaxies have been discovered: iPTF16geu
(Goobar et al., 2017), and SN Zwicky (Goobar et al., 2023). Both of these
glSNe, however, originated from compact systems with short time delays
($\lesssim 1$ day), which did not allow for precise fractional time-delay
measurements. In retrospect, this issue should have been expected for current
magnitude-limited surveys as they are biased towards systems with high
magnifications that tend to have shorter time-delays and smaller image
separations (Sainz de Murieta et al., 2023; Sagués Carracedo et al., 2024).
The Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST;
Collaboration et al. 2009) will increase the rate of galaxy-scale glSNe
discoveries substantially (Arendse et al., 2023; Goldstein et al., 2019; Sainz
de Murieta et al., 2023; Wojtak et al., 2019), with $\mathcal{O}(10^{2})$
forecast over the survey’s 10 year period. Roughly $40\%$ of discovered glSNe
will have time delays that are long enough to allow for a precise ($\lesssim$
10% precision) time-delay measurement (Sainz de Murieta et al., 2023).

Whilst strong lensing time delays are inversely proportional to $H_{0}$, the
constant of proportionality can only be inferred from strong lens modelling.
In order to obtain an accurate lens model and break degeneracies that come
with lens modelling, information about a lensed host is required (Suyu et al.,
2014; Wong et al., 2016; Bonvin et al., 2016; Rusu et al., 2020). In some
cases, the offset between the supernova and the bright core of its host galaxy
results in a glSN that does not appear to have a strongly-lensed host. In
other cases, some discoverable glSNe will appear in hosts that do not reach
the magnitude required to be detectable with wide-field optical surveys, and
will therefore appear hostless (Ryczanowski et al., 2020). Both of these types
of systems will limit our ability to obtain an accurate model of the lens and
therefore to obtain a precise estimate of $H_{0}$.

The major challenge of the large datasets provided by Rubin will be to
identify gravitationally lensed supernovae amongst the alert stream of
$\mathcal{O}$($10^{7}$) transient alerts per night. Due to the survey’s slow
cadence, additional follow-up will be needed in order to obtain precise time-
delay measurements (Huber et al., 2019). Several search strategies for glSNe
have been suggested, e.g. Quimby et al. 2014; Wojtak et al. 2019; Arendse et
al. 2023; Goldstein et al. 2019; Shu et al. 2018. For most of these methods,
the search for glSNe is hindered by our inability to distinguish them from
much more common contaminants (Magee et al., 2023), such as un-lensed SNe
hosted within potential lenses, AGN(i), faint uncatalogued variable stars,
TDEs, ANTs… Monitoring known lens systems sidesteps the contaminant problem
(apart from transients in the deflector), but at the cost of missing glSNe
that are not in already known lens systems.

The focus of this work is to investigate the frequency of glSNe hosts being
lensed and the prospects of watchlists of galaxy-scale lenses for glSNe
searches. We aim to answer two questions. 1) If imaging is deep enough, what
fraction of glSNe will have a multiply imaged host? 2) What fraction of glSNe
will occur in galaxy-galaxy lensing systems that could plausibly be discovered
beforehand in wide-field imaging surveys? To answer these questions, we
calculate how many glSNe will have a strongly lensed host galaxy, and how many
of these lensed hosts can be discovered with current and upcoming surveys. We
use this to compute rates for glSNe that will allow for an accurate lens
model. We also estimate how many more glSNe we can detect by targeting strong
lenses, accounting for the rate of false positives of this method. We simulate
a realistic population of gravitationally lensed supernovae and their lensed
host galaxies in Section 2, and we define our criteria for detecting our
sources. In Section 3 we study the rates of glSNe discoverable with this
method and the properties of their systems, as well as the implications of our
results for upcoming time-delay cosmography estimates. We discuss the
limitations of our method and potential challenges of this strategy in Section
4 and finalise with our conclusions in Section 5.

##  2 Simulations

###  2.1 SN and host population: the SCOTCH catalogue

We draw our source population from the Simulated Catalogue of Optical
Transients and Correlated Hosts (SCOTCH; Lokken et al. 2023). This simulated
catalogue is based on observational host galaxy properties from GHOST
(Gagliano et al., 2021), a database of 16,175 spectroscopically classified SNe
and the properties of their host galaxies that was designed to allow for
realistic host-transient simulations. The transients in the SCOTCH catalogue
are simulated following Kessler et al. (2019) and are each associated to
simulated galaxies from CosmoDC2 (Korytov et al., 2019) catalogues based on
their intrinsic properties. Similar catalogues have already shown great
agreement with data (e.g. Vincenzi et al. 2021).

The SCOTCH catalogue contains 11 transient classes. For this work we simulate
both SN Ia and several subtypes of core-collapse supernovae: SN II, SN IIn, SN
Ib and SN Ic. We group the last two in one class, SN Ibc, as both of these
types are hydrogen-poor “stripped-envelope” CCSNe (Modjaz et al. 2019) and are
both found in hydrogen-poor environments (Modjaz et al., 2019). Following
previous rate estimates (Wojtak et al., 2019; Goldstein et al., 2019), we
separate the SN IIn class, characterised by the presence of narrow spectral
line components. The remaining classes of core-collapse SNe are grouped in the
SN II class. We also simulate SLSN I (hereafter SLSN), which are found at much
lower rates than other subtypes but are extremely bright events, characterised
for reaching magnitudes beyond -21 mag in optical bands (Moriya et al., 2018).

The transients are simulated following the framework of the SNANA simulation
code (Kessler et al., 2009). The SNANA simulations include redshift-dependent
volumetric rates for each transient class. As detailed in Kessler et al.
(2019), these are measured from existing observations, meaning they become
increasingly uncertain with redshift. The SN Ia rate is defined as a piecewise
function that follows the rates determined by Dilday et al. (2008) at redshift
$z<1$ and the rate used by Hounsell et al. (2018) for redshift $z>1$, which
came from measurements by Rodney et al. (2014). The core-collapse SN rate is
given by Strolger et al. (2015), with the relative fractions of SN II, SN IIn,
and SN Ibc corresponding to Li et al. (2011). The rate of SLSNe appears to be
consistent with the cosmic star formation history (Prajs et al., 2016).
Therefore, the SLSN redshift dependent rate is modelled using the star
formation history calculated by Madau & Dickinson (2014), with $R(0)=2\times
10^{-8}$ yr-1Mpc-3. As well as redshift dependent rates, the SNANA simulation
code also includes spectral light curve templates for each transient class
from Pierel et al. (2018). Lokken et al. (2023) provides details about the
templates used for every transient class, which are mostly based on those
presented by Kessler et al. (2019).

Each of the supernovae simulated in the SCOTCH catalogue is then assigned a
host galaxy from the CosmoDC2 catalogue. In order to do this, host libraries
are tailored for different transient categories containing correlations in
host galaxy colour and magnitude. Due to the small sample size of some SN
classes in GHOST, the hosts of archival GHOST SNe are divided into three
broader categories: SN Ia hosts, SN II hosts, and H-poor SN hosts, which will
be used for SN Ibc and SLSNe. Using the properties of galaxies in the host
libraries, a weightmap was defined, which is a multidimensional grid of host
properties with values in the grid corresponding to the probability of a host
with those properties to host a SN of a certain subtype. In order to match the
host libraries to the observed distributions in GHOST, each galaxy was
assigned a weight by interpolating the weightmap values at the value
corresponding to its properties. Each SN is assigned a host galaxy based on
this weight.

After the host galaxy has been chosen, each system is assigned an offset
between the supernova and the centre of the host galaxy. The distribution of
these offsets is motivated by prior studies that find that SN surface density
profiles roughly trace the surface brightness distribution of the disks that
host them (Hakobyan et al., 2009). Despite the fact their offset prescriptions
were not fitted to observations, the results were compared with SNe Ia data
from the Young Supernova Experiment (YSE) Data Release 1 (Aleo et al., 2023),
see Lokken et al. (2023) Fig. 11. We note this comparison has not been done
for other subtypes, which might show statistically significant differences
from our catalogue. The distribution of the SN-host galaxy offsets for a
sample at redshift $z\approx 0.6$ are shown in Fig. 1.

####  2.1.1 Tidal Disruption Events and Kilonovae

We also simulate Tidal Disruption Events (TDEs), which occur the tidal forces
of a supermassive black hole (SMBH) tear apart an orbiting star (Hills, 1975);
and kilonovae (KNe), which are transient events due to neutron star mergers
(Li & Paczyński, 1998) in order to calculate lensed transient rates. These
rates are much less certain than for supernovae: following SCOTCH, we use the
rates of (Kochanek, 2016) and the Kessler et al. (2019) templates for TDEs and
the rates of (Scolnic et al., 2018) (modified according to Kessler et al.
2019) combined with the Kasen et al. (2017) templates for KNe.

We will not investigate the lensed hosts of these rare transients because
SCOTCH does not include realistic host-transient associations for these. The
majority of TDEs happen in the cores of the host, so if the TDE is lensed the
host most likely will be. Observations of kilonovae are scarce, leading to a
poor understanding of the kilonovae-host connection, especially at source
redshifts relevant to lensing.

Figure 1: Relative SN offsets from the host centre in units of the host half
light radius, for each source type at redshift $z=0.6\pm 0.05$.

###  2.2 Deflector population

We use the methods described by Collett (2015) to simulate a population of
lens galaxies. Each of our lenses is assigned a density profile corresponding
to a Singular Isothermal Ellipsoid (SIE; Kormann et al. 1994), which provides
a good approximation to the mass profile of elliptical galaxies. We assume the
lensing contribution of matter along the line of sight to be negligible. Thus,
our simple deflector model is described by the velocity dispersion of the
lens, $\sigma_{l}$; the ellipticity, q; and the lens redshift $z_{l}$. We use
the following distributions to sample these parameters, which are the same as
Sainz de Murieta et al. (2023).

The velocity dispersion distribution of galaxies is given by a modified
Schechter function (Sheth et al., 2003),

$\frac{dN}{dV}=dn=\phi_{*}\left(\frac{\sigma}{\sigma_{*}}\right)^{\alpha}\exp\left[-\left(\frac{\sigma}{\sigma_{*}}\right)^{\beta}\right]\frac{\beta}{\Gamma(\alpha/\beta)}\frac{d\sigma}{\sigma},$
(1)

where $\Gamma$ is the gamma function and $dn$ is the differential number of
lenses per unit volume. We use the parameter values derived by Bernardi et al.
(2010) from the SDSS DR6 data describing all galaxy types,
$\phi_{*}=2.099\times 10^{-2}$ $(h/0.7)^{3}$ Mpc-3, $\sigma_{*}=113.78$ km
s-1, $\alpha=0.94$, $\beta=1.85$. We draw velocity dispersions for each of the
deflectors from this distribution. We populate the sky with lenses uniformly
in comoving volume and assume no redshift evolution.

We follow Collett (2015) to describe the distribution of the ellipticity for a
given value of $\sigma_{v}$ using a Rayleigh distribution

$P(1-q|s)=\frac{1-q}{s^{2}}\exp\Big{[}-\frac{1}{2}\Big{(}\frac{1-q}{s}\Big{)}^{2}\Big{]},$
(2)

with $s=(A+B\ \sigma_{v})$, A = 0.38 and B = 5.7 $\times$ $10^{-4}$ km-1 s.
This distribution takes into account that more massive galaxies tend to be
closer to spherical. We truncate our distribution at $q=0.2$ to exclude
highly-flattened mass profiles.

The probability of a supernova at redshift $z_{s}$ to have a lens with
ellipticity $q$, redshift $z_{l}$ and velocity disperion $\sigma$ will be
given by the optical depth of the system. This is defined as the area of the
sky that is lensed by each lensing galaxy. It is given by

$\tau(q,zl,\sigma,z_{s})=\frac{1}{A_{\mathrm{sky}}}f(q)\pi\theta_{\mathrm{Ein}}^{2}(z_{l},\sigma,z_{s}),$
(3)

where $A_{\mathrm{sky}}$ is the area of the sky, $f(q)$ is a function that
describes the ratio between the optical depth of an elliptical lens of axis
ratio $q$ and a spherical lens with the same Einstein radius and
$\theta_{\mathrm{Ein}}$ is the Einstein radius of the system. For an SIE, this
is a function of the velocity dispersion of the galaxy, $\sigma_{v}$, and the
angular diameter distances between observer, lens and source:

$\theta_{\rm{Ein}}^{\mathrm{SIE}}=4\pi\left(\frac{\sigma_{v}}{c}\right)^{2}\frac{D_{\rm{ls}}}{D_{\rm{s}}}.$
(4)

Following Sainz de Murieta et al. (2023), we sample lens redshifts in the
range $0<z_{l}<1.5$, velocity dispersions in the range $100<\sigma<400$ km
s-1and axis ratios $0.2<q<1$ following these distributions in order to
generate a mock catalogue of 100 000 lens galaxies.

###  2.3 Unresolved glSNe lightcurves

In practice, most glSN lightcurves will appear unresolved, since typical image
separations are smaller than the scatter caused by atmospheric seeing
($\lesssim$ 1") . This means that in most cases, the individual images will
not be distinguishable, but will appear as a single object combining the flux
from all the individual images (Goldstein et al., 2019).

We first simulate our original unlensed lightcurves using the set of templates
from Kessler et al. (2019). We aim to replicate the magnitude distribution
from SCOTCH. In order to achieve this, we set our model to have the same
$z$-band magnitude at peak as the lightcurves in SCOTCH once correcting for
host-galaxy extinction111We draw our host $A_{V}$ values from the same
distribution but do not have the exact value for each lightcurve. Therefore
our simulated lightcurves will not be identical but will share similar
properties on a population level..

We then apply magnification and time-delay corrections to the unlensed
lightcurves, so that the combined unresolved glSN light curve is given by

$F_{\mathrm{total}}(z_{s},t)=\mu_{1}F_{1}(z_{s},t)+\sum_{i=2}^{N}\mu_{i}F_{i}(z_{s},t-\Delta
t_{i}),$ (5)

where $F_{i}\left(z_{s},t\right)$ is the flux derived from the spectral
template, which is dependent on the SN redshift and observer-frame time,
$\mu_{i}$ and $\Delta t_{i}$ are the corresponding magnifications and time
delays for each of the images, and $N$ is the total number of images for each
supernova.

We model the effects of reddening due to Milky Way extinction and dust in the
lens and host galaxy. For Milky Way extinction we simulate our glSNe at random
positions in the sky and calculate Milky Way extinction from galactic dust
maps from Schlegel et al. (1998). We model host galaxy extinction using the
same colour law used to describe Milky Way reddening (Fitzpatrick, 1999). We
draw $A_{V}$ from a "Galactic Line of Sight" distribution, following Eqn. 2 of
Wood-Vasey et al. (2007), and choose $R_{V}$ = 3.1. We follow Feindt et al.
(2019) to simulate dust at the redshift of the lens galaxies by drawing the
$E(B-V)$ value from an exponential distribution with $\lambda=0.11$.

Whilst the effects of microlensing by stars can add noise to individual glSN
lightcurves (Dobler & Keeton, 2006), it has been found to have a modest impact
on population level glSNe yields (Arendse et al., 2023; Sainz de Murieta et
al., 2023). For simplicity, we ignore the effects of microlensing in this
work.

###  2.4 Defining discoverable systems

We define different detectability criteria for time-varying and static
sources. The discoverability of transient sources will be dependent on the
observing strategy of the imaging survey used to detect them, while our
ability to identify galaxy-scale lensing systems is highly dependent on the
size of the system and our ability to distinguish lensing features (i.e arcs;
Rojas et al. 2023).

####  2.4.1 Discoverable glSNe with LSST

LSST is a 10 year optical survey of the Southern Hemisphere Sky carried out by
the Vera C. Rubin Observatory (Ivezić et al., 2019). It is set to begin
operations in early 2025. Roughly 90% of its observing time will be allocated
to the "Baseline Survey Strategy", a "wide-fast-deep" (WFD) survey covering
about half of the sky in ugrizy filters through 15-second exposures. The WFD
survey is expected to use a rolling cadence, which means some areas of its
footprint will be assigned more frequent visits. This will improve the light
curve sampling of the objects discovered in those fast-cadence areas. In
addition to this, up to $10\%$ of Rubin observing time will be dedicated to
other programs, including a "deep drilling" survey that will cover a smaller
area of the sky at a significantly higher cadence (Bianco et al.,
2021)222Specific details about the survey observing strategy are described in
https://pstn-055.lsst.io/. The median effective seeing in the ugrizy bands are
$(1.10,1.03,0.99,0.95,0.93,0.92)$ arcseconds respectively. The median single-
visit 5$\sigma$ depths for the WFD fields are
$(23.9,25.0,24.7,24.0,23.3,22.1)$ in the ugrizy bands.

We use simsurvey 333http://github.com/ZwickyTransientFacility/simsurvey
(Feindt et al., 2019) to simulate glSNe light curves as they would be observed
with LSST. We use the survey plan baseline_v3.0 from the Operations Simulator
(OpSim)444The 10-year survey plan simulations baseline_v3.0 can be downloaded
from http://astro-lsst-01.astro.washington.edu:8080/. This observing strategy
uses two filters at a time using a mix of uniform cadence and a half-sky
rolling cadence in the "deep drilling" survey region. Updated versions of this
strategy do not introduce significant changes that would impact our strategy
(Arendse et al., 2023). For the CCDs, we rescale the ZTF corners to account
for the field of view of LSST. We define different discoverability criteria:
one 3$\sigma$ detection, one 5$\sigma$ detection, and three 5$\sigma$
detections (all in any of the bands) to explore the potential of this method
to identify glSNe earlier in order to request follow-up observations. The
glSNe are distributed evenly across the LSST footprint. We simulate 10,000
glSNe per redshift bin of size $\Delta z=0.1$ and their lightcurves are
simulated as described in Section 2.3.

####  2.4.2 Discoverable lenses

A fraction of the glSNe will also have a magnified host; however, not all of
them will be detectable as strong lenses, either because the host will be
fainter than the magnitude limit of the survey, or due to lack of obvious
shearing. Most imaging-based strong lens searches require at least some level
of human input. We establish a set of requirements for a lens to be classified
as detectable. Firstly, and most importantly, we need the centre of the host
galaxy to be magnified, that is

$\theta_{\mathrm{Ein}}^{2}>x^{2}_{s}+y^{2}_{s},$ (6)

where $\theta_{\mathrm{Ein}}$ is the Einstein radius of the system and
$(x_{s},y_{s})$ are the unlensed source positions. Human lens finders will
look for typically arc-like features, the presence of a similar-colour counter
image, bluer features, or a morphology typical of a strongly lensed system. We
summarise this by requiring tangential shearing of the arcs to be observable,
so

$\mu_{\mathrm{TOT}}>3,$ (7)

where $\mu_{\mathrm{TOT}}$ is the total magnification of the source galaxy. We
also require there to be enough separation between the features in order for
us to distinguish between the different images

$\theta_{\mathrm{Ein}}>\mathrm{seeing}.$ (8)

Finally, we require the lensed galaxy to be brighter than some magnitude
threshold $m_{\mathrm{thres}}$ defined based on the survey,

$m_{\mathrm{host}}-2.5\mathrm{log}(\mu_{\mathrm{TOT}})<m_{\mathrm{thres}}.$
(9)

Due to these requirements, we consider two types of surveys. Upcoming ground
based surveys, such as LSST, will show an improvement in limiting magnitude,
but the seeing is impacted by atmospheric effects. Space-based surveys such as
Euclid avoid this issue, leading to better image resolution. These surveys
have the ability to distinguish lenses with smaller Einstein radii (Eqn. 8),
but will require brighter lenses and will survey a smaller fraction of the
sky. For these reasons, we set the LSST $i$-band threshold to 25 magnitudes
and the seeing cut at $0.8$ arcseconds. For Euclid, we set the magnitude
threshold at $m_{\mathrm{thres}}=22.5$ magnitudes in the $i$-band and a seeing
of 0.4 arcseconds.

The synergy between these surveys has already been suggested for supernova
studies, as combining data in the Euclid Deep Drilling Fields with the LSST
deep field observations can reduce systematics in measuring SN distances at
higher redshifts (Bailey et al., 2023). In our case we aim to maximize our
chances of identifying gravitational lenses by combining the strengths of both
survey types: increased sensitivity to faint lensing events and enhanced
resolution for resolving lensing configurations in more compact systems.

###  2.5 Sample weighting

To calculate rate estimates of glSNe for each class of SN, we calculate the
optical depth for each lensing system using Eqn. 3. Each glSN will then be
assigned a weight

$w_{i}=\mathrm{f}_{\mathrm{sky}}\frac{N_{\mathrm{glSNe}}}{N_{\mathrm{sim}}}\tau(q,z_{l},\sigma,z_{s})p_{\mathrm{det,i}}$
(10)

where $\mathrm{f}_{\mathrm{sky}}$ is the fraction of the sky we simulated
using simsurvey (which in our case is 75% of the sky), $N_{\mathrm{sim}}$ is
the size of our catalogue, $\tau$ is the optical depth and
$p_{\mathrm{det,i}}$ is the number of times the supernova was recovered with
simsurvey using the LSST observing strategy divided by the number of times it
was simulated. The total number of glSNe in the sky for each subtype,
$N_{\mathrm{glSNe}}$ is calculated by computing the integral

$N_{\mathrm{glSN}}=\int_{0}^{2}dz_{s}\frac{1}{1+z_{s}}\frac{dN}{dz_{s}}P_{\mathrm{SL}}(z_{s}),$
(11)

where $dN$/$dz_{s}$ is the volumetric rate of glSNe of a certain type and the
probability of a supernova being strongly lensed, $P_{\mathrm{SL}}(z_{s})$ is
given by

$P_{\mathrm{SL}}(z_{s})=\int_{0}^{1.5}dz_{l}\int_{0.2}^{1}dq\int_{100}^{400}d\sigma\tau(q,z_{l},\sigma,z_{s})p(q|\sigma)\frac{dN}{dz_{l}d\sigma}.$
(12)

The probability $p(q|\sigma)$ represents the distribution of axis ratio given
a velocity dispersion, and is derived from Eqn. 2. We calculate
${dN}/{dz_{l}d\sigma}$ from combining Eqn. 1 with the differential comoving
volume. The number of discovered glSNe of each subtype will then be given by
the sum of the weights

$\sum_{i}w_{i}=N_{\mathrm{glSNe,det}}/\ \mathrm{year}.$ (13)

To perform our simulations, we sample 10 000 SNe and hosts from the SCOTCH
catalogues for each SN subtype. Each of these SNe are then assigned 10 lens
galaxies from our deflector catalogue, drawn from a distribution proportional
to the optical depth of the lensing systems. Each SN is placed within the area
of strong lensing of the lens. We then use lenstronomy
555https://github.com/sibirrer/lenstronomy (Birrer & Amara, 2018; Birrer et
al., 2021) to compute the positions of their images, as well as the
magnification and time-delays for each of the images. We use the offset
between the SN and its host galaxy given in the SCOTCH catalogues to define
the host galaxy position, and use
LensPop666https://github.com/tcollett/LensPop (Collett, 2015) to solve the
lens equation and simulate images for the galaxies in our catalogue. A sample
of these simulated images can be found in Fig. 2. If the host is strongly
lensed, we recover its total magnification and magnitude in different bands.

Figure 2: Sample of simulated lenses from our catalogue. The first two columns
show the lenses without the supernova. Column 1 shows a mock Euclid $VIS$-band
image and Column 2 a mock LSST $gri$ composite. Column 3 shows a mock LSST
image including the supernova images at peak luminosity of the first image.
Column 4 shows what the system would look like with an orbit of HST follow up
($gri$, with 960s per band), including the supernova images. The lens in Row 1
would be discoverable in both LSST and Euclid. Row 2 shows a lens that would
be discoverable in Euclid but not LSST and Row 3 shows an example of a glSN
without a lensed host, similar to SN Zwicky.

##  3 Results and Discussion

###  3.1 Rates of lensed transients

We follow the method described in Section 2.5 to calculate the annual rates of
different SN subtypes discovered by LSST. These are shown in Table 1. We find
a total of 88 glSNe per year will be discoverable with LSST, out of which 39
will be glSNe Ia, 6 will be glSNe Ibc, 41 will be of Type II, 1.5 will be Type
IIn and 0.5 will be SLSNe. We also forecast 0.2 lensed TDEs and 10-3 lensed
kilonovae per year. These rates are sensitive to a range of assumptions, and
their uncertainties are discussed in detail in Section 4. The discoverable
glSNe are found in systems with median velocity dispersion $\sigma_{v}=197$ km
s-1, lens redshift $z_{l}=0.28$, source redshift $z_{s}=0.68$, and Einstein
radius $\theta_{\mathrm{Ein}}=0.53^{\prime\prime}$. The median longest777quads
have more than one-time delay time-delay in these systems is $\Delta
t_{\mathrm{max}}$ = 16 days, and the median maximum angular separation between
images of the same SN is $\theta_{\mathrm{max}}=1.02^{\prime\prime}$. This
means a large fraction of glSNe will have at least partially resolved images
and time-delays that are long enough to allow for precise time-delay
measurements ($\Delta$t>10 days; Pierel & Rodney 2019).

Class | Lensed Rate (yr-1)
---|---
SN Ia | 39
SN Ibc | 6
SN II | 41
SN IIn | 1.5
SLSN | 0.5
TDE | 0.2
Kilonovae | $10^{-3}$
Total |  $\approx$ 88
Table 1: Yearly rates of lensed transients by class

###  3.2 Lensed hosts of glSNe in LSST and Euclid

With our calculated rates of discovered glSNe, we now investigate which
systems will also feature strongly lensed host galaxies. Figure 3 shows the
distribution of offsets between the centre of the host galaxy and the centre
of the lens. We find that $\approx 54$ % of glSNe discoverable with LSST will
have a host that is located within the area of strong lensing. Dividing this
result by subtype, we find 51% of glSNe Ia will have a strongly lensed host,
42% of glSNe Ibc, 58% of glSNe II, 46% of glSNe IIn and 68% of glSLSNe. Only
33% of galaxies have a total magnification $\mu_{\mathrm{TOT,host}}>3$,
required for observing tangential shearing of the arcs as explained in Section
2.4. Table 2 shows a breakdown by SN subtype of these results.

Figure 3: Offset between the centre of the lens and the centre of the glSN host in units of the Einstein radius of the system. 54% of discoverable glSN host centres will be found in the unshaded region, and therefore will be strongly lensed. Class | $\mu_{\mathrm{TOT,host}}>2$ | $\mu_{\mathrm{TOT,host}}>3$
---|---|---
SN Ia | 51% | 30%
SN Ibc | 42% | 26%
SN II | 58% | 36%
SN IIn | 46% | 28%
SLSN | 68% | 37%
Total | 54% | 33%
Table 2: Fraction of glSNe with a strongly lensed host per subtype. The total
fraction is the average weighted by the relative rates of each subtype.

In reality, not all of these hosts will be discoverable with current and
upcoming surveys. As detailed in Ryczanowski et al. (2020), some of these will
be fainter than the detection threshold of the survey, making the glSNe appear
hostless. In order to define what fraction of these will be detectable with
LSST and Euclid, we need to specify our discovery criteria for each of these
surveys.

In Section 2.4 we discussed the conditions for discovering a lensed host. For
the case of LSST, we defined a magnitude threshold requiring objects to be
brighter than 25 mag in the $i$-band and an Einstein radius threshold
$\theta_{\mathrm{Ein}}>0.8"$. We find 40% of glSNe host galaxies have these
characteristics. For the Euclid survey, we require objects to have a VIS
magnitude888$\mathrm{VIS}=\frac{1}{3}(r+i+z)$, following Collett (2015)
brighter than 22.5 mag and Einstein radius $\theta_{\mathrm{Ein}}>0.4"$. We
find that 68% of the systems have these characteristics. Between the two
surveys, a total of 72% glSN hosts will be discoverable as lensing systems.
Table 3 shows these results for each SN subtype.

We note 25 % of lensed hosts are not discoverable with either of these
surveys. Upcoming deeper space-based surveys such as the Nancy Grace Roman
Space Telescope (Roman; Spergel et al. 2015) will be able to push the limiting
magnitude threshold required for detection, increasing the fraction of
discoverable glSN hosts and rates of glSNe (Pierel et al., 2021). The
combination of an Einstein radius threshold of $\theta_{\mathrm{Ein}}>0.4"$
and a magnitude threshold of 25 mag., merging both of the advantages of LSST
and Euclid, would allow for the discovery of 82% of lensed hosts. This means
18% of systems will still be too compact to be detectable, but these systems
are less likely to be interesting for time-delay cosmography and would pose
their own challenges for lens modelling.

It is worth noting that the rates cited require ideal lens discovery
conditions. Whilst substantial progress has been made in reducing billion
object surveys to samples of 1000s of lens candidates (More et al., 2016;
Jacobs et al., 2019; Li et al., 2021; Rojas et al., 2023) scaling strong lens
discovery systems up to Euclid and LSST has not yet been demonstrated. We
examine how much our values change for more pessimistic discovery prospects.
For the case of LSST, if we adjust the minimum Einstein radius threshold to
$\theta_{\mathrm{Ein}}>1"$, we find 25% of lensed glSN hosts are now
discoverable, as opposed to our previous 40 %. When we increase the Einstein
radius threshold for our Euclid lenses to $\theta_{\mathrm{Ein}}>0.5"$, the
loss is less severe, with 59% of the strongly lensed hosts still being
discoverable (c.f. 68% before). Alternatively, we can consider the case in
which we are able to meet the Einstein radius discovery threshold but do not
reach the discovery depth. For the case of LSST, reducing the magnitude
threshold to 24 still yields $\approx 40\%$ of lensed hosts, with a negligible
decrease in sample size. Euclid rates are much more sensitive to a loss of
depth: reducing the VIS threshold depth by a magnitude (to 21.5) loses a fifth
of the Euclid discoverable lens sample, with 52% of all lensed hosts now being
discoverable in Euclid data alone.

In the worst case scenario, with the more pessimistic thresholds in for both
Einstein radii and limiting magnitude, we find $54\%$ of lensed hosts are
discoverable, with $25\%$ coming from LSST and $46\%$ from Euclid, and an
overlap of 17% between both surveys. These results highlight that many of the
lensed SNe host systems will be close to detection thresholds as galaxy-galaxy
lenses. Combining surveys that reach fainter depths and better angular
resolution for galaxy-galaxy lens discovery will be important to mitigate risk
and maximise the discovery potential for galaxy-galaxy lensing.

| Total | SN Ia | SN Ibc | SN II | SNIIn | SLSN
---|---|---|---|---|---|---
Discoverable lenses (LSST + Euclid) | 75% | 72% | 81% | 76% | 73% | 43%
LSST | 40% | 38% | 44% | 42% | 39% | 30%
Euclid | 68% | 66% | 79% | 69% | 67% | 23%
Both | 34% | 32% | 43% | 35% | 33% | 10%
Total fraction with a lensed discoverable host | 25% | 22% | 21% | 27% | 20% | 16%
Table 3: Fraction of lensed hosts with $\mu_{\mathrm{TOT}}>3$ discoverable
with different imaging surveys Figure 4: Distribution of detectable lensed
host galaxies as a function of source i-band magnitude and Einstein radius.
The nominal discovery regions of LSST and Euclid are shown in red and orange
respectively. The left panel shows all glSNe host systems whereas the right
panel only shows golden glSNe for cosmology with long time delays and large
image separations. Since the golden sample is biased towards larger image
separations a larger fraction have a discoverable host.

We also consider whether the properties of the discoverable lenses differ from
those that are not discoverable, and explore how both of these populations
differ from the population of systems with unlensed hosts, shown in Fig. 5. As
expected, the main differences between the population of discoverable lensed
hosts and not discoverable lensed hosts are driven by our angular separation
requirement. The median Einstein radius of discoverable lenses is
$\theta_{\mathrm{Ein}}\approx 0.84"$, whereas the median Einstein radius of
not discoverable lenses is $\theta_{\mathrm{Ein}}\approx 0.67"$. The median
Einstein radius of unlensed systems peaks at lower Einstein radii than the
original population of discoverable glSNe, with $\theta_{\mathrm{Ein}}\approx
0.44"$. The median time-delay for unlensed hosts is 13 days, whereas for
discoverable lenses is 29 days and for not discoverable lenses it is 23 days.
The redshift distribution of sources peaks at $z_{s}\approx 0.69$ for
discoverable sources and $z_{s}\approx 0.72$ for not discoverable lenses, as
opposed to $z_{s}\approx 0.63$ for unlensed hosts. The properties of the
lenses in our sample are once again determined by our Einstein radius cut, as
slightly higher velocity dispersions are preferred for discoverable systems,
with mean $\sigma_{v}=229$ km s-1. Not discoverable lenses have
$\sigma_{v}=211$ km s-1 and unlensed hosts have a mean $\sigma_{v}=187$ km
s-1. These results show that the Einstein radius cut for discoverability has a
significant impact on the glSNe properties.

Figure 5: Properties of lensed SN systems, broken down into subpopulations.
Red shows glSNe where the host is not multiply imaged. For Yellow, the host is
multiply imaged and discoverable as a galaxy-galaxy lens in either Euclid or
LSST. Blue shows systems where the host is multiply imaged but not
identifiable as such with LSST or Euclid. Purple are the glSNe with lensed
hosts, time delay > 10 days and image separation > 0.8 arcseconds. Overall,
lensed systems have greater Einstein radii and longer time-delays.

Figure 6 shows the distribution of glSN offsets from their host galaxy centre
for strongly lensed hosts and unlensed hosts. We show that glSNe in galaxy-
galaxy lensing systems happen closer to the centre of their host galaxy, with
a median separation of 0.19 half-light radii whereas those in unlensed hosts
have a median separation of 0.35 half-light radii. This is expected, as the
closer the SN is to the host galaxy, the more likely it is for the host to
also be strongly lensed. For very compact systems where the supernova happens
close to its host, this will pose the additional challenge of being able to
resolve the supernova from the nucleus of the host galaxy.

Figure 6: SN offset from its host centre in units of host half-light radii.
SNe in strongly lensed hosts are located closer to the nucleus of their host
galaxy than those in unlensed hosts.

###  3.3 Implications for precision cosmology

A fraction of discoverable glSNe will belong to the "golden sample" for time-
delay cosmography, which is defined as the systems for which we believe we can
confidently measure a precise ($\sim$10% precision) time-delay measurement. We
define this as systems with image separations of $\theta_{\mathrm{max}}$> 0.8”
and $\Delta t_{\mathrm{max}}$ >10 days, due to the fact that we can measure
time-delays at the $\approx 1$d level (Pierel & Rodney 2019). So far these
requirements have been simply defined for glSNe Ia, partly due to their
standardisable nature and also due to the fact they are the only subtype that
has been detected so far in galaxy-scale lenses. For this reason, in this
subsection we will focus on glSNe Ia. It is worth noting that with the
increased detection threshold of LSST, more glSNe of all subtypes will be
observable, and time-delay measurements will be obtainable from other SN
subtypes, although our current ability to measure time-delays from other
subtypes will most likely require longer time delays (Huber et al., 2022;
Pierel et al., 2021).

As we have shown in Section 3.2, glSNe in strongly lensed systems have
slightly greater angular separations and time-delays than those systems in
which the host is not strongly lensed. This is because of a selection bias:
the physical separation between host and SN is independent of the lens, but
the larger the Einstein radius of the lens the greater the chance that both
host and SN are multiply imaged. Time delays and image separations scale with
the Einstein radius.

40% of glSNe Ia in the golden sample have a lensed host with
$\mu_{\mathrm{TOT,host}}>3$, compared to $30\%$ of the whole sample, and the
majority of these (91%) systems will be discoverable with either LSST or
Euclid. This is because the golden sample requires long time delays and image
separation which implies larger than typical Einstein radii. Fig. 5b shows the
Einstein radius of lenses and lensed source magnitudes for the golden sample.

47% of glSNe Ia in our simulations, that is $\approx 18$ per year, belong to
this "golden sample". According to this $\approx 7$ glSNe Ia per year will
have the properties suitable for time-delay cosmography and occur in a system
that could already have been identified as a galaxy-galaxy lens. Over the 10
years of the survey, this will add up to 70 glSNe and assume a $10\%$ error
$H_{0}$ measurement per system, this would imply an $\approx 1.3\%$ precise
measurement of the Hubble constant from glSNe detected through monitoring
galaxy-galaxy lenses.

###  3.4 Time of discovery and follow-up

Normally, multiple high signal-to-noise ratio (S/N) observations are required
before claiming that a new transient could be lensed. When monitoring known
lensing systems, if a new transient is discovered we can immediately classify
the transient as interesting since the chance of it being lensed is much
higher than for a random field. This would enable much faster spectroscopic
confirmation and follow-up observations to begin.

In this section, we explore how many more glSNe can be detected in
discoverable lenses if we reduce the number of detections required for
"discovery". We compare the results from our initial analysis, which requires
three 5$\sigma$ observations, with the more relaxed requirements of one
5$\sigma$ and one 3$\sigma$ observation. We first look at the change in the
fraction of simulated systems that were detected as a function of the glSN
redshift for each of our detection cuts, shown in Fig. 7. Relaxing our
detection cuts from three to just one 5$\sigma$ detection allows for an 8%
increase in rates. Setting our detection requirement to a single 3$\sigma$
detection not only increases rates by 30% with respect to our initial cuts,
but also increases the maximum redshift at which glSNe can be detected. This
is because normally higher redshift observations have lower S/N.

The success of time-delay cosmography also relies on our ability to detect
transients as soon as possible in order to request follow-up observations, as
other studies have shown that LSST alone does not have the cadence required
for a good enough sampling of light curves to do precision time-delay
measurements without follow-up from other telescopes (Huber et al., 2019). We
first consider at what time in the observer-frame glSNe can be detected. We
find 29% glSNe Ia are discovered 10 days before their B-band peak. Requiring
just one detection this fraction becomes 42%, and 45% for a 3$\sigma$
detection. Ideally, we would detect our glSNe 20 observer-frame days before
peak in order to maximise the chance of obtaining follow-up observations. For
the three 5$\sigma$ requirement only 4% ($\approx 2$ per year) satisfy this
requirement. For the one 5$\sigma$ requirement this fraction increases to 8%,
and further increases to 15% for the one 3$\sigma$ requirement. This is driven
by the fact the 3$\sigma$ detection allows us to detect glSNe at higher
redshift, which have a slower evolution as their lightcurve is stretched due
to time dilation.

We convert these times to rest-frame to find at what time in their evolution
glSNe will be detected. The distribution of discovery phases of the first
images occurring in discoverable lenses is shown in Fig. 8. For our 3
5$\sigma$ detection requirement, 14 % of glSNe Ia are discovered 10 days
before their B-band peak. Requiring just one 5$\sigma$ detection, this
fraction becomes 22% and 26% for a single 3$\sigma$ detection.

Figure 7: Rate of glSNe Ia discovered in LSST for different discovery
detection requirements. Reducing the S/N requirement substantially increases
the yield of glSNe past $z$=1 Figure 8: Probability density function for the
time of discovery of of glSNe Ia in LSST. The top (bottom) panel shows time
relative to peak in the rest-frame (observer-frame).

####  3.4.1 False positives - the rate of unlensed SNe in lensing galaxies

A lower S/N requirement could be beneficial for finding more glSNe, but also
comes at the cost of a higher false positive rate. When a transient is
discovered in a known lensing system it could be a glSN, but it could also be
an unlensed supernova in the lens galaxy. With a significant number of
observations in different colour bands, this second scenario can be discarded
by fitting observations to supernova templates at the redshift of the lens.
Robustly distinguishing this scenario however requires a considerable number
of observations to compare against the SN template. It cannot be done with a
single detection.

In order to estimate the rate of false positives, we first calculate the
observer-frame SN rate in the foreground galaxies. In order to do this we
follow the method described by Li et al. (2012) to obtain the rest-frame SN Ia
rate for each lens galaxy, which depends on both the stellar mass and star-
formation rate of the galaxy.

The stellar mass for each of the lenses is determined from their velocity
dispersion following Cannarozzo et al. (2020). We choose the expression for
the velocity dispersion ($\sigma_{v}$)-stellar mass($M_{\star}$) relation to
be

$\mathrm{log}\left(\frac{\sigma_{v}}{\mathrm{km}\mathrm{s^{-1}}}\right)\approx
2.21+0.18\ \mathrm{log}\left(\frac{M_{\star}}{10^{11}M_{\odot}}\right)+0.4\
\mathrm{log}(1+z).$ (14)

In their work, this expression is valid for a fiducial sample with redshift
$z<1$. We choose this due to the fact most of our lenses will be found up to
redshift $z$ $\approx$ 1 and LSST will not be able to discover unlensed SNe Ia
beyond that redshift. As expected, our lenses are massive galaxies, with a
median stellar mass log $M_{\star}/M_{\odot}$ = $11.2^{+0.6}_{-0.7}$. The
rest-frame yearly SN Ia rate ($\textrm{SNR}_{\textrm{Ia}}$) is given by a two-
component model (Scannapieco & Bildsten, 2005; Mannucci et al., 2005; Sullivan
et al., 2006)

$\left(\frac{\textrm{SNR}_{\textrm{Ia}}}{\textrm{yr}^{-1}}\right)=\hat{A}\cdot
10^{-10}\left(\frac{M_{\star}}{\textrm{M}_{\odot}}\right)^{\alpha}+\hat{B}\cdot
10^{-3}\left(\frac{\textrm{SFR}}{\textrm{M}_{\odot}\textrm{yr}^{-1}}\right).$
(15)

We define the parameters $\hat{A}=1.05\pm 0.16$ and $\alpha=0.68\pm 0.01$
following Smith et al. (2012). The parameter $\hat{B}$ can be related to the
$\textrm{SNR}_{\textrm{cc}}-\textrm{SFR}$ relation

$\hat{B}=k_{\textrm{cc}}\Theta,$ (16)

where $\Theta=\textrm{SNR}_{\textrm{Ia}}/\textrm{SNR}_{\textrm{cc}}$. At
redshift $z<1$, the ratio of
$\textrm{SNR}_{\textrm{Ia}}/\textrm{SNR}_{\textrm{cc}}$ approximately ranges
between $\Theta=1/2$ and $\Theta=1/4$, so we use $\Theta=1/4$ for the lower
limit on our rates and $\Theta=1/2$ for the upper limit. The parameter
$k_{\textrm{cc}}$ can be determined by measuring the ratio between
$\textrm{SNR}_{\textrm{cc}}$ and SFR. Using observational data this ratio has
been constrained to be $k_{\textrm{cc}}=7.5\pm 2.5$ (Scannapieco & Bildsten,
2005). We use $k_{\textrm{cc}}=10$ as the upper limit, and $k_{\textrm{cc}}=5$
as the lower limit.

As our lenses are red elliptical galaxies, we only focus on the rate component
that is dependent of the stellar mass of the system. We find our lenses have a
median SN Ia rate of $\approx$ 0.005 SN Ia yr-1. For each lens, we use the
observer-frame SN rate to estimate what fraction of SNe could actually be
discovered by simulating a sample of $10^{6}$ SNe Ia in the redshift range
$z\in[0,1]$ and determining the fraction of these that would be discoverable.
With our total detectable rate, given 10 000 lenses in LSST, the supernova
rate in these lenses is of 229 SNe yr-1. The rate of discoverable SNe as a
function of the lens Einstein radius is shown in Fig. 9.

Figure 9: Rate of unlensed SN Ia occurring in the lensing galaxy as a function
of Einstein radius of the system.

Following a similar method, we use Eqn. 15 to calculate the SFR for every host
galaxy. The host galaxies in our sample are less massive than our lenses, with
median stellar mass log $M_{\star}/M_{\odot}$ = $9.9^{+0.7}_{-0.9}$. These
younger populations do undergo star-formation, with median log SFR =
$0.0^{+0.6}_{-0.73}$. We find the median SN Ia rate per host galaxy to be
between 0.002-0.005 SN Ia per rest-frame year, from our lower and upper limits
respectively. We convert these to observer-frame rates and use our results
from simsurvey to estimate the rate of glSNe in every host galaxy that will be
discoverable.

We define the false positive rate (FPR) of this search method as

$\mathrm{FPR}=\frac{1}{N_{\mathrm{sim}}}\sum_{i}^{\mathrm{N_{sim}}}\frac{\mathrm{FP}_{i}}{\mathrm{FP}_{i}+\mathrm{TP}_{i}},$
(17)

where $N_{\mathrm{sim}}$ is the number of lenses in our simulation. For each
lensed system $i$, $\mathrm{FP}$ represents the observer-frame rate of
detectable unlensed SNe Ia in the lens and $\mathrm{TP}$ is the observer-frame
rate of detectable glSNe Ia in the source. We find FPR = 0.64 – 0.77 depending
on our rate parameters, which is consistent with the overall low probability
of having a lensing event. Most SNe detections around lenses will actually
correspond to unlensed transients in the lens galaxy. Relaxing our detection
requirements to one 3$\sigma$ detection increases this FPR further to 0.84.
Whilst such a false positive rate is likely tolerable given the high
scientific value of lensed SNe, further cuts could be introduced to reduce the
false positive rates, for example the position of the transient within the
lensed system (first images form close to or beyond the Einstein radius, where
the lens stellar density is lower) or the colour and magnitude at discovery.

##  4 Limitations

The forecasts presented in this work are based on Monte Carlo simulations
which, despite being well motivated by observations, are based on a series of
assumptions that could bias our predictions. These should therefore be seen as
a general estimate. In this subsection we comment on the main limitations of
our work and what effect this could have on our predictions.

###  4.1 Source population

Our simulations make use of the SCOTCH catalogue to establish correlations
between properties of SNe and their host galaxies. Lokken et al. (2023) list a
series of assumptions that were required to simulate the source catalogue.
While they discuss in detail the implications of these assumptions, the most
relevant one to this work is that the catalogue does not account for any kind
of redshift evolution of the host galaxy properties. This is due to the fact
SCOTCH itself is based on the GHOST catalogues (Gagliano et al., 2021), which
contain observational properties from ground-based surveys and are mostly
complete up to redshift $z$ < 0.6 (but not strictly volume or redshift-limited
to $z$ < 0.6). Thus, these catalogues probably do not accurately represent the
supernova-host relationships beyond redshift $z$ > 0.6. Recent work by Popovic
et al. (2024) shows evidence for evolving galaxy properties with redshift,
which had already been hinted at in the past (Nicolas et al., 2021) but was
hard to determine partly due to Malmquist bias (Popovic et al., 2021). The
discovery of larger samples of SNe at higher redshifts with upcoming deeper
photometric surveys will allow us to confirm whether there exists some
redshift evolution in host-transient correlations.

####  4.1.1 Rates and luminosity functions

Aside from the relationship between the host galaxy and the supernova, the
properties of the supernova itself will also have an impact on our results.
This becomes evident when considering our rate forecasts, which are roughly
consistent with previous rate estimates (e.g. Goldstein et al. 2019, Wojtak et
al. 2019, Arendse et al. 2023, Sainz de Murieta et al. 2023). One exception is
our rate estimates for SNe IIn. Previous rate estimates find that this subtype
of core-collapse supernovae will dominate glSN rates, due to their bright
luminosity functions that extend to peak absolute magnitudes brighter than
$M_{B}\gtrsim-19$. The rates of discoverable SNe rely heavily on a wide range
of assumptions about the source population, including the SED and probability
distribution of the intrinsic luminosity.

One of the key differences between our study and previous works is the choice
of volumetric rates in our simulations. A comparison between the different
rates used in different glSN yield studies are shown in Fig. 10. Core-collapse
supernovae rates follow the star formation history (SFH) (e.g. Madau &
Dickinson 2014). Throughout this work we used the same volumetric rates as
Kessler et al. (2019), which are given by Strolger et al. (2015). These
results are consistent with the SFR densities derived from dust-corrected UV
emission, but disagree with SN rates that follow the SFH evolution from Madau
& Dickinson (2014). Wojtak et al. (2019) use a model for core-collapse rates
that is directly proportional to the star formation history by Madau &
Dickinson (2014). Goldstein et al. (2019) use the core-collapse volumetric
rates by Li et al. (2011), which state an evolution of the rates
$\propto(1+z)^{3.6}$. These rates show a very different evolution with
redshift, and so a much larger number of core-collapse supernovae beyond
redshift $z\sim 0.8$.

There also exist some discrepancies in the choice of volumetric SN Ia rates.
While we again use the volumetric rates from Kessler et al. (2019), based on
two different observational studies of supernovae, Wojtak et al. (2019) used a
more theoretical approach based on the star formation rate and the delay-time
distribution. Goldstein et al. (2019) used the volumetric rates from Sullivan
et al. (2006). These rates, once again, show a higher number of SNe Ia in the
Universe. In practice, it is highly likely that none of our volumetric rate
estimates are a completely accurate representation of the volumetric rate of
supernovae in our Universe. The discovery of more glSNe at higher redshifts
will allow for a better understanding of the redshift evolution of the
volumetric rate.

Figure 10: Volumetric SN rates used in this work (from Kessler et al. 2019),
compared with those used in Wojtak et al. 2019 and Goldstein et al. 2019. The
rates in this worked are represented by the blue lines. The rates we have used
show a smaller number density of supernovae at high redshift. Figure 11:
Luminosity function used in this work for the peak magnitude of different
core-collapse SNe (based on Kessler et al. 2019) compared with those used by
Goldstein et al. 2019 and Wojtak et al. 2019.

Another difference is the choice of luminosity function, particularly for
core-collapse SNe. Figure 11 shows the luminosity function for the SN
templates in this work (based on Kessler et al. 2019), Wojtak et al. (2019)
and Goldstein et al. (2019). The luminosity distribution of SN II used in this
work has a narrower spread, which means that even if the median peak magnitude
of our population of SN II is brighter, the ones in other studies extend to
brighter absolute magnitudes. This is due to the templates used by Kessler et
al. (2019). Conversely, our populations of SNIbc and SNIIn are significantly
fainter, with differences reaching up to four magnitudes. These faint glSNe
would require extremely high magnifications in order to push them beyond the
detection threshold, which are more unlikely as they can only take place if
the SN is located in a region near the caustic. It is therefore unsurprising
that our rate estimates are much lower.

In order to quantify the impact of luminosity functions and volumetric rate
assumptions on rate forecasts, we change some of our initial assumptions and
study their impact on rates. For SNe Ia, which have better understood
luminosity functions, we simply use different volumetric rates. Using the
rates from Wojtak et al. (2019) results in a $\approx 15\%$ increase in glSNe
Ia rates. The use of volumetric rates from Goldstein et al. (2019) results in
a 5% decrease in yields. For SN II, which dominate our rate forecasts, we need
to take into account both the differences in rates and luminosity function. We
re-run our simulations with the different assumptions chosen by Wojtak et al.
(2019) and Goldstein et al. (2019). As expected given Fig. 10 and Fig. 11, the
use of these rates results in a large increase in rates. Using the assumptions
from Goldstein et al. (2019) results in a 65% increase in glSN II rates. The
assumptions in Wojtak et al. (2019) result in a 20% increase in rates. Despite
the differences in rate forecasts, their impact on our estimates for the
fraction of lensed hosts is negligible.

These differences highlight the intrinsic uncertainty in glSN rate forecasts,
and the importance of detailing the assumptions such forecasts are built upon.
Upcoming deeper surveys will allow us to increase our sample of core-collapse
supernovae, which will improve our understanding of their properties at higher
redshifts.

####  4.1.2 Scatter due to dust and microlensing

Our simulations accounted for extinction from the Milky Way and dust in both
the host and lens galaxies. The uncertainty coming from the distribution of
dust in galaxies is a source of significant scatter, especially for core-
collapse supernovae, which happen in star-forming and dusty environments
(Hatano et al., 1998). Along with their dimmer lightcurves, dust is believed
to have a big impact on the lower rates of core-collapse supernovae (Mannucci
et al., 2003). Our treatment of dust was the same for both SNe Ia and core-
collapse supernovae, therefore it is likely that we underestimated the effects
of dust, making the core-collapse rates an upper limit.

For simplicity, we did not account for the effects of microlensing. In most
scenarios, this will slightly demagnify the glSNe, with some glSNe being more
magnified (example histograms for the distribution of microlensing
magnifications can be found in Weisenbach et al. 2021, 2024). Some recent rate
forecasts (Arendse et al., 2023; Sainz de Murieta et al., 2023) study the
effect of microlensing on glSNe yields and find that for surveys with the
depth of LSST, microlensing demagnification pushes some supernovae below the
detection threshold, leading to an 8% decrease in rates. We introduce a simple
microlensing model in order to study whether this would have a significant
impact on our rates. We re-simulated our glSN lightcurves adding a 20% scatter
in the flux of the unlensed model. This means some supernovae will be up to
20% dimmer than originally simulated but others will be up to 20% brighter. We
find a 9% decrease in our glSN estimates, which is consistent with Arendse et
al. (2023); Sainz de Murieta et al. (2023). We also checked the effect of
microlensing on our lensed host fraction estimates and found that accounting
for microlensing had a negligible impact on our results.

###  4.2 Deflector population

Another source of uncertainty will come from the choice of lens model, which
will have an effect on the time-delays and magnifications of our supernovae.
Our lens population was simulated according to Collett (2015). We assumed SIE
profiles for all lens galaxies, which has been shown to be a good
approximation for massive elliptical galaxies (Auger et al., 2009). We did not
include any kind of evolution on the velocity dispersion with redshift. This
is consistent with results from Bezanson et al. (2011) and Shu et al. (2012),
at redshifts up to $z\leq 0.5$, where most of our lenses are found. Recent
work by Ferrami & Wyithe (2024) has shown the redshift dependence of the
velocity dispersion function of the deflector population can have an impact on
strong lensing estimates. If we instead use the velocity dispersion function
from Geng et al. (2021), which displays a negative evolution of the number of
deflectors with redshift, our forecast rate of glSNe Ia decreases by $33\%$.
However, because this velocity dispersion function removes high redshift
lenses, we actually see a slight increase in the fraction of glSNe Ia systems
with a multiply imaged host: the fraction discoverable with LSST and/or Euclid
increases by $8\%$. We conclude that the uncertainty on the $z>0.5$ velocity
dispersion function produces substantial uncertainty on the forecast rate of
glSNe, but it is less important for the rate of glSNe in discoverable lenses.

###  4.3 Discoverability criteria

Our actual ability to discover lenses, as well as our ability to identify
glSNe will also have an impact on our rate estimates. In Section 3.2 we have
already quantified the impact of our Einstein radius and magnitude cuts on the
rate of lensed glSN hosts. We showed a difference of $\approx$ 20% in the
number of lensed hosts that are discoverable between our more optimistic and
more pessimistic scenarios.

Finally, we made the assumption that all the glSNe can be resolved from both
the lens and the host, without loss of discovery depth. However, as shown in
Fig. 6, our findings indicate that discoverable glSNe in galaxy-galaxy lenses
generally occur closer to the centres of their hosts. As a result, some of
these glSNe may not be detected in the LSST difference image analysis
pipelines. Quantifying this effect is beyond the scope of this paper, but it
is likely that this will decrease the rate of glSNe discovered in LSST.

##  5 Conclusions

In this work we have investigated the hosts of lensed SNe. We built a
population of galaxy-scale lenses and put them in front of simulated SNe and
their host galaxies using the SCOTCH catalogue (Lokken et al., 2023). We did
this to answer two questions: how often are the hosts of glSNe also strongly
lensed, and how many glSNe should we expect to discover by monitoring the
population of galaxy-galaxy lenses that Euclid and LSST will soon deliver. We
simulated SNe type Ia, Ibc, II, IIn and SLSNe. Our model predicted annual
rates of 88 with LSST, with type II supernovae being the dominant subtype
followed closely by type Ia supernovae.

A fraction of these glSNe, however, will not have a strongly lensed host.
Thus, even if their time-delays are long enough to be measured precisely, it
will not be possible to use them to infer the value of $H_{0}$, as lens
modelling is needed to convert time delays into a value of the Hubble
constant. We determined that in $54\%$ of cases, both glSNe and their host
galaxy are lensed. $25\%$ of all glSNe should have a multiply imaged host that
is sufficiently bright, magnified and resolved that they would be detectable
as a galaxy-galaxy lens in at least one of the wide-field surveys of Euclid or
LSST. Because such systems have larger Einstein radii, they typically have
longer time delays and are more suited for cosmography.

Monitoring known galaxy-galaxy lenses can potentially allow us to be less
strict in discovery thresholds since any change in these systems is either a
lensed transient or a false positive transient in the lensing galaxy. Changing
from three 5$\sigma$ detections to a single 3$\sigma$ detection would increase
glSNe rates by 30% as well as allow us to discover them 4 days earlier on
average. Discovering glSNe early gives us more time to take follow-up with
other instruments and study the early evolution of glSNe. One drawback of this
approach is the false positive rate: 84% of $3\sigma$ changes will be due to
unlensed SNe in the lens galaxy.

Our simulations rely on assumptions regarding the mass profile of lenses, the
volumetric rates of unlensed transients, and supernova lightcurve templates.
We have quantified the effect of some of these assumptions and concluded that
they have a bigger impact on rate estimates than on the fractions of glSNe
with lensed hosts. The different choices likely explain why published forecast
glSNe rates can differ by a factor of two (Sainz de Murieta et al., 2023;
Goldstein et al., 2019; Arendse et al., 2023; Wojtak et al., 2019).

Monitoring known gravitational lenses for new transients is an established
field, but is currently limited by the number of known lenses. In the LSST
era, only a quarter of all glSNe will occur in galaxy-galaxy lenses that could
plausibly be discovered beforehand in Euclid and/or LSST, but this rises to a
third for the glSNe that are golden for cosmology. Finding and monitoring
large samples of galaxy-galaxy lenses can therefore give us a realistic
pathway to $H_{0}$ with glSNe even if searching for glSNe in the full LSST
alert stream proves challenging. We have found that 7 glSNe Ia suitable for
cosmology should be discoverable with LSST per year by monitoring the galaxy-
galaxy lenses we expect to discover in LSST and Euclid. Achieving this yield
should result in a 1.3% precision measurement of $H_{0}$ after 10 years,
assuming we have the galaxy-galaxy watch list in place at the start of the
LSST survey.

## Acknowledgements

We are grateful to Nikki Arendse, Suhail Dhawan and Anupreeta More for helpful
conversations that have enriched this work. This project has received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme (LensEra: grant agreement No 945536).

T. E. C. is funded by a Royal Society University Research Fellowship. This
work is supported by the National Science Foundation under Cooperative
Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and
Fundamental Interactions, http://iaifi.org/).

## Data Availability

All data are simulated: our simulations are available from the corresponding
author upon request.

## References

  * Aleo et al. (2023) Aleo P. D., et al., 2023, ApJS, 266, 9
  * Arendse et al. (2023) Arendse N., et al., 2023, arXiv e-prints, p. arXiv:2312.04621
  * Auger et al. (2009) Auger M. W., Treu T., Bolton A. S., Gavazzi R., Koopmans L. V. E., Marshall P. J., Bundy K., Moustakas L. A., 2009, ApJ, 705, 1099
  * Bailey et al. (2023) Bailey A. C., Vincenzi M., Scolnic D., Cuillandre J. C., Rhodes J., Hook I., Peterson E. R., Popovic B., 2023, MNRAS, 524, 5432
  * Bernardi et al. (2010) Bernardi M., Shankar F., Hyde J. B., Mei S., Marulli F., Sheth R. K., 2010, MNRAS, 404, 2087
  * Bezanson et al. (2011) Bezanson R., et al., 2011, ApJ, 737, L31
  * Bianco et al. (2021) Bianco F. B., et al., 2021, A&AS, 258, 1
  * Birrer & Amara (2018) Birrer S., Amara A., 2018, Physics of the Dark Universe, 22, 189
  * Birrer et al. (2021) Birrer S., et al., 2021, The Journal of Open Source Software, 6, 3283
  * Birrer et al. (2022) Birrer S., Millon M., Sluse D., Shajib A. J., Courbin F., Koopmans L. V. E., Suyu S. H., Treu T., 2022, arXiv e-prints, p. arXiv:2210.10833
  * Bonvin et al. (2016) Bonvin V., et al., 2016, Monthly Notices of the Royal Astronomical Society, 465, 4914
  * Cannarozzo et al. (2020) Cannarozzo C., Sonnenfeld A., Nipoti C., 2020, Monthly Notices of the Royal Astronomical Society, 498, 1101
  * Chen et al. (2022) Chen W., et al., 2022, Nature, 611, 256
  * Collaboration et al. (2009) Collaboration L. S., et al., 2009, LSST Science Book, Version 2.0 (arXiv:0912.0201)
  * Collett (2015) Collett T. E., 2015, ApJ, 811, 20
  * Di Valentino et al. (2021) Di Valentino E., et al., 2021, Classical and Quantum Gravity, 38, 153001
  * Dilday et al. (2008) Dilday B., et al., 2008, ApJ, 682, 262
  * Dobler & Keeton (2006) Dobler G., Keeton C. R., 2006, ApJ, 653, 1391
  * Feindt et al. (2019) Feindt U., Nordin J., Rigault M., Brinnel V., Dhawan S., Goobar A., Kowalski M., 2019, J. Cosmology Astropart. Phys., 2019, 005
  * Ferrami & Wyithe (2024) Ferrami G., Wyithe S., 2024, arXiv e-prints, p. arXiv:2404.03143
  * Fitzpatrick (1999) Fitzpatrick E. L., 1999, Publications of the Astronomical Society of the Pacific, 111, 63
  * Frye et al. (2023) Frye B., et al., 2023, Transient Name Server AstroNote, 96, 1
  * Gagliano et al. (2021) Gagliano A., Narayan G., Engel A., Carrasco Kind M., LSST Dark Energy Science Collaboration 2021, ApJ, 908, 170
  * Geng et al. (2021) Geng S., Cao S., Liu Y., Liu T., Biesiada M., Lian Y., 2021, MNRAS, 503, 1319
  * Goldstein et al. (2019) Goldstein D. A., Nugent P. E., Goobar A., 2019, ApJS, 243, 6
  * Goobar et al. (2017) Goobar A., et al., 2017, Science, 356, 291
  * Goobar et al. (2023) Goobar A., et al., 2023, Nature Astronomy, 7, 1098
  * Grillo et al. (2024) Grillo C., Pagano L., Rosati P., Suyu S. H., 2024, A&A, 684, L23
  * Hakobyan et al. (2009) Hakobyan A. A., Mamon G. A., Petrosian A. R., Kunth D., Turatto M., 2009, A&A, 508, 1259
  * Hatano et al. (1998) Hatano K., Branch D., Deaton J., 1998, The Astrophysical Journal, 502, 177
  * Hills (1975) Hills J. G., 1975, Nature, 254, 295
  * Hounsell et al. (2018) Hounsell R., et al., 2018, ApJ, 867, 23
  * Huber et al. (2019) Huber S., et al., 2019, A&A, 631, A161
  * Huber et al. (2022) Huber S., et al., 2022, A&A, 658, A157
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Jacobs et al. (2019) Jacobs C., et al., 2019, MNRAS, 484, 5330
  * Kasen et al. (2017) Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Nature, 551, 80
  * Kelly et al. (2015) Kelly P. L., et al., 2015, Science, 347, 1123
  * Kelly et al. (2022) Kelly P., et al., 2022, Transient Name Server AstroNote, 169, 1
  * Kelly et al. (2023) Kelly P. L., et al., 2023, Science
  * Kessler et al. (2009) Kessler R., et al., 2009, Publications of the Astronomical Society of the Pacific, 121, 1028
  * Kessler et al. (2019) Kessler R., et al., 2019, PASP, 131, 094501
  * Kochanek (2016) Kochanek C. S., 2016, MNRAS, 461, 371
  * Kormann et al. (1994) Kormann R., Schneider P., Bartelmann M., 1994, A&A, 284, 285
  * Korytov et al. (2019) Korytov D., et al., 2019, ApJS, 245, 26
  * Kundić et al. (1997) Kundić T., et al., 1997, ApJ, 482, 75
  * Li & Paczyński (1998) Li L.-X., Paczyński B., 1998, ApJ, 507, L59
  * Li et al. (2011) Li W., et al., 2011, MNRAS, 412, 1441
  * Li et al. (2012) Li X., Hjorth J., Richard J., 2012, J. Cosmology Astropart. Phys., 2012, 015
  * Li et al. (2021) Li R., et al., 2021, ApJ, 923, 16
  * Lokken et al. (2023) Lokken M., et al., 2023, MNRAS, 520, 2887
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Magee et al. (2023) Magee M. R., Sainz de Murieta A., Collett T. E., Enzi W., 2023, MNRAS, 525, 542
  * Mannucci et al. (2003) Mannucci F., et al., 2003, A&A, 401, 519
  * Mannucci et al. (2005) Mannucci F., Della Valle M., Panagia N., Cappellaro E., Cresci G., Maiolino R., Petrosian A., Turatto M., 2005, A&A, 433, 807
  * Millon et al. (2020) Millon M., et al., 2020, A&A, 640, A105
  * Modjaz et al. (2019) Modjaz M., Gutiérrez C. P., Arcavi I., 2019, Nature Astronomy, 3, 717
  * More et al. (2016) More A., et al., 2016, MNRAS, 455, 1191
  * Moriya et al. (2018) Moriya T. J., Sorokina E. I., Chevalier R. A., 2018, Space Sci. Rev., 214, 59
  * Nicolas et al. (2021) Nicolas N., et al., 2021, A&A, 649, A74
  * Pierel & Rodney (2019) Pierel J. D. R., Rodney S., 2019, ApJ, 876, 107
  * Pierel et al. (2018) Pierel J. D. R., et al., 2018, PASP, 130, 114504
  * Pierel et al. (2021) Pierel J. D. R., Rodney S., Vernardos G., Oguri M., Kessler R., Anguita T., 2021, ApJ, 908, 190
  * Pierel et al. (2024) Pierel J. D. R., et al., 2024, ApJ, 967, L37
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Popovic et al. (2021) Popovic B., Brout D., Kessler R., Scolnic D., Lu L., 2021, ApJ, 913, 49
  * Popovic et al. (2024) Popovic B., et al., 2024, arXiv e-prints, p. arXiv:2406.06215
  * Prajs et al. (2016) Prajs S., et al., 2016, Monthly Notices of the Royal Astronomical Society, 464, 3568
  * Quimby et al. (2014) Quimby R. M., et al., 2014, Science, 344, 396
  * Refsdal (1964) Refsdal S., 1964, MNRAS, 128, 307
  * Riess et al. (2022) Riess A. G., et al., 2022, ApJ, 934, L7
  * Rodney et al. (2014) Rodney S. A., et al., 2014, AJ, 148, 13
  * Rodney et al. (2021) Rodney S. A., Brammer G. B., Pierel J. D. R., Richard J., Toft S., O’Connor K. F., Akhshik M., Whitaker K. E., 2021, Nature, 5, 1118
  * Rojas et al. (2023) Rojas K., et al., 2023, Monthly Notices of the Royal Astronomical Society, 523, 4413
  * Rusu et al. (2020) Rusu C. E., et al., 2020, MNRAS, 498, 1440
  * Ryczanowski et al. (2020) Ryczanowski D., Smith G. P., Bianconi M., Massey R., Robertson A., Jauzac M., 2020, Monthly Notices of the Royal Astronomical Society, 495, 1666
  * Sagués Carracedo et al. (2024) Sagués Carracedo A., et al., 2024, arXiv e-prints, p. arXiv:2406.00052
  * Sainz de Murieta et al. (2023) Sainz de Murieta A., Collett T. E., Magee M. R., Weisenbach L., Krawczyk C. M., Enzi W., 2023, MNRAS, 526, 4296
  * Scannapieco & Bildsten (2005) Scannapieco E., Bildsten L., 2005, ApJ, 629, L85
  * Schechter et al. (1997) Schechter P. L., et al., 1997, ApJ, 475, L85
  * Schlegel et al. (1998) Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525
  * Scolnic et al. (2018) Scolnic D., et al., 2018, ApJ, 852, L3
  * Sheth et al. (2003) Sheth R. K., et al., 2003, ApJ, 594, 225
  * Shu et al. (2012) Shu Y., Bolton A. S., Schlegel D. J., Dawson K. S., Wake D. A., Brownstein J. R., Brinkmann J., Weaver B. A., 2012, AJ, 143, 90
  * Shu et al. (2018) Shu Y., Bolton A. S., Mao S., Kang X., Li G., Soraisam M., 2018, ApJ, 864, 91
  * Smith et al. (2012) Smith M., et al., 2012, ApJ, 755, 61
  * Spergel et al. (2015) Spergel D., et al., 2015, arXiv e-prints, p. arXiv:1503.03757
  * Strolger et al. (2015) Strolger L.-G., et al., 2015, ApJ, 813, 93
  * Sullivan et al. (2006) Sullivan M., et al., 2006, ApJ, 648, 868
  * Suyu et al. (2014) Suyu S. H., et al., 2014, ApJ, 788, L35
  * Treu (2010) Treu T., 2010, ARA&A, 48, 87
  * Treu et al. (2022) Treu T., Suyu S. H., Marshall P. J., 2022, A&ARv, 30, 8
  * Vincenzi et al. (2021) Vincenzi M., et al., 2021, Monthly Notices of the Royal Astronomical Society, 505, 2819
  * Walsh et al. (1979) Walsh D., Carswell R. F., Weymann R. J., 1979, Nature, 279, 381
  * Weisenbach et al. (2021) Weisenbach L., Schechter P. L., Pontula S., 2021, ApJ, 922, 70
  * Weisenbach et al. (2024) Weisenbach L., Collett T., Sainz de Murieta A., Krawczyk C., Vernardos G., Enzi W., Lundgren A., 2024, arXiv e-prints, p. arXiv:2403.03264
  * Wojtak et al. (2019) Wojtak R., Hjorth J., Gall C., 2019, MNRAS, 487, 3342
  * Wong et al. (2016) Wong K. C., et al., 2016, Monthly Notices of the Royal Astronomical Society, 465, 4895
  * Wong et al. (2020) Wong K. C., et al., 2020, MNRAS, 498, 1420
  * Wood-Vasey et al. (2007) Wood-Vasey W. M., et al., 2007, ApJ, 666, 694