© The Authors 2025

1. Introduction

The precise measurement of the difference in the time of arrival between different lensed components of a gravitational lens is the first step in time-delay cosmography (TDC), proposed by Refsdal (1964). Time-delay cosmography is enabled by lensing, when the light of a single source experiences path differences through cosmological distances, and opens the possibility of measuring the Hubble constant, H₀¹. For the measurement of this difference in the time of arrival to be possible, the distant source must be very bright and variable: supernovae fit this requirement, and so do quasars. In practice, the same photometric variations in the source can be observed at different times in the light curves of each lensed component, and the temporal shift required to align these light curves is called the time delay. The total time delay, Δt_tot, between the images depends on the geometric difference between the optical paths, Δt_geom, and on the difference between the lens potential values at the positions of the lensed images, Δt_grav. The total time delay is the sum of the geometric and gravitational delays and is directly proportional to the so-called time-delay distance, a ratio of the three angular-diameter distances, between the observer and the source, the observer and the lens, and the lens and the source (e.g. Refsdal 1964; Suyu et al. 2010).

Time-delay cosmography does not rely on cosmic microwave background (CMB) measurements or on any intermediate standard ruler or standard candles, and it does not involve the difficult construction of a multi-step ladder. The primary element in TDC is the time-delay measurement, but turning this delay into a time-delay distance, and finally into H₀, requires a mass model for the lens galaxy and a measurement of the contributions of nearby objects along the line of sight. We refer the reader to Birrer et al. (2024) and Treu & Shajib (2023) for reviews of the methodology of the modelling and line-of-sight mass contribution.

For a given lensed quasar system, the uncertainty in the inferred H₀ value can be divided into systematic and statistical errors. The systematic part of the error budget is typically dominated by the accuracy of the mass model and its inherent degeneracies, the most famous of which is the mass sheet degeneracy (MSD; e.g. Falco et al. 1985; Gorenstein et al. 1988; Kochanek 2002; Schneider & Sluse 2013, 2014; Blum & Teodori 2021). These can be mitigated by using spatially resolved kinematics of the lensing galaxy (e.g. Yıldırım et al. 2020; Shajib et al. 2023), which brings other but distinct degeneracies. On average, the statistical uncertainty in H₀ is roughly divided among contributions from lens modelling, line-of-sight characterisation, and time delay measurements.

TDCOSMO (Time Delay COSMOgraphy) is a collaborative project resulting from the fusion of COSMOGRAIL (Courbin et al. 2018), H0LiCOW (Suyu et al. 2017), STRIDES (Treu et al. 2018), and SHARP (Chen et al. 2019). The current objective of TDCOSMO is reaching a 1% precision on H₀, both by expanding the pool of ancillary-complete lensed quasars and understanding and controlling systematic errors (see, e.g., Birrer & Treu 2021).

This paper in the TDCOSMO series focuses on the primary ingredient needed for the method to work: time-delay measurements between multiple images of strongly lensed sources at cosmological distances, here quasars. Measuring lensed quasars time delays requires high-quality light curves with adequate temporal sampling over long periods of time, in an effort to avoid interruptions other than the inevitable seasonal gaps. Typically, a one-year baseline is required, but most monitoring campaigns carried out so far have obtained several years per object, sometimes up to two decades.

There are two main difficulties to overcome in order to achieve time-delay measurements at optical wavelengths. The first is that the lensed images are faint (20.2 r-mag on average for the present sample), and most often blended in ground-based images: astronomical seeing is thereby a major limitation. Given that the image separations in lensed quasars are typically not much larger than the average astronomical seeing, the data need to be processed with deblending techniques such as point spread function (PSF) fitting or more sophisticated forward modelling of the pixels (e.g. Millon et al. 2024; Michalewicz et al. 2023; Cantale et al. 2016; Magain et al. 1998). High signal-to-noise ratio (S/N) imaging is mandatory to make the deblending reliable, requiring fairly large telescopes (typically 2 m) and long exposure times (typically 30 min on-target). This achieves photometric precisions between a few milli-magnitudes to a few tens of milli-magnitudes for the brighter targets (∼19 r-mag).

The second difficulty is that microlensing by stars in the lensing galaxy introduces additional variations in the quasar images. These variations are uncorrelated in each quasar image and occur on timescales of the order of months to years. Notably, while this effect typically increases the uncertainty in time delays of lensed quasars, it also affects lensed type Ia supernovae and compromises their standard candle nature altogether (Foxley-Marrable et al. 2018; Weisenbach et al. 2024). For monitoring observations with temporal sampling of the order of one observation every 3−4 days and photometric precision of 10 milli-mag, it becomes necessary to monitor any given system for several years, averaging out microlensing over long periods of time. This has been the strategy adopted by most past monitoring campaigns (e.g. Burud et al. 2000, 2002a,b; Hjorth et al. 2002; Ullán et al. 2003; Kochanek et al. 2006; Shalyapin & Goicoechea 2017, 2019; Muñoz et al. 2022; Shalyapin et al. 2023), including the COSMOGRAIL program that has measured more than 30 time delays (Courbin et al. 2005; Eigenbrod et al. 2005; Bonvin et al. 2019; Millon et al. 2020a). The future Vera C. Rubin Observatory (Ivezić et al. 2019, Rubin-LSST) observations of lensed quasars will fall in this regime, with the drawback that years will be necessary to measure time delays in new objects, but with the advantage that Rubin-LSST will monitor all targets in the southern hemisphere, even those not identified as lensed quasars yet.

In pre-Rubin-LSST times, when objects are still monitored one by one, it has been shown that time delays can be measured within 1−2 years, provided daily high-S/N observations are possible (Courbin et al. 2018; Millon et al. 2020b). Such high-cadence and high-S/N data allow for the capture of low-amplitude and fast quasar variations occurring on shorter timescales than microlensing, and hence act as a natural frequency comb, discriminating between the two types of signals.

All known bright and well-separated quadruply lensed quasars in the southern hemisphere have been monitored for time delays in the past. This work tackles most of the known remaining such objects in the south, which are fainter and more narrowly separated (and thereby more difficult) than the ones monitored in earlier programs – roughly doubling the sample of southern time-delay lenses. The determination of time delays was achieved thanks to consistent, high-quality optical observations at two telescopes, and refinements of the data processing.

It should be noted that while this program was originally designed to observe each object with a single facility to achieve homogeneous data quality and photometric calibration, there were several interruptions due to the impact of the COVID-19 pandemic on the operations of the telescopes. As a consequence, we were affected by large gaps in our light curves and often had to observe objects with the two telescopes, as COVID-19 did not impact the two observatories in the same way.

We describe the modalities, observations, and statistics in Sect. 2. The methodology employed for extracting the light curves from the 21 582 exposures is detailed in Sect. 3. Section 4 lays out the methodology of the estimation of the time delays, the results of which are presented in Sect. 5. Next, in Sect. 6 we briefly discuss what the future holds for time-delay determination in light of the upcoming Rubin-LSST observations and the findings of this program. Finally, we draw conclusions in Sect. 7. The magnitudes quoted herein are calibrated on the Vega scale.

2. Observations and statistics

2.1. Targets

We monitored 22 lensed quasars visible from the southern hemisphere, including 19 quadruply lensed ones. This sample of 19 represents about half the known quadruply lensed quasars under 20° declination – the other half being brighter and having mostly been monitored for time delays already. We list co-ordinates, redshifts when available, and discovery papers in Table 1.

2.2. Observing facilities

The two instruments that observed our targets are multi-charged-coupled device imagers. For our purpose, a single charged-coupled device (CCD) provided a sufficient field of view, so we elected the one with the best read noise characteristics.

2.2.1. OmegaCAM/VLT Survey Telescope

The VLT Survey Telescope² (VST) located at the ESO Paranal Observatory is a wide-field survey telescope with a primary mirror diameter of 2.65 meters. It is currently the largest telescope in the world solely dedicated to sky surveys in visible light, with a 1 square-degree field of view. Its imager is the OmegaCAM, a large 268-megapixel camera, with 32 CCDs. Combining the 13 fields observed with this instrument yields a 13 square degrees survey. We used CCD #13 only for time delays – a small fraction of the total area covered by OmegaCAM.

2.2.2. WFI/MPG-ESO 2.2-meter telescope

The MPG/ESO 2.2-meter telescope³ (2p2), operational since 1984, is located at the La Silla Observatory. Its Wide Field Imager (WFI), a focal reducer-type camera of 68 megapixels at the Cassegrain focus of the telescope, provides a field of view of roughly 0.25 square degree. Of the eight CCDs constituting the WFI, we used CCD #4.

2.3. Data acquisition and statistics

The photometric data were obtained from the VST and 2p2 telescopes, via daily monitoring sessions from September 2018 to December 2023. Each session consisted of four dithered 320-second exposures through the SDSS-r (VST) or ESO-Rc (2p2) filters. Observations mostly occurred at an airmass below 1.5, with exceptions to extend the visibility window: this is crucial for longer time delays, for which the shifted seasonal light curves overlap only briefly. The median cadence was 1 day for all the 22 monitored targets, which is precisely the objective of this program. However, some sessions were missed due to bad weather, technical problems, or maintenance. Furthermore, the exposures were filtered in the final step of the data processing based on the seeing, sky level, normalisation of the flux in the image, and shape of the PSF. Specifically, at the VST, 8821 exposures were kept out of a total of 8944 captured. At the 2p2, 12 184 out of 12 637 were kept. Thus, the average effective cadences are slightly higher than the medians – ranging from 1.3 to 2.2 days at the VST, and from 1.3 to 1.8 day at the 2p2. The median seeings (measured from the images) were $0\stackrel{″}{.}91$ at the VST and $1\stackrel{″}{.}06$ at the 2p2. Their distribution is shown in Fig. 1. We provide a comprehensive overview of the key observation statistics of the monitored targets in Table C.1.

