# 𝐵⁢𝑉⁢𝑅⁢𝐼 photometry of QSO 0957+561A, B: Observations, new reduction method and time delay

## 1 Introduction

Since the discovery of the first gravitational lens, the Twin QSO 0957+561 (Walsh et al. 1979), the system has been subject to the most rigorous attempt to measure the time delay between its components. The specially propitious configuration of QSO 0957+561, two images separated by $`6\stackrel{}{\mathrm{.}}1`$ of a quasar ($`z=1.41`$) lensed by a galaxy ($`z=0.36`$) placed at the center of a cluster of galaxies, make it suitable for photometric monitorings. Although the time delay controversy has recently been solved, establishing a value for the time delay of $`420`$ days (Oscoz et al. 1997; Kundić et al. 1997; Pelt et al. 1998a), an on-going monitoring of the QSO components and comparison between the light-curves may yield important results, both for the study of the physical properties of quasars (Peterson 1993, Gould & Miralda-Escudé 1997) and for the detection of possible microlensing events (Gott 1981, Pelt et al. 1998b). Moreover, and although it will not lead to substantial changes in the value of the Hubble constant, a secure statement of the time delay is crucial for microlensing studies.

The main requirement to obtain useful information from the light-curves of the two components is a high level of photometric accuracy. However,the QSO 0957+561 is a very complicated system for two main reasons: i) the proximity of the point-like QSO components; and ii) the extended light distribution of the underlying lensing galaxy. Moreover, the whole scenario presents an additional complication due to the large amount of available data to reduce and analyze (Vanderriest et al. 1989; Press et al. 1992; Schild & Thomson 1995), hence automatic photometry codes become mandatory. Up to now, the only automated solution presented was developed by Colley & Schild (1999), who used HST data to subtract the lens galaxy and estimate the level of cross talk between the QSO components by selecting reference stars.

In this paper an alternative solution to this problem is presented: PSF fitting by means of DAOPHOT software. To check the feasibility of this new technique, it was applied to a sample of simulated data.

The dataset presented here is the result of three years monitoring, from 1996 February to 1998 July, a program which included 220 sessions of observation in the $`R`$ band, 62 in the $`B`$ band, 72 in the $`V`$ band and 68 in the $`I`$ band. The data acquisition and reduction processes are explained in detail in §2 and §3, respectively. The software environment used for the different reduction and analysis processes was IRAF<sup>1</sup><sup>1</sup>1IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. (Image Reduction and Analysis Facility, see http://www.noao.edu for more information), so any task or package referred to elsewhere is included in the IRAF environment. Section 4 is devoted to presenting the observed light curves and in §5 we discuss on the time delay obtained from these data. Finally, a brief summary of the results is given in §6.

## 2 Data Acquisition

Lens monitoring was performed in three consecutive seasons, 1996 February to June, 1996 October to 1997 July, and 1997 October to 1998 July (96, 97 and 98 season respectively hereby), using the CCD camera of the 82 cm IAC-80 telescope, sited at the Instituto de Astrofísica de Canarias’ Teide Observatory (Tenerife, Canary Islands, Spain). A Thomson 1024$`\times `$1024 chip was used, offering a field of nearly 7$`\stackrel{}{\mathrm{.}}`$5. Standard $`BVRI`$ broadband filters were used for the observations, corresponding fairly closely to the Landolt system (Landolt 1992). The IAC Time Allocation Committee awarded time for two kinds of observing runs: routine observation nights (RON nights hereafter) in which we could make use of 1200 s per night, and normal observation runs (NON nights hereafter) in which the telescope was available during the whole night for our project. The observational procedure was as follows:

* RON nights: on dark nights one image of 1200 s was taken, otherwise (moon nights) several short exposures, each of 300-400 s, were performed and then re-centered on selected field stars, and averaged to give the total exposure. The position of each individual field star was measured using the imexamine task and images were combined using the imcombine task.
* NON nights: under photometric conditions, $`BVRI`$ photometry of QSO0957+561 was performed. Landolt standard fields (Landolt 1992) were observed to provide the photometric calibration. When the nights were not of photometric quality several exposures of 1200 s each were performed in every filter to obtain a final deep exposure by averaging them.

The final data set is composed by 15 $`B`$, 14 $`V`$, 44 $`R`$, 19 $`I`$ brightness measurements in the 96 season; 13 $`B`$, 25 $`V`$, 72 $`R`$, 18 $`I`$ data points in the 97 season; and 34 $`B`$, 33 $`V`$, 104 $`R`$, 31 $`I`$ data points in the 98 season. High quality photometry in $`BVRI`$ was obtained on 30 nights during 1997 and 1998. Mean results for two reference stars (H and D, see Fig. 1) and QSO components are given in Table 1. It is important to mention that QSOB magnitudes were corrected from light of the lens galaxy following the procedure explained in §3.1.

## 3 Data reduction process

A remarkable characteristic of the photometric data presented here is their high degree of homogeneity; they were obtained using the same telescope and instrumentation over the entire monitoring campaign. Therefore, the reduction process can be the same for all the frames. In a first step, the data were reduced using the ccdred package. The overscan was subtracted from the images, which were then flat-fielded using very high signal-to-noise master flats, each of them taken from the mean of ten sky flat exposures made shortly before the beginning of the observations. These basic CCD reductions (bias, flat-field) are crucial when the noise must be kept as low as possible. However, to attempt the observation of quasar brightness fluctuations of $``$ 0.01 magnitudes -in order to detect short-timescale microlensing events- a high level of photometric accuracy is needed. To this end it is crucial to separate every source of error adopting specific solutions for each of them. There are two main sources of error in CCD photometry of QSO 0957+561 system:

1. Extinction errors: It is known that the main part of the variability of the observed target magnitude is explained in terms of atmospheric extinction and air-mass variability. Extinction errors are complicated by color terms when broad multi-band photometry is dealt with.
2. Aperture Photometry errors: Due to the special configuration of QSO 0957+561 system, there are some specific aperture photometry errors to take into account. As demonstrated in Colley & Schild (1999), these errors are driven by seeing variations, and can be separated in two parts as follows, 1) Influence of the lens galaxy: Since the core of the giant elliptical lens galaxy of $`R=18.3`$ is only separated by $`1\mathrm{}`$ from the B image, most of the galaxy’s light lies inside the image B aperture, but outside the image A aperture. This effect could introduce errors of order 1-2% in the final measured fluxes from images A and B (see Colley & Schild 1999). 2) Overlapping of images. The separation between the two images is $`6\stackrel{}{\mathrm{.}}1`$ and hence, when poor seeing conditions prevail, there is an important effect of cross-contamination of light between the two quasar images.

