# Simulating time-varying strong lenses

G. Vernardos
Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale
de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a self-consistent and versatile forward modelling software package
that can produce time series and pixel-level simulations of time-varying
strongly lensed systems. The time dimension, which needs to take into account
different physical mechanisms for variability such as microlensing, has been
missing from existing approaches and it is of direct relevance to time delay,
and consequently $H_{\mathrm{0}}$, measurements and caustic crossing event
predictions. Such experiments are becoming more streamlined, especially with
the advent of time domain surveys, and understanding their systematic and
statistical uncertainties in a model-aware and physics-driven way can help
improve their accuracy and precision. Here we demonstrate the software’s
capabilities by exploring the effect of measuring time delays from lensed
quasars and supernovae in many wavelengths and under different microlensing
and intrinsic variability assumptions. In this initial application, we find
that the cadence of the observations and combining information from different
wavelengths plays an important role in the correct recovery of the time
delays. The Mock Lenses in Time (MOLET) software package is available at:
https://github.com/gvernard/molet

###### keywords:

gravitational lensing: micro – accretion, accretion discs

††pubyear: 2020††pagerange: Simulating time-varying strong
lenses–LABEL:lastpage

##  1 Introduction

Gravitational lensing refers to the deflection of light rays from a background
source due to the presence of a massive object close to the line of sight. In
the case of strong galaxy-scale lensing, a galaxy-source (the source) will
appear multiply imaged, magnified, and deformed into arc-like shapes, or even
complete rings, around a galaxy-lens (the lens). The lensed images generally
remain static in time, unless the source happens to contain some time varying
and bright astrophysical object, namely, either an Active Galactic Nucleus
(AGN) that manifests as a quasar in the optical wavelengths or an exploding
supernova. In this case, gravitational lensing as a function of time has
implications for measuring the Hubble constant, $H_{\mathrm{0}}$, a major and
very timely cosmological application.

The current paradigm of the standard cosmological model, postulating a flat
Universe consisting mainly of dark matter and dark energy (Komatsu et al.,
2011), has been challenged by precise measurements of $H_{\mathrm{0}}$:
different probes - the cosmic microwave background (CMB) radiation in the
early Universe (Planck Collaboration et al., 2020) and the local cosmic
distance ladder (e.g. Riess et al., 2016) \- give values that disagree at the
3.4$\sigma$ level, which could imply the existence of new physics (Freedman,
2017). Strong lensing provides an entirely independent physical probe to
measure $H_{\mathrm{0}}$ under different assumptions and systematics (Treu &
Marshall, 2016; Oguri, 2019), potentially resolving the current crisis in
cosmology (de Grijs et al., 2017). Currently, there is an increased interest
in understanding the systematic uncertainties in measuring $H_{\mathrm{0}}$
through lensed quasars (Millon et al., 2020c) and supernovae (Suyu et al.,
2020). However, these efforts require robust measurements of time delays and
good knowledge of the lens mass distribution: from the large scales of the
order of several kpc, down to the distribution of stellar mass objects within
the lens. The former can be affected by several factors, like an incomplete
description of the main lens mass-model (e.g. Birrer & Treu, 2019; Millon et
al., 2020b), the presence of substructure (Keeton & Moustakas, 2009; Gilman et
al., 2020b), and the mass density along the line of sight (i.e. the mass-sheet
degeneracy, Rusu et al., 2017; Wong et al., 2017; Shajib et al., 2019). In the
latter case, the presence of compact stellar-mass objects near the line of
sight induces microlensing of point-like source light components (Schmidt &
Wambsganss, 2010), like the supernovae and quasar accretion discs of interest
here111We disregard the possible, yet unlikely, scenario of a very bright star
in the source crossing a caustic or a microcaustic., which can contaminate the
observed light curves by introducing an additional extrinsic time variability
(Liao, 2019) and time delay (Tie & Kochanek, 2018).

Apart from measuring $H_{\mathrm{0}}$, lensed quasars have two additional
notable astrophysical applications: performing tomography of supermassive
black hole immediate environments and detecting the presence and measuring the
mass of dark massive substructures ($\sim 10^{7}$ M⊙ or less) around the lens
(Dalal & Kochanek, 2002; Metcalf & Zhao, 2002) or along the line of sight. In
both cases, the time varying signal of the quasar and, most importantly,
superimposed microlensing modulations play a crucial role. In the first case,
it is the onset of those rare microlensing high magnification events, most
notably caustic crossings ($\approx$1 per decade per system, Mosquera &
Kochanek, 2011), that carry the desired, but still elusive, information on the
geometry and kinematics of quasar central regions under the effect of the
extreme gravity of the supermassive black hole (Padovani et al., 2017). Such
measurements would complement reverberation mapping (e.g. Fausnaugh et al.,
2017) and direct imaging by the Event Horizon Telescope222Which, however, will
be able to directly image only 2-3 more nearby quasars, like it did for M87.
(Akiyama et al., 2019) by a bigger sample of objects that are also at
cosmological distances (having implications on quasar formation, black hole
growth, galaxy evolution in general, etc). In the second case, microlensing
can mimic the effect of (de)magnification (Schechter & Wambsganss, 2002)
caused by dark substructures on individual quasar images (MacLeod et al.,
2013; Nierenberg et al., 2017; Stacey & McKean, 2018, e.g.) and produce
similarly ‘anomalous’ flux ratios that can be used to constrain dark matter
(Gilman et al., 2020a; Harvey et al., 2020; Chan et al., 2020). However,
because there is a chromatic dependence of microlensing through the effective
source size, its effect can be mitigated and/or disentangled by
multiwavelength flux ratios. Both of these applications are very timely
because of the upcoming monitoring campaign of the Rubin Observatory Legacy
Survey of Space and Time (LSST; LSST Science Collaborations, 2009), which is
expected to discover up to a hundred microlensing events per year (Neira et
al., 2020), and the Euclid space-based survey, which will provide high
resolution multiwavelength snapshots, rich in microlensing signal (Vernardos,
2019).

Here, we present the Mock Lenses in Time (MOLET) software package that
incorporates self consistently time variability of lensed quasars and
supernovae, and we explore the effect of such variability on measuring time
delays. In order to report high precision measurements of $H_{\mathrm{0}}$, it
is crucial to apply self-consistently our current understanding of the inter-
dependency and degeneracies of all the physical processes involved in
measuring the time delays. So far, microlensing effects have been approximated
to various degrees, or even ignored, due to their prohibitively large
parameter space and associated cost of the required ray-shooting simulations
(with a few exceptions, e.g. Bonvin et al., 2018; Chen et al., 2018; Foxley-
Marrable et al., 2018; Goldstein et al., 2018). In this work, we take
advantage of the GERLUMPH microlensing simulations and software tools
(Vernardos et al., 2014; Vernardos & Fluke, 2014) to create a fully self-
consistent pixel-level forward-simulation tool, which integrates all the
related aspects of lensing and observations in a coherent way. In Section 2,
we outline the various physical systems and model components of relevance and
in Section 3 the implemented forward-simulation framework. Section 4 showcases
a novel application of MOLET: measuring time delays and their uncertainties
from lensed quasars and supernovae with different observations. We present our
conclusions in Section 5.

##  2 Description of physical processes

Here we present the different physical processes that become relevant at
different space and time scales. Emphasis is given on their inter-dependencies
and the limitations and degeneracies of related models - we do not attempt to
provide a complete review on lens theory and applications. Clarifications and
specific references are provided whenever necessary, while the reader is
referred to the reviews by Keeton (2010) and Schmidt & Wambsganss (2010), the
lecture notes of Schneider et al. (2006), or any other textbook on
gravitational lensing for more details.

###  2.1 Time-independent lensing

The lens mass model, or macromodel, refers to the mass distribution of the
main deflector - a two-dimensional projected mass distribution on the plane of
the sky (under the thin lens approximation) - and it can explicitly include
nearby companion galaxies and substructures (in the same lens plane), and/or
more subtle tidal effects as external shear. The corresponding lens potential
determines a number of key physical properties, namely the deflection angles
(the first order derivatives of the potential), the magnification matrix (the
Jacobian of the potential), and the light arrival-time surface, which are all
a function of the position on the plane of the sky. The latter has a number of
extreme points (viz. a maximum, minima, and saddle-points) that define the
locations where multiple images of a point-source would appear (according to
Fermat’s principle) and their difference in arrival-time, i.e. the time
delays. The deflection angles create the overall appearance of the lensing
system by acting on the extended source brightness profile and creating lensed
arcs, rings, and multiple extended images - in the following we refer to these
features collectively as ‘extended lensed features’ and use the term ‘multiple
images’ to refer exclusively to the point-source. Finally, the magnification
matrix describes the change in shape and total area of a surface element from
the source to the lens plane, thus determining the magnification of each
multiple image.

