Improving the light curves of gravitationally lensed quasars with Gaia proper motion data

We show how to significantly improve difference image analysis (DIA) of gravitationally lensed quasars over long periods of time using Gaia proper motions. DIA requires the subtraction of a reference image from the individual images of a monitoring campaign, using stars in the field to align the images. Since the proper motion of the stars can be of the same order as the pixel size during a several-year campaign, we use Gaia DR3 proper motions to enable a correct image alignment. The proper motion corrected star positions can be aligned by the ISIS package. DIA is carried out using the HOTPAnTS package. We apply point spread function (PSF) photometry to obtain light curves and add a proper motion correction of the PSF star to GALFIT. We apply our method to the light curves of the three gravitationally lensed quasars HE1104-1805, HE2149-2745 and Q2237+0305 in the R and V band, respectively, obtained using 1 m telescopes of the Las Cumbres Observatory from 2014 to 2022. We show that the image alignment and the determination of the lensed quasar positions is significantly improved by this method. The light curves of individual quasar images display intrinsic quasar variations and are affected by chromatic microlensing.


Introduction
Multiply imaged quasars produced via the strong gravitational lensing effect are used for various applications in astrophysics and cosmology (see e.g. the recent review Shajib et al. 2022, and references therein).Since the first detection of such a system (Walsh et al. 1979), about 220 lensed quasars have been discovered to date 1 .Compact objects in the foreground lensing galaxy close to the line of sight to a quasar image can cause additional microlensing signals detectable as uncorrelated brightness variations in the light curves of the different quasar images (see Chang &Refsdal 1979 ande.g. reviews Schmidt &Wambsganss 2010;Vernardos et al. 2024).These microlensing signals can be used to test model predictions for the structure of quasars (Wambsganss & Paczynski 1991;Kochanek 2004;Anguita et al. 2008;Eigenbrod et al. 2008;Poindexter & Kochanek 2010;Morgan et al. 2018), such as the temperature profile of the accretion disk (Shakura & Sunyaev 1973).In this work, we present a method to obtain long-term quasar light curves using point spread function (PSF) photometry and difference imaging analysis (DIA) with a focus on improving the results by using Gaia proper motion data (Gaia Collaboration 2016, 2023).We note that it is also possible to determine proper motions directly from DIA (Skowron et al. 2014), but here, we use the Gaia values.We apply this method to three lensed quasars and determine their light curves using data taken at the Las Cumbres ⋆ R and V band light curves for all three sources are available at the CDS ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/683/A119 1 https://research.ast.cam.ac.uk/lensedquasars/Observatory (LCO2 , Brown et al. 2013) over 9 yr.In Sect.2, we summarise some properties of our targets before elaborating on our method in detail in Sect.3, present the resulting light curves in Sect. 4 and conclude in Sect. 5.

Targets
Here we state a few properties of the lensed quasars HE1104-1805, HE2149-2745 and Q2237+0305, which we analysed.

HE2149-2745
The quasar HE2149-2745, discovered by Wisotzki et al. (1996) in HES as well, with z Q = 2.033 and z L = 0.603 (Eigenbrod et al. 2007), also has two images that are separated by 1.70 arcsec.The time delay was determined to be ∆t AB = (103 ± 12) days by Burud et al. (2002), where image A leads.Eulaers & Magain (2011) reanalysed the same data, report an additional possible −9.3 days, which is the value we adopt.

Q2237+0305
Q2237+0305, also known as Einstein Cross, was discovered by Huchra et al. (1985).It is visible as four quasar images at z Q = 1.695 close to the center of a barred spiral galaxy known as Huchra's lens (see Fig. 1) at z L = 0.0394.The four images have angular separations of ≲ 1.8 arcsec.We use ∆t = 0 days for all time delays between the four images, since estimates (obtained with different techniques) for the multiple time delays are typically of the order of hours (see e.g.Schneider et al. 1988;Wambsganss & Paczynski 1994;Schmidt et al. 1998;Dai et al. 2003;Vakulik et al. 2006;Koptelova et al. 2006;Berdina & Tsvetkova 2018) in contrast to our light curves over 9 yr.This is the system in which quasar microlensing was first detected (Irwin et al. 1989;Corrigan et al. 1991).