Fig. 1. Astronomical seeing distribution in the images used to infer time delays in this data release.

3. Light curve extraction

We used the lightcurver Python package (Dux 2024) to extract the photometry of our targets. lightcurver provides the infrastructure to apply STARRED (Michalewicz et al. 2023; Millon et al. 2024) to time-series images. The method implements a type of generalised PSF fitting method for light curve extraction, whereby we take advantage of the hundreds of dithered observations available for each target to accurately constrain the astrometry of the lensed images, as well as the non-point-source components of the targets. The process is an extension of drizzling, where a high-resolution model represented on a grid of pixels is forward-fitted to all observations simultaneously – but, unlike drizzling, the PSF and point-source fluxes here can vary between observations. The method requires a precise zero-point calibration of the observed frames, which we detail in the paragraphs below and is implemented in lightcurver. It should be noted that when a target was observed with the two telescopes, the light curve extraction process was carried out once per telescope-dataset due to different pixel scales.

3.1. Pre-processing and bookkeeping

Initial steps on the individual exposures included flat-fielding using sky flats, and bias subtraction. Flat-fielding is particularly important to achieve the target photometric precision: a constant zero point is needed across the field, because stars are used as references for calculating the relative flux normalisation between epochs. Next, a conservative sky model was calculated and subtracted from each image using sep (Barbary 2016; Bertin & Arnouts 1996). The subtraction of the sky was straightforward, probably due to our targets all being in extragalactic fields. An astrometric solution was then found for each image using the Astrometry.net software (Lang et al. 2010). This allowed us to check that the pixel scale (0$\stackrel{″}{.}$213 and 0$\stackrel{″}{.}$237 per pixel for OmegaCAM and WFI, respectively) was constant across epochs, and also permitted the easy identification and elimination of bad pointings. Unlike what was done in previous monitoring campaigns (e.g., Courbin et al. 2018; Bonvin et al. 2019; Millon et al. 2020a,b), our frames were not interpolated onto a reference one, because interpolation leads to the loss of sub-pixel information – making the drizzling-like process mentioned above less effective.

3.2. Selection of calibration stars

The success of a light curve extraction depends enormously on the ability to precisely calculate the proper normalisation (zero point) of each image. In turn, the normalisation requires an excellent PSF model. Stars suitable for the construction of a PSF model near each lensed target were selected from the Gaia catalogues (Gaia Collaboration 2016, 2023). Specifically, we requested that each selected star has a g-mag between 16.8 and 19.5 (to avoid saturation given our set-up, while maintaining a high S/N), low variability (phot_g_mean_flux_over_error above 100), and is well fitted by a point-source solution in Gaia (astrometric_excess_noise_sig below 3). The (typically) ten closest such stars were elected both as PSF models and flux-normalisation references – depending on the field, within one to three arcminutes away from the target. Due to dithering and pointing irregularities, not all selected stars were always in the footprints of all frames, but this is not a problem.

3.3. Empirical PSF modelling

Cutouts of the stars and lens systems were extracted from each image. The noise map of each cutout was calculated accounting for Poisson noise (square root of the data in electrons) and read plus sky background (standard deviation of the noise in the background). Masks were implemented by strongly boosting the noise maps at the desired locations. Specifically, masks were needed for eliminating cosmic rays, identified using the L.A.Cosmic algorithm (van Dokkum 2001), and bad detector rows. Empirical PSF models were then constructed with STARRED: STARRED first fits a single elliptical Moffat (Moffat 1969) profile to the stars cutouts, then optimises the pixels in a two-times supersampled pixelated grid to fit the remaining residuals. The natural positional shifts of the stars (with respect to the grid of pixels of the CCD) allow for the recovery of sub-pixel information, which is why the model is best described on a supersampled grid of pixels. The overfitting of the noise is prevented by regularising with an isotropic family of wavelets (starlets). We note that STARRED builds a PSF model, but also and more importantly a kernel that transforms a two-pixel full width at half maximum (FWHM) Gaussian into the PSF. We used these kernels to forward-model the lenses, which were also represented on a twice-supersampled pixelated grid with a two-pixel resolution, such that convolving with the kernel would reproduce the PSF of the data (see Sect. 3.6 below). lightcurver and STARRED provide infrastructure for dealing with spatial distortion of the PSF, but our datasets did not require it: there were no significant residuals when fitting the PSF on the stars. Thus, we did not feel a need to account for colour differences between the stars either.

Thanks to the precautions taken when preparing star cutouts, the building of a good initial guess based on measured seeing, and the robust optimisation algorithm leveraged by STARRED (adabelief, Zhuang et al. 2020), the PSF generation procedure was extremely reliable, even with the occasional out-of-focus or trailed exposures (see Fig. 2). Only a few dozen images had to be discarded because of failed PSF modelling.

Fig. 2. Example simultaneous fit of a PSF model to ten stars, in a (rare) difficult case of unfocused observation. The frame chosen for this figure could still be used for photometry, not leading to an outlier in the light curves. The kernel, at the top left, transforms a two-pixel FWHM Gaussian into the PSF at the bottom left. Note that for our lens photometry we used the kernels and not the PSFs (see Sect. 3.6 and Fig. 3). The bottom row shows the residuals of the joint fit. Similar residuals were obtained in all conditions of seeing, focusing, and tracking.

3.4. Relative zero-point calibration

The relative zero points between frames were calculated by comparing the fluxes of the stars that had also served as PSF references. To avoid systematic errors, we used the same infrastructure (see Fig. 3 and paragraph about lens photometry below) to extract both star fluxes and lens fluxes. For the stars, however, the pixelated background was set to zero, making the process equivalent to PSF photometry. For each calibration star, we obtained a time series of flux measurements across all the frames in which it appeared. The raw flux time series for each star, F_i, j (where i indexes the stars and j indexes the frames), was normalised by its median:

Fig. 3. Overview of the final step of light curve extraction: the forward-modelling of the lens cutouts. Note that the colour coding is only used to differentiate point-source and pixelated components of the models, and does not indicate observations in multiple filters. Top: Components of the forward-modelling process of the zero-point-calibrated cutouts, starting with a high-resolution model containing a fully tunable grid of pixels and point sources of tunable flux (added to the grid of pixels as two-pixel FWHM Gaussians). The model is degraded to the resolution of the data, and the residuals are minimised through gradient-based optimisation of the high resolution models. Bottom: Exploring the result of the forward modelling for a selected epoch. A is the fitted high-resolution model. B is the kernel, which transforms a point source in A into the PSF of the data. C is the convolution of A and B, which is then downsampled down to the pixel scale of the data, yielding D. A, C, and D have the background component in red and the point-source component in white. E is the data cutout. F and G are again the data, but with the background and point-source components subtracted, respectively. Finally, H shows the residuals upon subtraction of the two components from the data. This method allows for the precise determination of the astrometry of the point sources, reveals the morphology of the non-point source components at high resolution, and provides precise photometry unbiased by flux leakage between point source and non-point source components.

$$\begin{array}{c}\hfill F_{i,j}^{\mathrm{norm}}=\frac{F_{i,j}}{\mathrm{median}(F_{i,:})}·\end{array}$$(1)

Next, we introduced a global scaling factor, c_i, for each star to minimise the intra-frame scatter among stars. These factors, c_i, were determined by minimising the scatter in the weighted mean flux for each frame, j, subject to the constraint that mean(c_i) = 1. The total cost function that the c_i were requested to minimise was

$$\begin{array}{c}\hfill \mathrm{Cost}={\displaystyle \sum _j}\frac{{\displaystyle {\sum }_i}{w}_{i,j}{(c_i{F}_{i,j}^{\mathrm{norm}}-{\mu }_j)}^2}{{\displaystyle {\sum }_i}{w}_{i,j}},\end{array}$$(2)

where the weighted mean for frame j is

$$\begin{array}{c}\hfill {\mu }_j=\frac{{\displaystyle {\sum }_i}{w}_{i,j}{c}_i{F}_{i,j}^{\mathrm{norm}}}{{\displaystyle {\sum }_i}{w}_{i,j}},w_{i,j}=\frac{1}{{\sigma }_{i,j}^2}·\end{array}$$(3)

After applying the optimised scaling factors, the adjusted normalised fluxes of all calibration stars within each frame were combined using a weighted average with 3-σ rejection, yielding the normalisation coefficients, N_j. The frames were brought to the same zero point by dividing their pixels by the N_j coefficients. This single coefficient per frame neglects the inevitable spatial zero-point variation, often caused by flat-fielding (multiplicative) or sky subtraction (additive) errors. The hope is that this effect is mitigated by the calibration stars being spread in a small region around the target. Nonetheless, it might introduce systematic normalisation errors not fully accounted for by the scatter of the zero points derived from different stars⁴. This scatter, σ_zp, was retained for later use in estimating the error bars on the fluxes of the lensed images.

3.5. Absolute zero-point calibration

The absolute magnitude calibration of each frame was made using Gaia magnitudes of our chosen calibration stars: the Gaia magnitudes were approximately converted to the Sloan Digital Sky Survey (SDSS) r filter using the G − r relation given by Evans et al. (2018, Table A.2). In the case of the 2p2 data, this introduces a further inaccuracy due to the slight mismatch between the SDSS-r and ESO-Rc filters, so, even though the relative night-to-night calibration is very precise, the absolute calibration is by nature approximate. Overall, the average absolute Vega zero points were 31.44 mag and 31.54 mag for the VST and 2p2 telescopes, respectively, with global scatters (due to observing at different airmasses, with different sky transparencies) of 0.05 mag and 0.15 mag, respectively.

3.6. Photometry of lensed images