The magnification matrix can be expressed as a function of the total shear
vector, $\boldsymbol{\gamma}$, and the total convergence, $\kappa$, fields,
with the latter being in fact the surface mass density. Convergence is split
between two components: one that is smooth across all scales, corresponding to
dark matter and gas within the lens, and a grainy component in the form of
compact stars and stellar remnants. Compact matter on the macroscopic (kpc)
scale can be described by a smooth field as well, however, at scales
comparable to the point-source size (spanning angles of $\approx 10^{-6}$
arcsec on the sky) it consists of distinct objects that act as microlenses
(point-mass lenses). Hence, these three fields, $\boldsymbol{\gamma}$,
$\kappa$, and the smooth-to-total matter fraction, $s$ (or the compact matter
convergence, $\kappa_{\star}$), fully describe, in a statistical sense, any
subsequent (time-dependent) microlensing variability of the multiple
images333The orientation of the shear vector determines the angle of elongated
microlensing caustics, increasingly present as the magnitude of the shear
grows. This direction is important with respect to the point-source velocity
because the source may happen to move along a trough, or valley, between high
magnification regions of clustered elongated caustics. Such troughs are devoid
of any microlensing caustics and therefore do not induce any high
magnification, rather demagnify the source..

Mass distributions along the line of sight to the lens and source and at
different redshifts can further act as lenses/sources. Such a scenario,
however, requires the close alignment of three or more massive astrophysical
objects along a given direction through the Universe, which, although less
frequent than the already rare usual single deflector lensing, has been
observationally confirmed (see Collett & Smith, 2020, for the recent discovery
of an impressive third source in the J0946$+$1006 lensing system). Here we
focus on the time-variability of the lenses, leaving substructure effects
along the line of sight (Despali et al., 2018) and multiplane lensing (Gavazzi
et al., 2008; Petkova et al., 2014) as future extensions to our method.

###  2.2 Time variability

The intrinsic physical mechanisms responsible for the time variability of
quasars and the light profile evolution of supernovae are quite complex and
sometimes poorly understood. Their time-variable observational signatures,
however, are well studied and can be described by empirical models. In the
following, we always consider a given intrinsic variability light curve of the
source at its rest-frame in the absence of lensing.

The collective effect of tens to thousands of microlenses on a point-source is
described by a magnification map. Such maps describe the total
increase/decrease in flux from the sum of the areas of the created micro-
images (unresolved multiple images of the point-source, separated by $\approx
10^{-6}$ arcsec, whose surface brightness is conserved) as a function of
source plane position, and consist of intertwined networks of caustics: lines
separating high from low magnification regions, corresponding to a higher or
lower number of micro-images (always by an increment of two, see Schneider et
al., 2006). If any part of the source, which can now have finite dimensions
when compared to the caustics, overlaps with any such line, there is an
increase in magnification; the higher the fraction of total source flux
contained by the part touching the caustic, the more dramatic the
magnification (with the extreme being a true point-source, which has an
infinite magnification on the caustic, see Fluke & Webster, 1999).

The number of microlenses per unit area in the lens plane is directly related
to the local value of $\kappa_{\star}$, and is either set to the sum of equal
mass microlenses or the integral (weighted sum) of a mass function444For all
practical purposes, microlensing is always taking place within length scales
far too small for $\kappa_{\star}$, $\kappa$, and $\boldsymbol{\gamma}$, to
diverge from their values at the location of each multiple image.. The
microlens positions are practically unobservable, hence each map results from
a random realization of a microlens ensemble. As long as the map is wide
enough, it can be considered as an accurate statistical representation -
Vernardos & Fluke (2013) find a value of 25 Einstein radii555The Einstein
radius is the typical lensing scale length, or caustic size, of a single lens,
depending only on the lens mass and the angular diameter distances to the
source and the lens. For an ensemble of lenses, the mean mass is used. to be
adequate. In practice, magnification maps are produced as high resolution
pixellated versions of caustic networks, using variations of the numerical
inverse ray-shooting technique (Kayser et al., 1986).

The position of the source on the magnification map is another unobservable,
random variable. Nonetheless, the source is not static: the relative
velocities of the observer, source, and lens666Although there is a velocity
vector associated with the proper motion of each microlens, which causes the
caustic pattern to move leading to “faster” variability (e.g. Kundic &
Wambsganss, 1993), its magnitude becomes significant only in very low redshift
systems. Poindexter & Kochanek (2010) explicitly include such a component in
their calculations, while Wyithe et al. (2000) have shown that this effect can
be approximated by a bulk velocity for all the microlenses based on the
lensing galaxy’s velocity dispersion., result in a motion of the source on a
trajectory across the magnification map defined by an effective velocity
vector on the source plane (Kayser et al., 1986; Kochanek, 2004; Neira et al.,
2020). In addition to the source displacement, there can be a variation in the
source effective size for reverberating quasars and expanding supernovae.
Hence, different parts of the source can undergo a varying degree of
magnification, which can be up to a few magnitudes over a period of months to
years, and produce a time-dependent signal as they cross or expand into
neighbouring regions, filled or devoid of caustics. This signal is
superimposed on the intrinsic source signal and independent between the
multiple images (extrinsic). The different possible such sources that can be
found in lensing systems are:

  1. (i)

Moving disc. This is the most commonly used quasar microlensing variability
model (e.g. Morgan et al., 2012; Morgan et al., 2018; Cornachione et al.,
2020). It consists of a fixed accretion disc brightness profile that is
crossing the magnification map at a given velocity. The intrinsic variability
is assumed to take place simultaneously across the entire disc, scaling its
brightness up or down.

  2. (ii)

Moving and reverberating disc. The quasar instrinsic variability and its light
profile are, in fact, connected through physical mechanisms and geometry. For
example, when a temperature variation is produced at the center of the
accretion disc, it propagates outwards with a finite speed, changing the
brightness of the regions it crosses with a time lag with respect to the
center. This is the so-called ‘lamp-post’ model (Cackett et al., 2007), which,
in combination with microlensing, has been found to introduce an additional
time delay (Tie & Kochanek, 2018). The brightness profile will be time
dependent but not with a simple overall scaling, as in the previous case,
rather with a time-lagged brightening of different regions of the accretion
disc, which are also differently microlensed. The net effect is a time varying
brightness profile that is crossing the magnification map.

  3. (iii)

Moving disc with un-microlensed component. This model ignores the response of
the innermost regions of the accretion disc to central variations but includes
a delayed emission coming from the much larger Broad Line emission Region
(BLR, e.g. Sluse & Tewes, 2014). The latter one is much less affected by
microlensing but the delayed intrinsic variability shapes the eventual signal
by introducing more small-scale structure in the light curves (Paic et al.
2021, in prep.).

  4. (iv)

Expanding supernova. This is the most commonly used supernova microlensing
variability model (e.g. Foxley-Marrable et al., 2018; Goldstein et al., 2018).
The main assumption is that the supernova expansion velocity ($\sim 10^{4}$
km/s) is much higher than the effective source plane velocity ($\sim 10^{2}$
km/s), therefore, the source center-of-light can be considered static. In this
case, microlensing variability arises due to different regions of the
supernova light profile overlapping with different magnification regions of
the map as it expands. The expansion velocity and light profile geometry can
vary as a function of time.

  5. (v)

Expanding and moving supernova. In some cases, the above assumption might not
be true: in systems having a maximum observer velocity (CMB dipole transverse
velocity, see figure 1 of Neira et al., 2020), a low redshift
($z_{\mathrm{L}}\lesssim 0.15$) and a high-mass lens, the mean effective
source plane velocity can be of the order of 10 per cent of the supernova
expansion velocity. This can be considerably higher if the lens and source
peculiar velocities happen to align parallel to the observer’s velocity as
well. Although less probable, such a system is plausible and requires a
simultaneous expansion and translation of the source across a magnification
map to be taken into account.

In all the above scenaria, the most crucial structural parameter of the
microlensed point-source is its size, expressed by the half-light radius,
$r_{\mathrm{1/2}}$. It has been shown that the detailed geometry of the source
brightness profile does not play any significant role in bulk microlensing
properties (Mortonson et al., 2005; Vernardos & Tsagkatakis, 2019, although
the opposite holds for short-lived high magnification events). If the source
is smaller than a single magnification map pixel then values from the map can
be used directly (subject to shot noise from the ray-shooting simulations).
Such a situation is not optimal for expanding or slow-moving sources though,
because if there is a caustic present, its fine structure is washed-out over
the entire pixel and hence a higher resolution map is required to correctly
capture the induced variability. In the case of the source being larger than a
map pixel, then a convolution of the map with the light profile is required
before extracting any magnification values.