Methods
We want to measure the light curves of the multiple images of three strongly lensed quasars HE1104-1805, HE2149-2745 and Q2237+0305 in the R and V bands, respectively, using data taken at LCO since 2014.LCO is a global network of robotic telescopes and we use data from the 1 m telescopes at four locations (mainly Cerro Tololo, Chile; but also: Sutherland, South Africa; Siding Springs, Australia; McDonald, USA).To extract the light curves from the data obtained, the four main steps (which we discuss in more detail in the following four subsections) are: Sect.3.1: align all images from LCO with respect to each other using the ISIS package 3 from Alard (2000), Sect.3.2: combine all images from the same night and produce a high S/N and small seeing reference image for the final step, tracting the reference image transformed with an appropriate convolution kernel from all combined images) using the HOTPAnTS software4 by Becker (2015) which implements the algorithm from Alard & Lupton (1998) and Alard (2000), and finally extract the light curves by again using GALFIT and therefore applying PSF photometry on the difference images.Overall, the method is similar to Giannini et al. (2017), but we now use Gaia data to improve the astrometry.We adapted the ISIS, GALFIT and HOTPAnTS codes to deal with specialities of our data such as low number of stars, the multiple images of our quasars, as well as to include improvements using Gaia proper motion data.All steps of our reduction are implemented in python or are python codes that run ISIS, GALFIT and HOT-PAnTS on our data.The modified codes (except for GALFIT) and the python codes are available on GitHub5 .

Improving the image alignment with Gaia data
In order to successfully combine the single images of each night and to apply difference imaging on these combined images with respect to a reference image, all the single images have to be aligned.To achieve this, ISIS determines the positions of N ∼ 500 to 1000 stars in all images of one quasar in one band and with that aligns the images relative to each other by applying small shifts, rotations and scalings to the single images in order to minimize the sum of the squared offsets.But since our images are taken up to ∼9 yr apart, the stars positions change due to their proper motion, which makes it difficult for ISIS to achieve proper image alignment.
Here we present our method to improve the alignment significantly.We want to emphasise that the method will be widely applicable and useful in time domain astronomy.
ISIS identifies as many stars as possible in an astrometric reference image and determines all their positions r ref i .It then does the same i.e. determining the star positions r i for all images and minimises all the star deviations to the reference image ∆r i = r i − r ref i for each image and corrects that image accordingly.The star positions relative to the reference image are (additionally to image misalignment) influenced by proper motions and the time difference between the current image and the reference image ∆t = t − t ref .Since we choose a high signal-to-noise (S/N) and low seeing image of the sample (typically ∼1.3 arcsec) as reference image, it is of better quality as most other images, which is why we use it to apply the following correction.
We determine the proper motions from Gaia data release 3 (Gaia Collaboration 2023, DR36 ) of stars in our fields with a Renormalised Unit Weight Error RUWE < 1.4 (Lindegren 2018) and proper motion values with S /N ≥ 5 to correct their positions in our reference image in the calculations of ISIS.As preparation, the stars in the neighbourhood of the quasar from the Gaia DR3 data are matched to the stars ISIS detects and intends to use for the alignment (see Fig. 1).This gives us a list of star positions as detected by ISIS, but additionally with their Gaia proper motions in the local tangent plane in right ascension µ α * ,i = µ α,i cos δ i and declination µ δ,i .For each image we determine ∆t and the list of proper motions of all stars and then propagate all star positions in the reference image taken at time t ref to the position where they should be at the time of the current image t with where the proper motions µ α * ,i and µ δ,i have to be converted to have units of pixel yr −1 with the pixel scale of 0.389 arcsec pixel −17 .ISIS then continues unchanged, simply using our corrected star positions r ref i (t) from Eq. ( 1) in the reference image to align the current image at time t to the reference image.
We applied this method to our datasets.To quantify the level of alignment of the images, in Fig. 2 we show the dispersions σ of the star position deviations with respect to the reference image ∆r i after ISIS aligned the images: where ∆r is the average star deviation vector that would ideally vanish.The dispersion σ for each image is calculated from σ x and σ y , the dispersions in x and y pixel-direction as they are determined by ISIS.In Fig. 2, we also compare the results for the dispersions from applying and not applying the correction.
The median dispersions are shown as well demonstrating that the level of alignment improved significantly using Gaia.