We performed the lens photometry using a forward modelling approach to accurately capture all the flux in the vicinity of the targets, without relying on high-resolution imaging for contaminant subtraction. Specifically, we constructed a high-resolution model that included point sources with adjustable flux (the quasar light curves) and a constant pixelated background component (absorbing extended photometric contaminants). This model was expressed on a twice-supersampled grid of pixels, which matched the supersampling of the PSF models. Each pixel in this supersampled grid was allowed to vary: we call this part the ‘pixelated background’. On top of it, the point sources were injected as 2-pixel FWHM Gaussians. Together, the two components form our ‘supersampled model’. The astrometry of the point sources and the pixelated background were common to all epochs⁵. The point-source fluxes, translations (implemented by interpolation of the supersampled model), and constant background terms were allowed to vary per epoch. To compare the high-resolution models to the data, we convolved the models with the kernels, which convert the 2-pixel FWHM Gaussians to the epoch-specific PSFs. The resulting images (one per epoch) were then downsampled twice to match the resolution of the data. These steps of degrading the supersampled models to the resolution of the data are illustrated in the lower part of Fig. 3. The model parameters were tuned to minimise the chi-squared, computed by summing the squared noise-normalised residuals between data and model-images, across epochs and pixels. The number of fitted parameters was typically quite high (∼10 000, see Table 2), but not an issue for gradient-based optimisation, which is enabled here by the automatic differentiability of jax, the framework in which STARRED is written. The top of Fig. 3 provides a schematic depiction of the process. The modelling fully explained the data for all targets (reduced chi-squared ∼ 1) and yielded a flux for each point source, p, in a given exposure, i: F_p, i.

3.7. Flux uncertainty estimation

Again thanks to the autodifferentiable nature of STARRED, we have access to an estimate of the uncertainties of each parameter through the Hessian matrix evaluated at the local minimum reached during optimisation. Specifically, we computed the diagonal of the Fisher information matrix, which, in the absence of deblending issues, yields error bars that are representative of the photon noise error. Thus, we denote this uncertainty as ${\sigma }_{p,i}^{\mathrm{photon}}$: the photon noise on the flux of a single point source. This is to be combined with the error on the relative zero point of the frame in question:

$$\begin{array}{c}\hfill {\sigma }_{p,i}=\sqrt{{\left({\sigma }_{p,i}^{\mathrm{photon}}\right)}^2+{\sigma }_{\mathrm{zp}}^2}.\end{array}$$(4)

The final light curves were obtained by combining the measurements (usually four frames, i = 1, 2, 3, 4) within a single night using a weighted average:

$$\begin{array}{c}\hfill F_p=\frac{1}{Z}{\displaystyle \sum _{i=1}^4}\frac{F_{p,i}}{{\sigma }_{p,i}},\end{array}$$(5)

with

$$\begin{array}{c}\hfill Z={\displaystyle \sum _{i=1}^4}\frac{1}{{\sigma }_{p,i}},{\sigma }_p=1/Z.\end{array}$$(6)

σ_p is the uncertainty on the flux of point source p on a given night.

We also kept track of the standard deviation within F_p, i, denoted hereafter by ΔF_p, as it provides a more empirical estimation of the uncertainty. We denote the same quantity, but in units of magnitude, ΔM_p. The median σ_p over the entire dataset is 0.017 mag, and 0.020 mag for ΔM_p. This slightly higher observed scatter compared to photon noise and normalisation errors can be attributed to a few very narrowly separated objects, in which flux leakage between point sources occurs. These are discussed on a case-by-case basis in Appendix A. Figure 4 shows what uncertainty can be expected in a night of monitoring, for a single lensed image of a given magnitude. For the purpose of time-delay determination, we took our final uncertainty on the magnitude of a point source on a given night as the average of ΔF_p and σ_p. This prescription qualitatively increased the robustness of the error bars compared to the local scatter in the extracted light curves. As an example, we show the extracted light curves of J1537−3010 in Fig. 5, while the others are delegated to Appendix B.

Fig. 4. Average nightly scatter of the measured lensed-image magnitudes, plotted against the mean magnitude of the lensed image at hand. The scatters are empirical, and thereby include photon noise, read noise, deblending errors, and normalisation errors between frames. The dashed line is the noise estimate given by the ESO/WFI exposure time calculator (ETC), for a point source in the r band in good seeing conditions (1″), average moon illumination, and 320 second exposures. The similarity between idealised and empirical errors indicates that our procedure has an excellent deblending performance. For reference, the line fit in blue indicates the measured photometric performance given the magnitude of a given target point source.

Fig. 5. Light curves extracted from the imaging data of J1537−3010 (labelling of the lensed images in Fig. 6), monitored at both observatories and spanning four seasons, with the 2020 and 2021 ones shortened by the COVID-19 pandemic. Note that the empty season gaps are cut from the plot. A magnitude offset was added to the individual curves as is indicated in the legend, for display purposes.

3.8. The result: Light curves and models

Besides the light curves themselves, the second products of our method are the high-resolution models fitted to all epochs simultaneously. These were fitted blindly together with the photometry, without prior knowledge of the astrometry or using existing Hubble Space Telescope (HST) images as a prior. Despite this, the fitted models accurately capture the morphology of Einstein rings, lensing galaxies, or nearby faint galaxies, while also recovering the HST astrometry with excellent precision (see Sect. 6.3). The models are shown in Fig. 6, in which we provide for each lens the labelling of the lensed images, a typical data frame from the telescope, the fitted high-resolution model, and an HST image when publicly available – for many, HST imaging only became available after the photometric extraction was conducted. The faithfulness of the high-resolution models, together with the reduced χ² of the fits (∼1 in all cases), shows that the relative zero points between epochs were calibrated with high precision, and that the PSF models were well estimated in all epochs.

Fig. 6. Imaging preview of all the monitored objects. For each lens, we include a representative exposure from the dataset with median seeing and low sky background, the high-resolution fitted model (see Sect. 3.6), and, when available, an HST/Wide Field Camera 3 (WFC3) image for comparison. The high-resolution models live on a grid twice-supersampled compared to the original images, with a two-pixel FWHM PSF. The resolution (FWHM of the PSF) in the models is thereby 0$\stackrel{″}{.}$213 for VST data, and 0$\stackrel{″}{.}$237 for 2p2 data. This is to be compared with the 0$\stackrel{″}{.}$07 resolution of the HST/WFC3 images. Note also that the F814W filter is much wider and redder than the r band, such that more details can be revealed in the redshifted objects of our fields. All cutouts are oriented with north facing up and east facing left.

4. Measuring time delays

Different methods have been proposed for the purpose of estimating time delays between light curves, such as autocorrelation (Press et al. 1992), or minimising the dispersion of curve differences (Pelt et al. 1996). These methods are fast, but they cannot properly account for complex extrinsic variations imprinted on the individual light curves. Extrinsic variations are always present in lensed quasar optical light curves, mainly due to the microlensing of the individual lensed images by stars of the lensing galaxy. Thus, superimposing light curves requires accounting for both the intrinsic variation of the quasar, common to all the lensed images, and extrinsic variations individual to each lensed image.

4.1. Matching curves with splines

A purely data-driven method doing exactly this, which was proven to work well in the 2014 time-delay challenge (Liao et al. 2014) as has been reported by Bonvin et al. (2016), is the simultaneous alignment of the curves using free-knot splines. In this method, each curve is given an extrinsic model of variation that modifies its magnitude in a smooth way. The modified curves are also shifted in time, until they all match a single free-knot spline representing the intrinsic variation in the quasar. All parameters – of the splines and the shifts – are optimised simultaneously, until a match is found. An example is shown in Fig. 7, in which three curves were matched to a single spline with timeshifts and modulations. This is the strategy we used herein, as is implemented in the PyCS3 (Millon et al. 2020c) toolbox, which uses polynomials or a free-knot B-spline (Molinari et al. 2004) to represent curves. The general methodology is given in the third figure of Millon et al. (2020a), with the exception being that we did not include the regression difference method in the final estimates. The splines were indeed found to be a more robust and accurate method than regression differences on simulated data (see Millon 2021, Sect. 2.4.4). The extrinsic variations were represented with polynomials or cubic splines, depending on the complexity of the light curves. In order of complexity, the extrinsic models, M_ext, included polynomials of degree one (a linear function) or three (a cubic polynomial), a cubic spline with one internal knot forced to the centre, or free-knot cubic spline with one or more internal knots.

Fig. 7. Curve matching example, in which a single spline (in black) matches all curves simultaneously. Due to microlensing in particular, which adds independent lower-frequency variations to each curve, the curves could not be matched in this way with time and magnitude shifts only. Thus, an extra modulation is applied to each separately. These modulations are represented by lower order splines, displayed at the bottom. Physically, the black spline models the intrinsic variations of the source, while the coloured splines undo the effects of microlensing. Linking with the notation of Sect. 4, Mext is here a cubic spline with a single internal knot of fixed position, and Nint = 9. To determine a time delay and its uncertainty, this matching was repeated thousands of times with artificially further-modulated curves (mimicking microlensing), other realisations of the noise, and different freedoms given to the splines. What this matching is like for the other targets can be found in a Jupyter notebook. This notebook is part of the repository containing the photometry extracted herein and the code to estimate the time delays.

https://github.com/duxfrederic/TDCOSMO_XVII_time_delays/blob/main/initial_guess.ipynb

However, not all extrinsic models were used on all sets of curves. Often, the simplest modulations would not have allowed for a proper matching of the curves: in such cases, strong systematic errors would have dominated the error budget due to the imperfect fit. Similarly, for curves with both little intrinsic and extrinsic variations, electing a complex extrinsic modulation would have resulted in a completely degenerate model, ultimately yielding overestimated uncertainties. Hence, for each dataset we started with the simplest extrinsic model that permitted a good alignment of the curves, and included two more levels of complexity. For example, if a set of curves could be well aligned (that is, matched within their noise envelope over their entire length) with splines of one internal knot, splines with two and three internal knots were added to the pool of possible extrinsic variation models as well.

