# Harnessing Unresolved Lensed Quasars: The Mathematical Foundation of the
Fluctuation Curve

Satadru Bag<EMAIL_ADDRESS>Korea Astronomy and Space Science Institute,
Daejeon 34055, Korea  Wuhyun Sohn Korea Astronomy and Space Science Institute,
Daejeon 34055, Korea  Arman Shafieloo<EMAIL_ADDRESS>Korea Astronomy and
Space Science Institute, Daejeon 34055, Korea  University of Science and
Technology, Daejeon 34113, Korea  Kai Liao School of Physics and Technology,
Wuhan University, Wuhan 430072, China

###### Abstract

Strong gravitational lensed quasars (QSOs) have emerged as powerful and novel
cosmic probes as they can deliver crucial cosmological information, such as a
measurement of the Hubble constant, independent of other probes. Although the
upcoming LSST survey is expected to discover $10^{3}-10^{4}$ lensed QSOs, a
large fraction will remain unresolved due to seeing. The stochastic nature of
the quasar intrinsic flux makes it challenging to identify lensed ones and
measure the time delays using unresolved light curve data only. In this
regard, Bag et al. (2022) introduced a data-driven technique based on the
minimization of the fluctuation in the reconstructed image light curves. In
this article, we delve deeper into the mathematical foundation of this
approach. We show that the lensing signal in the fluctuation curve is
dominated by the auto-correlation function (ACF) of the derivative of the
joint light curve. This explains why the fluctuation curve enables the
detection of the lensed QSOs only using the joint light curve, without making
assumptions about QSO flux variability, nor requiring any additional
information. We show that the ACF of the derivative of the joint light curve
is more reliable than the ACF of the joint light curve itself because
intrinsic quasar flux variability shows significant auto-correlation up to a
few hundred days (as they follow a red power spectrum). In addition, we show
that the minimization of fluctuation approach provides even better precision
and recall as compared to the ACF of the derivative of the joint light curve
when the data have significant observational noise.

## References

  *

##  1 Introduction

Strong gravitational lensed systems have emerged as a powerful and novel
cosmic probe (see, e.g.,Treu & Marshall (2016) for a review). They can deliver
cosmological information independent of other probes such as the Type Ia
supernovae (SNe), Baryon Acoustic Osculations (BAO), and Cosmic Microwave
Background (CMB). Time delay measurements, together with accurate lens
modelling, allow us to directly estimate the present epoch value of the cosmic
expansion rate, i.e. the Hubble constant ($H_{0}$) (Refsdal & Bondi, 1964;
Refsdal, 1964; Saha et al., 2006; Oguri, 2007; Bonvin et al., 2017; Wong et
al., 2020; Birrer et al., 2020; Birrer & Treu, 2021). Therefore, ‘time delay
cosmography’ can play a crucial role in elucidating the ongoing $H_{0}$
tension between the local measurements (Riess et al., 2022; Abdalla et al.,
2022) and early universe probes like the CMB (Aghanim et al., 2020). Other
applications of strong lensing in cosmology and astrophysics are summarised in
the review by Treu (2010).

For time delay measurements, one needs time-variable sources, such as quasars
(QSOs) and supernovae (SNe). Lensed SNe(Oguri, 2019; Liao et al., 2022; Suyu
et al., 2023) are extremely rare as only four with multiple images have been
discovered so far (Kelly et al., 2015; Goobar et al., 2017; Rodney et al.,
2021; Goobar et al., 2022). In comparison, lensed QSOs are more abundant and
thus remain to be the primary source for the time delay cosmography (however,
lensed SNe could be at the forefront of time delay cosmography in the next
decade (Suyu et al., 2020)). Although hundreds of lensed QSOs are known (Lemon
et al., 2022), only a few have been used for cosmology. For example, using
only six ‘good quality’ lensed QSOs the H0LiCOW team (Suyu et al., 2017)
measured the Hubble constant with $2.4\%$ uncertainty (Wong et al., 2020),
under standard assumptions about the mass density profile of the deflector,
and a seventh brings the precision to 2% (Shajib et al., 2020; Millon et al.,
2020). However, if one drops the assumptions and adopts density profiles that
are maximally degenerate with $H_{0}$ through the mass sheet degeneracy, the
uncertainty increases to $9\%$ (Birrer et al., 2020), highlighting the need
for substantially larger samples. While the uncertainty can also be reduced by
additional information per lens (especially stellar kinematics), a powerful
way of improving the $H_{0}$ precision is to increase the sample volume
significantly (Sonnenfeld, 2021). For example, observations of hundreds of
lensed systems will deliver sub-percent uncertainty and accuracy, regardless
of any assumption on their mass density profile (Birrer & Treu, 2021; Jee et
al., 2016).

The angular separation of images are typically of the order $\sim 1-2$ arcsec
for galaxy scale lenses (Narayan & Bartelmann, 1996; Treu, 2010). Therefore,
one is required to first resolve the images by using either a sufficiently
high-resolution (ground-based or space) telescope or through spectroscopy.
Then the individual light curves need to be monitored at sufficient resolution
for several years in order to measure the time delays (e.g., Tewes et al.,
2013; Liao et al., 2015; Millon et al., 2020). This can be difficult as these
observations are expensive. On the other hand, unresolved light curves may be
available “for free” from synoptic surveys. For example, we expect a lot of
lensed QSOs to be partially resolved or completely unresolved in the wide
field surveys, such as Zwicky Transient Facility (ZTF) (Bellm et al., 2019)
and Legacy Survey of Space and Time (LSST) (LSST Science Collaboration et al.,
2009, 2017), because their angular resolution is limited by seeing. In such
cases, we can observe the joint light curve of the unresolved system that is a
blend of the individual light curves. A robust method of detecting lensed QSOs
through unresolved light curves can take advantage of the more abundant
smaller telescopes. Thus, the importance of this approach cannot be overstated
for boosting the sample size of the observed lensed systems 111Although this
work focuses on the unresolved lensed QSOs, similar exercises for the
unresolved lensed supernovae have been pursued in the literature (Bag et al.,
2021; Denissenya et al., 2021; Denissenya & Linder, 2022)..

There are multiple other advantages in working with unresolved systems. Since
there is no need for resolving the images a priori, this approach can be
applied to the light curve data from ongoing time domain surveys such as ZTF
(Bellm et al., 2019), Pan-STARRS1 (Chambers et al., 2016). This will become
more crucial when the upcoming Vera Rubin Observatory starts the LSST (LSST
Science Collaboration et al., 2009, 2017). This approach also evades any
degeneracy between a binary QSO pair and a doubly lensed QSO that creates
confusion in lens detection using the resolved photometry (Peng et al., 1999;
Mortlock et al., 1999).

Recently, multiple different techniques have been proposed to identify the
lensed systems and to measure their time delays using the joint light curves
(Shu et al., 2021; Springer & Ofek, 2021a, b; Biggio et al., 2022). The
primary challenge in detecting the lensed cases using the joint light curves
is that the intrinsic QSO light curves are highly stochastic and show a vast
diversity in the flux variability. Therefore, any assumption on the flux
variability can lead to biased results with low precision and high false
positive detection rate when the real light curves are not well described by
the assumption. Therefore, it is extremely important to be model agnostic for
achieving higher recall and precision as well as for reliable unbiased
results.

In paper 1 (Bag et al., 2022) we introduced a novel data-driven method that
can detect the lensed QSOs and simultaneously measures the time delays only
using the joint light curves, most importantly neither assuming anything about
the quasar flux variability nor using any additional information. The
technique is based on the empirical observation that the reconstructed image
light curves corresponding to incorrect time delays exhibit more fluctuation
than the ones reconstructed using the correct time delay. Although Bag et al.
(2022) demonstrates that this approach is successful in the presence of
significant noise (e.g. ZTF-like noise) and on existing data quality, it lacks
a formal explanation as to how the minimization of the fluctuation works. This
article looks deeper into the mathematical formalism of this approach. We
attempt to understand the mathematical reasoning behind the empirical
observations which laid the foundation of this approach. A clear insight into
the mechanism should allow us to explore the strengths and possible
limitations of this approach.

The paper is organized as follows. Section 2 recapitulates the method
introduced by Bag et al. (2022). In Section 3, we provide a detailed
mathematical explanation for the method’s ability to detect lensed systems. We
show that the lensing signal in the fluctuation curve is dominated by the
auto-correlation function (ACF) of the derivative of the joint light curve.
Section 4 explains why the ACF of the derivative of the joint light curve
performs better in identifying lensed systems than the ACF of the joint light
curve itself. Finally, in Section 5 we compare the performance of ACF of the
derivative of the joint light curves against the full fluctuation curves. We
conclude our findings in the Section 6.

##  2 Reconstructing the underlying light curves and fluctuations in them

The joint light curve of an unresolved lensed QSO having $N_{I}$ images is the
sum of the image light curves,

$F(t)=\sum_{j=1}^{N_{I}}a_{j}\mathcal{F}(t-T_{j})\;,$ (1)

where the individual image light curves can be described by a common function
$\mathcal{F}(t)$ but with different magnifications ($a_{j}$) and time delays
($T_{j}$). To break the degeneracy among the images, we choose the ordering
such that $a_{1}\geq a_{2}\geq\cdots\geq a_{N_{I}}$ without loss of
generality. For simplicity, let us first consider double systems (two images)
with

$F(t)=f(t)+\widehat{\mu}f(t-\widehat{\Delta t})~{},{\rm
where}~{}~{}f(t)=a_{1}\mathcal{F}(t-T_{1})~{},~{}\widehat{\mu}=a_{2}/a_{1}~{}{\rm
and}~{}\widehat{\Delta t}=T_{2}-T_{1}\;.$ (2)

Here, $f(t)$ is the light curve of the brighter image which we take as a
reference. The other (fainter) image’s light curve is then given by
$\widehat{\mu}f(t-\widehat{\Delta t})$, where $\widehat{\mu}\leq 1$ and
$\widehat{\Delta t}$ correspond to the magnification ratio and time delay,
respectively, with respect to the reference (brighter) image. Note that
$\widehat{\Delta t}$ can be positive or negative; a positive (negative)
$\widehat{\Delta t}$ implies that the fainter image arrives later (earlier) in
time than the brighter image.

It is difficult to model the quasar flux variability due to its highly
stochastic nature. Instead, one can reconstruct the light curve of the
brighter image following Bag et al. (2022) (see also Geiger & Schneider, 1996)
as

$\displaystyle f_{\rm rec}(t;\mu,\Delta t)$ $\displaystyle=F(t)-\mu f_{\rm
rec}(t-\Delta t;\mu,\Delta t)\;,$ (3) $\displaystyle=F(t)-\mu F(t-\Delta
t)+\mu^{2}F(t-2\Delta t)$ $\displaystyle\qquad\qquad-\mu^{3}F(t-3\Delta
t)+\mu^{4}F(t-4\Delta t)-\ldots\;,$
$\displaystyle=\sum_{n=0}^{\infty}\left(-\mu\right)^{n}F(t-n\cdot\Delta t)\;,$
(4)

from the joint light curve $F(t)$, given $\Delta t$ and $\mu<1$ which need not
be equal to the true underlying values. In this article we denote the true
magnification ratio and time delay by $\\{\widehat{\mu},~{}\widehat{\Delta
t}\\}$ to avoid confusion with the free parameters of the reconstruction,
$\\{\mu,~{}\Delta t\\}$, which are same as $\\{\mu_{\rm try},~{}\Delta t_{\rm
try}\\}$ in the paper 1 (Bag et al., 2022).

Note that the higher order terms in Eq. (4) require $F(t)$ outside its
observed range. Since the quasar flux cannot be predicted, we assume that
$F(t)$ remains flat outside its observed range. As discussed in Geiger &
Schneider (1996); Bag et al. (2022), this convention has a negligible effect
on the reconstruction except near either of the boundaries (depending upon the
sign of $\Delta t$).