As explained above, the amount of archived data is so large (more that 1 thousand 1kx1k CCD images) that an automated photometry code is necessary. For extinction errors, the best and traditional method to work with is to measure differential photometry with several field stars close to the lens components (Kjeldsen & Frandsen 1992).

However, the solution for aperture photometry errors presents a higher level of difficulty. The only automated solution offered to date is explained in Colley & Schild (1999). These authors used HST data (Bernstein et al. 1997) of the lens galaxy for subtraction and reference stars to estimate the level of cross-talk between the images. After these corrections, they found that photometry is reliable to about 5.5 mmag (0.55 %) over three consecutive nights of real data. In this paper an alternative solution to the problem is proposed: PSF fitting using DAOPHOT software. A new, completely automatic IRAF task, pho2com, has been developed. Using a sample of simulated data it is demonstrated that the proposed scheme can reach high precision photometry; 0.5 % for B component and 0.2 % for A component. The following two sections are devoted to explaining each of the adopted solutions to eliminate CCD photometry errors.

### 3.1 PSF photometry: the pho2com IRAF task

It is well known that PSF fitting is the most precise method to carry out photometry of faint and/or crowded field stars, whereas aperture photometry is better for brighter, isolated stars. In order to benefit from these facts the pho2com task, written in the IRAF command language, combines aperture photometry (APPHOT package) and PSF fitting (DAOPHOT package) as explained below. Before applying the pho2com task it was necessary to select an image as a reference image and re-center all the frames, using accurate centroid determination from field stars, to the reference one. The pho2com photometry has two main iterations:

* Iteration 1: accurate sky background determination.

A precise determination of the sky background is extremely important for accurate photometry. There are mainly two different ways to find the sky background: global-sky or local-sky determination. Whereas in the local-sky method the sky value is calculated from pixels around objects, in the global-sky determination the sky is described by a simple, slowly changing function of the position in the field, e.g., a plane. This last method is the most precise, but uncrowded fields are necessary in order to prevent sky level changes from field stars. This is the case of TwQSO field where most pixels see a background sky value unperturbed by stars, so the global-sky option was used for sky determination. The main steps of current iteration are:

1. Reference stars and QSO components were removed from the frame using PSF fittings (allstar DAOPHOT task). This was done, as explained above, to prevent perturbations from these objects in the sky determination.
    2. The sky level was determined by means of a smooth surface fitting (imsurf task) to the frame. The resulting image of iteration 1 is a sky-subtracted frame.