Together with the extrinsic models, the intrinsic variation of the quasar was modelled with a free-knot cubic spline, with N_int internal knots⁶. A first value of N_int was selected such that the finest features of the best curve could be well captured by the spline. Next, N_int was decreased until a certain smoothing of the finest features would be noticed. In the interval between the two values, an additional N_int value was added to the pool of possible models. We call a pair of models (N_int,  M_ext) an estimator, and as such end up with a grid of nine estimators, with three different intrinsic splines of N_int internal knots, and three models of extrinsic variations, M_ext.

4.2. Microlensing bias and error estimation

To estimate the reliability and uncertainty of a time-delay estimation, the inference was performed repeatedly with each estimator on mock curves with distinct realisations of the observed noise. Importantly, red noise mimicking the effects of microlensing is injected into the mock curves: this essentially transforms a systematic error – arising from fitting microlensing, which distorts the curves in a way that can bias the time-delay estimates – into a statistical error. The main estimate and its uncertainty can then be read from the histogram of optimised time delays. The separate pools of delays measured from mocks (one per estimator) are then combined as is described in Sect. 3.3 of Millon et al. (2020a), with the tension parameter, τ (Bonvin et al. 2018), set to 0.5. This provides a way to combine estimates from different microlensing and intrinsic variation models, halfway between selecting the estimate with the best precision and marginalising over all estimators, and effectively eliminates outliers. This process can also lead to asymmetrical error bars, but on its own cannot account for the covariance between delay pairs.

4.3. Covariance between delays

A quadruply lensed quasar provides six independent pairs of curves from which a delay can be measured. However, the strategy of PyCS3 is maximising the available signal in each temporal bin by shifting all four curves together: implying only three independent time-shifts. As such, a covariance is bound to exist between different delay pairs, and can slightly change the overall precision in the application of TDC.

To compute our covariance matrices, we assume that the delays obtained from shifting mock light curves are normally distributed: this is a good approximation, but the mocks further than 4σ from the mean were discarded as to not artificially increase the error bars because of failed fits. Then, the covariance is read from its definition:

$$\begin{array}{c}\hfill \mathrm{cov}(p_1,p_2)=\mathbb{E}\left[(p_1-\mathbb{E}[p_1])(p_2-\mathbb{E}[p_2])\right].\end{array}$$(7)

This was computed on the pools of mocks that passed the τ = 0.5 selection mentioned above. We also made sure that the square root of the diagonal – the standard deviations, σ_d, on individual delays, neglecting covariance – were always within 20% of the confidence interval read directly from the histograms of mocks; that is,

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{d}}\sim \frac{P_{84}-P_{16}}{2},\end{array}$$(8)

where P_n are the n^th percentiles read from the histograms of mocks. This ensured that the standard deviations given by the diagonal of the covariance matrix are indeed representative of the actual widths of the distributions of the delays optimised on mock curves. Moreover, any bias in the median value (recovered versus injected) of the mock-optimised time-delays was added in quadrature to the diagonal of the covariance matrix, but this contribution was always negligible compared to the statistical uncertainties. We note that the resulting covariance matrix is ill-conditioned: for a quadruply lensed quasar, it is a 6 × 6 matrix of rank 3 (due to dealing with only three independent shifts). Nonetheless, for convenience, we provide all six delays together with their singular 6 × 6 covariance matrix. In practice, one should choose a reference lensed image and keep only the three delays relative to it, with the corresponding 3 × 3 covariance sub-matrix.

5. Results

Overall, this set of TDCOSMO lenses has proven more challenging than lenses of previous COSMOGRAIL publications, in particular due to the fainter (20.2 average magnitude versus 19.2 in Millon et al. 2020b) and more compact lensing configurations, all the brighter, well-separated targets having been monitored for time delays in the past. Nevertheless, at least one reliable time delay could be inferred for each monitored lens. Still, this shows that obtaining reliable photometry of even fainter and more compact lensed quasars will become more difficult, at least with seeing-limited observations. Not all the targets presented in this work have been monitored for the first time: four were already discussed in COSMOGRAIL publications; namely, HE J0230−2130, J0832+0404, WFI J2026−4536 by Millon et al. (2020a), and WG J0214−2105 as well as PS J1606−2333 by the previous data release of the present program (Millon et al. 2020b).

We delay the discussion of the challenges encountered in each target, as well as delay values and covariance matrices, to Appendix A. Here, we provide a summary of the relative uncertainties in each target in Table 3, and the labelling of the lensed images can be found in Fig. 6.

5.1. Uncertainty propagating to the Hubble constant

The time-delay part of the H₀ error budget is roughly directly given by the uncertainty of the time delay itself. This is a valid working assumption, but the situation is complicated when several delays are available. In the absence of correlation, we could argue that a quadruply lensed quasar provides three fully independent delays. These can then be combined as independent Gaussian estimates, with the resulting uncertainty scaling as

$$\begin{array}{c}\hfill \mathrm{relative}\mathrm{uncertainty}\sim \frac{1}{{\displaystyle {\sum }_i}\frac{{\tau }_i}{\mathrm{\Delta }{\tau }_i}},\end{array}$$(9)

where the τ_i and Δτ_i are three individual delay and uncertainty pairs. This is a benchmark of the best possible achievable precision, and works in the absence of correlation between delays. On the other hand, if the delays are closely correlated, then the relative uncertainty that maps to H₀ is roughly that of the most precise delay:

$$\begin{array}{c}\hfill \mathrm{relative}\mathrm{uncertainty}\sim {\displaystyle \underset{j}{min}}\left(\frac{\mathrm{\Delta }{\tau }_j}{{\tau }_j}\right),\end{array}$$(10)

where j denotes all possible delay pairs. It is likely that the actual uncertainty mapping to H₀ in a TDC analysis will fall somewhere between the two above estimates. We give both estimates and the most precise delay for each system in Table 3.

6. Discussion

6.1. A look back on intrinsic and extrinsic variability

This program aimed for a very high cadence of monitoring and high photometric precision per epoch so that extrinsic variations could be better disentangled from the faster intrinsic variation of the quasar. To verify whether this was realised, we provide an estimate of the intrinsic variability of each source quasar, and a comparison to that of the extrinsic variability of their lensed images. We turned to the structure function (SF) to measure the variability of a light curve:

$$\begin{array}{c}\hfill \mathrm{SF}(\mathrm{\Delta }\tau )=\sqrt{{\langle {(m_j-m_i)}^2-{\sigma }_i^2-{\sigma }_j^2\rangle }_{t_i-t_j\in [0,\mathrm{\Delta }\tau )}},\end{array}$$(11)

where i and j denote epochs. So, for each bin of lag time Δτ, the SF is the average squared differences of points of magnitude, m. We also subtracted the squared uncertainties in each point to avoid an artificial contribution to the SF from noise. However, given that all lensed images will show at least some extrinsic variability, one cannot access the true intrinsic variability of the quasar: we were instead forced to elect a reference lensed image, in which we assumed no extrinsic modulation. An initial working hypothesis is that the further a lensed image is from the halo of stars of the lensing galaxy, the less microlensing it will experience on average. Hence, assuming that the estimated delays are correct, microlensing curves were calculated by subtracting the furthest-image curve from shifted and interpolated versions of the other light curves. If the initial working hypothesis is correct, then we should observe a further dependence of the amount of variability in the microlensing curve as a function of distance to the lens. This was realised in our light curves, as is shown in Fig. 8: the maximum observed SF of extrinsic variations tends to be lower when the lensed image is further from the lensing galaxy.

Fig. 8. Maximum observed SF of the microlensing curves, plotted against the distance to the lens of the lensed image divided by the effective radius of the lensing galaxy.

Next, the quasars’ ‘intrinsic’ curves were built by aligning all the curves with their respective delay, and modulating them with the smoothest extrinsic variability model used in the time-delay estimation so that they fitted the curve of the reference image. The SFs were then calculated after Eq. (11); they are displayed in Fig. 9. We see that the intrinsic variations in the source quasar (in black) are mostly dominant at short lag times, a regime that can only be probed with high-frequency and -precision monitoring. However, extrinsic variation can also operate on short timescales. An example is HE J0230−2130, in which a sudden rise was observed in image A at around MJD 60 275 (Fig. B.5). This event is visible in the SFs as well, with that of A rising at much lower lag times than that of the other lensed images. Fast events observed in a single light curve only are common, and while they could enable source science by combining light curves in multiple colours and a mass model of the lens, they also complicate the time-delay estimation, especially in regimes with lower sampling frequencies.

Fig. 9. Structure functions in the observed frame, of the best estimate of intrinsic variability of the source quasars (in black), as well as that of the extrinsic variations mainly due to microlensing. The absence of markers at certain lag times indicate that any potential variation was below the noise level. The error bars were estimated by bootstrapping. The SF of the extrinsic variations only catches up with the intrinsic SF at later time lags, which this program of high frequency monitoring takes advantage of by being able to probe lag times below 10 days. The SF is cut at 200 days, the characteristic duration of our monitoring seasons.

6.2. Comparison with LSST

In terms of depth, this program has remarkable similarities to the upcoming LSST. The 5-σ depth of a single LSST visit (two back-to-back 15 seconds exposures) is expected to be about 24.5 mag in the r band (Bianco et al. 2022), similar to the empirical depth measured on single visits of this program (24.5 mag, in 4 × 320 seconds exposures). Overall, this indicates that the per-visit photometric precision and performance in deblending can be expected to be identical. Next, the total number of visits including a given target is expected to be about 180 in the r band for LSST over the 10-years operation period. The mean total number of visits per target of this program is 227, ranging from 85 (WFI J2026−4536) to 435 (GRAL J1131−4419). These are again comparable figures, with LSST falling slightly short but making up for it by providing this number of visits multiple times in different filters. The main difference, however, is cadence, with visits in a given filter spaced by about 20 days on average for LSST, compared to the much faster one-day cadence of this program. Looking at the SFs of Fig. 9, we see that even with 20-day sampling, the intrinsic variability of our pool of lensed quasars is still much higher than that of microlensing in about half the sample, likely permitting the capture of large and fast variations in multiple images.