We emphasise that by restricting ourselves to $\mu<1$, that ensures the
convergence of the sum in Eq. (4), we are essentially reconstructing the light
curve of the brighter image without any loss of generality. Moreover, for any
choice of $\\{\mu<1,~{}\Delta t\\}$, Eq. (4) gives a unique solution for the
brightest image light curve; the corresponding fainter image light curve would
be $\mu f_{\rm rec}(t-\Delta t;\mu,\Delta t)$. Only when $\mu=\widehat{\mu}$
and $\Delta t=\widehat{\Delta t}$, however, Eq. (4) recovers the true image
light curve, $f_{\rm rec}(t)=f(t)$. Still, the equation $F(t)=f_{\rm
rec}(t;\mu,\Delta t)+\mu f_{\rm rec}(t-\Delta t;\mu,\Delta t)$ is exactly
satisfied for any choice of $\\{\mu<1,~{}\Delta t\\}$ by construction. This
demonstrates the mathematical degeneracy present in this lensed detection
problem as discussed in Geiger & Schneider (1996); Bag et al. (2022); without
any prior assumptions on $f(t)$, any choice of $\\{\mu,~{}\Delta t\\}$ can
yield a (unique) lensing solution to the joint light curve $F(t)$.

In paper 1 (Bag et al., 2022), we introduce a completely data-driven technique
to break the degeneracy in the time delay by minimizing the fluctuation in
reconstructed image light curves. To quantify the amount of fluctuation in a
reconstructed image light curve $f_{\rm rec}(t;\mu,\Delta t)$, Bag et al.
(2022) uses the simple metric,

$\epsilon(\Delta t;\mu)\equiv\sum_{i=0}^{N_{D}-1}\left[f_{\rm
rec}(t_{i+1};\mu,\Delta t)-f_{\rm rec}(t_{i};\mu,\Delta t)\right]^{2}\;.$ (5)

Fixing the trial magnification ratio ($\mu$) to an arbitrary value (less than
unity), one reconstructs the brighter image light curve for a number of
$\Delta t$ and finally looks for minima in the fluctuation curve,
$\epsilon(\Delta t)$.

Figure 1: The fluctuation curves for an example of a double system simulated
(in perfect condition with negligible observation noise) using damped random
walk (DRW) process with time delay $\widehat{\Delta t}=22.0$ days and the
magnification ratio $\widehat{\mu}=0.86$. The amount of fluctuation
($\epsilon$ calculated from Eq. (5)) in the reconstructed image light curves
is plotted against the trial $\Delta t$. The three panels correspond to three
arbitrary choices of trial magnification ratio, $\mu=0.3,0.5,0.7$. The dashed
vertical lines in each plot mark the true time delay, $\Delta
t=\pm\widehat{\Delta t}$, where we find prominent pairs of secondary minima
that can be used to detect the lensed system and measure the time delay. This
figure is qualitatively similar to Fig. 2 of Bag et al. (2022).

Fig. 1 demonstrates how one can detect the unresolved lensed QSOs through the
fluctuation curves by considering an example of a double (2-image) system that
is simulated using the damped random walk (DRW) process in perfect conditions
(marginal observational noise and 1 day of cadence). The true magnification
ratio and time delay are set to: $\widehat{\mu}=0.86$ and $\widehat{\Delta
t}=22.0$ day. The three panels show $\epsilon(\Delta t)$ for different
(arbitrary) choices of $\mu$. We make the following observations from Fig. 1.

  * •

The fluctuation curve $\epsilon(\Delta t)$ is highly symmetric with respect to
$\Delta t=0$.

  * •

There exists a global minimum at $\Delta t=0$; the amount of fluctuations in
the reconstructed light curve is minimised when we assume that the system is
not lensed. This is a generic feature of all fluctuation curves irrespective
of the system being lensed or not, and hence can be ignored 222While trying to
detect the lensed QSOs through the minimal fluctuation, as proposed by Bag et
al. (2022), this unlensed solution always remains to be a feasible lensing
solution to the unresolved problem in the absence of any assumption on the
flux variability..

  * •

Strikingly, when $\Delta t=\pm\widehat{\Delta t}$ (shown by the dashed
vertical lines in each panel), we observe a pair of prominent secondary minima
that correctly identify the system as lensed and simultaneously estimate the
time delay accurately. On the other hand, if the system is unlensed (i.e.
$\widehat{\mu}=0$ or $\widehat{\Delta t}=0$) we do not find any prominent
secondary minimum in $\epsilon(\Delta t)$ as demonstrated in Bag et al.
(2022).

  * •

All of the observations above are somewhat insensitive to the choice of $\mu$.
Therefore, the true magnification ratio cannot be estimated in this approach.
Nevertheless, one can accurately measure the time delays from the location of
the pair of secondary minima in the fluctuation curves.

In the next section, we provide the mathematical reasoning behind all the
above empirical observations.

##  3 Why is fluctuation minimised for correct time delay?

Substituting Eq. (4) into Eq. (5), one gets

$\displaystyle\epsilon(\Delta t)$
$\displaystyle=\sum_{i}^{N_{D}-1}\left[\sum_{n}^{\infty}\\{\left(-\mu\right)^{n}F(t_{i+1}-n\Delta
t)\\}-\sum_{m}^{\infty}\\{\left(-\mu\right)^{m}F(t_{i}-m\Delta
t)\\}\right]^{2},$ (6)
$\displaystyle=\sum_{i}^{N_{D}-1}\left[\\{F(t_{i+1})-F(t_{i})\\}-\mu\\{F(t_{i+1}-\Delta
t)-F(t_{i}-\Delta t)\\}\right.$
$\displaystyle\qquad\left.+\mu^{2}\\{F(t_{i+1}-2\Delta t)-F(t_{i}-2\Delta
t)\\}-\mu^{3}\\{F(t_{i+1}-3\Delta t)-F(t_{i}-3\Delta t)\\}+\ldots\right]^{2}.$
(7)

For convenience, let us define the difference in the successive observed joint
flux as a separate time series,

$H(t_{i})\equiv F(t_{i+1})-F(t_{i})\;.$ (8)

Using Eq. (8) one can simplify Eq. (7) to

$\displaystyle\epsilon(\Delta t)$ $\displaystyle=\sum_{i}\left[H(t_{i})-\mu
H(t_{i}-\Delta t)+\mu^{2}H(t_{i}-2\Delta t)-\mu^{3}H(t_{i}-3\Delta
t)+\ldots\right]^{2}=\sum_{i}\left[\sum_{n}^{\infty}\left(-\mu\right)^{n}H(t_{i}-n\cdot\Delta
t)\right]^{2}.$ (9)

Further expanding the squared term, we get

$\displaystyle\epsilon(\Delta t)$
$\displaystyle=\sum_{i}H(t_{i})^{2}-2\mu\sum_{i}H(t_{i})H(t_{i}-\Delta
t)+\mu^{2}\left[\sum_{i}H(t_{i}-\Delta t)^{2}+2\sum_{i}H(t_{i})H(t_{i}-2\Delta
t)\right]$ (10)
$\displaystyle\qquad-2\mu^{3}\left[\sum_{i}H(t_{i})H(t_{i}-3\Delta
t)+\sum_{i}H(t_{i}-\Delta t)H(t_{i}-2\Delta t)\right]+\ldots\;.$

Thus we can express $\epsilon(\Delta t)$ as a power series in $\mu$,

$\epsilon(\Delta t)=\epsilon_{0}+\mu\epsilon_{1}(\Delta
t)+\mu^{2}\epsilon_{2}(\Delta t)+\mu^{3}\epsilon_{3}(\Delta
t)+\mu^{4}\epsilon_{4}(\Delta t)+\ldots\;,$ (11)

where

$\displaystyle\epsilon_{0}$ $\displaystyle=\sum_{i}H(t_{i})^{2}\;,$ (12)
$\displaystyle\epsilon_{1}(\Delta t)$
$\displaystyle=-2\sum_{i}H(t_{i})H(t_{i}-\Delta t)\;,$ (13)
$\displaystyle\epsilon_{2}(\Delta t)$ $\displaystyle=\sum_{i}H(t_{i}-\Delta
t)^{2}+2\sum_{i}H(t_{i})H(t_{i}-2\Delta t)\;,$ (14)
$\displaystyle\epsilon_{3}(\Delta t)$
$\displaystyle=-2\sum_{i}H(t_{i})H(t_{i}-3\Delta t)-2\sum_{i}H(t_{i}-\Delta
t)H(t_{i}-2\Delta t)\;,$ (15)

and so on. Therefore, the fluctuation curve can be written in the following
closed form,

$\epsilon(\Delta
t)=\sum_{n=0}^{\infty}\left(-\mu\right)^{n}\sum_{m=0}^{\lfloor
n/2\rfloor}\left(2-\delta_{m,n-m}\right)\sum_{i}\left[H(t_{i}-m\Delta
t)H(t_{i}-(n-m)\Delta t)\right]\;,$ (16)

where the floor function $\lfloor n/2\rfloor$ returns the highest integer
equal to or below $n/2$. Since $\mu<1$, the contribution from the higher order
terms gets suppressed rapidly.

Let us now take a closer look at the difference series $H(t)$ (defined in Eq.
(8)), which can be recast as

$\displaystyle H(t_{i})$
$\displaystyle=\left[f(t_{i+1})-f(t_{i})\right]+\widehat{\mu}\left[f(t_{i+1}-\widehat{\Delta
t})-f(t_{i}-\widehat{\Delta t})\right]\;,$ (17)
$\displaystyle=h(t_{i})+\widehat{\mu}~{}h(t_{i}-\widehat{\Delta t})\;,$ (18)

where

$h(t_{i})\equiv f(t_{i+1})-f(t_{i})$ (19)

can be regarded as another time series. For a uniformly sampled data, $H(t)$
and $h(t)$ are proportional to the derivative of $F(t)$ and $f(t)$,
respectively. In this article, we refer to $H(t)$ and $h(t)$ as the
derivatives of $F(t)$ and $f(t)$, respectively. However, for non-uniformly
sampled data, the former two are just difference series and the mathematical
arguments remain intact. From Eq. (18) note that $H(t)$ follows similar
lensing equation as $F(t)$ in Eq. (2) but with $h(t)$ as the underlying time
series.

The QSO intrinsic light curve, $f(t)$, follows a stochastic process that is
‘wide-sense stationary’; the mean and covariance properties remain constant
over the time. It is easy to check that the time series $F(t)$, $h(t)$ and
$H(t)$ are also wide-sense stationary. For a generic wide-sense stationary
time series $X(t)$,

$\displaystyle\langle X(t)\rangle_{t}\equiv\frac{1}{N_{D}}\sum_{i}X(t_{i})$
$\displaystyle\approx\langle X(t+t_{c})\rangle_{t}\;,$ (20)
$\displaystyle\sigma^{2}_{X(t)}\equiv\langle X(t)^{2}\rangle_{t}-\langle
X(t)\rangle_{t}^{2}$ $\displaystyle\approx\langle
X(t+t_{c})^{2}\rangle_{t}-\langle
X(t+t_{c})\rangle_{t}^{2}\equiv\sigma^{2}_{X(t+t_{c})}$ (21)

and

$\langle X(t+t_{c_{1}})X(t+t_{c_{2}})\rangle_{t}\approx\langle
X(t)X(t+\left[t_{c_{2}}-t_{c_{1}}\right])\rangle_{t}\approx\text{function
of}~{}|t_{c_{2}}-t_{c_{1}}|~{}\text{only}\;,$ (22)

for any shift in time much smaller than the observed time range
($t_{c},~{}t_{c_{1}},~{}t_{c_{2}}\ll t_{\rm range}$). Therefore,

$\langle H(t)\rangle_{t}=\langle F(t_{i+1})\rangle_{t}-\langle
F(t_{i})\rangle_{t}\approx 0$ (23)

since we assume that the statistical properties of the joint light curve
remain invariant under translations, i.e. $\langle
F(t_{i+1})\rangle_{t}\approx\langle F(t_{i})\rangle_{t}$.

Now that we assembled all the necessary tools, we proceed to explain the
characteristics of the fluctuation curve that are described in Section 2 and
illustrated in Fig. 1.

  * •

First, notice that Eq. (10) is invariant under $\Delta t\to-\Delta t$ as long
as $|\Delta t|\ll t_{\rm range}$, since $H(t)$ also follows Eqs. (20) – (22).
This explains the symmetry in the $\epsilon(\Delta t)$ curve with respect to
$\Delta t=0$.

  * •