##  3 Simulation technique and implementation

Each macromodel determines the values of $\kappa,\gamma,s$ at the locations of
the multiple images, which, in turn, define the statistical properties of a
magnification map. Computing a single map requires ray-tracing simulations
with deflections due to hundreds to tens of thousands of microlenses - a
computationally demanding task (Bate & Fluke, 2012). Although there are
several forward-simulation techniques to generate static lenses (see the
review by Plazas, 2020), none of these include the simultaneous self-
consistent production of microlensing observables, viz. flux ratios and light
curves, precisely due to the prohibitive associated computational cost.

The novelty of the approach presented here is the self-consistent inclusion of
microlensing through pre-computed magnification maps and software
acceleration. The GERLUMPH777http://gerlumph.swin.edu.au parameter survey
provides a systematic coverage of the entire $\kappa,\gamma,s$ parameter space
of interest with high-quality readily available magnification maps (Vernardos
et al., 2014; Vernardos & Fluke, 2014). These maps have a sufficiently large
width (25 Einstein radii), high resolution (400 pixels per Einstein radius),
accurate statistics ($>10^{9}$ light rays shot per map), and are scalable with
the mean microlens mass888GERLUMPH maps have been created using a fixed mass
for the microlenses. Although distributing the baryonic mass to an ensemble of
stellar mass objects using a given Initial Mass Function prescription is
physically justified, it has been shown to have a minimal effect (Wambsganss
et al., 1990; Lewis & Irwin, 1995; Wyithe & Turner, 2001).. Software
acceleration by Graphics Processing Units, which enabled GERLUMPH, is also
critical in performing convolutions between maps and finite-sized point-
sources, leading to speed-ups of $\times 10$ (Vernardos et al., 2015). This is
particularly relevant in the case of expanding sources, which require a
convolution of a differently sized profile at each time step (as opposed to a
fixed source moving across a map that requires only one convolution). The
software library GERLUMPHPP999https://github.com/gvernard/gerlumphpp, provides
an easy interface to harness GPU acceleration (but still compatible with
CPUs), and has already been used in a number of microlensing applications
(Foxley-Marrable et al., 2018; Neira et al., 2020).

Our simulation technique does not depend on any specific assumption on lens or
source populations and instrument specifications. Theoretical or empirical
relations between different physical models/properties can be provided
independently, e.g. one could relate the mass with the light profile of a
fiducial elliptical lens galaxy using the fundamental plane, make further
assumptions on the mass-to-light ratio or the total and dark matter mass
profile to determine $\kappa_{\star}$, assume different point-source light
profile, geometry, effective velocity, intrinsic variability, etc. A basic
choice of specific physical models is provided, e.g. mass, light profile, etc,
however, our method can easily be extended by additional user-provided
functions and data. A flowchart describing the simulations’ pipeline is shown
in Fig. 1; below we summarize its various components:

Figure 1: Flowchart of the MOLET software package.

  1. (i)

Cosmology. We calculate the angular diameter distances to the lens, source,
and between them, based on a given cosmological model. The value of
$H_{\mathrm{0}}$ needs to be provided here.

  2. (ii)

Lens. The following properties of the deflector are required:

     * •

Mass model: it can consist of a single, or several analytical models (see
Keeton, 2001, for a catalogue), or pre-computed lens potential values,
$\psi(x,y)$, which, for example, can be based on numerical N-body
hydrodynamical simulations (e.g. Mukherjee et al., 2018). An external shear
vector can be included as well. We currently provide the Singular Isothermal
Ellipsoid (SIE; Kormann et al., 1994) and softened power-law (Barkana, 1998)
models.

     * •

Light profile: it can have one or more analytical components [here we provide
Sérsic Sersic (1963) and Gaussian profiles], or be an actual image of a
galaxy. In the latter case, it should have a resolution higher than the
simulated instrument to avoid introducing artifacts during the ray shooting
and binning to the final pixel size.

     * •

Compact mass model: in the case of a point-source, such a model is required to
calculate $\kappa_{\star}$ at the location of its multiple images. It can be
an analytical profile (we provide Sérsic and Gaussian), or a mass-to-light
ratio value to convert the light profile directly into compact stellar mass.

  3. (iii)

Extended source. Any light profile, similarly to the lens light profile, can
be used for the source. Its location with respect to the macromodel caustics
determines the final appearance of the lensing system, i.e. the number of
multiple extended images and if they are merged into arcs or a ring.

  4. (iv)

Point source. The location of the point-source with respect to the macromodel
caustics determines the number of multiple images and their configuration,
i.e. cross, cusp, or fold. The half-light radius, $r_{\mathrm{1/2}}$, is the
main parameter of the point-source’s actual light profile, while its detailed
shape becomes relevant only with respect to high magnification events,
observed with a very fine temporal resolution ($\sim$day). A mean mass value
of the compact objects acting as microlenses is required to set the Einstein
radius and consequently the microlensing magnification map length scale. Here
we use GERLUMPH maps, but any magnification map could be used, e.g. generated
according to a microlens mass function, as long as it is converted to the
compatible format first. Regarding the evolution of the point-source with
respect to the microlensing magnification map caustics (translational movement
or expansion), we can have the following cases:

     * •

Moving disc. We use the effective source plane velocity model from Neira et
al. (2020), which consists of: a fixed velocity vector for the observer based
on the CMB velocity dipole that is a function of the system’s sky coordinates,
a randomly directed vector of the combined peculiar velocities of the lens and
the source with magnitude drawn from a corresponding Gaussian distribution
(eq. 5 of Neira et al., 2020), and a random velocity for the microlenses equal
in magnitude to the velocity dispersion in the lens.

     * •

Moving and reverberating disc. We use the ‘lamp-post’ quasar variability model
(Cackett et al., 2007) to relate a central temperature variation to a
brightness profile variation that propagates outwards through the accretion
disc. For an analytic light profile, one needs to compute its first derivative
and relate it to the central variability via a time lag (see Tie & Kochanek,
2018). Hence, at each time step a different light profile is produced, which
is then microlensed as it moves across a magnification map. Here we provide a
lamp-post variability model for the thin disc (Shakura & Sunyaev, 1973).

     * •

Moving disc with un-microlensed component. We use an additive non-microlensed
intrinsic flux for each multiple image (see eq. 1 here and 4 of Sluse & Tewes,
2014).

     * •

Expanding supernova. We use a fixed expansion velocity and shape for the
supernova light profile, but our method can be easily extended to include time
varying velocities and profile shapes (also for the moving model below).

     * •

Expanding and moving supernova. We use the same effective velocity model as a
moving disc and combine it with a supernova expanding light profile.

  5. (v)

Instrument. The following characteristics of existing or hypothetical
instruments need to be specified:

     * •

Resolution: the pixel size of the detector’s camera.

     * •

Field of view: the dimensions of the simulated patch of sky.

     * •

PSF: the point-spread function.

     * •

Observing strategy: a sequence of dates when the simulated system was
observed. This can be related to an assumed RA and DEC of the simulated
system, combined with a specific facility’s observing schedule (e.g. see call
for white papers on LSST cadence optimization101010ls.st/doc-28382).

     * •

Noise: here we provide a simple uniform Gaussian noise model calculated from a
given signal-to-noise ratio. Existing simulators of varying complexity can
also be used (e.g Lage, 2019, for LSST’s camera).

To produce the extended lensed images, a grid of pixels is assumed on the lens
plane and then deflected backwards on to the source plane using the deflection
angles computed from the lens mass model. Once the location of the deflected
rays on the source plane is determined, the lensed brightness values are
computed by evaluating an analytical source profile or interpolating between
the associated adjacent pixels of a numerical one. At the same time, the
critical and caustic lines of the mass model are computed, which are used in
validating the multiple image configuration.

To find the multiple images of a point source at a given location, we cover
the lens plane completely with a coarse grid of triangles, which we then
deflect on to the source plane using the lens mass model. Finding which
triangles contain the point source gives an approximate location of its
multiple images. We then redefine a finer grid of triangles around these
approximate locations and repeat the process until the size of the grid
becomes smaller than 1/100 of the final simulated instrument’s resolution. The
final output is the positions of the multiple source images and an associated
uncertainty, which is below 1 per cent.