The other half will fall into a regime similar to that of the original COSMOGRAIL strategy, in which the contributions of extrinsic and intrinsic variation can only be disentangled thanks to the longer baseline. LSST will still have the advantage of observing in multiple filters: with the prior knowledge of how intrinsic variation transforms from band to band (shift and distortion, see, e.g., Chan et al. 2020), one might be able to disentangle the features due to actual intrinsic variation (to be matched when shifting curves), and those due to extrinsic variation (to be ignored for the purpose of time-delay estimation). However, this would introduce a dependence on our understanding of the geometry of the source quasar to hand. It would also require additional work on the side of time-delay estimation and modelling and would not be as straightforward as the shifts permitted by the high-frequency monitoring of this program.

6.3. Astrometric precision and lens models

Currently, high-quality time series ground imaging is a lot more expensive than high-resolution imaging. This will not hold true any longer once LSST starts operations – in the southern sky at least. Thus, being able to vet lensed quasars from ground-based time series imaging – that is, evaluate their suitability for cosmography or other applications – would be useful. For this, we need (i) to detect lensing features that can be used to constrain a lens model, which we showed in Fig. 6, and (ii) accurate astrometry for preliminary lens modelling. We compare the astrometry of the lensed images derived from our high-resolution model fitting with measurements from HST imaging data reported in Table 3 of Schmidt et al. (2023). We find a root-mean-squared deviation of 11 milliarcseconds (mas), which roughly matches the uncertainties reported by Schmidt et al. (2023) (6 mas per axis). In this sample of lenses, only one problematic case was found: PS J0030−1525, for which the position of image D was off by 100 mas. This was caused by the narrow separation between the bright image A (r ∼ 19.3) and faint image D (r ∼ 23.2), as well as contamination by the bright lensing galaxy. Thus, our method achieves a level of precision comparable to HST imaging in non-pathological cases – the pathological ones being easily identifiable.

We also suggest the possibility of fitting lens models directly on the ground imaging time series data. This is achievable by connecting, for example, STARRED and Herculens (Galan et al. 2022): given a lensing model, Herculens would predict the lens-plane image, while STARRED would provide a chi-squared to minimise by comparison with the dithered time series of imaging, given the PSF of each frame. Because both code bases are autodifferentiable, the gradient of the lensing parameters can be back-propagated from the chi-squared to the lensing parameters. This strategy will be the subject of future work.

6.4. Identifying lensed quasars through matched variability

The infrastructure of light curve extraction presented in this work is largely automated and can scale up to a large number of blended targets. It automatically selects calibration stars from Gaia, adapts to those available in each epoch, and accounts for frame rotation and PSF distortion, making it ready for the future LSST. Thus, for lensed quasar searches where spectroscopic follow-up is the limiting factor, the present infrastructure can leverage imaging time series to create pure samples of lensed quasars by selecting candidates exhibiting identical (modulo time delay) variations in all their point-like components.

7. Conclusions

We presented new time delays in 22 lensed quasars, a sample constituting the faintest and most narrowly separated targets for which time delays could be measured to date. Of the 22 systems, 8 are TDC-ready with a time-delay error budget below 5%, and 9 will likely contribute errors comparable to those of the other components of a TDC analysis. Regarding infrastructure, the improved photometric pipeline enables an extraction of light curves limited by photon noise only, for separations between lensed images down to 0$\stackrel{″}{.}$5, which is an improvement over the work presented in Millon et al. (2020a) in which deblending errors still play a role for lower-separation objects. Previously measured delays from shorter light curves have also been confirmed, with us finding compatible results even in unlucky cases in which the first season was strongly affected by strong microlensing. This demonstrates the robustness of our time-delay measurements and error estimation techniques. We have also shown how the careful matching of the zero points of the wide-field frames and forward-fitting the pixels can yield high-resolution models of the targets, compatible with HST observations in morphology and astrometry. In the times of LSST, the capacity for high-resolution follow-up imaging will be dwarfed by the number of unveiled targets, but each will have been imaged from the ground hundreds of times. Thus, the ability to reliably convert a deep, seeing-limited survey to a high-resolution image will become very valuable.

Data availability

Light curve tables are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/697/A139

Acknowledgments

The access to the 2.2m telescope was made possible through an agreement with the Max Planck Institute for Astronomy (MPIA): we would like to thank Thomas Henning, Hans-Walter Rix, and Roland Gredel for their help in making the agreement possible. We thank the European Southern Observatory (ESO) for accommodating the use of the VLT Survey Telescope (VST) amid the COVID-19 pandemic constraints. COSMOGRAIL is supported by the Swiss National Science Foundation (SNSF) and by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (COSMICLENS: grant agreement No 787886). This project has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Sklodovska-Curie grant agreement No 101105725. Paul Schechter generously allowed us to add his weekly-spaced points of J0924+0219 captured at the 2p2 to our time-delay estimation. The data analysis was made accessible by the Python ecosystem, relying heavily on libraries such as Numpy (Harris et al. 2020), Scipy (Virtanen et al. 2020), Astropy (Astropy Collaboration 2013, 2018), and Matplotlib (Hunter 2007). Gaia: This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. TA acknowledges support from ANID-FONDECYT Regular Project 1240105, ANID Millennium Science Initiative AIM23-0001 and ANID BASAL project FB210003. Support for MR is provided by the Dirección de Investigación of the Universidad Católica de la Santísima Concepción with the project DIREG 10/2023. MM acknowledges support by the SNSF (Swiss National Science Foundation) through mobility grant P500PT_203114 and return CH grant P5R5PT_225598. CDF acknowledges support for this work from the National Science Foundation under Grant No. AST-2407278 and the UC Davis College of Letters and Sciences. V.M. acknowledges support from ANID FONDECYT Regular grant number 1231418, Millennium Science Initiative, AIM23-0001, and Centro de Astrofísica de Valparaíso CIDI N21.

References

Appendix A: Notes about individual objects and time delays

Below we comment on the data, difficulties and results regarding each target. Note the convention of our time delays: a positive delay implies the second image arrives first. For example, A B = 10 days means that a given feature will first be seen in B, and will arrive 10 days later in A. We also estimate how the uncertainties in the obtained delays will map into H₀ by providing a benchmark precision: that is, we take all the delays involving one lensed image, and combine their relative precision assuming Gaussian uncertainties. The ideal scenario is reaching a below 5% combined precision, as this is where the time-delay errors become subdominant compared to that of the modelling and line-of-sight components of TDC (Suyu et al. 2014, 2017; Wong et al. 2020). Note that this only is a benchmark estimation of the precision that will map to H₀, while in reality the covariance between delays will also play a role in the cosmographic inference (see Sect. 5).

DES J0029−3814 This quadruply lensed quasar (Schechter et al. in prep.) has source redshift z_s ∼ 2.81. It was monitored at the 2p2 telescope from May 2019 to January 2020. The resulting one-season-long light curves, visible in Fig. B.1, have well defined long term variations which allow a relatively precise alignment. The low frequency nature of the observed variation makes it degenerate with potential microlensing, such that different freedoms afforded to the extrinsic modulations yield slightly different solutions. Some higher frequency variations are also observed, but their amplitude close to the noise level makes them unable to significantly break the microlensing degeneracy. This is reflected in the uncertainties we claim on the different delays. Overall, as expected because of the symmetry axis of the system, we observe a very small delay between B and C, and a larger delay between D and the trio A, B and C. The best precision on a single delay is 7% for B D, and the combined precision estimate that will map to H₀ is 3%, combining the delays relative to D. We provide the delay values and covariance matrix in Table A.1, below.

PS J0030−1525 This z_s = 3.36 quadruply lensed quasar was discovered by Lemon et al. (2018). We collected one season worth of imaging data from June 2021 to January 2022 with the 2p2 telescope, and even though it looks like a doubly imaged quasar only in the ground-based imaging, the forward modelling allows for reliable deblending of the lensing galaxy and the four images, revealing even the faint galactic satellite to the East of the lens. Due to the low separation in the folding pair however, we can only perform reasonable photometry of the sum of their fluxes. This is not a problem for TDC as their difference in time of arrival is expected to be very small. The data results in high-quality light curves with sharp features (Fig. B.2), even though the quality of the C curve is degraded by the photon noise of the lensing galaxy. A robust determination of the time delay between AD, the folding pair, and B and C could be conducted thanks to the two relatively high frequency peaks observed in all the curves. With a 50 days timescale of the intrinsic variations, the time delays are insensitive to the freedom given to the extrinsic variation models of each curve. Overall, the time-delay uncertainties are dominated by the photon noise of the curves. The best precision on a single delay is 12% for B C, and the combined precision estimate that will map to H₀ is 7%, combining the delays relative to C.

DES J0053−2012 Discovered by Lemon et al. (2020), this is a quadruply lensed quasar with source redshift z_s ∼ 3.8. Four light curves were extracted from a full season at the 2p2, starting in June 2021 and ending in January 2022, resulting in the light curves shown in Fig. B.3. The three brightest images, A, B and C, show enough high-frequency variation to reliably constrain time delays, all insensitive to how free the extrinsic modulation is. The D curve, however, is extremely noisy due to the faintness of image D, whose r-magnitude lies beyond 21.9. From prior knowledge that the A D delay should be quite large and negative, it appears that the best match is around A D ∼ −90 days. No strong microlensing is observed in any of the images as the curves can be matched within their noise with only time and magnitude shifts, but D is both the noisiest and the closest lensed image to the lensing galaxy: the most prone to undetectable extrinsic variation. Thus, even though more flexible microlensing models were included in the estimation of the error bars of the delays, those involving D should still be taken with care. So, excluding D, the best precision on a single delay is 10% for A B, and the combined precision estimate that will map to H₀ is 6%, combining the delays relative to A.