The first term in Eq. (10) or (11), $\epsilon_{0}$, is constant (independent
of $\Delta t$) and hence can be ignored. The next leading order term
$\epsilon_{1}(\Delta t)$ is proportional to the correlation coefficient
between $H(t)$ and $H(t-\Delta t)$ defined as (also known as the auto-
correlation function (ACF), see Eq. (A1))

${\rm ACF}(H;\Delta t)=\frac{\sum_{i}\left[H(t)-\langle
H(t)\rangle_{t}\right]\left[H(t-\Delta t)-\langle H(t-\Delta
t)\rangle_{t}\right]}{\sqrt{\sum_{i}\left[H(t)-\langle
H(t)\rangle_{t}\right]^{2}\cdot\sum_{i}\left[H(t-\Delta t)-\langle H(t-\Delta
t)\rangle_{t}\right]^{2}}}\approx\frac{1}{\epsilon_{0}}\sum_{i}H(t_{i})H(t_{i}-\Delta
t)\;,$ (24)

where one can identify that $\langle H(t)\rangle_{t}\approx\langle H(t-\Delta
t)\rangle_{t}\approx 0$ and the denominator simply reduces to
$\sum_{i}H(t_{i})^{2}=\epsilon_{0}$. Therefore, Eq. (10) can be approximated
as

$\epsilon(\Delta t)\approx\epsilon_{0}\left[1-2\mu{\rm ACF}(H;\Delta
t)\right]+\mathcal{O}(\mu^{2})\;.$ (25)

  * •

The auto-correlation function, ${\rm ACF}(H;\Delta t)$, is maximised to unity
at $\Delta t=0$. This in turn minimizes the $\epsilon_{1}(\Delta t)$ term that
dominates the $\Delta t$ dependence of the fluctuation curve in Eq. (11) at
the leading order in $\mu$ 333This could also be obtained from the inequality,
$\epsilon_{1}(\Delta t)=-2\sum_{i}H(t_{i})H(t_{i}-\Delta
t)\geq-\sum_{i}H(t_{i})^{2}-\sum_{i}H(t_{i}-\Delta
t)^{2}\approx-2\sum_{i}H(t_{i})^{2}=-2\epsilon_{0}$ (26) following the fact
that $2xy\leq x^{2}+y^{2}$. The equality holds when $H(t_{i})=H(t_{i}-\Delta
t)$ for all $i$ requiring $\Delta t=0$ unless $H(t)$ is a constant function.
Therefore, $\epsilon_{1}(\Delta t)$ is minimized at $\Delta t=0$. Note that
this argument does not require $\langle H(t)\rangle_{t}\approx 0$. .
Therefore, one always observes a global minimum in the $\epsilon(\Delta t)$
curve at $\Delta t=0$ irrespective of the system being lensed or not:
$\epsilon(0)\approx\epsilon_{0}-2\mu\epsilon_{0}+\mathcal{O}(\mu^{2})$. In
fact, it is easy to calculate the exact value
$\epsilon(0)=\epsilon_{0}/(1+\mu)^{2}$ as given in Appendix B.1.

  * •

When the system is lensed, ${\rm ACF}(H;\Delta t)$ is also maximised locally
at $\Delta t=\pm\widehat{\Delta t}$ as the correlation function picks up
excess power due to the matching of two sets of same intrinsic features that
are separated by the true time delay in the joint light curve (see Appendix A
for detailed derivations). Thus the fluctuation curve, $\epsilon(\Delta t)$,
shows a pair of secondary minima at $\Delta t=\pm\widehat{\Delta t}$ following
Eq. (25). When $h$ is pure white noise, ${\rm ACF}(H;\Delta
t=\pm\widehat{\Delta t})\approx\widehat{\mu}/(1+\widehat{\mu}^{2})$ as evident
from Eq. (A6) for long observation ranges. The exact height of the secondary
minima in $\epsilon(\Delta t)$ is calculated in Appendix B.2 under this
approximation. On the other hand, when the system is unlensed
($\widehat{\mu}=0$), we get no such secondary minima. Therefore, by detecting
this pair of secondary minima, which are the dominant part of the lensing
features in the fluctuation curve, one can identify the system as lensed.
Simultaneously, the location of this minima pair allows us to estimate the
time delay of the system.

  * •

As the $\epsilon_{i}$ terms in Eq. (11), which contain the $\Delta t$
variable, are independent of the choice of $\mu$, the locations of the lensing
features (local minima and maxima due to lensing) in the fluctuation curve
$\epsilon(\Delta t)$ is very much insensitive to $\mu$. (However, note that
for a higher value of $\mu$ higher order terms in Eq. (11) also contribute
which makes the curve more fluctuating as evident from the three panels of
Fig. 1.) Therefore, we cannot determine the true magnification ratio in this
approach, although the time delay can be recovered accurately 444In an ideal
condition where $h$ is uncorrelated in time, the observation range is long,
data have negligible noise etc, one can in principle determine the
magnification ratio from the height of the secondary minima at
$\pm\widehat{\Delta t}$ using Eq. (B5) (or from the height of the secondary
peaks at ${\rm ACF}(H;\pm\widehat{\Delta t})$). But this won’t be reliable for
the realistic QSO light curves with unknown correlation and especially in the
presence of significant observational noise..

Figure 2: The top-left panel shows the joint light curve (in perfect condition
with negligible observation noise) of the double system considered in Fig. 1:
$\widehat{\Delta t}=22.0$ days and $\widehat{\mu}=0.86$. The difference series
for this system, $H(t)$ as defined in Eq. (8) (which is simply the derivative
of $F(t)$ when it is sampled uniformly), is shown in the bottom-left panel.
Its auto-correlation function, ${\rm ACF}(H;\Delta t)$ shown in the top-right
panel, exhibits a prominent and sharp pair of secondary maxima at $\Delta
t=\pm\widehat{\Delta t}$ through which one can easily detect the lensed
system. Finally, the bottom-right panel compares the full fluctuation curve
($\epsilon(\Delta t)$, same as in the left panel of Fig. 1) using the trial
$\mu=0.3$ with the contribution from the ${\rm ACF}(H;\Delta t)$ to it.

Fig. 2 demonstrates the above arguments for the same example considered in
Fig. 1. The top-left panel shows the joint light curve (in perfect conditions
with negligible noise and one day of cadence) of the doubly lensed system
simulated using the damped random walk (DRW) template with $\widehat{\Delta
t}=22.0$ days and $\widehat{\mu}=0.86$. The derivative of the joint light
curve, $H(t)$, is shown in the bottom-left panel. The top-right panel displays
the auto-correlation function of $H(t)$ where we can clearly find the pair of
secondary maxima at $\Delta t=\pm\widehat{\Delta t}=\pm 22.0$ days. Finally,
we compare the full fluctuation curve ($\epsilon(\Delta t)$) with its main
contributing term (linear in $\mu$) stemming from ${\rm ACF}(H;\Delta t)$ for
an arbitrary trial $\mu=0.3$. We notice that the features in the fluctuation
curve are dominated by the ${\rm ACF}(H;\Delta t)$, which is also responsible
for the secondary minima in $\epsilon(\Delta t)$ appearing at $\Delta
t=\pm\widehat{\Delta t}$. The system is therefore identified as lensed. Note
that the auto-correlation of the joint light curve itself, i.e. ${\rm
ACF}(F;\Delta t)$, does not show the lensing peaks in this case; hence it is
not reliable for lens detection as explained in Section 4 in more detail.

###  3.1 Contribution from the higher order terms

For an arbitrarily long time series ($t_{\rm range}\to\infty$), one can keep
substituting

$\sum_{i}^{N_{D}}H(t+t_{c})H(t+t_{c}-n\Delta t)\approx{\rm ACF}(H;n\Delta
t)\cdot\epsilon_{0}\;,$ (27)

for any $n$ ignoring the boundary effect. Here $t_{c}$ is any constant shift
in time and $\epsilon_{0}$ is given by (12). Therefore, Eqs. (13) – (16) can
be simplified as

$\displaystyle\epsilon_{1}(\Delta t)$ $\displaystyle=-\epsilon_{0}\left[2{\rm
ACF}(H;\Delta t)\right]\;,$ (28) $\displaystyle\epsilon_{2}(\Delta t)$
$\displaystyle=\epsilon_{0}\left[2{\rm ACF}(H;2\Delta t)+1\right]\;,$
$\displaystyle\epsilon_{3}(\Delta t)$ $\displaystyle=-\epsilon_{0}\left[2{\rm
ACF}(H;3\Delta t)+2{\rm ACF}(H;\Delta t)\right]\;,$
$\displaystyle\epsilon_{4}(\Delta t)$ $\displaystyle=\epsilon_{0}\left[2{\rm
ACF}(H;4\Delta t)+2{\rm ACF}(H;2\Delta t)+1\right]\;,$ $\displaystyle\ldots$
$\displaystyle\epsilon_{n}(\Delta t)$
$\displaystyle=\begin{cases}\left(-1\right)^{n}\epsilon_{0}\left[2{\rm
ACF}(H;n\Delta t)+2{\rm ACF}(H;(n-2)\Delta t)+\ldots+2{\rm ACF}(H;\Delta
t)\right]&\text{if $n$ is odd}\;.\\\
\left(-1\right)^{n}\epsilon_{0}\left[2{\rm ACF}(H;n\Delta t)+2{\rm
ACF}(H;(n-2)\Delta t)+\ldots+2{\rm ACF}(H;2\Delta t)+1\right]&\text{if $n$ is
even}\;.\end{cases}$ (29)

Note that ${\rm ACF}(H;n\Delta t)$ exhibits a pair of peaks at $\Delta
t=\pm\widehat{\Delta t}/n$. Thus we get pairs of local minima (maxima) in the
fluctuation curve at $\pm\widehat{\Delta t}/n$ for each odd (even) $n$ from
Eq. (11). We refer to these features as the lensing signal in the fluctuation
curve. Remarkably, from each odd term, we get a pair of minima at $\Delta
t=\pm\widehat{\Delta t}$ which, despite being suppressed by a factor of
$\mu^{n}$, enhances our target pair of secondary minima and increases its
detectability. On the other hand, we find local maxima in the fluctuation
curve at $\Delta t=\pm\widehat{\Delta t}/2,\widehat{\Delta t}/4,\ldots$ from
the even terms. Therefore, as more and more terms contribute to the
$\epsilon(\Delta t)$ curve, we get slight enhancement in the lensing minima
(as compared to the lensing peaks in ${\rm ACF}(H;\Delta t)$) but at the
expense of more extrema that increase the overall oscillations in the
fluctuation curve. Note that for a finite time series, we still expect similar
higher-order features in the fluctuation curve up to the order $n\ll t_{\rm
range}/|\widehat{\Delta t}|$.

###  3.2 Quad systems

Figure 3: The top panel shows the auto-correlation function of the derivative
of the joint light curve ($H(t)$, defined in Eq. (8)) for a quad system (4
images in reality). The bottom panel shows the fluctuation curve for the same
system, analysed assuming just two images. The six dashed vertical lines on
either side of $\Delta t=0$ in each panel represent the six relative time
delays between the four actual images.

It is possible to generalise the image reconstruction, Eq. (4), for more than
two image systems. However, for a generic system with $N_{I}$ images, ${\rm
ACF}(H;\Delta t)$ exhibits $N_{I}(N_{I}-1)/2$ pairs of secondary maxima, one
for each relative time delay $|T_{i}-T_{j}|$ ($i\neq j$). Note that if the
time delays between two pairs of images are too close to each other, or more
explicitly, if the difference between two time delays is smaller than the
observation cadence or the smoothing time scale (used to deal with noisy
data), then the two corresponding peaks will merge into one in ${\rm
ACF}(H;\Delta t)$. For example, ACF$(H;\Delta t)$ for a quad system (four
images blended together in reality) can show up to six pairs of secondary
maxima at the time delays $\Delta t=\pm\Delta t_{i,j\neq i}$ given that
$\Delta t_{i,j\neq i}$ are well separated from each other. As ${\rm
ACF}(H;\Delta t)$ dominates the lensing signal in the fluctuation curve, the
latter also shows the same number of secondary minima pairs in this case even
if one assumes only two images in the reconstruction analysis, i.e. following
Eq. (4) with two images only. Thus, one can identify the quad systems using
the method presented in paper 1 (Bag et al., 2022) by detecting multiple (up
to six) pairs of minima in the $\epsilon(\Delta t)$ curve.