The intrinsically variable flux of the point-source needs to be provided as an
input (through simulations, e.g. see Sartori et al., 2019; Kessler et al.,
2009, for quasars and supernovae respectively), and has to span at least
$t_{\mathrm{obs}}+\Delta t_{\mathrm{max}}$, where $t_{\mathrm{obs}}$ is the
duration of the observing period and $\Delta t_{\mathrm{max}}$ the maximum
time delay. External microlensing variability is calculated separately and
consists of a varying multiplicative magnification factor for the time period
$t_{\mathrm{obs}}$. The final point-source signal for each multiple image q is
then:

$F_{\mathrm{q}}(t)=M_{\mathrm{q}}\mu_{\mathrm{q}}(t)F_{\mathrm{in}}(t-\Delta
t_{\mathrm{q}})+M_{\mathrm{q}}F_{\mathrm{rev}}(t-\Delta t_{\mathrm{q}}),$ (1)

where $M$ is the constant macro-magnification, $\mu(t)$ the microlensing
magnification, $\Delta t$ the time delay, and $F_{\mathrm{in}}$ the intrinsic
source flux. The additive term is only valid in the presence of an un-
microlensed flux component (more pertinent to quasar sources),
$F_{\mathrm{rev}}$, that can be modelled by a response of the BLR to the
central variability:

$F_{\mathrm{rev}}(t)=f_{\mathrm{BLR}}\left[\Psi(t,\tau)\circledast
F_{\mathrm{in}}(t)\right],$ (2)

where $\Psi$ is a transfer function, $\tau$ a time lag usually corresponding
to the light travel time from the central regions to the BLR, and
$f_{\mathrm{BLR}}$ the fraction of this flux that is eventually re-emitted in
the same wavelength.

The three simulated light components, viz. extended lensed features, point-
source, and lens light, are created first in a super-resolved grid, having
$\times 10$ smaller pixels than the final simulated instrument’s resolution.
Convolution with the PSF takes place in this grid, over-sampling the PSF to
match the smaller pixel size if necessary. For the point-source light, we
simply scale an image of the PSF (convolution with a delta function) according
to the time-varying signal from eq. (1) for each multiple image. The images
are then binned to the instrument’s resolution and noise is added111111We note
that the convolution and binning operations are not necessarily permutable and
reversing their order may overly smooth out the lensed images or, conversely,
introduce spurious pixel edge effects..

The most critical observational parameter with respect to time-variable lenses
is organizing their monitoring, i.e. the frequency and distribution of the
observations over time. Neira et al. (2020) have shown the effect of different
observing (LSST-specific) strategies in relation with detecting high
magnification events in lensed quasars and Goldstein et al. (2018) with
respect to lensed supernovae from ZTF and LSST. While factors like
instrumental effects, noise, and artifacts, do affect the resulting simulated
images as every other astronomical observation, modelling them correctly,
which depends on various physical mechanisms affecting incoming photons from
the lensed system reaching Earth’s atmoshpere (or Earth’s orbit) and finally
being recorded by an instrument’s detector, is of secondary importance for the
purposes examined here. In any case, our method can be easily extended by
additional instrumental modules, consisting of a model for the PSF, noise, and
artifacts.

The essential final products of the above procedure are the continuously
(daily) sampled light curves of each multiple image. Based on this, one can
study the effect of different observing strategies for each science
application. In addition, our method produces pixel-level simulations of
cutout images of lens systems at different time intervals. Such images are
necessary to calibrate methods that analyze the actual observational pixel
data, from deblending the lensed images to lens mass modelling.

###  3.1 The MOLET software package

The independent simulation components, viz. extended lensed features, point-
source image locations, microlensing, etc, are reflected in the modular design
of the MOLET software package outlined in Fig. 1. Each procedural component is
a C++ program written using the GERLUMPHPP and
VKLLIB121212https://github.com/gvernard/vkl_lib libraries for microlensing and
strong lensing respectively. Minimal shell scripts are used to manage some
operating system specific tasks (e.g. checking and creating files), images are
written in .fits format, and data in .json files, easily shared between the
different modules. The code is open source and each module can be replaced by
any other program that performs the same scientific computations, as long as
the output data format is adhered to. All models presented above are
implemented in this first version of the software, except a moving and
reverberating disc and an expanding and moving supernova (both cases of a
moving and variable in size source). These are planned to be included in a
future version, together with a self-consistent way to get a reverberated flux
component from the BLR (eq. 2) and possibly a basic intrinsic variability
model (like the damped random walk, e.g. MacLeod et al., 2010). All third-
party software required to run MOLET is provided in software containers,
easily deployable in any operating system or supercomputer. The source code
and detailed instructions can be found at the project’s
homepage131313https://github.com/gvernard/molet.

Figure 2: Mock image of RXJ 1131-1231 and associated magnification maps at the
locations of the quasar multiple images. The maps are consisted with the mass
model, i.e. have the corresponding $\kappa,\gamma,s$ (see Table 1), as well as
being aligned with the local total shear vector. A trial trajectory for the
quasar point source is shown with an effective velocity drawn from the
distribution given in fig. 3 of Neira et al. (2020). Its length has been
multiplied by a factor of 5 to better visualize that the trajectories in each
map are parallel - constituting, in fact, just a single trajectory on the
source plane where the different map realizations are projected. The
corresponding microlensing light curves are shown in Fig. 3. Figure 3: Top:
assumed intrinsic variability of the quasar point source of RXJ 1131-1231. The
blue-shaded area shows the $\approx$10-year time-window of the assumed
monitoring campaign and the orange-shaded one the time window into the past
provided by the maximum time delay of the system (corresponding to the
unobservable maximum image that has a delay of $\approx 150$ days). Bottom:
microlensing magnification for each image, and corresponding zoomed-in region
of the magnification maps with the source trajectory indicated by the arrow
(same as the trajectories in Fig. 2). The horizontal dashed lines correspond
to the macro-magnification of each image, listed in Table 1. The shaded areas
indicate the Poisson error of the magnification in the GERLUMPH map pixels,
resulting from the numerical inverse ray-shooting technique (Kayser et al.,
1986) that was used to generate them.

The novelty of MOLET is employing microlensing simulations in real time and
self-consistently. This is enabled by the $>75,000$ pre-computed GERLUMPH
magnification maps that cover extensively the $\kappa,\gamma,s$ parameter
space of relevance. Running MOLET locally requires to have copies of the
required magnification maps; when such copies are not present a link to
download them is provided automatically.

Generating microlensing light curves is the most computationally intensive
task. In particular, if several convolutions between maps and source profiles
are required, as is the case for expanding or time-varying source light
profiles, the computational cost on the CPU can become prohibitive. However,
GPUs can accelerate this task by a factor of 20 (Vernardos & Fluke, 2014) and
MOLET can use the GPU-compiled version of GERLUMPHPP to achieve this.

##  4 Results

To demonstrate the capabilities of MOLET, we present simulated monitoring of
the quasar RXJ 1131-1231(Sluse et al., 2006) over a period of $\approx$10
years. This particular object is interesting because it is known to undergo
microlensing, it is being currently monitored by the COSMOGRAIL program
(Millon et al., 2020b), and it is extensively used to derive $H_{\mathrm{0}}$
from time delay studies (e.g. Suyu et al., 2013; Wong et al., 2020). Here, we
use the results of Chen et al. (2016) for the parameters of the lensing
potential, the lens light profile, the quasar point source location, the
reconstructed light profile of the extended source, and the reconstructed PSF
(see their tab. 1, fig. 9, and fig. 10 respectively). Consequently, the
simulated data are for the Keck Adaptive Optics K′ band, which has a pixel-
scale of 0.04 arcsec (see fig. 1 of Chen et al., 2016). A mock image for the
system is shown in Fig. 2.

Figure 4: Example of continuous (lines) and sampled (points) mock light curves
and corresponding mock Keck K-band cutouts of RXJ 1131-1231. The sampling of
the light curves is the same as the one from the COSMOGRAIL observations. The
cutouts are shown for illustration purposes only and have not been used to
extract the light curves. Here we have scaled the lens and source brightness
and skewed the color coding of the pixels to better show the variability of
the multiple images, e.g. a clear brightening of image A due to microlensing
is observed around 2500 days.