* Iteration 2: object photometry.

As commented above, the pho2com tasks uses aperture photometry for reference stars and PSF fitting for TwQSO components. Following Stetson (1987), the PSF is defined from a small sample of isolated stars (G, H, E, D stars in our case). The PSF fit has two components: an analytic and an empirical one. For the 2D analytic function the user can select between an elliptical Gaussian, an elliptical moffat function, an elliptical Lorentzian and a Penny function consisting of an elliptical Gaussian core and Lorentzian wings. These functions were applied to each frame, selecting the one which yields the smaller scatter in the fit. For the PSF empirical component a linear variation with position in the image proved to give the best results. The main steps in iteration 2 are:

1. Applying aperture photometry with a variable aperture of radius=2xFWHM (the FWHM was measured from reference stars) the reference stars fluxes were extracted. It is important to remember that the frames resulting from iteration 1 are sky-subtracted, and therefore the sky background value was forced to zero in the aperture intensity extractions.
    2. PSF fit photometry, with a variable aperture of radius=FWHM, was applied to all the objects.
    3. Aperture corrections were computed from the previous data to compare the QSO component fluxes with reference stars (aperture correction will transform data with radius=FWHM to radius=2xFWHM) and standard stars (aperture correction will transform data with radius=FWHM to photometric standard star radius, normally 4xFWHM).

A sample of simulated astronomical data was chosen in order to test the performance of the pho2com task. Simulations were made with the artdata package. Each simulated frame included the lens galaxy, the A and B quasar components and the D and H reference stars (see Table 1 for photometric data). The lens galaxy was created with a de Vacouleurs (elliptical) light distribution, $`I(r)=exp\{7.67[(r/R_e)^{(1/4)}1]\}`$ with $`R_e=4\stackrel{}{\mathrm{.}}5`$, taking into account published HST data (Bernstein et al. 1997) and ground-based photometry (Schild & Weekes 1984, Bernstein et al. 1993). The accurate position of each object was also defined using HST astrometry. Finally 200 simulated images were created with the mkobjects task. The only free parameter (see Table 2) was the atmospheric seeing, which was simulated with values between $`0\stackrel{}{\mathrm{.}}9`$-$`2\stackrel{}{\mathrm{.}}7`$ (see Figs. 2 and 3). Effects of pixellation and noise were included (for more details see mknoise task). Noise effects were considered by adding a Gaussian and Poisson noise to the images, which have a constant background (for each filter a mean sky value is deduced from real data). This kind of ideal photometry is not, of course, a full noise description. In any case the main error sources (lens galaxy light contamination and cross-talk between components) were included in the simulated images so the final estimated errors should be considered first order ones, where high order corrections (faint neighboring stars or galaxies, basic CCD reductions.., see Gilliland et al. 1991) are neglected. Aperture (with a fixed radius of $`3\mathrm{}`$) and pho2com photometry was applied to the simulated images. Differential light curves are plotted in Figs. 2 and 3. Correlations with seeing variations are clear. Although the seeing profile is the same for the two reference stars and the QSO 0957+561 components, fixed aperture photometry has final mean errors of $`1.5\%`$, $`2.2\%`$ for A and B components respectively.

PSF fitting photometry improves aperture photometry magnitudes, but subtraction of the lens galaxy is still not perfect and some of its light is present in the final B component magnitude; therefore the final QSOB magnitudes are over-estimated. To correct B magnitudes from underlying galaxy light linear relations between seeing and magnitude errors were calculated by means of simulated data. Figures 4 and 5 are plots of magnitude errors for A and B components versus seeing for $`BVRI`$ filters. After correcting data for these errors, final errors of $`0.2\%`$, $`0.5\%`$ were obtained, for the A and B simulated components, respectively. Two main conclusions can be deduced: 1) As explained above, the B component presents higher errors than A, due to its proximity to the lens galaxy; 2) Because the lens galaxy is extremely red (Schild & Weekes 1984), QSOB magnitude errors are larger in the red colors. Real data were also corrected for underlying galaxy light using the linear correlations of Figs. 4 and 5.

### 3.2 Differential photometry

The basic technique of differential photometry is very simple, and consists in determining the difference, in terms of magnitude, of the A and B images to selected field stars.

The transformation equations used to obtain the standard magnitudes are the following:

$`b=`$ $`B+B_0+B_1(BV)+B_2X`$ (1)   
$`v=`$ $`V+V_0+V_1(BV)+V_2X`$   
$`r=`$ $`R+R_0+R_1(VR)+R_2X`$   
$`i=`$ $`I+I_0+I_1(RI)+I_2X,`$   
where $`BVRI`$ are the standard magnitudes; $`bvri`$ are the instrumental magnitudes (i.e. $`r=2.5log[F_r]`$, where $`F_r`$ is the object flux through a predefined aperture); $`X`$ is the airmass; and $`(B_0,V_0,R_0,I_0)`$, $`(B_1,V_1,R_1,I_1)`$, $`(B_2,V_2,R_2,I_2)`$ are the zero-point constants, the color term coefficients and the extinction coefficients, respectively, determined from observations of standard stars. For a given object, the main source of magnitude variability can be explained in terms of atmospheric extinction and air mass variability. The usual way to remove this error is to use a comparison star observed at the same time under the same conditions (this is one of the main advantages of CCD observations). Under this assumption, the differential magnitude, for instance $`R`$, is then found as

$$r_or_s=(R_oR_s)+R_1[(VR)_o(VR)_s],$$
(2)   
where subindices $`o`$,$`s`$ represent the object and comparison star respectively. The term $`R_1[(VR)_o(VR)_s]`$ is very important and is null only if the color term of the system is equal to zero, $`R_1=0`$, or the target object and the companion star have similar colors, $`(VR)_o=(VR)_s`$. In $`BVRI`$ photometry, color terms are not zero and, to decrease errors it is necessary to have similar spectra for the object and the comparison star. In this case it is possible to approximate $`R_o=R_s+(r_or_s)`$.

Figure 1 shows the field of QSO 0957+561 in the R band obtained as a combination of all the individual images taken during the three seasons. The total equivalent exposure time is 51.46 hours and the limiting R magnitude 25. The set of potential comparison stars, F, G, H, E, D, were examined differentially in sets of 4 versus one star. This allowed us to establish the stability of each comparison star. After careful analysis, only two stars -D and H- were selected as reference stars for differential photometry. Photometric errors were calculated using the statistical error analysis developed by Howell et al. (1988), which uses the rms of the differential photometry of comparison stars (H-D in our case) to deduce the photometric errors of QSO components A and B. In initial rms calculations the derived values are higher than expected so Eq. 2 was considered which, for selected reference stars, can be written as

$$r_Hr_D=(R_HR_D)+R_1colVR_{HD},$$
(3)   
where $`colVR_{HD}=(VR)_H(VR)_D`$ which, taking into account the data in Table 1, is equal to 0.08. The color terms $`B_1,V_1,R_1,I_1`$ are not normally expected to change during the course of a night, as they are due to the mismatch between the instrumental bandpasses and the standard Johnson $`BVRI`$ bandpasses. However, instrumental bandpasses are derived as the convolution of the mirror reflectivities, the filter transmissions, and the chip response, so significant changes are indeed expected in the course of a season. Under this assumption, Eq. 3 can be formulated as

$$r_Hr_D=(R_HR_D)+f_R(JD)colVR_{HD},$$
(4)   
where $`f_R(JD)=R_1`$ is a smooth function of Julian Day which fits the possible time changes of the color term $`R_1`$. This equation is demonstrated in Fig. 6, where we plotted the color term $`R_1`$ derived from Landolt standard stars and the same term derived from Eq. 4 using observational data from reference stars H,D. The curve is a parabolic fitting to reference star data which has, due to error propagation, large errorbars ($`0.1`$).If it is assumed that parabolic fitting represents real data without noise it is clear that the smooth variations in the differential light curves of reference stars H,D are mainly due to changes in color terms. To correct $`R`$ data of color term variations (the process is equivalent for the other filters) the following steps were taken: 1) from the differential light curve of reference stars the $`f_R(JD)`$ function was calculated by means of a parabolic fitting; 2) for reference star data the term $`f_R(JD)colVR_{HD}`$ was directly subtracted from $`r_Hr_D`$ observational data obtaining the differential magnitude values $`R_HR_D`$; and 3) for QSO data it was necessary to assume mean constant values for $`colVR_{AD}=(VR)_A(VR)_D=0.09`$ and $`colVR_{BD}=(VR)_B(VR)_D=0.16`$, and the final corrected R magnitudes are