Combining images
Typically, three or four images of the same object where taken per night in the same filter with an exposure time of 180 s each.
These are combined into one image per night (but we do nor combine data from different telescopes) to improve the S/N.In the same step, an error image is calculated, since a high quality error image is needed in order for GALFIT and HOTPAnTS to work properly and reliably in the following steps.The images are combined, and the associated error image is calculated, using Eqs.( 3) and (4) below, which describe a weighted mean with the possibility of rejecting the N rej /2 minimum and N rej /2 maximum values in every pixel position of the N photometrically scaled single images (min-max-rejection) in order to get rid of bad pixels from e.g.cosmic rays (typically we use N rej = 2).The photometric scales s i measure the relative general brightness and exposure of the images.They are of the order of one and are determined by aperture photometry of multiple (∼20) reference stars (checked for non-variability over time with respect to each other) in the field via s i = median F j 1 /F j i , where F j i is the flux of reference star j in image i.The weights ω i of the single images are chosen to be the inverse of the photometric scales s i 8 .The combined image z is given by where the z i are the background b i reduced single images d i , i.e. z i = d i − b i , i counts over the sorted values of s i z i where the min-max-rejection is already applied and (x, y) are the individual pixel positions.By using Gaussian error propagation on Eq. ( 3) and Poisson statistics i.e. ∆d i (x, y) = d i (x, y), we arrive at the error image given by where the correction factor of the min-max-rejection f rej , as well as the background flux error ∆b i and the scale error ∆s i are determined statistically.For the correction factor, we use f rej = noise of background in combined image median background of error image with f rej =1 , which typically is of the order of one.