This has been illustrated in Fig. 3 where we consider an example of a quad
unresolved lensed system simulated using DRW process with time delays:
$\widehat{\Delta t}=11.0,31.0,68.0$ days with respect to the first image
(again in perfect condition). The top and bottom panels show the auto-
correlation function for the derivative, ACF$(H;\Delta t)$ and the fluctuation
curve $\epsilon(\Delta t)$ (using two image analysis), respectively. In both
panels, we find six pairs of prominent extrema at $\Delta t=\pm 11,~{}\pm
31,~{}\pm 68,~{}\pm 20,~{}\pm 57,~{}\pm 37$ days. Therefore, one can detect
the quad systems using both approaches. In a blind analysis, if we detect one
pair of secondary minima in the $\epsilon(\Delta t)$ (or secondary maxima in
${\rm ACF}(H;\Delta t)$) curve, we can identify the system as a doubly lensed
QSO. On the other hand, if we find multiple such pairs (up to six), we can
detect it as a quad system since 3-image systems do not exist in reality.

Thus, in this section, we explain the mathematical origin of the
characteristics of the fluctuation curves which play the pivotal role in
developing the method introduced in paper 1 (Bag et al., 2022). However, we
are left with two additional questions: (i) why can one detect the lensed QSOs
more reliably using ${\rm ACF}(H;\Delta t)$ as compared to ${\rm ACF}(F;\Delta
t)$, and (ii) what are the advantages of using the fluctuation curve over the
simpler ${\rm ACF}(H;\Delta t)$ in detecting lensed QSOs? These questions are
addressed below in Sections 4 and 5, respectively.

##  4 ${\rm ACF}(F;\Delta t)$ vs ${\rm ACF}(H;\Delta t)$

Figure 4: The top-left panel shows auto-correlation in the intrinsic QSO light
curve, ${\rm ACF}(f;\Delta t)$, for 5 out of 1000 realisations (randomly
chosen) simulated using the DRW template. The expectation value (ensemble
average) and $68\%$ percentile are marked by the dashed black curve and the
shaded region respectively. The top-right panel shows ${\rm ACF}(h;\Delta t)$
for the same five realisations along with the expected value (black dashed
curve). The bottom two panels compare ${\rm ACF}(F;\Delta t)$ and ${\rm
ACF}(H;\Delta t)$ for the double systems constructed using these five
realisations. The dashed vertical lines in the bottom two panels represent the
true time delay: $\Delta t=\pm 50$ days for these systems. Since ${\rm
ACF}(h;\Delta t)$ is typically much narrower than ${\rm ACF}(f;\Delta t)$, we
expect to find the lensing peaks in ${\rm ACF}(H;\Delta t)$ more reliably than
in ${\rm ACF}(F;\Delta t)$, as evident from the comparison of the bottom two
panels.

One can in principle reliably detect the lensed systems using the auto-
correlation function of the joint light curve, ${\rm ACF}(F;\Delta t)$, if the
intrinsic light curve $f(t)$ is uncorrelated in time (white noise) as
explained in Appendix A. However, the QSO light curves can have long time
scale correlations that violate Eq. (A4). In this case, the existence of the
lensing peaks in ${\rm ACF}(F;\Delta t)$ depends on the characteristics of the
auto-correlation function of the intrinsic light curve, ${\rm ACF}(f;\Delta
t)$. Even when Eq. (A4) is not strictly valid, we can expect excess power in
${\rm ACF}(F;\Delta t)$ at $\Delta t=\pm\widehat{\Delta t}$ from Eq. (A3) if
${\rm ACF}(f;\Delta t)$ decays reasonably sharply (from its peak at $\Delta
t=0$), i.e. if $f(t)$ and $f(t\pm\Delta t)$ are uncorrelated for $\lvert\Delta
t\rvert\gtrsim\lvert\widehat{\Delta t}\rvert$. In other words, one can detect
the lensing peaks in ${\rm ACF}(F;\Delta t)$ only if ${\rm ACF}(f;\Delta t)$
is narrowly peaked around $\Delta t=0$ even in the ideal condition with
negligible observation noise. This is discussed in Appendix A.1 in detail (see
Eq. (A7) for the precise condition). Since $F(t)$ and its derivative $H(t)$
follow the same lensing equation (compare Eqs. (2) and (18)), the above
criterion is also applicable to ${\rm ACF}(h;\Delta t)$ in order to find the
lensing peaks in ${\rm ACF}(H;\Delta t)$.

However, QSO flux variability typically shows temporal correlation till a few
hundred to even thousand days (Kelly et al., 2009; MacLeod et al., 2010). The
expectation value of the auto-correlation function of light curves generated
using damped random walk (DRW) process (Kelly et al., 2009; MacLeod et al.,
2010; Zu et al., 2013) decays exponentially,

$\langle{\rm ACF}(f;\Delta t)\rangle_{E}=\exp(-\lvert\Delta t\rvert/\tau)\;,$
(30)

but not sharply since the decay time scale, $\tau$, is typically
$\mathcal{O}(10^{2}-10^{3})$ days. Here, $\langle\cdot\rangle_{E}$ denotes an
ensemble average over all possible realizations of $f(t)$ with the same
stochastic properties.

Therefore, it can be difficult to detect the secondary peaks in ${\rm
ACF}(F;\Delta t)$, and one can also have a significant number of false
positive detections (Geiger & Schneider, 1996; Shu et al., 2021). This has
been demonstrated in Fig. 4 using numerical simulations. Here we simulate
$1000$ realisations of the intrinsic QSO light curves ($f(t)$) using the DRW
template. We fix the correlation time scale to $\tau=10^{2.5}\approx 316$ days
which is consistent with the findings of Kelly et al. (2009); MacLeod et al.
(2010) (see also Dobler et al. (2015)) throughout the paper for our
illustration purpose. Then we construct the joint light curves ($F(t)$)
separately for each realisation following Eq. (2) with $\widehat{\mu}=0.86$
and $\widehat{\Delta t}=50$ days, kept same across the realisations. For
simplicity, we consider the perfect condition with marginal noise in the data.
The solid curves in the top-left panel show ${\rm ACF}(f;\Delta t)$ for five
randomly selected samples, whereas the dashed back curve and the shaded region
represent the ensemble average of ${\rm ACF}(f;\Delta t)$ and the $68\%$
quantile, respectively 555Eq. (30) is true only if the observation range is
much larger than correlation scale, i.e. $t_{\rm range}\gg\tau$ ($\tau$ has
been set to $10^{2.5}$ days for these simulations) so that Ergodicity is
observed. So the black dashed curve in the top-left panel of Fig. 4 coincides
with Eq. (30) for much longer observation range.. It is evident that there
exists significant correlation in $f(t)$ till a few hundred days as the ${\rm
ACF}(f;\Delta t)$ curves decay slowly with $|\Delta t|$. Also, notice that the
$68\%$ quantile region of ${\rm ACF}(f;\Delta t)$ expands with $\lvert\Delta
t\rvert$, thus some realisations of ${\rm ACF}(f;\Delta t)$ can exhibit local
maxima as in the case for the red and purple curves in the top-left panel.
Thus, using ${\rm ACF}(F;\Delta t)$ one can get a substantial number of false
lensed detections in true unlensed cases because of these maxima.

The ${\rm ACF}(F;\Delta t)$ for the same five realisations have been shown in
the bottom-left panel, the dashed vertical lines represent the true time
delay, $\Delta t=\pm 50$ days, for these systems. Although for some
realisations (purple curve) one can find excess power in ${\rm ACF}(F;\Delta
t)$ at $\Delta t=\pm\widehat{\Delta t}$, for others (the green, red and blue
curves) this is not true.

On the other hand, since DRW behaves like a random walk at small time scales
($t\ll\tau$), its derivative $h(t)$ behaves like white noise obeying Eq. (A4)
at this limit. This is evident from the top-right panel of Fig. 4 where we
show ${\rm ACF}(h;\Delta t)$ for the same five realisations with different
colours. For all realisations we find that ${\rm ACF}(h;\Delta t)$ becomes
very small for $\Delta t\neq 0$. In fact the ensemble average, shown by the
dashed black curve, follows $\langle{\rm ACF}(h;\Delta
t)\rangle\approx\delta_{\Delta t,0}$. Hence, we get a prominent pair of
secondary maxima in ${\rm ACF}(H;\Delta t)$ at $\Delta t=\pm\widehat{\Delta
t}$ for all the realisations as evident from the bottom-right panel.

Lastly, let us consider time series with different correlation time scales
$\tau$ in Eq. (30), even if they do not describe the QSO light curves
accurately. A smaller (larger) $\tau$ leads to a narrower (broader) peak in
$\langle{\rm ACF}(f;\Delta t)\rangle$, which in turn increases (decreases) the
probability of ${\rm ACF}(F;\Delta t)$ showing the lensing peaks, according to
Appendix A.1. However, in both limits, ${\rm ACF}(h;\Delta t)$ has a sharp
peak at $\Delta t=0$ so that one can always find the lensing peaks in ${\rm
ACF}(H;\Delta t)$.

This exercise using the DRW template thus demonstrates that ${\rm
ACF}(H;\Delta t)$ outperforms ${\rm ACF}(F;\Delta t)$ in terms of
detectability of the lensing peaks. However, this conclusion is not restricted
to DRW and unbound Random walk templates. It stands valid for any ‘red-type’
power spectrum, as argued below in Section 4.1 and in Appendix C more
explicitly.

###  4.1 Connection to the power spectrum

Let us define the power spectrum, $P_{f}(\omega)$, of a time series as the
two-point correlation function in the Fourier space;

$\langle\tilde{f}(\omega)\tilde{f}(\omega^{\prime})^{*}\rangle_{E}=P_{f}(\omega)\delta(\omega-\omega^{\prime})\;,$
(31)

where $\tilde{f}(\omega)$ is the Fourier transform of the time series $f(t)$.
We assume that the Fourier modes are uncorrelated. To be precise, we assume
that the QSO intrinsic light curve $f(t)$ is ‘wide-sense stationary’ as its
mean and covariance properties do not vary over time (see (20) – (22)). Under
these assumptions, the Wiener-Khinchin theorem states that the (expected)
auto-correlation function of $f(t)$ is given by the Fourier transform of the
power spectrum (Wiener, 1930; Khintchine, 1934; Einstein, 1914). Therefore,
one can determine the auto-correlation function of a time series by studying
its power spectrum. This is especially useful since a derivative in the time
domain corresponds to a multiplication by $i\omega$ in the Fourier domain, and
hence $P_{h}(\omega)=\omega^{2}P_{f}(\omega)$.

The simplest example is when the underlying signal consists purely of white
noise: $f_{\textrm{white}}(t)$ with the power spectrum given by
${P_{f_{\textrm{white}}}(\omega)=\sigma_{\textrm{white}}^{2}}=\text{constant}$.
Taking the Fourier transform yields the expected value of the auto-correlation
function ACF$(f_{\textrm{white}};\Delta t)\propto\delta(\Delta t)$. Since this
is sharply peaked around $\Delta t=0$, we can accurately retrieve the lensing
peaks in the ${\rm ACF}(F;\Delta t)$ curve. Further details for this white
noise scenario are discussed in Appendix A.

For the damped random walk (DRW) templates often used for describing the
intrinsic QSO light curves, the power spectrum takes the form

${P_{f_{\textrm{DRW}}}(\omega)\propto\tau/(1+\tau^{2}\omega^{2})}\;,$ (32)

which is the Fourier transform of Eq. (30) up to a normalisation factor. Note
that power spectra like $P_{f_{\textrm{DRW}}}(\omega)$ which decay with
$|\omega|$ give rise to red noise and hence can be classified as ‘red-type’
power spectra. The decay time scale $\tau$ parameterises how long it takes for
the time series to forget the fluctuations that happened in the past. Since
$P_{f_{\mathrm{DRW}}}(\omega)\approx$ constant for $\tau^{2}\omega^{2}\ll 1$,
the DRW is mostly independent of its past and behaves like a white noise when
$\tau$ is small.