For the quasar point-source, we use the coordinates of its recovered position
on the source plane141414G.C.F. Chen, private communication. to find the
corresponding multiple images created by the lens potential. MOLET currently
requires the intrinsic variability of the quasar to be provided as input,
which we get as a realization of a Damped Random Walk (DRW, see MacLeod et
al., 2010, for a justification of this choice) with the characteristic
timescale and variance151515We use the values of 80 days and 13 for these two
parameters - D. Sluse, private communication. parameters derived from the
observed light curve of RXJ 1131-1231(Sluse & Tewes, 2014). A realization of
an intrinsic light curve is shown in the top panel of Fig. 3. Finally, there
is a secondary point source visible in the Keck AO data, attributed to a star
forming region of the quasar host galaxy, which, however, we omit, as
monitoring of such faint point sources is unrealistic with current
instruments. Nevertheless, dealing self-consistently with (the rare case of)
multiple point sources and their associated microlensing can be an easily
implemented extension of MOLET.

The local values of $\kappa$ and $\gamma$ for the multiple quasar images
result directly from the lens potential. For $s$, one needs to provide a
distribution for the projected compact (stellar) mass of the lens. Usually,
the mass is assumed to follow the lensing galaxy’s light distribution with
some mass-to-light ratio. Proceeding in this way would require the calibration
of the fluxes reported in Chen et al. (2016) to standard filters. Instead, and
because our goal is not to study the actual distribution of stars within the
lens, we assume a compact mass profile identical to the main Sérsic light
component of the lens but scale its amplitude using a reasonable value for the
mass density in stars, in this case $\rho_{\rm eff}=1.0$ kg/m2. The
magnification maps at the locations of the multiple images, shown in Fig. 2,
are directly retrieved from GERLUMPH, with their $\kappa,\gamma$, and $s$
values matching the ones resulting from the lens macromodel as closely as
possible (see Table 1). These square maps are 25 (10000) Einstein radii
(pixels) wide. We set a uniform mass of 1 M⊙ for the microlenses.

Table 1: Properties of the multiple images located at the coordinates $x$ and $y$ on the image plane (in arcsec). The shear angle, $\phi_{\gamma}$, is given in degrees east-of-north and the time delay, $\Delta t$, is given in days with respect to image B. $\mu$ is the macro-magnification. The last three rows show the closest $\kappa,\gamma,s$ values from GERLUMPH. | A | B | C | D
---|---|---|---|---
| Macromodel
x | 2.27 | 2.23 | 2.76 | 5.38
y | 3.38 | 4.6 | 2.45 | 4.31
$\kappa$ | 0.49 | 0.44 | 0.46 | 0.93
$\gamma$ | 0.55 | 0.49 | 0.46 | 0.92
$\phi_{\gamma}$ | -164.44 | -11.75 | -139.72 | -166.47
$s$ | 0.72 | 0.69 | 0.71 | 0.76
$\mu$ | -23.07 | 13.19 | 13.36 | -1.19
$\Delta t$ | 0.8 | 0 | 0.2 | 94.6
| GERLUMPH
$\kappa$ | 0.49 | 0.44 | 0.46 | 0.95
$\gamma$ | 0.55 | 0.49 | 0.46 | 0.91
$s$ | 0.7 | 0.7 | 0.7 | 0.8

For the accretion disc and its effective velocity across the source plane, we
adopt the same models as Neira et al. (2020), whose parameter values are given
in their table 1. The half-light radius of the accretion disc as a function of
wavelength is given by a simple power-law model with a slope equal to $4/3$,
i.e. the same as the thin-disc model; for an observing wavelength of 2190 nm
(the midpoint of the standard photometric K filter) the resulting half-light
radius of the disc is $7.47\times 10^{14}$cm, or 0.016 times the value of the
microlenses’ Einstein radius on the source plane. The light profile of the
disc is assumed to have the shape of a normal distribution viewed face-on. The
fiducial effective velocity that produces the trajectories on the
magnification maps shown in Fig. 2 has a magnitude of 294 km/s, which is very
close to the expected mean, and a $-126$ deg direction east-of-north (see fig.
3 of Neira et al., 2020). Except image D, for which it is rare to lie close to
a caustic (see the corresponding magnification map in Fig. 2), all the other
images have light curves that cross dense caustic regions on the maps,
producing the rich in features light curves shown in Fig. 3.

For the time sequence of the observations we choose 600 observing dates from
the COSMOGRAIL light curve of RXJ 1131-1231, which have intervals of roughly
once a week and season gaps, with an arbitrary starting point (only the
relative time intervals are relevant). Photometric errors can be added
directly to the light curves according to instrument specifications, depth,
etc, as we will perform below, or calculated from the simulated cutouts
through deblending and light curve extraction techniques. An example of
resulting simulated light curves taking into account all the variability
effects and the macroscopic time delays is shown in Fig. 4.

We can also simulate multiwavelength light curves for RXJ 1131-1231, like the
ones expected in the 6 LSST filters. After choosing a specific observing
strategy (see Neira et al., 2020, who examined several observational
strategies for this system), the only input to determine the dates of each
visit at the patch of sky containing the lens are its coordinates. The
accretion disc model is the same as the one used above, leading to similarly
small sizes: in the u band the source is smaller than a magnification map
pixel and in the y band the source is still $\approx 3$ times smaller than the
one assumed in Figs. 3 and 4. For simplicity, we keep the same intrinsic light
curve across the 6 filters, although smoother (sharper) intrinsic variability
is expected for longer (shorter) wavelengths. The differences in disc size for
our model are of the order of a light day, too small to produce significant
variations in the intrinsic light curves anyway. An example of resulting
multi-wavelength light curves is shown in Fig. 5, where we do observe the
expected microlensing behaviour: the bluer (redder) the wavelength the sharper
(smoother) the features.

Figure 5: Example of LSST-observed light curves corresponding to those shown
in Fig. 4. Here, the smaller accretion disc gives rise to sharper high
magnification events. The colossus_2664 LSST observing strategy was used in
this particular example to determine the time intervals between observations
in each filter (see Neira et al., 2020, for more examples). The numbers in the
legend indicate the number of data points in each band. The error bars of
0.005 mag are smaller than the size of the data points. Figure 6: Example of
continuous (lines) and sampled (points) mock light curves and selected cutouts
in the LSST i band, for a hypothetical system having exactly the same mass and
light configuration as RXJ 1131-1231, but having replaced the quasar source
with an exploding supernova. The observing dates result from the the
colossus_2664 LSST observing strategy. Here we have scaled the brightness of
the supernova host galaxy and skewed the color coding of the pixels to better
show the variability of the multiple images. Image D becomes noticeable around
JD61165. The first cutout shows the lens before the supernova explosion
(without any point-source light).

Finally, we can also use MOLET to produce mock lensed supernovae light curves.
As an example, we keep the same lens configuration but change the nature of
the source from a variable quasar to an exploding supernova. We simulate the
observations with the same LSST observing strategy as in the previous case. We
assume a brighter than typical Type Ia supernova that peaks at 16 mag in the i
band, having the geometry of a uniform face-on disc that expands with a
constant velocity of $10^{5}$ km/s and with the same size in each wavelength.
This corresponds to a size range from $<10^{15}$ cm ($<0.02$ Einstein radii)
on day one, about 8 map pixels, to $1.6\times 10^{17}$ cm (4 Einstein radii)
after 186 days, when we assume that the phenomenon stops161616Long term
microlensing effects in supernovae have been very little explored so far.. We
generate the intrinsic supernova LSST light curves from the software package
SNCosmo 171717https://sncosmo.readthedocs.io/en/stable/, and assume an
arbitrary date for the onset of the explosion. Fig. 6 shows an example of such
simulated light curves in the LSST i band, for a single microlensing
realization.

###  4.1 Measuring time delays

Measuring time delays in each of the above observational setups has its own
challenges and requirements; here, we assess the effect of different intrinsic
and microlensing variability realizations. We use resolved light curves, not
cutouts, therefore the only relevant instrumental parameters are: for quasars,
the wavelength of the observations that sets the corresponding accretion disc
size, and for supernovae, the cadence of the observations. We keep the same
lens potential, point source position, and time delay distance (1935 Mpc, Chen
et al., 2016), therefore the same time delays, listed in Table 1. For a flat
Lambda Cold Dark Matter cosmological model with $\Omega_{\mathrm{m}}=0.3$, the
time delay distance corresponds to a value of 84 km/s/Mpc for
$H_{\mathrm{0}}$. Our focus here is to assess the uncertainty and bias of
techniques that measure the time delays, not $H_{\mathrm{0}}$, hence such a
high value is not a problem (although not too surprising either for this
specific lens, e.g. Wong et al., 2020). The measured COSMOGRAIL time delay
between the leading image B and image D - the most meaningful to examine since
the other delays are too small, between 1 and 2 days - is $93.7^{+2}_{-2}$
days (Millon et al., 2020c), practically the same as the $94.6$ days we get
under the assumptions of our model.