WG J0214−2105 Discovered by Spiniello et al. (2019), this quadruply imaged quasar has source redshift z_s ∼ 3.24. This object was monitored at the 2p2 for two seasons (June 2018 – February 2019, August 2019 – February 2020), and then for one more season at the VST (August 2021 – February 2022). The two latter seasons were shorter than the first due to the COVID-19 pandemic. The 2p2 data were already presented in Millon et al. (2020b), but the curves were extracted again with the new method. These and the new ones are shown in Fig. B.4. We provide a re-analysis together with the new VST curves, improving the precision and changing the value of the A B delay slightly. The A D delay is also changed (although by less than 2-σ upon quadratically combining error bars), possibly due to the more flexible microlensing models allowed for in the present estimation. Due to combining three seasons worth of data, the microlensing effects become very noticeable, such that splines with several internal knots (up to 10 for curve D, whose shape is qualitatively different from the other curves) modulating each curve were needed to properly align the light curves. Nevertheless, the many high-frequency and S/N features observed in the light curves permitted a precise alignment. The best precision on a single delay is 5% for B C, and the combined precision estimate that will map to H₀ is 2%, combining the delays relative to C.

HE 0230−2130 Discovered by Wisotzki et al. (1999), this is a quadruply lensed quasar with source redshift z_s ∼ 2.162 and lens redshift z_l ∼ 0.523 (Eigenbrod et al. 2006). It was monitored for a short season (started late due to the pandemic), from December 2020 to February 2021. Next, another season was attempted, starting in August 2021 but quickly interrupted the next month, again due to the pandemic. The resulting curves are visible in Fig. B.5. The first short season could constrain the time delays quite well, but it was clear that additional data would help: thus, it was scheduled again at the 2p2 starting in June 2023 and ending in January 2024. During this last 2p2 season, a clear microlensing event affected the A image, making the estimation of the time delays involving this image more difficult: the alignment of the curves required very flexible extrinsic modulations of A. This is reflected in the uncertainties of the A B, A C and A D delays. This object has also been the subject of a two- and a three-years monitoring campaigns at the Euler Swiss telescope, the results of which were presented in Millon et al. (2020a). Thanks to the higher S/N provided by the 2p2 telescope compared to the 1.2 m Euler Swiss Telescope, we were this time able to reliably deblend the A and B images, and to also provide the delays involving the fourth image, D. Our time-delay estimations are in slight tension with the measurement presented in Millon et al. (2020a): 1.2σ and 1.5σ for the A C and B C delays respectively (both referred to as A C in Millon et al. (2020a), as their A is the sum of our A and B). A re-analysis of the light curves presented in Millon et al. (2020a) by Donnan et al. (2021) provides much tighter error bars on the delays, and is in even more tension with the measurements we present herein from the new, deeper data. In particular, we find a tension of 2- and 3-σ respectively when comparing our A C and B C delays to their (17.7 ± 2.1) days A C value. Their A D and C D values are, also, roughly 2- and 3-σ away from ours. We explain this difference by the higher quality data used herein, with uncertainties that account for the possible degeneracies due to extrinsic modulation. Overall, with the delays of this work, the best precision on a single delay is 10% for C D, and the combined precision estimate that will map to H₀ is 4%, combining the delays relative to D.

SDSS J0248+1913 This quadruply imaged quasar has source redshift z_s ∼ 2.424, and was discovered both with the technique presented in Ostrovski et al. (2017) and with HST follow-up presented in Shajib et al. (2019). The system was also discovered independently by Delchambre et al. (2019). A monitoring season was started at the VST at the end of September 2021, but the pandemic effects quickly thwarted the attempt at the end of December. No time delay could be extracted from the resulting light curves, which we do not present due to the poor sampling and short duration. Nevertheless, the existence of 8 hours of OmegaCAM exposure time on this field is worth mentioning.

WISE J0259−1635 Discovered by Schechter et al. (2018), this quadruply imaged quasar has source redshift z_s ∼ 2.16. A first monitoring season was captured at the VST, starting in August 2021 and ending in March 2022. A second season was then taken at the 2p2, starting in June 2023 and ending in January 2024. High-precision time delays could be obtained thanks to the high S/Ns of the extracted curves, shown in Fig. B.6. The recorded features have timescales of 5-15 days, which is well separated from the ∼100 days observed extrinsic modulations attributable to microlensing. Thus, the time delays are insensitive to how flexible the extrinsic modulation models are, such that the uncertainty budget is mostly due to photon noise. The best precision on a single delay is 9% for B C, and the combined precision estimate that will map to H₀ is 4%, combining the delays relative to C.

J0420−4037 This quadruply lensed quasar has source redshift z_s ∼ 2.4. It was discovered (Ostrovski et al. in prep.) with the method described in Ostrovski et al. (2017). It was monitored at the VST for one season from October 2020 to March 2021, resulting in the four curves provided in Fig. B.7. Two of these show high S/N variations (A and B), while the other two are noisier, albeit with visible features still (C and D). The short delays in this lens, however, prevent us from reaching a high relative precision in any of the pairs of images: the only definitely incompatible with zero delays are those involving the D lensed image. The best precision on a single delay is 33% for A C, and the combined precision estimate that will map to H₀ is 11%, combining the delays relative to C.

DES J0602−4335 This doubly imaged quasar, discovered by Dawes et al. (2023), was monitored at the 2p2 from the end of November 2020 to the end of March 2021. The redshift of the source is z_s ∼ 2.92(1), determined from archival ESO EFOSC2 data.⁷ The lens is an edge-on galaxy, which the forward-modelling could readily deblend from image A. Even though the extracted light curves (Fig. B.8) are only 120 days long, two rises with timescale ∼20 days were captured in both curves. These can hardly be degenerate with extrinsic variations given their short timescale, and as such, the measured delays are not sensitive to the freedom given to the extrinsic modulation. The resulting time delay is therefore very precise, its uncertainty mostly being due to the photon noise of the curves. The precision is here 6% on the single delay, A B.

J0607−2152 Discovered by Stern et al. (2021) and Lemon et al. (2023), this quadruply imaged quasar has source redshift z_s ∼ 1.302. This object was first monitored for a short season at the VST, starting in October 2020 and ending in March 2021. Monitoring was then resumed for an additional season, from August 2021 to April 2022, resulting in the light curves displayed in Fig. B.9. The light curve of the B image was unfortunately too noisy to produce a reliable estimation of a time delay, due to its faintness (22.1 r-mag), but also to the presence of a very bright star in the field of view (YY Leporis), adding to the photon noise and degrading the quality of the background subtraction. A, C and D however, all show sharp features that could be aligned without degeneracies with potential extrinsic modulations. Overall, the best precision on a single delay is 12% for C D, and the combined precision estimate that will map to H₀ is 7%, combining the delays relative to D. Additionally, the high-resolution model shows a bright Einstein ring with visible structure, making this target a good candidate for deep, sharp follow-up imaging.⁸ Moreover, Lemon et al. (2023) suggest that a single galaxy model is not sufficient to explain the lensing configuration – indeed our fitted model shows a potential second galaxy very close to image B that supports their suggestion. We can also confirm the image configuration suggested by Lemon et al. (2023) over that of Stern et al. (2021) who suggested image B was between A and D.

J0659+1629 Discovered by Delchambre et al. (2019) and Lemon et al. (2023), this is a quadruply lensed quasar with source redshift z_s ∼ 3.09 and lens redshift z_l ∼ 0.766 (Stern et al. 2021). This object was observed for two seasons at the VST: a short one from November 2020 to March 2021, and another from October 2021 to April 2022. We also include archival data from the Las Cumbres Observatory⁹, with a mean cadence of 2 days, ranging from December 2019 to May 2020. The resulting light curves (Fig. B.10) can be aligned without strong extrinsic modulations, except for that of image C: this is not surprising considering its proximity to a second galaxy, which both the HST imaging and high-resolution model fitted on VST data reveal. Once the required freedom is given to the C curve, the curves can be aligned within their noise envelope with a preferred solution involving D ahead by about 270 days. However, it is clear that a 330 days shift of D could fit the curves almost as well, albeit involving stronger extrinsic modulations. Moreover, our data cannot constrain potential delays beyond 400 days. We provide in Table A.10 delays and uncertainties around the best solution, but one should keep in mind the second possible solution (Table A.11), and the possibility that the delays are much bigger. To use this lens in a TDC study, one should carefully confront the mass model predictions of the time delays with the different solutions we find herein, and only move forward if one of the possibilities is favoured beyond all doubts by the mass model. Because J0659+1629 will land in the LSST footprint, there is hope for additional constraints on the delay in the long run.

SDSS J0832+0404 This doubly imaged quasar was discovered by Oguri et al. (2008). It has source redshift z_s ∼ 1.115. It was monitored at the 2p2 from the end of November 2017 to the beginning of June 2018. The extracted curves (Fig. B.11) show plenty of high-frequency features. However, due to the long delay, the overlap is only ∼70 days. This makes the time delay sensitive to the local slope of the extrinsic variation in each lensed image, which is reflected in the uncertainty on the time delay. The obtained value is, however, perfectly compatible with that reported by Millon et al. (2020a) at the Euler Swiss telescope. It is also compatible with the re-analysis of the Swiss telescope curves by Donnan et al. (2021). Combining the estimation of Millon et al. (2020a) with ours as independent measurements, assuming Gaussian uncertainties, we obtain a precision of 6%.

RX J0911+0551 Discovered by Bade et al. (1997), this is a quadruply imaged quasar with source redshift z_s ∼ 2.763 and lens redshift z_l ∼ 0.769 (Kneib et al. 2000). It was monitored at the 2p2 for two seasons, from October 2020 to March 2021 and from November 2021 to June 2022. The resulting curves are displayed in Fig. B.12. The ABCD delay is known to be 150 ± 6 or 146 ± 4 days (Eulaers & Magain 2011): a precision sufficient for TDC even with the potential multi-modality. Two high cadence monitoring seasons should have been enough to much better constrain the D delay (due to a season of D overlapping both the first and second season), but the pandemic interruption of the first season leaves us with only one, short overlap. This raises a similar problem to the case of SDSS J0832+0404, making the individual delays imprecise. Nevertheless, we do have three delays involving D thanks to the A, B and C trio being reliably deblended by the forward-modelling. Thus, the best precision on a single delay is 12% for A D, and the combined precision estimate that will map to H₀ is 4%, combining the delays relative to D.