To describe realistic quasar light curves, $\tau$ should typically be of the
order $10^{2}$ – $10^{3}$ days. The DRW then behaves more similarly to an
(undamped) Gaussian random walk with a power spectrum
$P_{f}(\omega)\sim\omega^{-2}$ up to a time scale smaller than $\tau$, or
$\omega^{2}\gg 1/\tau^{2}$. The Fourier transform of such a power spectrum has
a broad peak at zero unlike a sharp delta-function-like peak of the white
noise one. Therefore, ${\rm ACF}(F;\Delta t)$ is less reliable for finding
lensed cases in this case.

Meanwhile, the derivative time series $h(t)$ has a much flatter power spectrum
with $P_{h}(\omega)\propto\omega^{2}/(1+\tau^{2}\omega^{2})$. This is nearly
constant for large $\omega$; $h(t)$ behaves like white noise. Thus, ${\rm
ACF}(h;\Delta t)$ tends to be much narrower and is expected to decay down
quickly as compared to ${\rm ACF}(f;\Delta t)$. This is consistent with our
previous findings that ${\rm ACF}(H;\Delta t)$ is more reliable than ${\rm
ACF}(F;\Delta t)$ for detecting the lensed systems. Indeed, a larger $\tau$
makes the $P_{f}(\omega)$ steeper (or redder) in Eq. (32) and leads to a
flatter $\langle{\rm ACF}(f;\Delta t)\rangle$ in Eq. (30). This in turn
reduces the possibility of the lensing peaks appearing in ${\rm ACF}(F;\Delta
t)$ and thus enables ${\rm ACF}(H;\Delta t)$ to outperform ${\rm ACF}(F;\Delta
t)$ even to a greater extent. In fact, this is true for all red type power
spectra ($P_{f}(\omega)$ decreases with increasing $|\omega|$) as illustrated
in Appendix C. Since the lensing features in the fluctuation curve
$\epsilon(\Delta t)$ predominantly arise from ${\rm ACF}(H;\Delta t)$ we
expect similar performance from both these approaches.

##  5 Fluctuation in reconstruction vs auto-correlation of derivative

Figure 5: The fluctuation curve and ${\rm ACF}(H;\Delta t)$ are compared after
being normalized using Eq. (33) for the same system as in Fig. 1 (i.e.
simulated with DRW template and with $\widehat{\Delta t}=22.0$ days and
$\widehat{\mu}=0.86$). We find secondary extrema of similar strengths in both
curves as the ${\rm ACF}(H;\Delta t)$ dominates the lensing signal in the
fluctuation curve. However, the latter still has more lensing signal due to
contributions from the higher order terms in Eq. (11) as illustrated by the
two inset plots which zoom into the two lensing/secondary minima.

Since the lensing signal in the fluctuation curve is dominated by the auto-
correlation of the derivative of the joint light curve ($H(t)$ defined in Eq.
(8)), it is interesting to test if ${\rm ACF}(H;\Delta t)$ can similarly be
used to detect the lensed QSOs. Note that ${\rm ACF}(H;\Delta t)$ would be
less computationally expensive and has a better physical interpretation.

Let us begin by comparing the lensing signals (the prominence of the secondary
maxima that can be used for lensed detection) in ${\rm ACF}(H;\Delta t)$ and
the fluctuation curve ($\epsilon(\Delta t)$) for the example presented in
Figs. 1 (left panel) and 2. We normalize both curves using the generic
transformation,

$\Sigma_{X(\Delta t)}\equiv\frac{X(\Delta t)-\langle X(\Delta
t)\rangle_{\Delta t}}{\sigma_{X(\Delta t)}}\;,$ (33)

where $\langle X(\Delta t)\rangle_{\Delta t}$ and $\sigma_{X(\Delta t)}$ are
the average and standard deviation taken over $\Delta t$ of a one-dimensional
function $X(\Delta t)$. (Thus $\Sigma_{X(\Delta t)}$ simply measures $X(\Delta
t)$ in the units of its standard deviation.) Fig. 5 compares
$\Sigma_{\epsilon(\Delta t)}$ (with trial $\mu=0.3$) and $-\Sigma_{{\rm
ACF}(H;\Delta t)}$ for that example system (simulated using DRW template with
$\widehat{\Delta t}=22.0$ days, $\widehat{\mu}=0.86$ and marginal observation
noise). As expected, both curves show prominent pair of secondary minima at
$\Delta t=\pm\widehat{\Delta t}=\pm 22.0$ days (marked by dashed vertical
lines) deeper than $\Sigma=-2$. However, we still find that the secondary
minima in the $\Sigma_{\epsilon(\Delta t)}$ are slightly deeper than that of
the $-\Sigma_{{\rm ACF}(H;\Delta t)}$ curve. This indicates that
$\epsilon(\Delta t)$ contains a more enhanced lensing signal than ${\rm
ACF}(H;\Delta t)$ in this case.

Figure 6: The normalised fluctuation curve (blue curve, from Bag et al.
(2022)) and ${\rm ACF}(H;\Delta t)$ (shown by the orange curve) are compared
for the COSMOGRAIL system SDSS J1226-0006.

Next, we test the approach based on ${\rm ACF}(H;\Delta t)$ on the same
validation and blind sets used in Bag et al. (2022) (simulated using the DRW
template). As usual, we first consider the perfect conditions where the
observational noise is marginal compared to the time variation in the
intrinsic QSO light curves. Using the same selection criteria, ${\rm
ACF}(H;\Delta t)$ can detect 19 out of 20 true lensed cases correctly. Among
the 20 unlensed cases, however, it gives one false positive case while
identifying the rest 19 unlensed systems correctly. In comparison, the
minimization of the fluctuation (in the reconstructed image light curves)
approach produces 1 false negative case but zero false positive cases for the
same data sets. We also notice that the signal in $\Sigma_{{\rm ACF}(H;\Delta
t)}$ is slightly diminished as compared to $\Sigma_{\epsilon(\Delta t)}$.

When we consider uncertainty in the observed light curve data, the ${\rm
ACF}(H;\Delta t)$ approach suffers more from the added noise than the
technique based on the fluctuation curve. We follow the same prescription
given in Bag et al. (2022) for both methods to handle noisy data; we smooth
the joint light curve for multiple smoothing scales and combine the
fluctuation curves obtained from each smoothed light curve. For the same
datasets with ZTF-like noise, ${\rm ACF}(H;\Delta t)$ approach detects only 8
out of the 20 lensed systems correctly. However, for 2 other lensed cases, it
detects the lensing nature based on peaks at completely wrong time delays
(hence these two should be counted as false positives). Furthermore, it
detects 15 out of 20 unlensed systems correctly but gives the rest 5 as false
positives. Thus, in combination, it produces a precision of $8/15$ (slightly
higher than $50\%$) and a recall of $8/20$. These numbers are significantly
inferior to that of the fluctuation curve method which produces a precision of
$12/13$ and a recall of $12/20$ on the same datasets (Bag et al., 2022). The
recall and precision values for these two approaches have been summarized in
Table 1. As before, we again notice that the lensing signal in $\Sigma_{{\rm
ACF}(H;\Delta t)}$ is typically diminished as compared to
$\Sigma_{\epsilon(\Delta t)}$.

Datasets from | Recall | Precision
---|---|---
${\rm ACF}(H;\Delta t)$ | Fluctuation curve | ${\rm ACF}(H;\Delta t)$ | Fluctuation curve
(Bag et al., 2022) |  | (Bag et al., 2022) |  | (Bag et al., 2022)
With marginal noise | $95\%$ | $95\%$ | $95\%$ | $100\%$
With ZTF-like noise | $40\%$ | $60\%$ | $53.33\%$ | $92.31\%$
Table 1: We compare the results from the two approaches – ${\rm ACF}(H;\Delta
t)$ vs the full fluctuation curve – in terms of recall (completeness) and
precision (purity) for the validation and blind sets considered in Bag et al.
(2022). The results corresponding to marginal noise and ZTF-like noise in the
data are presented in the top and bottom rows respectively.

We also test ${\rm ACF}(H;\Delta t)$ on the COSMOGRAIL system SDSS J1226-0006
which was used in Bag et al. (2022) as an example. The time delay estimated
using the resolved light curves by the COSMOGRAIL team is $33.7\pm 2.7$ day
for this system (Millon et al., 2020). Fig. 6 compares the fluctuation curve
(blue curve, same as in Fig. 21 of Bag et al. (2022)) and ${\rm ACF}(H;\Delta
t)$ (orange curve) after normalisation. Like the fluctuation curve,
$-\Sigma_{{\rm ACF}(H;\Delta t)}$ shows a pair of prominent secondary minima
symmetrically placed around $\Delta t=0$. Specifically, the secondary minima
in $-\Sigma_{{\rm ACF}(H;\Delta t)}$ occurs at $\Delta t=-28.9,~{}26.9$ days
with the depths $\Sigma=-1.40$ and $-1.50$ respectively leading to the final
time delay estimation of $\Delta t_{\rm est}=27.9$ days. In comparison,
$\Sigma_{\epsilon(\Delta t)}$ curve exhibits slightly overall deeper minima
($\Sigma=-1.35,-1.68$) at $\Delta t=-28.7,~{}30.5$ days that give rise to
$\Delta t_{\rm est}=29.6$ days (Bag et al., 2022) which is a better agreement
with the COSMOGRAIL results.

In conclusion, it is evident that the fluctuation curve approach as introduced
by Bag et al. (2022) performs better than ${\rm ACF}(H;\Delta t)$ although the
former method is dominated by the latter. This is because the secondary minima
(i.e. the lensing signal) are more prominent in the fluctuation curves,
$\epsilon(\Delta t)$, than the secondary maxima in the corresponding ${\rm
ACF}(H;\Delta t)$ curves. This is due to the fact that in $\epsilon(\Delta t)$
curves, the secondary minima at $\Delta t=\pm\widehat{\Delta t}$ are further
enhanced by all the odd terms in Eq. (11) as compared to ${\rm ACF}(H;\Delta
t)$ curves. Therefore, the minimization of fluctuation approach can be more
useful for marginal detection of the unresolved lensed QSOs. However, ${\rm
ACF}(H;\Delta t)$ can be very useful and quick crosscheck as it is
computationally inexpensive.

##  6 Conclusion

Bag et al. (2022) introduces a data-driven technique for detecting lensed QSOs
and for measuring their time delays only using the unresolved joint light
curve data by minimizing the fluctuations in the reconstructed image light
curves. In this article, we provide the proof as to how this method works. We
showed that the lensing signal in the simple fluctuation estimator given by
Eq. (5) is dominated by the auto-correlation of the derivative (difference
series in general for non-uniformly sampled data) of the joint light curve.
This observation explains all the characteristics of the fluctuation curve
that were used as the foundation of the technique proposed by Bag et al.
(2022). Above all, ${\rm ACF}(H;\Delta t)$ is locally maximized at $\Delta
t=\pm\widehat{\Delta t}$ and these lensing peaks manifest themselves as the
secondary minima in the fluctuation curve, $\epsilon(\Delta t)$, which have
been used to detect the lensed cases in Bag et al. (2022). Other interesting
results are summarized below.

  * •

We also showed that ${\rm ACF}(H;\Delta t)$ is more reliable than the auto-
correlation function of the joint light curve itself ${\rm ACF}(F;\Delta t)$,
because the intrinsic flux variability of QSOs is correlated in the time
domain, or in other words, the power spectra of the intrinsic quasar light
curves are of red type. Nevertheless, even if $P_{f}(w)$ is flat or of blue
type, ${\rm ACF}(H;\Delta t)$ can also find the lensed cases. Since the
primary contribution to the lensing signal in $\epsilon(\Delta t)$ comes from
the ${\rm ACF}(H;\Delta t)$, the minimization of the fluctuation approach
would be similarly successful in these scenarios.

  * •

However, the approach based on the fluctuation curve provides better recall
and precision over the ${\rm ACF}(H;\Delta t)$ when one considers significant
amount of noise in the joint light curve data. This is due to the higher-order
terms contributing in Eq. (11) that further enhance the pair of secondary
minima in the $\epsilon(\Delta t)$ curve (as compared to the lensing peaks in
${\rm ACF}(H;\Delta t)$).

  * •