$`R_A=R_D+(r_Ar_D)+f_R(JD)colVR_{AD}=15.163+(r_Ar_D)f_R(JD)0.09`$ (5)   
$`R_B=R_D+(r_Br_D)+f_R(JD)colVR_{BD}=15.163+(r_Br_D)f_R(JD)0.16.`$ (6)   
For the current system, the red spectra of the D reference star and those of QSO components are similar, so the derived color term correction values are rather small, $`0.5\%`$ for the $`R`$ and $`I`$ filters. On the contrary the QSO 0957+561 is bluer than the D star, and in this case color term errors become as high as $`2\%`$, $`5\%`$ for $`V`$,$`B`$ colors respectively. The final mean errors for reference stars and A,B component light curves are presented in Table 3.

## 4 $`BVRI`$ Light Curves

The results of our monitoring program are shown in Figs. 7, 8, 9, and 10 (the photometric data are available at URL http://www.iac.es/lent). In these figures we show the light-curves and error-bars for components A (black circles) and B (red squares) of Q0957+561 in $`R`$, $`B`$, $`V`$, and $`I`$ band, where the data for the B component are shifted by 425 days (the time delay estimate in this paper, see §5). Final magnitudes were calculated using the pho2com task and finally corrected for (1) the influence of the lens galaxy (see §3.1) and, (2) color term variations (see §3.2). Note the similar behavior of the curves for both components (especially in Fig. 7, corresponding to the $`R`$ band).

The robustness of the proposed photometry method can be assessed by comparing the magnitudes of the QSO 0957+561 A and B components deduced from monitoring light curves (averaged values) and Landolt standard star calibrations (see Table 1). The calculated values are presented in Table 4. The global agreement between both sets of magnitude values is clear.

The photometric data presented in Table 1 also needs discussion. In principle the colors of QSOA and QSOB, averaged over the monitoring campaign, should be essentially the same if sight-line-dependent extinction is ignored. A slight reddening is present in component A, although the significance of this excess, $`E(VR)=0.07\pm 0.08`$, is questionable. In order to verify the significance of the previous result we have plotted, in Figure 11, the $`V`$$``$$`R`$ color difference between components A and B, with B shifted by 425 days so that the emission time is the same for both components over the monitoring campaign averaged every 20 days. The “bluing” of component A is now clear and we may try to understand its origin:

1)A lens galaxy absorption effect would have produced a redder, and not a bluer, ($`V`$$``$$`R`$) color for image B.

2)The most likely explanation, proposed by Michalitsianos et al. (1997), is that the ray paths of lensed components intercept different regions of a galactic disk associated to the host galaxy of the source that is viewed pole on and situated in the quasar rest frame.

A preliminary analysis of the the data obtained in the $`R`$ filter has yielded an important result: component B is brighter than component A. The $`R`$ data have been averaged every 10 and 20 days and then the B component light curve has been shifted by 425 days. The average difference between components A and B is $`m_Bm_A=0.06`$ mag for both the 10- and 20-day average. Moreover, the averaged B/A magnification ratio is 1.06, varying between 1 and 1.12, in perfect agreement with the results described in Press & Rybicki (1998), indicating the prolongation of the long-timescale microlensing event during 1997 and 1998. At any rate, an exhaustive analysis of the long-timescale microlensing in the whole dataset is being conducted and will be presented in a future paper. This study will also include a comparison between the short-timescale microlensing during an epoch of calmness (96/97 seasons) and the rapid microlensing at a relatively active (but non-violent) epoch (the 97/98 seasons). The consequences for the population of dark-matter objects in the lensing galaxy and quasar properties will be also discussed and put into perspective.

## 5 Time Delay

Today, the historical controversy regarding the value of the time delay of Q0957+561A, B is almost solved. After twenty years of monitoring, recent data establish this value at around 420 days. There is, however, a small controversy between two values, $`\mathrm{\Delta }\tau _{\mathrm{BA}}`$ = 417 days (Kundić et al. 1997; Pelt et al. 1998a) and $`\mathrm{\Delta }\tau _{\mathrm{BA}}`$ = 424 days (Pelt et al. 1996; Oscoz et al. 1997; Pijpers 1997; Goicoechea et al. 1999). The difference (one week) is irrelevant in the Hubble constant calculations, but it may be crucial in order to detect microlensing events.