SDSS J0924+0219 Discovered by Inada et al. (2003), this is a quadruply lensed quasar with source redshift z_s ∼ 1.685 and lens redshift z_l ∼ 0.393. Monitoring started with a daily cadence in October 2020 at the VST, until March 2021. Then, three additional weekly cadenced seasons were acquired for microlensing studies. Two at the VST still, from December 2020 to April 2021 and from November 2021 to June 2022, and the third one at the 2p2 from November 2022 to June 2023. Shown in Fig. B.13, these weekly cadenced seasons propitiously coincided with fluctuations in the A and B images, providing a fast enough sampling to capture the variations appropriately. The C curve is much noisier, but its first season shows a definite oscillation that can be matched to A and B. The present A B estimation is more precise than, but also compatible with, the one presented in Millon et al. (2020b). A C on the other hand, was not given in Millon et al. (2020b) due to insufficient S/N, but was in the reanalysis by Donnan et al. (2021). Their value (-30.9 days) is in slight (1.5σ) tension with the one presented herein. The best precision on a single delay is 16% for B C, and the combined precision estimate that will map to H₀ is 7%, combining the delays relative to C.

GRAL J1131−4419 Discovered by Krone-Martins et al. (2018), this is a quadruply imaged quasar with source redshift z_s = 1.090 (Wertz et al. 2019). It was first monitored at the 2p2 for a short season starting in December 2019, but quickly interrupted by the COVID-19 pandemic in March 2020. Next, another short season was acquired at the VST starting in November 2020 and ending in April 2021. Finally, two additional seasons could be fully completed: a first one at the VST starting in November 2021 and ending in July 2022, and the second at the 2p2 starting in November 2022 and ending in June 2023. The narrow separation (∼0$\stackrel{″}{.}$45) of the A-B pair made it impossible to reliably deblend in most seeing conditions, so we only present the delays within the trio AB, C, and D, with AB the sum of the A and B fluxes. The resulting curves, AB, C and D are provided in Fig. B.14. Thanks to rapid, high-amplitude intrinsic variations and the availability of four monitoring seasons, excellent precision could be achieved in the delays. The best precision on a single delay is 7% for C D, and the combined precision estimate that will map to H₀ is 4%, combining the delays relative to D.

2M 1310−1714 Discovered by Lucey et al. (2018), this is a quadruply lensed quasar with source redshift z_s ∼ 1.975 and double-galaxy lens redshift z_l ∼ 0.293. This object benefited from three monitoring seasons at the VST. The first two, starting in December 2019 and 2020 respectively, both ended early in March of the following year because of the COVID-19 pandemic. Luckily, a third season, complete this time, was acquired starting in early January 2022 and ending in July 2022. High-S/N light curves (Fig. B.15) were extracted from the three brightest images, A, B and C. The light curve of the D image became noisier as monitoring progressed, as it got dimmer with time, reaching 21.6 mag by the end of the last VST season. Despite this, it still shows features that can be used to match it with the other three curves. The reached precision in the time delays is excellent thanks to the sharp observed variations, availability of three seasons, and moderate extrinsic modulation needed to align the light curves. The best precision on a single delay is 3% for A B, and the combined precision estimate that will map to H₀ is 1%, combining the delays relative to A.

J1537−3010 Discovered by Lemon et al. (2019) and Delchambre et al. (2019), this quadruply imaged quasar has source redshift z_s ∼ 1.721. This object was first monitored at the 2p2 before the COVID-19 pandemic, for a full season, from February to October 2019. Next, two seasons were begun, at the 2p2 in February 2020 and at the VST in February 2021, but both were quickly interrupted due to the pandemic. Finally, it was monitored for a final full season at the 2p2 again in 2022, starting in January and ending in September. The extracted light curves are those shown above, in Fig. 5. The combined ‘two and a half’ seasons provide the most stringent time-delay constraints obtained in this work, thanks to the source quasar’s large and fast intrinsic variations which are easy to separate from the microlensing signal. The best precision on a single delay is 2% for B C, and the combined precision estimate that will map to H₀ is 1%, combining the delays relative to B.

PS J1606−2333 Discovered by Lemon et al. (2018), this quadruply lensed quasar has source redshift z_s ∼ 1.69. It was first monitored for a season at the 2p2 from February to September 2018. This first season was already presented in Millon et al. (2020b), but did not allow the reliable determination of the time delays involving the D lensed image. Thus, another season was planned at the VST, starting in late February 2022 until late September 2022. Note that the curves were extracted from the 2p2 data again with the new method – the combined result is displayed in Fig. B.16. Together, the 2p2 and VST curves can this time constrain the D delay thanks to the higher frequency variations observed in the VST data. The delays within the A, B and C trio determined with the 2p2 data by Millon et al. (2020b) are confirmed and made more precise thanks to this added season of data and re-analysis. The best precision on a single delay is 5% for A C, and the combined precision estimate that will map to H₀ is 2%, combining the delays relative to A.

WGD J2021−4115 This doubly imaged quasar, discovered by Agnello et al. (2018) has source redshift z_s ∼ 1.390 and lens redshift z_l ∼ 0.335. Its monitoring at the 2p2 telescope from April to September 2019 resulted in two well-defined light curves (Fig. B.17), albeit with slow variations somewhat degenerate with eventual microlensing, mainly due to the short overlaps between curves resulting from the relatively long -91 days delay. The delay uncertainty would thereby dominate the error budget of an H₀ estimation with TDC. However, the high-resolution fitted model does reveal signs of an Einstein ring, essential for constraining a precise mass model of the lens. This could motivate the push for obtaining deep, sharp imaging of the object despite it being only doubly imaged.¹⁰ Furthermore, the long delay, the brightness of the two images (r-mag 19.5 and 19.8) and their wide separation (2$\stackrel{″}{.}$7), all constitute ingredients that will likely permit a time-delay determination with LSST. Overall, we expect this object to become a compelling system for TDC after a few years of LSST operation, despite the relatively imprecise delay obtained herein. The precision on the single delay, A B, is 11%.

WFI J2026−4536 Discovered by Morgan et al. (2004), this quadruply imaged quasar has source redshift z_s ∼ 2.237. It was monitored at the VST for a short season between April and September 2022. The two brightest images, A and D, are too narrowly separated (0$\stackrel{″}{.}$3) to be reliably deblended. Their fluxes were thereby summed. Even though high-S/N features are present in the resulting set of light curves (Fig. B.18), patches of bad weather made the time-delay estimation degenerate: possible solutions include the overlapping of the features of one curve with the empty regions of another. This is reflected in the high uncertainty in the delays, for example the AD C delay: determined with a relative precision of more than 50%, it is nonetheless perfectly compatible with the value reported by Millon et al. (2020a). The best precision on a single delay is 19% for AD B, and the combined precision estimate that will map to H₀ is 12%, combining the delays relative to B.

WG J2038−4008 Discovered by Agnello et al. (2018), this quadruply lensed quasar has source redshift z_s ∼ 0.777 and lens redshift z_l ∼ 0.228 (Stern et al. 2021). It was monitored at the 2p2 telescope for three seasons: from April to December 2017, from June to December 2021, and from April to December 2022. It was also monitored for five seasons at the Euler Swiss Telescope, a 1.2m facility at the same observatory, with a two-day cadence. When double coverage occurred, only the 2p2 data were used in the time-delay estimation due to its much higher quality and cadence. Of note, the curves resulting from this dataset were already presented in Wong et al. (2024); a TDC analysis of the system. Not only is this lens on the faint side of the targets monitored during this campaign, it also is the one whose source quasar varied the least on the monitoring baseline. Fortunately, the combined datasets yield curves (Fig. B.19) that provide meaningful constraining power, thanks to the image D lagging behind the others by a longer time (about 30 days). An additional problem was encountered with the A C delay however, for which two solutions were possible: one positive (7.9 days), and the other negative (-5.3 days). The negative one was favoured by the data, and, more importantly, switching solution changes the ordering of arrival of the images. The ordering of arrival is very well constrained by mass models of the lens, and allowed us to further prefer the negative solution. Overall, the best precision on a single delay is 18% for A D, and the combined precision estimate that would map to H₀ is 9%, combining the delays relative to D. However, the total error budget in the TDC analysis by Wong et al. (2024) was roughly 25–30%, showing how correlation between delays might have played a role in degrading the precision of the final estimate.

WG J2100−4452 Discovered by Agnello & Spiniello (2019), this is a quadruply lensed quasar with source redshift z_s ∼ 0.92 and lens redshift z_l ∼ 0.203 (Spiniello et al. 2019). It was first monitored at the VST for a season, between April and November 2019. Next, it was observed simultaneously by both the 2p2 and VST for an additional season, with combined coverage starting in April 2022 and ending in December of the same year. The double coverage occurred due to scheduling constraints and the uncertainty due to the COVID-19 pandemic, and proved beneficial given that the 2p2 monitoring was interrupted for almost two months due to a snowstorm at La Silla observatory. The extracted light curves are displayed in Fig. B.20. The VST and 2p2 points are differentiated by the colour of their error bars: gray for VST, matching colour for 2p2. This double coverage also offers a reassuring test of the photometric extraction, as curves produced by the two facilities match within their noise after a slight magnitude offset due to the mismatch in filter. However, the two brightest and closest images A and B were too narrowly separated (0$\stackrel{″}{.}$53) to allow for a reliable deblending. Thus, their fluxes were summed, yielding the high-S/N AB light curve. C and D are much weaker, with light curves consequently much noisier. Fortunately, the oscillations starting at MJD 58700 are still well visible in both B and C as visible in Fig. B.20. These have timescales way below those we allow for our microlensing modulation, as such the estimated delays are not overly sensitive to the freedom granted to the extrinsic modulation. Thus, mainly due to the low S/N of the C and D curves, we obtain delays whose uncertainty budget is likely dominated by photon noise. The best precision on a single delay is 12% for C D. Combining all delays relative to C, the best possible precision mapping to H₀ would be 8%.