For a generic lensed system having $N_{I}$ images, ${\rm ACF}(H;\Delta t)$
displays $N_{I}(N_{I}-1)/2$ pairs of lensing local maxima. Likewise, the
fluctuation minimization method can be used to detect multiple imaged lensed
systems even if we assume two images in the reconstruction analysis a priori.
For example, by detecting one pair of prominent secondary minima in
$\epsilon(\Delta t)$ one can identify a double system, whereas if there exist
multiple such lensing minima pairs (up to six), the system must be a quad
(having 4 images).

Although we choose time delays of the order of tens of days as examples in
this work for demonstration purposes, one can in principle detect the lensing
minima pair in the fluctuation curve for any arbitrarily small time delay as
long as it is sufficiently larger than the observation cadence and the signal
to noise ratio is sufficiently high. However, in reality, the cadence could
vary from a few days to $\mathcal{O}(10)$ days depending upon the observation
conditions or the observing strategy; this puts a limitation on the
sensitivity of the method as the time delay needs to be larger than the
cadence for the lensing minima pair to emerge in the fluctuation curve.

The fact that the primary contribution to the lensing signal in the
fluctuation curves stems from ${\rm ACF}(H;\Delta t)$ also informs us about
some key benefits of the fluctuation minimization approach. The fluctuation
minimization approach should be able to handle microlensing up to a certain
limit as the auto-correlation can withstand moderate microlensing effects. We
plan to comprehensively investigate the effect of microlensing on the
performance of this method in the follow up work. Note that microlensing can
significantly alter the time delay measurements from the resolved image light
curves, up to a few days (Tie & Kochanek, 2018; Liao, 2020). It would be
interesting to see how this affects the results of our method based on the
unresolved fluxes.

To discern another crucial advantage, recall that the selection criteria for
lens detection using the fluctuation curve are so far based on only the pair
of minima at $\Delta t=\pm\widehat{\Delta t}$. The higher order terms further
put predictable features in the fluctuation curve at certain values of $\Delta
t$, e.g. a pair of minima (maxima) at $\Delta t=\pm\widehat{\Delta t}/n$ for
every odd (even) $n$. Although the higher order features are suppressed by
$\mu^{n}$, the first few of these features can nevertheless be exploited to
improve the selection criteria, potentially using deep learning.

Finally, we emphasise that the enhancement of fluctuations in the image light
curves reconstructed using wrong time delays is a fundamental characteristic
of the fluctuation curve approach. However, this article is restricted to the
simple metric Eq. (5) for quantifying the fluctuations. The lensing signal in
this estimator $\epsilon(\Delta t)$ is found to be dominated simply by ${\rm
ACF}(H;\Delta t)$. Nevertheless, there might be a better metric for
quantifying the fluctuations in the reconstructions that delivers better
results in terms of recall and precision. This remains an open question and is
worth investigating further.

As stressed out in Bag et al. (2022), to estimate the error in the time delay
measurements in this non-parametric approach one needs to statistically
analyse a large number of unresolved cases simulated in a variety of
observational conditions; e.g. considering different flux variations, cadence
distributions, noise levels, many microlensing realisations etc. This exercise
forms the focus of the follow up work.

## Acknowledgement

We thank Tommaso Treu for his crucial inputs to this work. SB also thanks Eric
V. Linder and Alex G. Kim for useful discussions. The Seondeok cluster at
Korea Astronomy and Space Science Institute has been used to carry out a part
of the analysis. A.S. would like to acknowledge the support by the National
Research Foundation of Korea NRF-2021M3F7A1082053 and the support of the Korea
Institute for Advanced Study (KIAS) grant funded by the government of Korea.
KL was supported by the National Natural Science Foundation of China (NSFC)
under Grant Nos. 12222302, 11973034 and Wuhan University talent research
start-up funds. A.S. and K.L. also acknowledge the support and hospitality
received from Beijing Normal University.

##  Appendix A Auto-correlation function

The auto-correlation function of a generic time series $x(t)$ is defined as

${\rm ACF}(x;\Delta t)=\frac{\sum_{i}^{N_{D}}\left[x(t_{i})-\langle
x(t)\rangle_{t}\right]\left[x(t_{i}-\Delta t)-\langle x(t-\Delta
t)\rangle_{t}\right]}{\sqrt{\sum_{i}^{N_{D}}\left[x(t_{i})-\langle
x(t)\rangle_{t}\right]^{2}\cdot\sum_{i}^{N_{D}}\left[x(t_{i}-\Delta t)-\langle
x(t-\Delta t)\rangle_{t}\right]^{2}}}\;.$ (A1)

For a wide-sense stationary (bounded and long) time series $\langle
x(t)\rangle_{t}\approx\langle x(t-\Delta t)\rangle_{t}$ which reduces the
denominator to $N_{D}\sigma^{2}_{x}$. Let us define the joint time series
$X(t)$ following the lensing equation (similar to $F(t)$ in Eq. (2) and $H(t)$
in Eq. (18))

$X(t)=x(t)+\widehat{\mu}x(t-\widehat{\Delta t})\;,$ (A2)

where $\widehat{\mu}$ and $\widehat{\Delta t}$ represent the (true)
magnification ratio and relative time delay as in Eq. (2) and (18). Thus,
$\\{x(t),X(t)\\}$ here are proxies for $\\{f(t),F(t)\\}$ or their derivatives
$\\{h(t),H(t)\\}$. Using Eq. (A1) it is easy to find that

${\rm ACF}(X;\Delta t)=\frac{{\rm ACF}(x;\Delta
t)+\frac{\widehat{\mu}}{1+\widehat{\mu}^{2}}\left[{\rm ACF}(x;\Delta
t-\widehat{\Delta t})+{\rm ACF}(x;\Delta t+\widehat{\Delta
t})\right]}{1+\frac{2\widehat{\mu}}{1+\widehat{\mu}^{2}}{\rm
ACF}(x;\widehat{\Delta t})}\;,$ (A3)

where the denominator is just a normalisation constant (independent of $\Delta
t$).

If the intrinsic time series $x(t)$ is uncorrelated in time (white noise),
sufficiently long and hence obeys

$\langle x(t)x(t-\Delta t)\rangle_{t}=\begin{cases}\langle
x(t)^{2}\rangle_{t}&\text{if }\Delta t=0\;,\\\ \langle
x(t)\rangle_{t}^{2}&\text{if }\Delta t\neq
0\;,\end{cases}\text{~{}therefore,~{}}{\rm ACF}(x;\Delta
t)=\begin{cases}1&\text{if }\Delta t=0\;,\\\ 0&\text{for all }\Delta t\neq
0\;,\end{cases}$ (A4)

the denominator of Eq. (A3) reduces to unity (in the generic cases with
$\widehat{\Delta t}\neq 0$) leading to

${\rm ACF}(X;\Delta t)={\rm ACF}(x;\Delta
t)+\frac{\widehat{\mu}}{1+\widehat{\mu}^{2}}\left[{\rm ACF}(x;\Delta
t-\widehat{\Delta t})+{\rm ACF}(x;\Delta t+\widehat{\Delta t})\right]\;,$ (A5)

and we can conclude the followings.

  * •

When $\Delta t=0$, Eq. (A5) trivially reduces to unity as only the first term
contributes. Naturally, the auto-correlation is always maximized at unity for
no shift in time ($\Delta t=0$).

  * •

On the other hand, for $\Delta t\neq 0$ in general, Eq. (A4) ensures that all
three terms in Eq. (A5) reduce to zero and ${\rm ACF}(X;\Delta t)\approx 0$.

  * •

Interestingly, ${\rm ACF}(H;\Delta t)$ is locally maximized at $\Delta
t=\pm\widehat{\Delta t}$. In the view of Eq. (A4), when $\Delta
t=\widehat{\Delta t}$ or $-\widehat{\Delta t}$ only the first or the second
term in the square bracket of Eq. (A5) contributes, that leads to

${\rm ACF}(X;\Delta t=\pm\widehat{\Delta
t})\approx\frac{\widehat{\mu}}{1+\widehat{\mu}^{2}}\leq\frac{1}{2}$ (A6)

for both cases. In fact, the set of intrinsic features in the underlying time
series $x(t)$ appears twice in $X(t)$, separated by $\widehat{\Delta t}$ (and
scaled by $\widehat{\mu}$, see Eq. (A2)). Since ${\rm ACF}(X;\Delta t)$
measures the correlation between two copies of $X(t)$ with one shifted by
$\Delta t$ in the time domain, it gains an excess power when $\Delta
t=\pm\widehat{\Delta t}$ due to matching of the two sets of the same features
($\pm$ sign accounts for the shift in either direction). Although this
argument applies to any generic $x(t)$ with a substantial amount of features,
we emphasise that the above equation is valid strictly when the $x(t)$ is pure
white noise following Eq. (A4).

  * •

For unlensed cases, when $\widehat{\mu}=0$ or $\widehat{\Delta t}=0$, this
pair of secondary maxima vanishes.

In summary, ${\rm ACF}(X;\Delta t)$ shows a pair of secondary maxima at
$\Delta t=\pm\widehat{\Delta t}$ (the lensing peaks) and remain zero at all
other $\Delta t\neq 0$ as long as the intrinsic time series, $x(t)$, is white
noise (no temporal correlation) and satisfies Eq. (A4). Note that the whole
analysis is applicable to both $F(t)$ and $H(t)$ with the underlying functions
$f(t)$ and $h(t)$ respectively.

###  A.1 When to expect lensing peaks in ${\rm ACF}(F;\Delta t)$ even if
$f(t)$ is not white noise?

Let us discuss the interesting case when Eq. (A4) is not valid strictly. This
is important for quasars since the intrinsic light curves, $f(t)$, typically
possess correlations up to a time scale of a few hundred days (Kelly et al.,
2009; MacLeod et al., 2010). In such cases, the auto-correlation function
${\rm ACF}(f;\Delta t)$ would show a broad peak at $\Delta t=0$, unlike a
delta function. ${\rm ACF}(f;\Delta t)$ is still expected to be symmetric
around $\Delta t=0$ and monotonically decaying with $|\Delta t|$ (the decay
rate depends on the correlation time scale). For instance, see the top-left
panel of Fig. 4 for typical examples of ${\rm ACF}(f;\Delta t)$ where $f(t)$
is generated from damped random walk (DRW) with a correlation scale of $\sim
300$ days.

The lensing peaks in ${\rm ACF}(F;\Delta t)$ arise from the terms in the
square bracket in the numerator of Eq. (A3) (after replacing $\\{x(t),X(t)\\}$
by $\\{f(t),F(t)\\}$). Thus, as $|\Delta t|$ approaches $|\widehat{\Delta t}|$
from the below, the lensing peaks could emerge only if the change in the
second term in the numerator of Eq. (A3) dominates over the change in the
first term. For all practical purposes (i.e. with the expected ${\rm
ACF}(f;\Delta t)$ being symmetric around $\Delta t=0$ and monotonically
decreasing with $|\Delta t|$), this requirement boils down to the condition,

$\left|{\rm ACF}^{\prime}(f;\Delta t\sim\pm\widehat{\Delta
t})\right|<\frac{\widehat{\mu}}{1+\widehat{\mu}^{2}}\left|{\rm
ACF}^{\prime}(f;\Delta t\sim 0)\right|\;,$ (A7)

where ${\rm ACF}^{\prime}(f;\Delta t)$ is the derivative of ${\rm
ACF}(f;\Delta t)$ with respect to $\Delta t$. The stronger the inequality is,
the steeper (more prominent) the lensing maxima pair in ${\rm ACF}(F;\Delta
t)$ would be. For our purpose, the above condition requires that ${\rm
ACF}(f;\Delta t)$ must have sufficiently narrow peak at $\Delta t=0$ so that
it stabilizes at $\lvert\Delta t\rvert\gtrsim\lvert\widehat{\Delta t}\rvert$
by decaying down sufficiently fast. We again emphasise that all the arguments
made in this section stand valid if we replace $F(t)$ and $f(t)$ by $H(t)$ and
$h(t)$ respectively. Note that here we assumed the best-case scenario when the
noise is negligible; the inclusion of observation noise brings in additional
complexities.

##  Appendix B Exact values of the global and secondary minima in the
fluctuation curve