First, we consider RXJ 1131-1231 observed in the r band of the ECAM
($600-725$nm, see Millon et al., 2020b, for details). This results in an
accretion disc $\approx 6$ times smaller than the one in the Keck K band used
in the examples shown in Figs. 3 and 4. Because of this smaller size, we
expect more pronounced short-range microlensing features to be present. The
cadence of the observations is the same as in the example shown in Fig. 4 and
photometric errors are added by hand; we use the relative errors from the
publicly available COSMOGRAIL light curve for RXJ 1131-1231 and assign them to
the simulated data points of each image in the same order in time but with an
arbitrary starting point. For the intrinsic variability, we assume the same
DRW model and parameters as before and create 10 realizations with a mean
apparent magnitude of 18, roughly matching the Euler telescope images in the r
band. For microlensing, we also create 10 realizations by sampling the
effective velocity magnitude and direction. We combine these using just the
first term in eq. (1) to create a set of 100 light curves to measure the time
delays from. A second set of another 100 simulations is produced, this time
including an unmicrolensed intrinsic flux component from the BLR (the second
term in eq. 1). This flux is calculated from the same intrinsic light curve
each time, assuming a top-hat transfer function, a time lag of 17 days from
the size of the BLR (Mosquera & Kochanek, 2011), and an eventual 20 per cent
of the reverberated unmicrolensed flux making it into our r band filter (see
Paic, 2021, in prep., for a similar application to the microlensed quasar Q
J0158-4325).

To measure the delays, we use the
PyCS181818https://cosmograil.gitlab.io/PyCS3/ software (Tewes et al., 2013) to
fit the light curves assuming spline polynomials to describe the intrinsic and
microlensing signal and get a point estimate of the delays for each mock (see
Millon et al., 2020a, for a description of the different fitting methods
available). In order to avoid convergence of the fit to some local minimum, we
always initialize the fitting algorithm with the true delays rounded to the
nearest day. Finally, we assume a small number of knots in the splines to
account for short timescale and high amplitude variability - as expected due
to the small accretion disc and reverberated flux. In Fig. 7 we show the
measured time delays for each image with respect to the leading image B. The
measurements are quite accurate, nicely clustering around the true values even
for the short time delay images A and C, but not precise, having a spread of
up to $\pm 5$ days. The mean measurements and 1$\sigma$ confidence intervals
are listed in Table 2, where we can see that the presence of an unmicrolensed
reverberated component decreases the precision for image D by roughly 1.5
days.

We replace the quasar with a supernova and create a set of 100 LSST light
curves with different microlensing (we keep the same intrinsic variability).
To facilitate the convergence of the fitting algorithm, we have increased the
LSST photometric errors to the fixed value of 0.005 mag. We use the
SNTD191919https://sntd.readthedocs.io/en/latest/ software package to measure
the time delays in 4 (r,i,z,y) of the 6 LSST bands, because u and g end up not
having enough data points in the given time period. A simultaneous fit across
these bands is performed using the SALT-2 (Guy et al., 2007) supernova
variability model (here we ignore the fitted parameters of the intrinsic light
curve and only keep the time delays). The measured time delays for a sample
with 100 different microlensing realizations are shown in Fig. 7 and the mean
and $1\sigma$ confidence intervals listed in Table 2. For image D, the spread
of these measurements is very similar to the one for a microlensed quasar with
an unmicrolensed flux component. However, it is important to note that LSST
will not have an adequate cadence to constrain time delays with supernovae,
hence the uncertainty in the measurements shown in Fig. 7 is expected to
reduce to 1-2 days when using dedicated follow-up (Huber et al., 2019).

Figure 7: Probability density of measured time delays for different samples of simulated light curves. The two samples of 100 simulated quasar light curves each, one with intrinsic and microlensing variability only (orange) and one including an unmicrolensed reverberated flux component (turqoise), give a very similar distribution of measurements, albeit slightly broader for the case with reverberation. A sample of 100 simulated multiband supernova light curves, including microlensing, has a somewhat wider spread than the quasar measurements for image D (the other images having too small a delay of 1-2 days). The means and $1\sigma$ confidence intervals of the distributions are listed in Table 2, while the true values of the delays are listed in Table 1 \- indicated here by the vertical lines. Figure 8: Probability density of time delays measured from a sample of 100 simulated light curves in each LSST filter. The number of data points in each light curve in each filter is shown in the legend of Fig. 5. The means and $1\sigma$ confidence intervals of the distributions are listed in Table 2, while the true values of the delays are listed in Table 1 \- indicated here by the vertical lines. Table 2: Measured mean time delays (in days) and $1\sigma$ confidence intervals from different datasets of 100 light curve realizations each. mock dataset | $\Delta t_{\mathrm{A}}$ | $\Delta t_{\mathrm{C}}$ | $\Delta t_{\mathrm{D}}$
---|---|---|---
Quasar+micro | $0.8_{-0.3}^{+0.4}$ | $0.1_{-0.4}^{+0.3}$ | $94.9_{-1.0}^{+1.5}$
Quasar+micro+rev. | $-0.1_{-1.2}^{+0.9}$ | $-0.2_{-1.4}^{+1.1}$ | $94.8_{-1.7}^{+2.2}$
Supernova | $1.4_{-2.7}^{+2.5}$ | $-0.2_{-1.8}^{+2.6}$ | $94.3_{-2.3}^{+2.1}$
LSST u | $-1.1_{-7.0}^{+3.5}$ | $-0.0_{-2.9}^{+1.7}$ | $96.5_{-3.7}^{+4.7}$
LSST g | $1.7_{-5.4}^{+1.8}$ | $0.5_{-2.0}^{+0.9}$ | $94.5_{-3.6}^{+2.9}$
LSST r | $0.7_{-0.7}^{+0.6}$ | $0.1_{-0.4}^{+0.6}$ | $99.2_{-2.8}^{+2.2}$
LSST i | $0.3_{-0.3}^{+0.2}$ | $0.0_{-0.1}^{+0.1}$ | $95.6_{-0.9}^{+0.8}$
LSST z | $0.7_{-0.4}^{+0.3}$ | $0.2_{-0.2}^{+0.2}$ | $98.0_{-1.7}^{+4.4}$
LSST y | $0.5_{-0.2}^{+0.4}$ | $0.1_{-0.1}^{+0.3}$ | $96.8_{-0.8}^{+1.0}$

Finally, we repeat the experiment for the lensed quasar RXJ 1131-1231 using
100 LSST light curve realizations and measure the time delays separately in
each filter with the same setup for PyCS. We show our results in Fig. 8, where
we observe that the accuracy and precision vary considerably per filter. For
images A and C, filters u and g have the largest spread, which could be
explained by sparser observations and increased short-scale variability due to
the smaller accretion disc, while the delays are well recovered in the other
filters. For the longest and most interesting delay of image D, the i and y
bands have the smallest uncertainty, albeit with a slight, 1-2 day systematic
offset towards longer delays. This result is harder to explain because it does
not correlate neither with the size of the accretion disc, nor with the number
of observations. The mean measurements and 1$\sigma$ intervals are shown in
Table 2.

##  5 Discussion and future prospects

We have presented a self-consistent and versatile forward modelling software
package that can produce time series and pixel-level simulations of time-
varying strongly lensed systems. In summary, the main input required is a mass
model for the lensing galaxy, including some assumption on its compact
(stellar) matter distribution, and a variability model for the point-source
(intrinsic and microlensing). This, together with the observing wavelength and
cadence, is sufficient to generate time series. For simulated images, in
addition to the above, the light profile of the lensing galaxy, the extended
light from the source, and the observing instruments’ resolution, PSF, and
noise properties are also required. The resulting images and time series can
be used independently of each other - in the latter case, error bars will need
to be provided separately by the user. Finally, we note that images for static
lenses can be generated as well by switching off the time-dependent part,
however, the originality and wider applicability of the method stem from the
latter.

There are no assumptions on specific populations or scaling relations inherent
in the method, meaning that a very broad combination of various models and
parameters is possible. Although this allows for unrealistic combinations, and
the responsibility of avoiding this is left to the user, it also allows for
maximum flexibility. For example, there is no constraint (or safeguard) on
matching an extended source’s brightness to its redshift, or the intrinsic
mean brightness of a quasar to its accretion disc size. Specifically, in this
first version of the software all the intrinsic variability light curves and
their relations (e.g. via transfer functions in case of reverberation, or
spectral features of supernovae across wavelengths) need to be computed
externally. However, all the physical mechanisms that are connected within the
context of strong lensing and microlensing are modelled self-consistently,
e.g. the underlying cosmological model, the convergence and shear fields
required for microlensing, conservation of extended source flux, etc.