PSF photometry and extracting quasar image positions improved with Gaia data
Next, we use PSF photometry with GALFIT on the combined images, not to extract the light curves, but to determine the position of quasar image A and the PSF for every night.They will be necessary when applying PSF photometry with GALFIT to the difference images in the final step (where positions and PSFs cannot be determined reliably).This procedure leads to significantly improved light curves.The PSF is a 30 × 30 pixel cutout of a star in the neighbourhood of the quasar we select for GAL-FIT.We modified GALFIT to include our multiple quasar model, which consists of copies of the PSF arranged in the orientation and fixed relative distances (to quasar image A).The relative coordinates of the images (see  (Falco et al. 2001).GALFIT then fits our 8 The reason for choosing ω i = s −1 i is that inverse variance weighting results in the smallest variance of all possible weights of a weighted average.Then, using Poisson statistics, i.e. that the variance in each pixel is given by the flux itself, we find this expression, where we approximate the flux for all pixels by the overall photometric scale of the image, as it is a measure of the typical brightness of the individual images. 9https://lweb.cfa.harvard.edu/castles/A119, page 3 of 9  2)) i.e. the standard deviation of all star position distances in mas between an aligned image and the reference image for all aligned images without (blue) and with (red) our new method using Gaia proper motion data.The improvement of the median dispersion over all images is shown as well.
quasar model to the combined images by minimizing the sum of squared differences of model flux and data (χ 2 red ) returning the fluxes of the images and the quasar position of image A (see also Giannini et al. 2017).
As an example, in Fig. 3 we show PSF photometry results of HE1104-1805 with GALFIT.The image separation is ∼3.2 arcsec (see Sect. 2.1) and the lensing galaxy is not visible in our data (but faintly in e.g.HST images).This works analogously with the quasar HE2149-2745: Again, two images are visible, they have an angular separation of only ∼1.7 arcsec (see Sect. 2.2) and again the lens galaxy is orders of magnitudes fainter and is not accounted for in our GALFIT model to determine the quasar position and PSF.The third quasar in our dataset, Q2237+0305, i.e. the Einstein Cross, is visible as four images of the quasar on a ring with a radius of ∼0.9 arcsec (see Sect. 2.3).For this object, the quasar images are barely visible in our data because of the bright lensing galaxy core (see Fig. 1), but well accessible to difference photometry, as we show below.
To adapt the GALFIT model to the Einstein Cross we have to fit four instead of two quasar images and also include the galaxy  model described in Giannini (2017), which consists of (i) a de Vaucouleurs profile (i.e. a Sersic profile with fixed sersic index of n = 4) for the galactic bulge and (ii) an exponential profile (n = 1) for the galactic disk, which are both given by the Sersic profile where I e = I(r e ) and b n ≈ 2n − 1/3 are chosen such that r e is the half-light radius.We used the values from Giannini (2017) to fix the scale radii, the galaxy orientation and the axis ratio, but the relative positions and brightnesses are not fixed.For the average magnitude difference of the galaxy components we find ∼0.9 mag as in Giannini (2017).In Fig. 4, we show the positions of image A of the Einstein Cross as a function of time from 2014 to 2023 as determined with GALFIT in our R data as data points (for the x and y directions respectively).These quasar image A positions r QSO (t) follow lines with non-vanishing slopes.But since we aligned our images (and showed that our method leads to well aligned images), the quasar position should not change with time in our images since the quasar is situated at cosmic distances.The linear trend in the positions comes from the proper motion of the PSF star, which moves inside its fixed 30 by 30 pixel box over the years and therefore the position of this box (fitted to the non-moving quasar) changes with time.The plotted lines are no fits to the data points, but lines with their slopes being calculated from the Gaia DR3 (Gaia Collaboration 2023) proper motion value for the PSF star µ PSF .The absolute position offset r0 is fixed with the median of The lines are no fits to these positions, but simply lines rQSO (t) with the proper motion µ PSF of the PSF star from Gaia as slope and the median of all proper motion corrected positions as overall offset r0 (see Eq. ( 6)).We did this for all six datasets and then used the lines as fixed quasar image A positions when applying GALFIT to the difference images to determine the light curves.
all proper motion corrected positions and we use these lines, calculated from Gaia proper motion data of the PSF star, as quasar image A positions: where τ(t) = t − 2014.0 yr.We use these improved quasar image A positions rQSO (t), instead of the data points r QSO (t), since they match the data perfectly.Additionally, this takes the median of the positions, giving a better value, since the (PSF proper motion corrected) quasar position must be constant.This was done for all six datasets (three quasars in the R band and the V band) and the improved quasar positions are kept fixed in the final determination of the light curves using difference images.

Difference imaging
We also construct a best seeing reference image from ∼10 low seeing images (∼1.3 arcsec) in the same way, as we constructed daily images above (see Appendix A for more information on the reference images).The difference images are the combined images minus the reference image.However, since the reference image has better seeing than the combined images, HOTPAnTS convolves the reference image with a convolution kernel estimated from the PSFs of 5 × 5 so-called stamps, which are 31 × 31 pixel boxes centered on stars covering the whole image10 , as described in Alard & Lupton (1998) and Alard (2000).This method requires the specification of the σ values for the three Gaussian components of the kernel function.We use 0.8 pixel, 2.4 pixel, 4.0 pixel, respectively, adapted to the range of seeing values obtained with the LCO telescopes.Moreover, each Gaussian component is multiplied with polynomials of order 4, 3, and 2, respectively.The variation of the kernel function over the 5 × 5 regions of our images is modelled with a spatial polynomial of order 2. Finally, the background is also allowed to vary with a spatial polynomial of order 2. The subtraction removes all constant sources in the image and especially removes all the light from the lens galaxy (see also Woźniak et al. 2000a,b).In Fig. 5, we show the difference images and the reference image of HE1104-1805 in the R band.
The remaining variation is due to quasar image brightness variations with respect to the reference image.We measure these with GALFIT (using only the quasar model since the galaxy light has been removed) with the improved quasar positions from the previous step (which are kept fixed), as well as the PSFs.The results are light curves of our quasar images.
One more improvement of the final light curves of the Einstein Cross was achieved by splitting the data into two intervals overlapping in 2018.We only work with reference images of the respective interval and carry out the difference imaging separately.This leads to two light curves for all images that overlap in 2018, which is where we link the brightnesses from times ≥2018 to those <2018.The calibration of the quasar fluxes is described in Sect. 4. It can be done consistently with either interval.For the Einstein Cross we used the first interval.The light curves of the Einstein Cross are improved since the stars do not move too far away from their positions in the reference image during these shorter intervals.The difference images from HOTPAnTS then contain less residuals.

Light curves of the lensed quasars HE1104-1805, HE2149-2745 and Q2237+0305
Applying the methods described in Sect.3, we reduced our six datasets of LCO data.The number of images in or after the different steps and the final number of epochs used for the light curves are given in Table 2.The main loss of images of about 14 % from all LCO images to the aligned images comes from the visual inspection of the images and subsequent excluding defocused or otherwise damaged images.Typically there are three to four aligned images per night that are then combined into one image as described and for each combined image there is an corresponding difference image (plus one or two reference images for each of our six datasets).In the final light curves, a small number of data points with very large uncertainties were removed (those with σ > 0.05 mag in HE1104-1805 and σ > 0.04 mag in HE2149-2745).The outliers were identified in a logarithmic histogram of the measured uncertainties.
Figure 6 shows the light curves i.e. the apparent magnitudes with time of the multiple quasar images extracted from the difference images with GALFIT as described in Sect.3.4.The quasar brightness in the reference images (for each quasar and each band) are determined by running GALFIT (using the multiple quasar model and additionally including the galaxy model only for the Einstein cross).The zero points are determined from multiple (∼20) reference stars with known brightnesses.The apparent magnitudes of the reference stars used to determine the zero points are calculated from the Gaia magnitudes G, G bp and G rp using the conversion formulas to the R and V band from Table C.2 in Riello et al. (2021).
Finally, the light curves of the images of the double quasars HE1104-1805 and HE2149-2745 are time delay corrected (see  Sects.2.1 and 2.2), as well as shifted with magnitude offsets for displaying reasons to separate the two bands.They show the intrinsic variations of the quasar as well as deviations between the images indicating microlensing events.We can compare the LCO light curve for HE2149-2745 with the R band light curve from Millon et al. (2020) because the observing intervals overlap in the years from 2014 to 2016.They are in good agreement, although the fainter image B displays larger scatter in the LCO data.
The light curves of the four images of Q2237+0305 have time delays ≲ 1 day (see Sect. 2.3) and are presented here with ∆t = 0 days and no relative magnitude offsets.The brightness variations in the images of the Einstein Cross appear mostly different and uncorrelated (except e.g. for a magnitude decline in 2019) indicating microlensing affects all images all the time (see e.g.Udalski et al. 2006;Millon et al. 2020, and references therein).
We have used Gaia proper motion data to improve the image alignment (see Sect. 3.1 and Fig. 2) and the quasar position (see Sect. 3.3 and Fig. 4).The effect of these improvements on the quasar light curves is illustrated in Fig. 7.The left panel and right panel compare the light curves of Q2237+0305 in the V band without and with these two Gaia steps.We note that we concentrate on the interval from 2018 to 2022 with the relevant second reference image.The effect is stronger for the fainter images (especially image C which is the faintest in that period of time), but also for images A and B the signal was improved.

Conclusion
Using data taken at LCO from 2014 to 2022, we have demonstrated our method to carry out DIA on long-term monitoring data of the strongly lensed quasars HE1104-1805, HE2149-2745 and Q2237+0305 in the R and V bands.The method makes use of Gaia DR3 proper motion data of field stars to improve (a) the image alignment as described in Sect.3.1 and (b) the quasar position (which needs to be fixed for PSF photometry using an observed PSF on the difference images) as described in Sect.3.3.
Additionally, we find that (depending on the object) it can be helpful to split the data into two or more intervals with separate reference images to conduct the difference imaging.The resulting photometry can be combined afterwards where the intervals overlap (in our case: two 5 yr intervals with 1 yr overlap for the data of the Einstein Cross).This works only if there are enough high quality images to create reference images for the different intervals.When is this step necessary?One has to consider that our method aligns the images based on the hypothetical star positions at the epoch of the reference image.HOTPAnTS, however, can only work with the observed data, and so the convolution kernel will be affected by this.In the case of Q2237+0305, we find that using two intervals leads to a consistent photometry over the whole observing campaign.
There are many current applications for difference imaging over long periods of time.A prominent example is the upcoming survey of the Vera C. Rubin Observatory.Whenever stars are used as a reference to align images from different nights/epochs the method introduced here to use Gaia proper motion data for improved difference imaging can be useful.Especially whenever PSF photometry is carried out on difference imaging data where the position of the objects needs to be fixed, one has to take care to include the proper motion of the PSF star.

Fig. 1 .
Fig. 1.LCO image of the Einstein Cross (Q2237+0305) in the R band in inverted grey scale (north is down, east is right and the field size is ∼12.4 arcmin × 4.2 arcmin).Overlayed as blue arrows are the proper motions µ i of stars in the field from Gaia DR3 data (as reference: the high eastward proper motion star on the top right has µ ≈ 55.4 mas yr −1 ) that we use to improve the image alignment with ISIS (see Sect. 3.1).The PSF star is marked.In the lower right, a zoom into the central part of the image (∼25.3 arcsec × 13.6 arcsec), where the Einstein Cross is situated, is shown (image names A-D are indicated).Here we use the reference image for difference imaging in the R band.

3
http://www2.iap.fr/users/alard/package.html Sect.3.3: extract the quasar image A positions and the images PSFs using a modified version of the GALFIT software version 2.0.3 from Peng et al. (2002), and Sect.3.4: apply difference image analysis (i.e.essentially sub- Fig. 2. Results of aligning all the images of HE1104-1805 (first row), HE2149-2745 (second row) and Q2237+0305 (third row) in the R (left column) and V (right column) band.We show the dispersion σ (see Eq. (2)) i.e. the standard deviation of all star position distances in mas between an aligned image and the reference image for all aligned images without (blue) and with (red) our new method using Gaia proper motion data.The improvement of the median dispersion over all images is shown as well.

Fig. 3 .
Fig. 3. GALFIT example of HE1104-1805 in the R band from July 3, 2019: (a) the combined LCO data of that night, (b) the quasar model made by GALFIT with two copies of the PSF fitted to the two quasar images in the orientation and distance fixed with HST data and (c) the residuals between data and model.

Fig. 4 .
Fig. 4. Q2237+0305 image A x and y pixel positions r QSO (t) as a function of time (in the R band) determined by GALFIT on the combined images shown as upper (blue) and lower (orange) data points respectively.The lines are no fits to these positions, but simply lines rQSO (t) with the proper motion µ PSF of the PSF star from Gaia as slope and the median of all proper motion corrected positions as overall offset r0 (see Eq. (6)).We did this for all six datasets and then used the lines as fixed quasar image A positions when applying GALFIT to the difference images to determine the light curves.

Fig. 5 .
Fig. 5. Difference images of HE1104-1805 in the R band from mid 2014 in the upper left to mid 2022 in the lower right (second last image).Shown is the flux in each pixel relative to reference image linearly in grey scale (where white means brighter and black darker), while the last image in the lower right shows the flux in the corresponding reference image logarithmically coloured.

Fig. 6 .
Fig. 6.Final light curves.Upper row: time delay corrected light curves of the images of HE1104-1805 and HE2149-2745 in the R and V bands together respectively, including offsets of the apparent magnitudes for the individual images for displaying reasons.Lower row: light curves of the four images of Q2237+0305 separated into R and V bands with no time delays or magnitude offsets included.

Fig. 7 .
Fig. 7. Light curves of Q2237+0305 in the V band without (left panel) and with (right panel) the combined effects of the Gaia proper motion data.We added magnitude offsets to the different images (but the same for both plots) to separate the light curves.

Table 2 .
Number of images after each reduction step.