###  B.1 Global minima in the fluctuation curve

For $\Delta t=0$, Eq. (3) reduces to $f_{\rm rec}(t)=\frac{1}{1+\mu}F(t)$
which reduces Eq. (5) to

$\epsilon(\Delta
t=0)=\frac{1}{(1+\mu)^{2}}\sum_{i}^{N_{D}}\left[F(t_{i+1})-F(t_{i})\right]^{2}=\frac{\epsilon_{0}}{\left(1+\mu\right)^{2}}\;.$
(B1)

Note that the above expression is valid even for the unlensed cases.

One can get the same expression from Eq. 11 as follows,

$\displaystyle\epsilon(\Delta t=0)$
$\displaystyle=\epsilon_{0}\left[1-2\mu+3\mu^{2}-4\mu^{3}+5\mu^{4}-6\mu^{5}+\ldots\right]\;,$
$\displaystyle=\epsilon_{0}\frac{\partial}{\partial\mu}\left[\mu-\mu^{2}+\mu^{3}-\mu^{4}+\mu^{5}-\mu^{6}++\ldots\right]\;,$
$\displaystyle=\epsilon_{0}\frac{\partial}{\partial\mu}\left[\frac{\mu}{1+\mu}\right]=\frac{\epsilon_{0}}{\left(1+\mu\right)^{2}}\;.$
(B2)

Here we identify that for $\Delta t=0$ all $\epsilon_{n}$ terms become
proportional to $\epsilon_{0}$, and the sum converges as $\mu<1$.

###  B.2 Secondary minima when $h(t)$ is white noise

If $h(t)$ is pure white noise (or in other words, if it follows Eq. (A4)) we
can analytically calculate the height of the secondary minima. Following Eq.
(28), we get

$\displaystyle\epsilon(\Delta t=\pm\widehat{\Delta t})$
$\displaystyle=\epsilon_{0}\left[1-2{\rm ACF}(H;\pm\widehat{\Delta
t})\mu+\mu^{2}-2{\rm ACF}(H;\pm\widehat{\Delta t})\mu^{3}+\mu^{4}-2{\rm
ACF}(H;\pm\widehat{\Delta t})\mu^{5}+\mu^{6}-\ldots\right]\;,$ (B3)
$\displaystyle=\epsilon_{0}[\left(1+\mu^{2}+\mu^{4}+\mu^{6}+\ldots\right)-2{\rm
ACF}(H;\pm\widehat{\Delta t})\left(\mu+\mu^{3}+\mu^{5}+\ldots\right)]\;,$ (B4)
$\displaystyle=\epsilon_{0}\left[\frac{1-2\mu{\rm ACF}(H;\pm\widehat{\Delta
t})}{1-\mu^{2}}\right]\;,$ (B5)

where we use ${\rm ACF}(H;n\cdot\widehat{\Delta t})\approx 0$ for $n\neq
0,~{}\pm 1$. Note from Eq. (A6) that ${\rm ACF}(H;\pm\widehat{\Delta
t})\approx\widehat{\mu}/(1+\widehat{\mu}^{2})\leq 1/2$ in this limit. It is
also easily seen that for unlensed cases, when ${\rm ACF}(H;\Delta t\neq
0)=0$, $\epsilon(\Delta t)=\epsilon_{0}/(1-\mu^{2})=\text{constant}$ for any
$\Delta t\neq 0$.

By comparing with Eq. (B2), one can trivially show that

$\frac{\epsilon(\Delta t=\pm\widehat{\Delta t})}{\epsilon(\Delta
t=0)}=\left(\frac{1+\mu}{1-\mu}\right)\left(1-2\mu{\rm
ACF}(H;\pm\widehat{\Delta t})\right)\approx 1+2\mu\left(1-{\rm
ACF}(H;\pm\widehat{\Delta t})\right)+\mathcal{O}(\mu^{2})\;,$ (B6)

and owing to ${\rm ACF}(H;\pm\widehat{\Delta t})\leq 1/2$ it is readily seen
that $\epsilon(\Delta t=0)<\epsilon(\Delta t=\pm\widehat{\Delta t})$. Thus the
central minimum at $\Delta t=0$ will always be deeper than the pair of lensing
minima.

##  Appendix C ${\rm ACF}(H;\Delta t)$ is more reliable than ${\rm
ACF}(F;\Delta t)$ for any red power spectrum

Figure 7: The left panel shows the power spectra ($P_{f}(\omega)$ in log-
scale) given by Eq. C3 for three values of $\gamma=1.8,1.0$ and $1/2$ by the
blue, orange and green solid curves respectively. Here we set
$\omega_{0}=10^{-6}$ so that the power spectra behave like
$|\omega|^{-\gamma}$ for the most part. The dashed curves represent the power
spectra of the derivative:
$P_{h}(\omega)=\omega^{2}P_{f}(\omega)\approx|\omega|^{2-\gamma}$. Each of the
three panels on the right displays an example of the intrinsic flux $f(t)$ (in
linear scale and arbitrary unit) generated from these power spectra. Note that
observational noise is marginal here. The horizontal dashed line in each right
panel represents the mean of the light curve. Figure 8: Comparing ${\rm
ACF}(f;\Delta t)$, ${\rm ACF}(F;\Delta t)$, ${\rm ACF}(h;\Delta t)$ and ${\rm
ACF}(H;\Delta t)$ for five random realisations out of 1000 for three red type
power-law power spectra:
$P(\omega)=(1+\lvert\omega/\omega_{0}\rvert)^{-\gamma}$ with
$\gamma=1.8~{}\text{(left panels)},1~{}\text{(middle panels)}$ and
$1/2~{}\text{(right panels)}$. Here, we set $\omega_{0}=10^{-6}$ day-1 so that
$P_{f}(\omega)$ behaves like a power-law for most part of it.

We have seen in Section 4 that the quasar intrinsic flux variability generated
from a random damped walk (DRW) follows a red power spectrum given by Eq.
(32). We found that ${\rm ACF}(H;\Delta t)$ is more likely to exhibit the
lensing peaks as compared to ${\rm ACF}(F;\Delta t)$ in this case due to the
presence of temporal correlation in $f(t)$ up to a few hundred days
realistically. Moreover, we discussed that a larger correlation scale ($\tau$)
leads to a steeper (or redder) $P_{f}(\omega)$ in Eq. (32) and a flatter
(expected) ${\rm ACF}(f;\Delta t)$ in Eq. (30), which enables ${\rm
ACF}(H;\Delta t)$ to outperform ${\rm ACF}(F;\Delta t)$ even more profoundly.
In this appendix, we argue that this is not restricted to just the DRW
template; for any time series with a red-type power spectrum, ${\rm
ACF}(H;\Delta t)$ is more reliable than ${\rm ACF}(F;\Delta t)$ for finding
the lensed cases.

Recall that the expected auto-correlation function of a ‘wide-sense
stationary’ time series is given by the Fourier transform of the power
spectrum. Since $\langle{\rm ACF}(f;\Delta t)\rangle_{E}$ is symmetric around
its peak at $\Delta t=0$, the power spectrum should also be symmetric around
$\omega=0$; $P_{f}(\omega)$ must be a function of $|\omega|$ only. Recall
further from Eqs. (A3) and (A7) that ${\rm ACF}(F;\Delta t)$ can show the
lensing peaks only if ${\rm ACF}(f;\Delta t)$ decays sharply from its peak at
$\Delta t=0$ and plateaus by $|\Delta t|\gtrsim|\widehat{\Delta t}|$. In other
words, the narrower the (central) peak in ${\rm ACF}(f;\Delta t)$, the better
the chance of finding the lensing peaks in ${\rm ACF}(F;\Delta t)$. The same
criterion applies to ${\rm ACF}(h;\Delta t)$ for the lensing peaks to appear
in ${\rm ACF}(H;\Delta t)$. However, for all practical purposes, a redder (or
steeper) power spectrum, $P_{f}(\omega)$, leads to a flatter $\langle{\rm
ACF}(f;\Delta t)\rangle_{E}$ which diminishes the probability of finding the
lensing peaks in ${\rm ACF}(F;\Delta t)$. A completely general proof of this
statement is challenging, so here we argue using some concrete examples.

First, consider a Gaussian power spectrum;

$P_{f}(\omega)\propto\exp(-a~{}\omega^{2}/2)\;.$ (C1)

Its Fourier transform is then also Gaussian with $\langle{\rm ACF}(f;\Delta
t)\rangle_{E}\propto\exp[{-\Delta t^{2}/(2a)}]$. Thus, having a large $a$ in
Eq. (C1) gives a narrow $P_{f}(\omega)$ but a broad $\langle{\rm ACF}(f;\Delta
t)\rangle_{E}$. This argument can be made more general; when we scale a
function by $t\rightarrow at$, its Fourier transform gets inversely scaled,

$\text{F.T}\left[g(at)\right]=\frac{1}{|a|}G(\omega/a)\;,$ (C2)

where $G(\omega)$ is the Fourier transform of a generic function $g(t)$. Thus,
if a scaling makes $P_{f}(\omega)$ steeper, $\langle{\rm ACF}(f;\Delta
t)\rangle_{E}$ would be flatter and vice-versa. On the other hand, the
derivative series $h(t)$ has always a bluer power spectrum than that of $f(t)$
as $P_{h}(\omega)=\omega^{2}P_{f}(\omega)$. Thus ${\rm ACF}(h;\Delta t)$ falls
sharper from its peak at $\Delta t=0$ as compared to ${\rm ACF}(f;\Delta t)$
resulting in more prominent lensing peaks in ${\rm ACF}(H;\Delta t)$.

Let us consider another example of red type power spectrum,

$P_{f}(\omega)\propto\left(1+|\omega/\omega_{0}|\right)^{-\gamma}\;,~{}\gamma>0\;,$
(C3)

which behaves like a power-law for $\omega\gg\omega_{0}$ but avoids blowing up
at $\omega\to 0$ by remaining stable at $|\omega|\ll\omega_{0}$. As the
analytical solution for the Fourier transform of Eq. (C3) does not exist for a
generic $\gamma$, we carry out a numerical analysis. We set
$\omega_{0}=10^{-6}$ day-1 so that $P_{f}(\omega)$ behaves like
$|\omega|^{-\gamma}$ in most part (except $\omega\approx 0$). In the left
panel of Fig. 7 we show three such power spectra with $\gamma=1.8,1.0$ and
$1/2$ by the solid blue, orange and green curves respectively. The power
spectrum of the derivative,
$P_{h}(\omega)=\omega^{2}P_{f}(\omega)\approx|\omega|^{2-\gamma}$, is shown by
the dashed curve with respective colour. The three panels on the right show
examples of the intrinsic light curves, $f(t)$, generated from the power
spectra with the three values of $\gamma$ in Eq. (C3); the dashed horizontal
line in each right-panel represents the mean of $f(t)$. It is evident that as
$P_{f}(\omega)$ becomes redder with larger $\gamma$, $f(t)$ shows correlation
till a longer time scale (i.e. two nearby points are more likely to be either
above or below the mean unless not separated sufficiently in time).

Similar to the example illustrated in Fig. 4 for the DRW template, we now
simulate 1000 realisations for the intrinsic flux variability separately for
each of these three power spectra; the corresponding results are arranged in
the three columns of Fig. 8. Then for each realisation, we construct a double
lensed system using Eq. (2) with $\widehat{\mu}=0.86$ and $\widehat{\Delta
t}=50.0$ days, again we consider the perfect condition with marginal noise in
the data for simplicity. The four rows in Fig. 8 show ${\rm ACF}(f;\Delta t)$,
${\rm ACF}(F;\Delta t)$, ${\rm ACF}(h;\Delta t)$, ${\rm ACF}(H;\Delta t)$
respectively from the top. Five random realisations have been shown by the
solid curves in each panel. The dashed black curve and the shaded region in
the first and third row (from the top) panels represent the ensemble average
and $68\%$ quantile around it. From the top-row panels, it is clearly evident
that the redder the power spectrum (i.e. with larger $\gamma$) is, the slower
the ${\rm ACF}(f;\Delta t)$ decays from unity and the broader the $68\%$
quantile is. Therefore, a redder power spectrum in turn reduces the
probability of detecting the lensing peaks in the ${\rm ACF}(F;\Delta t)$ as
shown in the second row (from the top) panels. Also, statistically, we tend to
get more false positive cases, e.g. the blue/red curve in the top-left panel.