Here we used a detailed mass model of RXJ 1131-1231 from the literature (Chen
et al., 2016) as a workbench to measure time delays with different methods
and/or data. Although we restricted the extent of the parameter space explored
- this is a demonstrative study anyway - useful conclusions can still be
drawn. We find evidence that more regular monitoring of lensed quasars, like
in the COSMOGRAIL program ($\approx$60 observations per year), performs better
in measuring time delays than LSST (approx. between 6 to 20 observations per
year depending on the filter), as shown in Table 2. Cadence, however, is not
the main limiting factor of the LSST measurements (see Fig. 8) and the role of
other factors is yet to be explored. Different components of the variable
signal coming from the quasar (e.g. different relative power of high and low
frequencies) seem to play a role in the precision of the measured time delays,
as indicated by the broader measurements recovered here while assigning an
additional 20 per cent of the intrinsic flux to an unmicrolensed component
reverberated by the BLR. Replacing the quasar with a supernova, leads to
similar measurements, but in this case this is driven by the low number of
data points ($\approx 15$ per filter, see Fig. 6), which should not be the
case if an early alert for a supernova transient is provided.

To further understand the impact of different physical mechanisms and
observing parameters on such time delay measuring experiments, we propose the
following systematic explorations of the parameter space:

  1. (i)

Investigate the effects of different intrinsic variability, cadence, and
length of the observing period of lensed quasars.

  2. (ii)

Athough it seems that the presence of a high frequency reverberated flux
component does not bias the delays by much, different realizations of such
variability will need to be explored to better understand its effect on time
delays.

  3. (iii)

The effect of uncertainties in the mass model parameters of the lensing galaxy
is reflected in microlensing through the $\kappa,\gamma$. Similarly, $s$ is
accompanied by assumptions on the stellar mass distribution and mass-to-light
ratio. These uncertainties can be incorporated by sampling a distribution of
these parameters when generating mock sets of light curves. Most importantly,
this will allow to estimate the covariance between, for example, the mass
model and the accuracy and precision of the measured time delays.

For all these tests, independent information and physically motivated priors
and covariances between parameters can and should be assumed or tested.
Finally, in the specific case of LSST, a time delay measuring algorithm
fitting simultaneously data in many filters could increase the performance.

In our study, we have used only mock time series data, to which we arbitrarily
introduced some observationally motivated error bars according to the assumed
instrument, and we did not make any use of the mock cutout images other than
in illustrative examples. However, using directly the cutouts, e.g. within
light curve extraction pipelines, could simultaneously allow for
testing/improving these pipelines and for having real instrument error bars.
This could be achieved by, for example, injecting mock systems in telescope
images at the early stages of data reduction and performing the same
procedures on them as on the real lenses. As a next step, high-resolution
space-based images could be used to derive lens models, which in turn would be
used to create several mock datasets in order to achieve model-aware source
deblending and light curve extraction. This is particularly useful in cases
where multiple images are too close (e.g. merging pairs in fold configuration
lenses, or low mass lenses), or when the extended source is bright enough to
contribute a significant flux component at the locations of the multiple
images.

We have demonstrated how our software can be combined with other methods, in
this case time-delay measurement codes, in a streamlined fashion.
Specifically, we envisage our open-source modelling software to be used in
connection with: i) analytic or numerical population models or catalogs (e.g.
Oguri & Marshall, 2010; Collett, 2015), ii) time delay measurement methods,
especially linked to $H_{\mathrm{0}}$ measuring pipelines, where covariance
between the delays and the mass model parameters can be provided, iii) machine
learning methods that require light curves to study the accretion disc size
(Vernardos & Tsagkatakis, 2019), microlensing events, and time delays, or
static images to infer substructure (Vernardos et al., 2020). We advocate that
the MOLET software package can become an indispensable tool for data intensive
future strong lensing science.

## Data availability

The data that support the findings of this study are openly available in
github at https://github.com/gvernard/molet

## Acknowledgements

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodovska-Curie grant
agreement No 897124. The author has also received support from the Schmidt
Futures foundation. I would like to thank Chih-Fan Chen and Sherry Suyu for
sharing their results and useful suggestions. Thanks also to Eric Paic and
Dominique Sluse for generating intrinsic light curve realizations of RXJ
1131-1231, and to Martin Millon and Justin Pierel for setting up the fits with
PyCS and SNTD respectively.