One of the “classical” ways of obtaining the time delay between components A and B of Q0957+561 is the computation of the $`A`$$``$$`B`$ cross-correlation (see Oscoz et al. 1997, and references therein). In the standard procedure, the maximum of the CCF (cross-correlation function) is identified with the time delay. However, the delay-peak generally has an irregular shape, and this fact causes a bias in the measurement of the time delay between the two components of the system. In this way, two different datasets could lead to two different estimates of the time delay that are in appreciable disagreement. The problem was considered by some authors in the past. Lehár et al. (1992) made a parabolic fit around the maximum of the cross-correlation function, whereas Haarsma et al. (1997) used a cubic polynomial fit to the delay-peak. Lehár et al. (1992) suggested that the delay-peak of the cross-correlation function should be closely traced by the central peak (around $`\tau =0`$) of the autocorrelation function. Moreover, other features of the cross-correlation function around lags $`\tau _1`$, $`\tau _2`$,… will be closely reproduced in the autocorrelation function around lags $`\tau _1\mathrm{\Delta }\tau _{\mathrm{BA}}`$, $`\tau _2\mathrm{\Delta }\tau _{\mathrm{BA}}`$,…, respectively.

In this paper we make use of the similarity between the discrete autocorrelation function (DAC) of the light curve of one of the components (B, for example) and the $`AB`$ discrete cross-correlation function (DCC) to improve the estimation of the time delay. The same origin of the A and B curves guarantees the fulfilment of the relationship DCC($`\tau `$) $``$ DAC($`\tau \mathrm{\Delta }\tau _{\mathrm{BA}}`$) in the absence of strong microlensing masking the QSO’s intrinsic variability. However, several questions such as the impossibility of observing the system during certain months of the year and the necessary lack of suitable edges, are additional drawbacks. So, the comparison between the DAC and the DCC from real data should be done by previously selecting a “clean” dataset, i.e., a homogeneous monitoring of both images during two active and clear (free from large gaps and microlensing) epochs separated by $``$420 days (the rough estimate of the time delay). Therefore, from the DAC and DCC functions, one can define the following function for every fixed value $`\theta `$ (days):

$$\delta ^2(\theta )=\left(\frac{1}{N}\right)\underset{i=1}{\overset{N}{}}S_i\left[\mathrm{DCC}(\tau _i)\mathrm{DAC}(\tau _i\theta )\right]^2,$$
(7)   
where $`S_i=1`$ when both the DCC and DAC are defined at $`\tau _i`$ and $`\tau _i\theta `$, respectively, and 0 otherwise. Equation (7) can be minimized to obtain $`\theta _0=\mathrm{\Delta }\tau _{\mathrm{BA}}`$, the most probable value for the time delay. This least squares comparison ($`\delta ^2`$-test) of the auto and cross-correlation functions enables the time delay to be determined by comparing two discrete series, DCC and DAC, which should, in general, have the same shape.

### 5.1 Simulated data

Prior to computing the time delay from real data, the $`\delta ^2`$-test was applied to some simulated datasets to verify its reliability when dealing with discrete and irregularly sampled datasets. Several sets of artificial photometric data with similar magnitudes, error-bars and time distribution to that of the observations collected at Teide Observatory were created. In this section we will use the same terminology as with real data; that is, the $`y`$-axis will be considered as magnitude, the $`x`$-axis as truncated JD, and the delay between both curves as time delay.

A program was developed to generate sets of dates, $`x_i`$, between 1800 and 2000 (JD), approximately, with a pseudo-random separation taken from a uniform distribution between zero and five days, obtained with the G05CAF NAG function. The time data were then alternately separated in two different subsets, corresponding to A and B light curves. A first value of the magnitude was then calculated with the equation $`y_i=F(x_i)`$, where $`F(x_i)`$ is a selected function of the dates $`x_i`$. The probability of measuring a value $`y`$ for each $`x_i`$, due to several “observational effects”, is proportional to $`\mathrm{e}^{[(yy_i)^2/2\sigma _i^2]}`$, and hence characterized by $`\sigma _i`$, or, equivalently, the variable $`d=yy_i`$ is distributed as $`\mathrm{e}^{[d^2/2\sigma _i^2]}`$. A $`\sigma _i`$ taking pseudo-random values between 0.01 and 0.03, obtained with the G05CAF NAG function, was generated for each $`x_i`$. From here the quantities $`d_i`$, pseudo-random numbers obtained from a normal Gaussian distribution with zero mean and standard deviation $`\sigma _i`$, were calculated with the G05DDF NAG function, allowing them to adopt positive or negative values. Finally, the magnitude was generated from the equation $`y_o=F(x_i)+d_i=y_i+d_i`$, with an error-bar of $`\sigma _i`$. The A component was forced to be brighter by adding 0.1 to the magnitudes of the B component (although this situation is not realistic, it may be illustrative); moreover, 420 days were subtracted from the JD of the A dataset to simulate the existence of a time delay. The result was two pseudo-random sampled functions with pseudo-random noise, a true delay of 420 days, and the B component 0.1 mag fainter than component A. The first two selected functions were:

$$\mathrm{F1}:y=17.17+0.5\mathrm{e}^{0.5f}\mathrm{sin}(f),\mathrm{where}f=\frac{\left(x1800\right)}{20}$$
(8)   
$$\mathrm{F2}:y=17.2+0.1\mathrm{sin}(f)\mathrm{sin}(4f),\mathrm{where}f=\frac{x}{40}$$
(9)   
An additional function, consistent with the actual variability of Q0957+561, was created. The raw observational data, with none of the modifications explained in this paper, from 97–98 seasons were selected. The light curves were then fitted by the function

$$\mathrm{F3}:y=17.070.16\mathrm{e}^f,\mathrm{where}f=\frac{(x15.8m)^2}{2(10+s)^2}$$
(10)   
$`m`$ being the mean of the JD in the selected range and $`s`$ its standard deviation. The resulting simulated data show a lower variability to that obtained from F1 and F2.

To calculate the DAC and the DCC functions the procedure described in Edelson & Krolik (1988, see also Oscoz et al. 1997) was followed. For two discrete data trains, $`a_i`$ and $`b_j`$, the formula corresponding to the DCC is

$$DCC(\tau )=\frac{1}{M}\frac{\left(a_i\overline{a}\right)\left(b_j\overline{b}\right)}{\sqrt{\left(\sigma _a^2e_a^2\right)\left(\sigma _b^2e_b^2\right)}},$$
(11)   
averaging over the $`M`$ pairs for which $`\tau \alpha \mathrm{\Delta }t_{ij}<\tau +\alpha `$, $`\alpha `$ and $`e_k`$ being the bin semi-size and the measurement error associated with the data set $`k`$, respectively. The expression for the DAC can be obtained in a straightforward manner from Eq. (11), while the expression for $`\delta ^2`$ is given by Eq. (7). Finally, to calculate the uncertainty in the estimation of the time delay a Monte Carlo algorithm with 1000 iterations was applied to the simulations (see Efron & Tibshirany 1986).

The three simulated clean datasets are shown in Fig. 12. Open circles correspond to the A component, while red filled squares correspond to the B component shifted by 420 days and with an offset in magnitude. As can be seen, the two first sets of simulated data (Fig. 12, a-b) could represent violent epochs in the source quasar, with episodes where the variability is as much as 0.2–0.3 mag in only 20–30 days. The last set (Fig. 12, c) represents an epoch with less variability than the observational one reported in Fig. 7. The $`\delta ^2`$-test was applied to each clean dataset in three different cases: (i) DAC obtained from the A component; (ii) DAC obtained from the B component; and (iii) similar to (i) but this time with a large gap in the light curve B (32 days for F1, 30 days for F2, and 30 days for F3). The resulting values for the time delay and the corresponding error (1$`\sigma `$) in days, see Table 5, clearly indicate that the $`\delta ^2`$-test offers good estimates in all the simulations, even considering the large error-bar generated for each point, the existence of “periodic” trends, and the presence of some gaps in some light curves. From Table 5 one can see that the maximum difference between the real and the central value of the derived time delay is of seven days (for a relatively inactive source) and the 1$`\sigma `$ intervals always include the true delay. An example of the performance of the $`\delta ^2`$-test has been plotted in Fig. 13. The DAC (open circles) for the A component shifted by 420 days versus the DCC (red squares) for F2 appear in the upper panel. There is a very good correspondence between both curves. Possible values of the time delay ($`\theta `$) versus the associated values $`\delta ^2(\theta )`$, normalized by its minimum value, have also been represented in the lower panel.

### 5.2 Real data

The success of the calculation of the time delay from simulated data, as shown in §5.1, made it reasonable to apply the $`\delta ^2`$-test to real data. The observations, collected at Teide Observatory, covered three consecutive seasons (1996, 1997, and 1998), with 220 different points in the $`R`$ band. Some points are affected by strong systematic effects and show a strong and simultaneous variation in both components. Once these points were discarded, their total number was reduced to 197. Taking into account the presence of two main gaps in the data—JD 2450242 to 2450347 and JD 2450637 to 2450729—roughly corresponding to the summer months, two different datasets (free from large gaps and edges) can be selected: DSI, corresponding to the 96–97 seasons, with 28 points for the A component and 27 points for the B component; and DSII, corresponding to the 97–98 seasons, with 44 points for the A component and 86 points for the B component. Both DSI and DSII have been represented in Fig. 14, where the B light curves have been shifted by 420 days and +0.06 mag. As can be seen, DSI corresponds to an epoch of significant calmness in the activity of the quasar which, together with the relatively small number of points, made it problematical for time-delay calculations. This fact was stated after some tests. On the contrary, DSII (the 97 and 98 seasons) shows some level of activity (although not as strong as in A95/B96) and, moreover, contains an appreciable number of points. Neither is there any clear evidence for any microlensing event, a fundamental requirement for selecting a clean dataset. So, DSII was finally used to perform time-delay calculations, i.e., DSII is our clean dataset.