Fig. A.1. Light curves of DES J0029−3814

J2205−3727 Discovered by Lemon et al. (2023), this is a quadruply imaged quasar with source redshift z_s ∼ 1.848. This is the faintest lens monitored by this program: image C in particular, with an average r-magnitude of 22.6, is the faintest quasar lensed image for which a time delay could be estimated. J2205−3727 was monitored for two seasons at the VST, the first one started in April 2019 and ended in December 2019, and the second from May to August 2022. The resulting light curves are shown in Fig. B.21. The curves of images A, B, and D all exhibit short-timescale features that can be used for alignment, but the flux of image C was often compatible with zero, especially on nights with poor seeing. Consequently, data points of the C curve from nights with seeing values exceeding 1.3 were discarded. After this filtering was made, the C curve was fortunately still well sampled enough to reveal the patterns seen in the other three curves, but a part of its constraining power was lost. This is well reflected in our estimated time-delay uncertainties, which are much larger in the delays involving it. Overall, we find time delays well compatible with zero within the A, B, C trio, but longer delays in the pairs involving D, the image across from the lensing galaxy. The error budget in the delays is mostly dominated by photon noise, as the estimated delays were insensitive to the complexity of the allowed extrinsic modulations. The best precision on a single delay is 15% for A D, and the combined precision estimate that will map to H₀ is 7%, combining the delays relative to D.

Appendix B: Light curves

Here we list the light curves of the monitored targets.

Fig. B.1. Light curves of PS J0030−1525.

Fig. B.2. Light curves of DES J0053−2012.

Fig. B.3. Light curves of WG J0214−2105.

Fig. B.4. Light curves of HE J0230−2130.

Fig. B.5. Light curves of ATLAS J0259−1635.

Fig. B.6. Light curves of DES J0420−4037.

Fig. B.7. Light curves of DES J0602−4335.

Fig. B.8. Light curves of J0607−2152.

Fig. B.9. Light curves of J0659+1629.

Fig. B.10. Light curves of J0832+0404.

Fig. B.11. Light curves of RX J0911+0551.

Fig. B.12. Light curves of J0924+0219.

Fig. B.13. Light curves of GRAL J1131−4419.

Fig. B.14. Light curves of 2M J1310−1714

Fig. B.15. Light curves of PS J1606−2333.

Fig. B.16. Light curves of WGD J2021−4115.

Fig. B.17. Light curves of WFI J2026−4536.

Fig. B.18. Light curves of WG J2038−4008.

Fig. B.19. Light curves of WG J2100−4452 This lens is a special case, it was monitored at the same time both at the 2p2 and at the VST due to scheduling constraints. The VST points are differentiated with gray error bars.

Fig. B.20. Light curves of J2205−3727.

Appendix C: Observational statistics

All Tables

All Figures

	Fig. 1. Astronomical seeing distribution in the images used to infer time delays in this data release.
In the text

Fig. 2. Example simultaneous fit of a PSF model to ten stars, in a (rare) difficult case of unfocused observation. The frame chosen for this figure could still be used for photometry, not leading to an outlier in the light curves. The kernel, at the top left, transforms a two-pixel FWHM Gaussian into the PSF at the bottom left. Note that for our lens photometry we used the kernels and not the PSFs (see Sect. 3.6 and Fig. 3). The bottom row shows the residuals of the joint fit. Similar residuals were obtained in all conditions of seeing, focusing, and tracking.

In the text

Fig. 3. Overview of the final step of light curve extraction: the forward-modelling of the lens cutouts. Note that the colour coding is only used to differentiate point-source and pixelated components of the models, and does not indicate observations in multiple filters. Top: Components of the forward-modelling process of the zero-point-calibrated cutouts, starting with a high-resolution model containing a fully tunable grid of pixels and point sources of tunable flux (added to the grid of pixels as two-pixel FWHM Gaussians). The model is degraded to the resolution of the data, and the residuals are minimised through gradient-based optimisation of the high resolution models. Bottom: Exploring the result of the forward modelling for a selected epoch. A is the fitted high-resolution model. B is the kernel, which transforms a point source in A into the PSF of the data. C is the convolution of A and B, which is then downsampled down to the pixel scale of the data, yielding D. A, C, and D have the background component in red and the point-source component in white. E is the data cutout. F and G are again the data, but with the background and point-source components subtracted, respectively. Finally, H shows the residuals upon subtraction of the two components from the data. This method allows for the precise determination of the astrometry of the point sources, reveals the morphology of the non-point source components at high resolution, and provides precise photometry unbiased by flux leakage between point source and non-point source components.

In the text

Fig. 4. Average nightly scatter of the measured lensed-image magnitudes, plotted against the mean magnitude of the lensed image at hand. The scatters are empirical, and thereby include photon noise, read noise, deblending errors, and normalisation errors between frames. The dashed line is the noise estimate given by the ESO/WFI exposure time calculator (ETC), for a point source in the r band in good seeing conditions (1″), average moon illumination, and 320 second exposures. The similarity between idealised and empirical errors indicates that our procedure has an excellent deblending performance. For reference, the line fit in blue indicates the measured photometric performance given the magnitude of a given target point source.

In the text

	Fig. 5. Light curves extracted from the imaging data of J1537−3010 (labelling of the lensed images in Fig. 6), monitored at both observatories and spanning four seasons, with the 2020 and 2021 ones shortened by the COVID-19 pandemic. Note that the empty season gaps are cut from the plot. A magnitude offset was added to the individual curves as is indicated in the legend, for display purposes.
In the text

Fig. 6. Imaging preview of all the monitored objects. For each lens, we include a representative exposure from the dataset with median seeing and low sky background, the high-resolution fitted model (see Sect. 3.6), and, when available, an HST/Wide Field Camera 3 (WFC3) image for comparison. The high-resolution models live on a grid twice-supersampled compared to the original images, with a two-pixel FWHM PSF. The resolution (FWHM of the PSF) in the models is thereby 0$\stackrel{″}{.}$213 for VST data, and 0$\stackrel{″}{.}$237 for 2p2 data. This is to be compared with the 0$\stackrel{″}{.}$07 resolution of the HST/WFC3 images. Note also that the F814W filter is much wider and redder than the r band, such that more details can be revealed in the redshifted objects of our fields. All cutouts are oriented with north facing up and east facing left.

In the text

Fig. 7. Curve matching example, in which a single spline (in black) matches all curves simultaneously. Due to microlensing in particular, which adds independent lower-frequency variations to each curve, the curves could not be matched in this way with time and magnitude shifts only. Thus, an extra modulation is applied to each separately. These modulations are represented by lower order splines, displayed at the bottom. Physically, the black spline models the intrinsic variations of the source, while the coloured splines undo the effects of microlensing. Linking with the notation of Sect. 4, Mext is here a cubic spline with a single internal knot of fixed position, and Nint = 9. To determine a time delay and its uncertainty, this matching was repeated thousands of times with artificially further-modulated curves (mimicking microlensing), other realisations of the noise, and different freedoms given to the splines. What this matching is like for the other targets can be found in a Jupyter notebook. This notebook is part of the repository containing the photometry extracted herein and the code to estimate the time delays.

In the text

	Fig. 8. Maximum observed SF of the microlensing curves, plotted against the distance to the lens of the lensed image divided by the effective radius of the lensing galaxy.
In the text

Fig. 9. Structure functions in the observed frame, of the best estimate of intrinsic variability of the source quasars (in black), as well as that of the extrinsic variations mainly due to microlensing. The absence of markers at certain lag times indicate that any potential variation was below the noise level. The error bars were estimated by bootstrapping. The SF of the extrinsic variations only catches up with the intrinsic SF at later time lags, which this program of high frequency monitoring takes advantage of by being able to probe lag times below 10 days. The SF is cut at 200 days, the characteristic duration of our monitoring seasons.

In the text

	Fig. A.1. Light curves of DES J0029−3814
In the text

	Fig. B.1. Light curves of PS J0030−1525.
In the text

	Fig. B.2. Light curves of DES J0053−2012.
In the text

	Fig. B.3. Light curves of WG J0214−2105.
In the text

	Fig. B.4. Light curves of HE J0230−2130.
In the text

	Fig. B.5. Light curves of ATLAS J0259−1635.
In the text

	Fig. B.6. Light curves of DES J0420−4037.
In the text

	Fig. B.7. Light curves of DES J0602−4335.
In the text

	Fig. B.8. Light curves of J0607−2152.
In the text

	Fig. B.9. Light curves of J0659+1629.
In the text

	Fig. B.10. Light curves of J0832+0404.
In the text

	Fig. B.11. Light curves of RX J0911+0551.
In the text

	Fig. B.12. Light curves of J0924+0219.
In the text

	Fig. B.13. Light curves of GRAL J1131−4419.
In the text

	Fig. B.14. Light curves of 2M J1310−1714
In the text

	Fig. B.15. Light curves of PS J1606−2333.
In the text

	Fig. B.16. Light curves of WGD J2021−4115.
In the text

	Fig. B.17. Light curves of WFI J2026−4536.
In the text

	Fig. B.18. Light curves of WG J2038−4008.
In the text

	Fig. B.19. Light curves of WG J2100−4452 This lens is a special case, it was monitored at the same time both at the 2p2 and at the VST due to scheduling constraints. The VST points are differentiated with gray error bars.
In the text

	Fig. B.20. Light curves of J2205−3727.
In the text