## References

  * Akiyama et al. (2019) Akiyama K., et al., 2019, The Astrophysical Journal, 875, L6
  * Barkana (1998) Barkana R., 1998, The Astrophysical Journal, 502, 531
  * Bate & Fluke (2012) Bate N. F., Fluke C. J., 2012, The Astrophysical Journal, 744, 90
  * Birrer & Treu (2019) Birrer S., Treu T., 2019, Monthly Notices of the Royal Astronomical Society, 489, 2097
  * Bonvin et al. (2018) Bonvin V., Chan J. H., Millon M., Rojas K., Courbin F., Chen G. C., Fassnacht C. D., Paic E., 2018, Astronomy and Astrophysics, 616, 1
  * Cackett et al. (2007) Cackett E. M., Horne K., Winkler H., 2007, Monthly Notices of the Royal Astronomical Society, 380, 669
  * Chan et al. (2020) Chan J. H. H., Schive H.-Y., Wong S.-K., Chiueh T., 2020, preprint (astro-ph/2002.10473), pp 1–6
  * Chen et al. (2016) Chen G. C. F., et al., 2016, Monthly Notices of the Royal Astronomical Society, 462, 3457
  * Chen et al. (2018) Chen G. C., et al., 2018, Monthly Notices of the Royal Astronomical Society, 481, 1115
  * Collett (2015) Collett T., 2015, The Astrophysical Journal, 811, 20
  * Collett & Smith (2020) Collett T. E., Smith R. J., 2020, Monthly Notices of the Royal Astronomical Society, 497, 1654
  * Cornachione et al. (2020) Cornachione M. A., Morgan C. W., Millon M., Bentz M. C., Courbin F., Bonvin V., Falco E. E., 2020, The Astrophysical Journal, 895, 125
  * Dalal & Kochanek (2002) Dalal N., Kochanek C. S., 2002, The Astrophysical Journal, 572, 25
  * Despali et al. (2018) Despali G., Vegetti S., White S. D., Giocoli C., van den Bosch F. C., 2018, Monthly Notices of the Royal Astronomical Society, 475, 5424
  * Fausnaugh et al. (2017) Fausnaugh M. M., Peterson B. M., Starkey D. A., Horne K., AGN Storm Collaboration 2017, Frontiers in Astronomy and Space Sciences, 4, 55
  * Fluke & Webster (1999) Fluke C. J., Webster R. L., 1999, Monthly Notices of the Royal Astronomical Society, 302, 68
  * Foxley-Marrable et al. (2018) Foxley-Marrable M., Collett T. E., Vernardos G., Goldstein D. A., Bacon D., 2018, Monthly Notices of the Royal Astronomical Society, 478, 5081
  * Freedman (2017) Freedman W. L., 2017, Nature Astronomy, 1, 121
  * Gavazzi et al. (2008) Gavazzi R., Treu T., Koopmans L. V. E., Bolton A. S., Moustakas L. A., Burles S., Marshall P. J., 2008, The Astrophysical Journal, 677, 1046
  * Gilman et al. (2020a) Gilman D., Birrer S., Nierenberg A., Treu T., Du X., Benson A., 2020a, Monthly Notices of the Royal Astronomical Society, 491, 6077
  * Gilman et al. (2020b) Gilman D., Birrer S., Treu T., 2020b, Astronomy & Astrophysics, 642, A194
  * Goldstein et al. (2018) Goldstein D. A., Nugent P. E., Kasen D. N., Collett T. E., 2018, The Astrophysical Journal, 855, 22
  * Guy et al. (2007) Guy J., et al., 2007, Astronomy and Astrophysics, 466, 11
  * Harvey et al. (2020) Harvey D., Valkenburg W., Tamone A., Boyarsky A., Courbin F., Lovell M., 2020, Monthly Notices of the Royal Astronomical Society, 491, 4247
  * Huber et al. (2019) Huber S., et al., 2019, Astronomy & Astrophysics, 631, A161
  * Kayser et al. (1986) Kayser R., Refsdal S., Stabell R., 1986, Astronomy & Astrophysics, 166, 36
  * Keeton (2001) Keeton C. R., 2001, preprint (astro-ph/0102341)
  * Keeton (2010) Keeton C. R., 2010, General Relativity and Gravitation, 42, 2151
  * Keeton & Moustakas (2009) Keeton C. R., Moustakas L. A., 2009, Astrophysical Journal, 699, 1720
  * Kessler et al. (2009) Kessler R., et al., 2009, Publications of the Astronomical Society of the Pacific, 121, 1028
  * Kochanek (2004) Kochanek C. S., 2004, The Astrophysical Journal, 605, 58
  * Komatsu et al. (2011) Komatsu E., et al., 2011, Astrophysical Journal, Supplement Series, 192, 18
  * Kormann et al. (1994) Kormann R., Schneider P., Bartelmann M., 1994, Astronomy & Astrophysics, 284, 285
  * Kundic & Wambsganss (1993) Kundic T., Wambsganss J., 1993, The Astrophysical Journal, 404, 455
  * LSST Science Collaborations (2009) LSST Science Collaborations 2009, preprint (astro-ph/0912.0201)
  * Lage (2019) Lage C., 2019, preprint (astro-ph/1911.09577)
  * Lewis & Irwin (1995) Lewis G. F., Irwin M. J., 1995, Monthly Notices of the Royal Astronomical Society, 276, 103
  * Liao (2019) Liao K., 2019, The Astrophysical Journal, 871, 113
  * MacLeod et al. (2010) MacLeod C. L., et al., 2010, The Astrophysical Journal, 721, 1014
  * MacLeod et al. (2013) MacLeod C. L., Jones R., Agol E., Kochanek C. S., 2013, The Astrophysical Journal, 773, 35
  * Metcalf & Zhao (2002) Metcalf R. B., Zhao H., 2002, The Astrophysical Journal Letters, 567, L5
  * Millon et al. (2020a) Millon M., Tewes M., Bonvin V., Lengen B., Courbin F., 2020a, Journal of Open Source Software, 5, 2654
  * Millon et al. (2020b) Millon M., et al., 2020b, Astronomy and Astrophysics, 639, 1
  * Millon et al. (2020c) Millon M., et al., 2020c, Astronomy & Astrophysics, 640, A105
  * Morgan et al. (2012) Morgan C. W., et al., 2012, The Astrophysical Journal, 756, 52
  * Morgan et al. (2018) Morgan C. W., Hyer G. E., Bonvin V., Mosquera A. M., Cornachione M., Courbin F., Kochanek C. S., Falco E. E., 2018, The Astrophysical Journal, 869, 106
  * Mortonson et al. (2005) Mortonson M. J., Schechter P. L., Wambsganss J., 2005, The Astrophysical Journal, 628, 594
  * Mosquera & Kochanek (2011) Mosquera A. M., Kochanek C. S., 2011, The Astrophysical Journal, 738, 96
  * Mukherjee et al. (2018) Mukherjee S., et al., 2018, Monthly Notices of the Royal Astronomical Society, 479, 4108
  * Neira et al. (2020) Neira F., Anguita T., Vernardos G., 2020, Monthly Notices of the Royal Astronomical Society, 495, 544
  * Nierenberg et al. (2017) Nierenberg A. M., et al., 2017, Monthly Notices of the Royal Astronomical Society, 471, 2224
  * Oguri (2019) Oguri M., 2019, Reports on Progress in Physics, 82
  * Oguri & Marshall (2010) Oguri M., Marshall P. J., 2010, Monthly Notices of the Royal Astronomical Society, 405, 2579
  * Padovani et al. (2017) Padovani P., et al., 2017, Astronomy and Astrophysics Review, 25
  * Paic (2021) Paic E., 2021, in preparation
  * Petkova et al. (2014) Petkova M., Metcalf R. B., Giocoli C., 2014, Monthly Notices of the Royal Astronomical Society, 445, 1954
  * Planck Collaboration et al. (2020) Planck Collaboration P., et al., 2020, Astronomy & Astrophysics, 641, A6
  * Plazas (2020) Plazas A. A., 2020, Symmetry, 12
  * Poindexter & Kochanek (2010) Poindexter S., Kochanek C. S., 2010, The Astrophysical Journal, 712, 658
  * Riess et al. (2016) Riess A. G., et al., 2016, The Astrophysical Journal, 826, 56
  * Rusu et al. (2017) Rusu C. E., et al., 2017, Monthly Notices of the Royal Astronomical Society, 467, 4220
  * Sartori et al. (2019) Sartori L. F., Trakhtenbrot B., Schawinski K., Caplar N., Treister E., Zhang C., 2019, The Astrophysical Journal, 883, 139
  * Schechter & Wambsganss (2002) Schechter P. L., Wambsganss J., 2002, The Astrophysical Journal, 580, 685
  * Schmidt & Wambsganss (2010) Schmidt R. W., Wambsganss J., 2010, Gen. Relativ. and Gravit., 42, 2127
  * Schneider et al. (2006) Schneider P., Kochanek C. S., Wambsganss J., 2006, in Meylan G., Jetzer P., North P., eds, Saas-Fee Advanced Course vol. 33. Springer, Berlin
  * Sersic (1963) Sersic J. L., 1963, Boletín de la Asociación Argentina de Astronomía, 6, 41
  * Shajib et al. (2019) Shajib A. J., et al., 2019, Monthly Notices of the Royal Astronomical Society, 483, 5649
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, Astronomy & Astrophysics, 24, 337
  * Sluse & Tewes (2014) Sluse D., Tewes M., 2014, Astronomy and Astrophysics, 571, 1
  * Sluse et al. (2006) Sluse D., Claeskens J.-F., Hutsemékers D., Surdej J., 2006, Astronomy & Astrophysics, 449, 539
  * Stacey & McKean (2018) Stacey H. R., McKean J. P., 2018, Monthly Notices of the Royal Astronomical Society: Letters, 481, L40
  * Suyu et al. (2013) Suyu S. H., et al., 2013, The Astrophysical Journal, 766, 70
  * Suyu et al. (2020) Suyu S. H., Huber S., Cañameras R., Kromer M., Schuldt S., Taubenberger S., Yıldırım A., Bonvin V., 2020, Astronomy & Astrophysics, 644, A162
  * Tewes et al. (2013) Tewes M., Courbin F., Meylan G., 2013, Astronomy & Astrophysics, 553, 120
  * Tie & Kochanek (2018) Tie S. S., Kochanek C. S., 2018, Monthly Notices of the Royal Astronomical Society, 473, 80
  * Treu & Marshall (2016) Treu T., Marshall P. J., 2016, Astronomy and Astrophysics Review, 24
  * Vernardos (2019) Vernardos G., 2019, Monthly Notices of the Royal Astronomical Society, 483, 5583
  * Vernardos & Fluke (2013) Vernardos G., Fluke C. J., 2013, Monthly Notices of the Royal Astronomical Society, 434, 832
  * Vernardos & Fluke (2014) Vernardos G., Fluke C. J., 2014, Astronomy & Computing, 6, 1
  * Vernardos & Tsagkatakis (2019) Vernardos G., Tsagkatakis G., 2019, Monthly Notices of the Royal Astronomical Society, 486, 1944
  * Vernardos et al. (2014) Vernardos G., Fluke C. J., Bate N. F., Croton D. J., 2014, The Astrophysical Journal Supplement Series, 211, 16
  * Vernardos et al. (2015) Vernardos G., Fluke C. J., Bate N. F., Croton D. J., Vohl D., 2015, The Astrophysical Journal Supplement Series, 217, 23
  * Vernardos et al. (2020) Vernardos G., Tsagkatakis G., Pantazis Y., 2020, Monthly Notices of the Royal Astronomical Society, 499, 5641
  * Wambsganss et al. (1990) Wambsganss J., Paczynski B., Schneider P., 1990, The Astrophysical Journal, 358, L33
  * Wong et al. (2017) Wong K. C., et al., 2017, Monthly Notices of the Royal Astronomical Society, 465, 4895
  * Wong et al. (2020) Wong K. C., Suyu S. H., Chen G. C.-F., Rusu C. E., Millon M., Sluse D., Bonvin V., Fassnacht C. D., 2020, Monthly Notices of the Royal Astronomical Society, 498, 1420
  * Wyithe & Turner (2001) Wyithe J. S. B., Turner E. L., 2001, Monthly Notices of the Royal Astronomical Society, 320, 21
  * Wyithe et al. (2000) Wyithe J. S. B., Webster R. L., Turner E. L., 2000, Monthly Notices of the Royal Astronomical Society, 312, 843
  * de Grijs et al. (2017) de Grijs R., Courbin F., Martínez-Vázquez C. E., Monelli M., Oguri M., Suyu S. H., 2017, Space Science Reviews, 212, 1743