The DAC and DCC functions were obtained with the same procedure as in §5.1, taking into account that the better monitoring of the B light curve as compared to that for the A light curve (see Fig. 14) made it more suitable for the DAC calculations. The application of the $`\delta ^2`$-test to the DAC and DCC curves appear in Fig. 15 (normalized as in §5.1), where the minimum of the $`\delta ^2`$-curve appears at 425 days, corresponding to the best delay. The uncertainty in our estimate of the time delay was obtained by using a Monte Carlo algorithm. A random-number generator added a variable to each point of DSII to simulate the effects of the observational errors (see Efron & Tibshirany 1986), standard bootstrap samples being thereby obtained. The $`\delta ^2`$-test was applied to the bootstrap samples to get the time delay in each case, repeating the process 10000 times, a number large enough for the results to be treated statistically. The use of the Monte Carlo algorithm led to a final value of $`425\pm 4`$ days ($`1\sigma `$). The uncertainty with the $`\delta ^2`$-test is better than the uncertainties obtained with the same clean dataset with other alternative methods, like the dispersion spectra and the discrete cross-correlation techniques, $`426\pm 12`$ days and $`428\pm 9`$ days, respectively (see Oscoz et al. 1997 and references therein). On the other hand, the $`\delta ^2`$-test with the DAC obtained from the A light curve gives a time delay of $`425\pm 5`$ days. The resulting DCC (filled squares) and DAC (for the A component, open circles) curves are presented in Fig. 16, where a bin semi-size of $`\alpha =`$ 20 days was used. The DAC has been shifted by 417 (upper panel) and 425 (lower panel) days. The disagreement between both curves is evident in the former case. Our study indicates that the time delay between components A and B of Q0957+561 must be in the interval 420–430 days and is therefore slightly different from the “standard” typical value of 417 days.

## 6 Conclusions

CCD observations of the gravitational lens system Q0957+561A,B in the $`BVRI`$ bands are presented in this paper. The observations, taken with the 82 cm IAC-80 telescope, at Teide Observatory, Spain, were made from the beginning of 1996 February to 1998 July, as part of an on-going lens monitoring program. An alternative method to obtain accurate multi-band CCD photometry of this object is presented. A new, completely automatic IRAF task, pho2com, has been developed. This code yields accurate photometry by simultaneously fitting a stellar two-dimensional profile to each QSO component by means of DAOPHOT software. Using a sample of simulated data, it is demonstrated that the proposed method can achieve high precision photometry, 0.5 % for B component and 0.2 % for A component. In this paper we show that it is also necessary to correct $`BVRI`$ photometry for color term variations during a season, and a possible procedure is presented. Although PSF fitting photometry improves aperture photometry errors, the subtraction of the lens galaxy is still not perfect and some of its light is present in the final B component magnitude, therefore the final QSOB magnitudes are over-estimated. To correct B magnitudes from underlying galaxy light, linear relations between seeing and magnitude errors are deduced by means of simulated data. A remarkable characteristic of the final presented light-curves is their high degree of homogeneity; they have been obtained using the same telescope and instrumentation during the three years of monitoring campaign.

A calculation of the time delay between both components by using a clean dataset has been performed. The resulting delay, obtained with a new test, the $`\delta ^2`$-test, is of 425 $`\pm `$ 4 days, slightly higher than the value previously accepted (417 days), but concordant with the results obtained by Pelt et al. (1996); Oscoz et al. (1997); Pijpers (1997) and Goicoechea et al. (1999).

We are especially grateful to E. E. Falco for advising us on the possible presence of strange points in our datasets and to F. Atrio for their help in understanding some hidden aspects in our statistical treatment. The authors would like to thank Dr. Jesús Jiménez and Dr. Francisco Garzón for making the telescope readily available to us. This work was supported by the P6/88 project of the Instituto de Astrofísica de Canarias (IAC), Universidad de Cantabria funds, and DGESIC (Spain) grant PB97-0220-C02.             

.

.

.

.

.

.

.

.

.