On the other hand, for smaller $\gamma$ (flatter power spectrum), ${\rm
ACF}(f;\Delta t)$ tends to decay faster that increases the possibility of
detecting lensed systems through ${\rm ACF}(F;\Delta t)$ (indeed, in the limit
$\gamma\to 0$, $f(t)$ becomes white noise and ${\rm ACF}(f;\Delta
t)\propto\delta(0)$). Nevertheless, the peak prominence is still inferior to
that of ${\rm ACF}(H;\Delta t)$ (bottom panels). In contrast, for all the
$\gamma$’s, ${\rm ACF}(h;\Delta t)$ decays sharply from its peak at $\Delta
t=0$ as evident from the dashed black curves in the third (from the top) row
panels. 666When $\gamma=2$, $P_{h}(\omega)\approx|\omega|^{2-\gamma}=$
constant and $h(t)$ becomes white noise. Hence for $\gamma<2$, $h(t)$ possess
(negative) correlation restricted to only the adjacent point that explains the
local minima in ${\rm ACF}(h;\Delta t)$ just next to $\Delta t=0$ on either
side. But, this does not affect the detection efficiency of ${\rm
ACF}(H;\Delta t)$. Thus, ${\rm ACF}(H;\Delta t)$ shows the sharp lensing peaks
at $\Delta t=\pm\widehat{\Delta t}$ (marked by the dashed vertical lines) for
all the realisations and for all the power spectra considered. Therefore, Fig.
8 (along with the Fig. 4 for DRW process) demonstrates that for any red type
power spectrum, ${\rm ACF}(H;\Delta t)$ is more reliable than ${\rm
ACF}(F;\Delta t)$ for finding the lensed systems and the advantage of using
${\rm ACF}(H;\Delta t)$ over ${\rm ACF}(F;\Delta t)$ increases for a redder
power spectrum 777For completeness, let us discuss if the power spectrum is
flat or blue ($\gamma>0$) although it does not describe QSO light curves well.
In these cases, ${\rm ACF}(F;\Delta t)$ shows sharp lensing peaks and so does
${\rm ACF}(H;\Delta t)$ since the correlation in $h(t)$ is limited to the
neighbouring point(s). Thus both perform well in these scenarios. This is also
illustrated in Appendix C of Bag et al. (2022). Since the lensing signal in
the $\epsilon(\Delta t)$ curve is dominated by ACF$(H;\Delta t)$, the
fluctuation statistics introduced in Bag et al. (2022) is also very successful
in all red power spectra scenarios.

## References

  * Abdalla et al. (2022) Abdalla, E., Abellán, G. F., Aboubrahim, A., et al. 2022, Journal of High Energy Astrophysics, 34, 49, doi: 10.1016/j.jheap.2022.04.002
  * Aghanim et al. (2020) Aghanim, N., et al. 2020, Astron. Astrophys., 641, A6, doi: 10.1051/0004-6361/201833910
  * Bag et al. (2021) Bag, S., Kim, A. G., Linder, E. V., & Shafieloo, A. 2021, Astrophys. J., 910, 65, doi: 10.3847/1538-4357/abe238
  * Bag et al. (2022) Bag, S., Shafieloo, A., Liao, K., & Treu, T. 2022, Astrophys. J., 927, 191, doi: 10.3847/1538-4357/ac51cb
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, Publ. Astron. Soc. Pac, 131, 018002, doi: 10.1088/1538-3873/aaecbe
  * Biggio et al. (2022) Biggio, L., Domi, A., Tosi, S., et al. 2022, Mon. Not. Roy. Astron. Soc., 515, 5665, doi: 10.1093/mnras/stac2034
  * Birrer & Treu (2021) Birrer, S., & Treu, T. 2021, Astron. Astrophys., 649, A61, doi: 10.1051/0004-6361/202039179
  * Birrer et al. (2020) Birrer, S., Shajib, A. J., Galan, A., et al. 2020, Astron. Astrophys., 643, A165, doi: 10.1051/0004-6361/202038861
  * Bonvin et al. (2017) Bonvin, V., Courbin, F., Suyu, S. H., et al. 2017, Mon. Not. Roy. Astron. Soc., 465, 4914, doi: 10.1093/mnras/stw3006
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560.  https://arxiv.org/abs/1612.05560
  * Denissenya et al. (2021) Denissenya, M., Bag, S., Kim, A. G., Linder, E. V., & Shafieloo, A. 2021, arXiv e-prints.  https://arxiv.org/abs/2109.13282
  * Denissenya & Linder (2022) Denissenya, M., & Linder, E. V. 2022, Mon. Not. Roy. Astron. Soc., 515, 977, doi: 10.1093/mnras/stac1726
  * Dobler et al. (2015) Dobler, G., Fassnacht, C. D., Treu, T., et al. 2015, Astrophys. J., 799, 168, doi: 10.1088/0004-637X/799/2/168
  * Einstein (1914) Einstein, A. 1914, Archives des Sciences, 37, 254
  * Geiger & Schneider (1996) Geiger, B., & Schneider, P. 1996, Mon. Not. Roy. Astron. Soc., 282, 530, doi: 10.1093/mnras/282.2.530
  * Goobar et al. (2017) Goobar, A., et al. 2017, Science, 356, 291, doi: 10.1126/science.aal2729
  * Goobar et al. (2022) Goobar, A., Johansson, J., Schulze, S., et al. 2022, arXiv e-prints, arXiv:2211.00656, doi: 10.48550/arXiv.2211.00656
  * Jee et al. (2016) Jee, I., Komatsu, E., Suyu, S. H., & Huterer, D. 2016, J. Cosmology Astropart. Phys, 2016, 031, doi: 10.1088/1475-7516/2016/04/031
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, Astrophys. J., 698, 895, doi: 10.1088/0004-637X/698/1/895
  * Kelly et al. (2015) Kelly, P. L., et al. 2015, Science, 347, 1123, doi: 10.1126/science.aaa3350
  * Khintchine (1934) Khintchine, A. 1934, Mathematische Annalen, 109, 604.  http://eudml.org/doc/159698
  * Lemon et al. (2022) Lemon, C., Anguita, T., Auger, M., et al. 2022, arXiv e-prints, arXiv:2206.07714.  https://arxiv.org/abs/2206.07714
  * Liao (2020) Liao, K. 2020, Astrophys. J., 899, L33, doi: 10.3847/2041-8213/abadfd
  * Liao et al. (2022) Liao, K., Biesiada, M., & Zhu, Z.-H. 2022, Chinese Physics Letters, 39, 119801, doi: 10.1088/0256-307X/39/11/119801
  * Liao et al. (2015) Liao, K., Treu, T., Marshall, P., et al. 2015, Astrophys. J., 800, 11, doi: 10.1088/0004-637X/800/1/11
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201.  https://arxiv.org/abs/0912.0201
  * LSST Science Collaboration et al. (2017) LSST Science Collaboration, Marshall, P., Anguita, T., et al. 2017, arXiv e-prints, arXiv:1708.04058.  https://arxiv.org/abs/1708.04058
  * MacLeod et al. (2010) MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, Astrophys. J., 721, 1014, doi: 10.1088/0004-637X/721/2/1014
  * Millon et al. (2020) Millon, M., Galan, A., Courbin, F., et al. 2020, Astron. Astrophys., 639, A101, doi: 10.1051/0004-6361/201937351
  * Millon et al. (2020) Millon, M., et al. 2020, Astron. Astrophys., 640, A105, doi: 10.1051/0004-6361/202037740
  * Mortlock et al. (1999) Mortlock, D. J., Webster, R. L., & Francis, P. J. 1999, Monthly Notices of the Royal Astronomical Society, 309, 836, doi: 10.1046/j.1365-8711.1999.02872.x
  * Narayan & Bartelmann (1996) Narayan, R., & Bartelmann, M. 1996, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/9606001
  * Oguri (2007) Oguri, M. 2007, The Astrophysical Journal, 660, 1, doi: 10.1086/513093
  * Oguri (2019) Oguri, M. 2019, Reports on Progress in Physics, 82, 126901, doi: 10.1088/1361-6633/ab4fc5
  * Peng et al. (1999) Peng, C. Y., Impey, C. D., Falco, E. E., et al. 1999, The Astrophysical Journal, 524, 572, doi: 10.1086/307860
  * Refsdal (1964) Refsdal, S. 1964, Monthly Notices of the Royal Astronomical Society, 128, 307, doi: 10.1093/mnras/128.4.307
  * Refsdal & Bondi (1964) Refsdal, S., & Bondi, H. 1964, Monthly Notices of the Royal Astronomical Society, 128, 295, doi: 10.1093/mnras/128.4.295
  * Riess et al. (2022) Riess, A. G., et al. 2022, Astrophys. J. Lett., 934, L7, doi: 10.3847/2041-8213/ac5c5b
  * Rodney et al. (2021) Rodney, S. A., Brammer, G. B., Pierel, J. D. R., et al. 2021, Nature Astronomy, doi: 10.1038/s41550-021-01450-9
  * Saha et al. (2006) Saha, P., Coles, J., Macciò, A. V., & Williams, L. L. R. 2006, Astrophys. J., 650, L17, doi: 10.1086/507583
  * Shajib et al. (2020) Shajib, A. J., Birrer, S., Treu, T., et al. 2020, Mon. Not. Roy. Astron. Soc., 494, 6072, doi: 10.1093/mnras/staa828
  * Shu et al. (2021) Shu, Y., Belokurov, V., & Evans, N. W. 2021, Mon. Not. Roy. Astron. Soc., 502, 2912, doi: 10.1093/mnras/stab241
  * Sonnenfeld (2021) Sonnenfeld, A. 2021, Astron. Astrophys., 656, A153, doi: 10.1051/0004-6361/202142062
  * Springer & Ofek (2021a) Springer, O. M., & Ofek, E. O. 2021a, Mon. Not. Roy. Astron. Soc., 508, 3166, doi: 10.1093/mnras/stab2432
  * Springer & Ofek (2021b) —. 2021b, Mon. Not. Roy. Astron. Soc., 506, 864, doi: 10.1093/mnras/stab1600
  * Suyu et al. (2023) Suyu, S. H., Goobar, A., Collett, T., More, A., & Vernardos, G. 2023, arXiv e-prints, arXiv:2301.07729, doi: 10.48550/arXiv.2301.07729
  * Suyu et al. (2017) Suyu, S. H., et al. 2017, Mon. Not. Roy. Astron. Soc., 468, 2590, doi: 10.1093/mnras/stx483
  * Suyu et al. (2020) Suyu, S. H., Huber, S., Cañameras, R., et al. 2020, Astron. Astrophys., 644, A162, doi: 10.1051/0004-6361/202037757
  * Tewes et al. (2013) Tewes, M., Courbin, F., & Meylan, G. 2013, Astron. Astrophys., 553, A120, doi: 10.1051/0004-6361/201220123
  * Tie & Kochanek (2018) Tie, S. S., & Kochanek, C. S. 2018, Mon. Not. Roy. Astron. Soc., 473, 80, doi: 10.1093/mnras/stx2348
  * Treu (2010) Treu, T. 2010, Ann. Rev. Astron. Astrophys., 48, 87, doi: 10.1146/annurev-astro-081309-130924
  * Treu & Marshall (2016) Treu, T., & Marshall, P. J. 2016, Astron. Astrophys. Rev., 24, 11, doi: 10.1007/s00159-016-0096-8
  * Wiener (1930) Wiener, N. 1930, Acta mathematica, 55, 117
  * Wong et al. (2020) Wong, K. C., et al. 2020, Mon. Not. Roy. Astron. Soc., 498, 1420, doi: 10.1093/mnras/stz3094
  * Wong et al. (2020) Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, Mon. Not. Roy. Astron. Soc., 498, 1420, doi: 10.1093/mnras/stz3094
  * Zu et al. (2013) Zu, Y., Kochanek, C. S., Kozłowski, S., & Udalski, A. 2013, Astrophys. J., 765, 106, doi: 10.1088/0004-637X/765/2/106