1 - Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern, 0315 Oslo, Norway
2 - Centre for Advanced Signal Processing, Copenhagen, Denmark
3 - Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

Abstract
We present a re-reduction of archival CCD frames of the doubly imaged quasar 0957+561 using a new photometry code. Aperture photometry with corrections for both cross contamination between the quasar images and galaxy contamination is performed on about 2650 R-band images from a five year period (1992-1997). From the brightness data a time delay of $424.9 ~\pm 1.2$ days is derived using two different statistical techniques. The amount of gravitational microlensing in the quasar light curves is briefly investigated, and we find unambiguous evidence of both long term and short term microlensing. We also note the unusual circumstance regarding time delay estimates for this gravitational lens. Estimates by different observers from different data sets or even with the same data sets give lag estimates differing by typically 8 days, and error bars of only a day or two. This probably indicates several complexities where the result of each estimate depends upon the details of the calculation.

Key words: gravitational lensing - quasars: individual: QSO 0957+561 - techniques: photometric - methods: data analysis

The first reported example of gravitational lensing, the twin quasar QSO 0957+561, was discovered in 1979 by Walsh et al. (1979). It is one of the most studied objects in modern cosmology, and the research and monitoring campaigns have mainly been fueled by the desire to measure the time delay, and thereby, to get an independent and direct estimate of the Hubble parameter (Refsdal 1964). In addition, several groups have tried to analyze the extrinsic variability in the light curves. This variability is assumed to be caused by gravitational microlensing (ML) by stars or MACHOs in the lensing galaxy (as predicted by Chang & Refsdal 1979).

QSO 0957+561 is a doubly imaged quasar at a redshift z=1.41, with the components A and B separated by about $6\hbox{$.\!\!^{\prime\prime}$ }2$ . A massive, elliptical cD galaxy (named G1) at z=0.36, located only $\simeq$ 1'' from the center of the B image, seems to be the principal lensing object.

The closely juxtaposed quasar images and the extended brightness profile of the lens galaxy make accurate photometry a challenge. During the 1980s and mid-1990s, standard aperture photometry was performed without any corrections for the light contamination between the quasar components (crosstalk) or from the lens galaxy, see e.g. Schild & Cholfin (1986); Schild (1990); Kundic et al. (1995). Later reduction schemes have tried to address the above-mentioned problems in order to reduce the chance of correlated (seeing-dependent) brightness variations in the light curves; e.g. Colley & Schild (1999, 2000), Serra-Ricart et al. (1999). Precise photometry is necessary for time delay determinations and for investigations of possible microlens-induced fluctuations in the brightness records.

In spite of extensive observations by several groups, the time delay ( $\tau$ ) between the two quasar images has proved hard to determine. Complicating factors include heterogeneous data sets, large temporal gaps in the data sets, and additional variability in one or both of the quasar images (microlensing). Even more than 15 years after the discovery, the time delay was not determined. However, there were two favored candidates; $\sim$ 540 days and $\sim$ 415 days. The results of the different investigations prior to 1997 are summarized in Haarsma et al. (1997). Kundic et al. (1997) convincingly settled the long-standing controversy in favor of the lower value, finding $\tau = 417 \pm 3$ days. Since 1995 different groups have reported values of $\tau$ in the range 416-425 days. The results seem again to concentrate around two values; 417 days (Kundic et al. 1997; Pelt et al. 1998; Colley & Schild 2000) and 424 days (Pelt et al. 1996; Oscoz et al. 1997; Pijpers 1997; Serra-Ricart et al. 1999; Oscoz et al. 2001).

QSO 0957+561 was the first system to provide strong indications of microlensing effects; uncorrelated brightness variations between the A and B images were found by Vanderriest et al. (1989). Several researchers have reported microlens-induced variability in the quasar light curves. Pelt et al. (1998) found unambiguous evidence of long time scale (order of several years) microlensing in the "difference light curve'' (DLC; A light curve minus time-shifted B curve). Results are ambiguous when it comes to the short time scale (lasting a few months) and rapid (less than a few weeks) microlensing events. Schild & Thomson (1995a), Schild (1996) and Colley & Schild (2000) have reported interesting high-frequency features in the brightness record, having amplitudes of only 0.03-0.05 mag and time scales of months and even weeks. Goicoechea et al. (1998) also found fluctuations which could be associated with microlensing events. However, Schmidt & Wambsganss (1998), Wambsganss et al. (2000) and Gil-Merino et al. (2001) all found DLCs with no clear microlensing signature, and notably no short time scale events with $\vert\Delta m\vert > 0.05~\textrm{mag}$ were observed. Gil-Merino et al. actually showed that the fluctuations in their DLC could be due to (several) observational noise processes. To reveal any rapid fluctuations caused by microlensing, high quality images with good temporal sampling are required.

This paper is mainly a summary of some results from a Master's thesis project by Ovaldsen (2002) undertaken at the Institute of Theoretical Astrophysics, University of Oslo, Norway. We shall here concentrate on the aperture photometry scheme, the time delay estimation and microlensing investigation. The data set consists of some 2650 archival CCD images of QSO 0957+561 covering a period of nearly five years (June 1992-April 1997). This data set has previously been reduced by one of the authors (RES), but with cruder corrections for crosstalk and galaxy contamination. In the next two sections we discuss the data set and briefly present the main principles of our photometry scheme. Then, from the final A and B light curves (Sect. 4), the time delay is determined using two different statistical techniques (Sect. 5). In Sect. 6 we briefly investigate the microlensing residual. The results are summarized and discussed in Sect. 7. Our new photometry gives, among other things, a time delay that differs significantly from the result we obtain when employing the same method on the old RES brightness data.

2 Data set

Although taken with the same telescope, the quality of the frames varies considerably. Cosmic rays, and especially bad pixels and bad columns, occur frequently. Quite a few frames exhibit varying background levels, not only in the form of a gradient across the image, but as bright or dark "patches'' at certain locations. Such frames may not have been properly calibrated. (Other artifacts from poor pre-processing are also seen).

The image headers do not contain all the desired information, e.g. the pixel size and the gain factor are often missing. The gain is fixed to 2.3 e^-/ADU. The pixel size is computed empirically for each frame, using the calculated positions of typically 6 field stars and the astrometry presented in the Guide Star Catalog II (GSC-II), see Table 1 .

Two different CCDs are employed. The scale of the first one is approximately 0.65''/pixel (binned mode) and 0.32''/pixel (unbinned), and that of the second one is 0.70''/pixel. The range of seeing values (FWHM) is approximately 1''-5'', with a mean value of around 2''. The global background is mainly between 100 and 2000 ADU (92% of the frames). The stellar images are typically non-circular, the PSF having a mean ellipticity of 0.09, equivalent to an axis ratio of 1.1. We also note that the PSF often departs from elliptical symmetry. The coma-like appearance is probably due to tracking errors and astigmatism in the camera optics.

The sampling of the observations must be regarded as very good. Besides the gaps in the summer months, the one day interval dominates. More than 90% of all time separations between consecutive observation runs are less than eight days.

A very different and more homogeneous data set, comprising some 1000 R- and V-band frames obtained over four consecutive nights, is discussed in a forthcoming paper (Ovaldsen et al. 2003).

3 Aperture photometry scheme

The software used to reduce and analyze the CCD frames was developed by JT and JEO. The entire package was written almost from scratch in the Interactive Data Language (IDL); it is specially adopted to the 0957+561 twin quasar system. Several sub-routines had been implemented by JT when working on other quasar lens systems. Many of these remained unchanged. All steps are automated; from detection and localization of objects, via field star photometry and calibration, to the actual quasar photometry. We will only describe the main features of our photometry scheme. The automatic source detection program, background determination, centering algorithms etc. are not discussed here. We refer to Ovaldsen (2002) for a complete and detailed treatment of the entire package.

3.1 Field star photometry and calibration

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F1.ps} \end{figure}$

Figure 1: The field surrounding QSO 0957+561A,B. North is up, east is left. The seven field stars (F, G, H, E, D, X and R) as well as the quasar components (A and B) are indicated. The lens galaxy, not visible here, is located approximately 1'' from B, slightly to the left of the vector from B to A. The frame is 4.7 arcmin on each side. (This particular image is taken with the 2.4 m Hiltner telescope, Michigan-Dartmouth-MIT observatory, and is of better quality than the frames in our data set both in terms of seeing and spatial resolution; pixel size = 0.275''/pixel, seeing $\approx 1''$ ).

Open with DEXTER

The separation between the two quasar images of $\sim$ 6'' motivates the use of 3'' radius apertures. Using the same aperture size for the comparison stars as for the target objects makes it easier to transform the quasar intensities into standard magnitudes. We use the stars F, G, H, E, D, X and R for the stellar photometry, see Fig. 1 and Table 1. Of course, the use of such small apertures is only valid if they collect the same fraction of the total light for all point sources, equivalent to the assumption that all point sources on a frame have the same PSF. We assume that this is approximately the case. Preferably, one should observe reference stars with the same spectral distribution as the primary targets; this would reduce the error when calculating the magnitudes of the target objects. In the case of 0957+561, the two quasar images are bluer than the field stars. However, since we only have single band observations, we are not able to correct for any color effects.

The large pixel size (mostly $\sim$ 0.65''/pixel) combined with the relatively small apertures ( $\rm radius=3''$ ) give a quite irregular polygon on the pixel array. To simulate a "perfect'' circular aperture we apply a weighting scheme for the pixels which lie on the border. The value of a partial pixel is calculated as the original pixel value multiplied by the ratio of the partial pixel area to the total (square) pixel area. We also tried to quantify the implication of a non-zero brightness gradient across the aperture border. This second-order correction, however, proved insignificant.

The local background level is calculated from 20 small apertures arranged in a circle of radius 20'' around the object of interest. The apertures containing cosmic rays, bad columns, sources etc. are automatically discarded.

Table 1: Reference star astrometry from GSC-II. Both right ascension, $\alpha$ , and declination, $\delta$ , are given in degrees. The rightmost column quotes the names of the stars as they appear in GSC-II.

We use seven comparison stars (Fig. 1) to determine the calibration level needed to put the quasar magnitudes on the standard system. The instrumental intensities are compared to the reference values, and any (5 $\sigma$ ) outliers are registered. Some frames contain fewer than seven comparison stars, but we require a minimum of three to proceed. If there is more than one outlier, the frame is simply discarded. Our fundamental assumption is thus that the intensities of all the present stars (except one possible outlier) should be consistent with the reference values. The measurement errors are taken into account. With this procedure we make sure that the calibration constant is calculated from "well-behaved'' stars and, consequently, that the quasar magnitudes will be as accurate as possible.

All measurement uncertainties (i.e. the standard deviation of the aperture intensity $\hat{I}$ ) are calculated using the formula

$\begin{displaymath} \sigma (\hat{I}) = \sqrt{I/g + n_{\rm a}\left(1+\frac{n_{\rm a}}{n_{\rm b}}\right)v_{\rm b}} \ . \end{displaymath}$

(1)

3.2 Photometry of the two quasar components

As previously mentioned, we use 3'' radius apertures centered on each quasar image. These apertures will be subject to seeing- (and ellipticity-) dependent light contamination from the neighboring quasar component and the underlying G1 lens galaxy.

Table 2: Statistics of the field star aperture photometry. For each star the mean magnitude and standard deviation are tabulated, along with the mean of the individual formal errors, and reference magnitude.

3.2.1 Galaxy subtraction

The position of the center and the orientation of the semi-major axis are calculated from the positions of the quasar images and from the relative astrometry of Bernstein et al. (1997). To synthesize the galaxy we start by oversampling the pixels four times. The value of each sub-pixel is computed by interpolating the brightness profile quadratically. The ellipticity and position angle are taken into account. When determining the calibration (or zero) level for the galaxy, the small 3'' apertures do not suffice. G1 is an extended object whose profile is much broader and totally different from that of the stars. For this reason it is important to calibrate a resolved object like G1 with the "total'' light from the comparison point sources, here taken to be the flux in apertures of radius 12''.

The model is finally "smeared out'' in accordance with the seeing on each particular image. This is done by convolving the synthesized galaxy with the image PSF. Having performed the proper scaling and positioning, the convolved galaxy image is simply subtracted from the frame.

3.2.2 Crosstalk correction

Several methods to minimize cross contamination between the A and B images were explored, some of which were similar to the procedure in Colley & Schild (1999). However, because we have to calculate the PSF for each frame (used in the modeling of the synthetic lens galaxy), we decided to utilize one of the characterizing features of the PSF-fitting technique. The A and B images are cleaned from the frame in an iterative fashion, thereby allowing aperture photometry to be performed on each quasar image after the galaxy is subtracted and after the neighboring twin is cleaned from the frame. The cleaning works well for a wide range of seeing conditions, and this way of eliminating the crosstalk between the A and B images proved to be significantly more robust than the other methods (it is, for instance, less sensitive to bad columns, bad pixels etc.).

4 Photometric results

4.1 Stellar photometry

2486 frames were analyzed with respect to reference star photometry, and the calibration procedure (see Sect. 3.1) accepted 2028 frames. Figure 2 shows the magnitudes of the seven reference stars for all accepted frames, as a function of Julian Day (J.D.). Table 2 shows some statistics for each star. We remark that the stars F, G, H, E, D, X and R were saturated on 280, 289, 97, 0, 2, 719 and 0 frames, respectively.

From each frame we only calculated the magnitudes of the comparison stars which were unsaturated and had intensities consistent with the reference values. The calibration constant was computed from the same stars. As we see from the plots, there are a few "outliers''. Some of these can be identified as points with large error bars, but some are true outliers in the sense that they should have been discarded. The calibration constant is not necessarily biased by such outliers in all cases, because it is derived from several (3-7) stars.

As expected, the scatter in magnitudes increases for fainter stars. The R star is by far the faintest, and consequently has the largest standard deviation for the calculated magnitudes, i.e. 13.5 mmag. Given that its brightness is $\sim$ 0.3 mag greater than the quasar images, this should indicate the minimum general dispersion to be expected in the quasar light curves due to photometric noise. After all, the photometry of the two quasar images is much more complicated than that of an isolated star; the galaxy subtraction and the crosstalk correction obviously increase the error budget of images A and B.

4.2 Quasar photometry

4.2.1 Binning and censoring

The number of data points (or accepted frames) per night, n, varies from one to six. The observation were made within a small time interval, thus it seems appropriate to combine all data points for a particular night, and only quote "binned'' magnitude values. It is almost impossible to perform a rigorous statistical analysis with so few data points (and an unknown number of outliers due to erratic photometry of the twin images). We decided to address this issue in a simple and transparent way. From the image A magnitudes and the image B magnitudes on each night, we computed the corresponding median values. With this approach, single outliers do not bias the results too much, at least for $n \geq 3$ . Obviously, for n = 1or 2, the median-filtering will not throw away possible outliers. However, without any a priori knowledge, this is about the best we can do. We quote error bars which correspond to the median points.

Before the median filter was applied we censored the data, accepting only images with seeing $\leq$ 3'' FWHM, background level $\leq$ 3000 ADU, PSF axis ratio $\leq$ 1.3 and AB-separation $6\hbox{$.\!\!^{\prime\prime}$ }175 \pm 0\hbox{$.\!\!^{\prime\prime}$ }05$ . The final data set was then reduced to 1720 images. The subsequent "binning'' yielded 422 data points for each quasar image.

4.2.2 Crosstalk and galaxy contamination

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F2.eps} \end{figure}$	Figure 2: Aperture photometry of the field stars around QSO 0957+561. Plotted are the R-band magnitudes of the reference stars from all accepted data frames, as a function of Julian Day. Dashed red line is the mean value, dotted green lines are the $1\sigma$ limit. The stars are arranged in order of decreasing brightness. Note that the scaling for the R star is different from the others.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F3.eps} \end{figure}$	Figure 3: The percent change in the A and B aperture flux when correcting for crosstalk, as a function of seeing. The crosstalk from the B (A) image into the A (B) aperture is marked with open squares (triangles). The corrections are virtually the same for the two apertures. A best fit quadratic curve is overplotted.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F4.eps} \end{figure}$	Figure 4: Light contamination (in percent) in the A and B apertures due to the G1 lens galaxy, as a function of seeing. Best fit quadratic curves are also plotted.
Open with DEXTER

Aperture photometry on images A and B was performed both with and without the corrections for crosstalk and galaxy contamination (see Sect. 3.2), so that we could check the performance. We now make a few comments on the results from analyzing the 1720 accepted images.

Figure 3 displays how the light from the two quasar images affects each aperture. We emphasize that this effect is not only a function of seeing; the ellipticity is certainly also a factor. In particular, consider an image where the PSF is highly elliptical and has a semi-major axis parallel to the line joining the centers of A and B. The crosstalk would be larger here compared to the case where the semi-major axis is perpendicular to the AB vector. (In fact, the scatter of the corrections can be significantly reduced by imposing stricter limits on the ellipticity.)

For completeness, a least squares second order polynomial fit was computed using all the data points (see also Colley & Schild 2000). The formula for the intensity corrections, $\delta I$ , reads

We also measured the flux in the quasar apertures before and after the galaxy model had been subtracted. Figure 4 shows the light contribution from the galaxy to the A and B apertures as a function of seeing. Note that crosstalk between A and B has already been "eliminated''. As seeing deteriorates, galaxy light systematically seeps out of the B aperture, but into the A aperture. Best fit quadratic curves are overplotted to guide the eye. Some of the scatter in the plots is due to the fact that A and B itself varied during the observational period (the G1 contribution is compared to the A and B fluxes on each particular frame). The contribution to the A image aperture is between 3 and 5%, and has a moderate scatter. The B image aperture has corrections of roughly 20-18% due to the G1 galaxy. Here, the scatter is rather large, having a "full amplitude'' of $\pm$ 2%. It does not seem to increase with seeing. This probably indicates that the calibration/zero level for the galaxy is determined equally well (or poor!) for the whole range of seeing conditions. Comparable corrections for subtracting the lens galaxy contribution were discussed in detail by Colley & Schild (2000), whereas the original RES reductions incorporated subtraction of a fixed 18.34 magnitude correction for the lens galaxy contribution to the B aperture flux.

4.2.3 Quasar light curves

Figure 5 displays the light curves of QSO 0957+561A,B corresponding to the period June 1992 to April 1997. We note that the light curves show variability on both short (order of weeks) and long (order of years) time scales. For some periods there is also an apparent "zero lag'' correlation between A and B. This is best seen for J.D.-2448000 $\gtrsim$ 1700. Such a correlation should in general not exist, because the signal from image B lags A by some $\tau\sim 420$ days. We have not been able to identify the cause of this frame-by-frame correlation. It has also been reported by Colley et al. (2003), and is always presumed to be some systematic effect caused by errors in the photometry; however the amplitude exceeds any known error source.

Although our "binning'' scheme only uses the nightly median magnitude value for each quasar image, we can still estimate the dispersion for nights with two or more accepted frames. The mean standard deviations of the magnitudes on each night are 12 mmag and 11 mmag for A and B, respectively. We note that the mean of the formal error bars, as seen in Fig. 5, is 17 mmag for both quasar images. The formal error bars are rather conservative, as they include the formal errors from Poisson statistics, galaxy subtraction and calibration (see Ovaldsen 2002 for details).

We also made a rough and simple estimate of the day-to-day dispersion within the A and B brightness data: First each light curve was smoothed with a 7-point median filter (making sure not to filter over gaps greater than three days). Then the original data (A or B) was subtracted from the corresponding median-filtered curve. We allowed a maximum time gap of 1.5 days between two data points to be subtracted. The residuals should thus probe fluctuations in the A and B light curves on this time scale. We assume that this very short time scale variability is not dominated by microlensing effects. The standard deviations of the residuals are $\sigma^{\rm resid}_{\rm A} = 11$ mmag, and $\sigma^{\rm resid}_{\rm B} = 10$ mmag. These values are quite consistent with similar estimates for the image set made by Schild & Thomson (1995b), who found 9.5 and 12.0 mmag; however their reductions lack the corrections for aperture crosstalk and galaxy subtraction that are strictly functions of seeing, and are more susceptible to systematic errors on time scales relevant to seeing changes.

$\begin{figure} \par\includegraphics[width=17cm,clip]{H4015F5.eps} \end{figure}$	Figure 5: Results from aperture photometry of QSO 0957+561A,B. There are 422 binned data points for each quasar image. The B data is shifted by -0.15 mag. Error bars are $1\sigma$ limits.
Open with DEXTER

5 Time delay

5.1 Dispersion estimates

The algorithm for the Dispersion estimation technique is included in the ISDA (Irregularly Spaced Data Analysis) package, designed by J. Pelt to perform various tasks on irregularly spaced time series. It is discussed extensively in Pelt et al. (1994, 1996), so we provide here only a short review. The principle of the Dispersion method is simply to measure the dispersion, D², between the A and B image light curves for different time shift values, $\tau$ . The true value should show up as a minimum in the dispersion spectrum, $D^2(\tau)$ . By dispersion we here mean the sum of the squared differences between the A and the B image magnitudes (see below).

$\begin{eqnarray*}A_i & = & q(t_i) + \epsilon_{\rm A}(t_i)\ , \ i=1,\ldots,N_{\rm... ...) + l(t_j) + \epsilon_{\rm B}(t_j)\ , \ j=1,\ldots,N_{\rm B} \ , \end{eqnarray*}$

The combined light curve (denoted in the formulae as C) is constructed by taking the A values as they are and "correcting'' the B data by l(t) and shifting them by $\tau$ :

$\begin{displaymath}C_k(t_k) = \left\{ \begin{array}{ll} A_i & {\rm if}\ t_k=t_i \\ B_j-l(t_j) & {\rm if}\ t_k=t_j+\tau, \end{array}\right. \end{displaymath}$

(3)

The accuracy of the observations is taken into account by using the statistical weights, $W_i=1/\epsilon_{\rm A}(t_i)$ and $W_j=1/\epsilon_{\rm B}(t_j)$ . The squared difference between two data points in the CLC (see estimates below) must be multiplied with the combined statistical weights $W_k = W_{i,j} = \frac{W_i W_j}{W_i+W_j}$ .

With $l(t)={\rm constant}$ , the A and B curves are considered to be unaffected by microlensing variability, differing only by the unknown ratio of the amplification factors of gravitational (macro-) lensing. We shall, however, also compute spectra where we account for slowly varying microlensing effects in one or both of the light curves. In these cases we set l(t) equal to a polynomial (typically of order two to eight). The B data is thus "modified'' by the perturbing polynomial, into B_j + l(t_j), and the coefficients of the polynomial are determined in such a way as to minimize the dispersion between A and B data.

In this analysis we use two different methods to estimate dispersions. The simplest one is

$\begin{displaymath} D^2_3 = \min_{l(t)}\displaystyle\frac{\sum_{k=1}^{N-1} ~S_k ... ...~G_k~ (C_{k+1}-C_k)^2} {2\sum_{k=1}^{N-1}~ S_k ~W_k ~G_k} \ , \end{displaymath}$

(5)

$\begin{displaymath}S_k = \left\{ \begin{array}{ll} 1 & {\rm if}\ \vert t_{k+1}... ...m if}\ \vert t_{k+1}-t_k\vert > \delta \ , \end{array}\right. \end{displaymath}$

(6)

$\begin{displaymath} D^2_{4,h} = \min_{l(t)} \frac{\sum_{n=1}^{N-1}~\sum_{m=n+1}^... ...n=1}^{N-1}~\sum_{m=n+1}^{N}~ S_{n,m}^{(h)}~W_{n,m}~ G_{n,m}} , \end{displaymath}$

(7)

$\begin{displaymath} S_{n,m}^{(2)} = \left\{ \begin{array}{ll} 1-\frac{\vert t_... ... {\rm if}\ \vert t_n-t_m\vert > \delta \ . \end{array}\right. \end{displaymath}$

(8)

ISDA contains a simple bootstrap procedure for calculating the error bars for the minima in the dispersion spectra. The CLC is smoothed and the corresponding residuals (data points minus "smoothed curve values'') are re-shuffled 1000 times to create bootstrap samples. Smoothing is performed using a 7-point median filter, with an upper limit (typically a few days) on the time separation between two successive data points.

The time delay value from a particular dispersion estimate is taken to be the mean of the time delay distribution (an example is given in Fig. 11). The standard deviation of the distribution gives the estimated error. We quote the 1 $\sigma$ errors.

5.1.1 Complete data set

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F6.ps} \end{figure}$	Figure 6: Window functions for the complete data set when employing the estimates D²_4,2 ( $\delta =5$ days; upper curve, $\delta =2$ days; middle curve) and D²₃, $\delta =5$ days (lower curve).
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F7.ps} \end{figure}$	Figure 7: D²_4,2 dispersion spectra, $\delta =3$ . The thin curve corresponds to l(t)=l₀, while the thick curve is computed by also accounting for additional fluctuations (l(t) is a 5th order polynomial).
Open with DEXTER

From the complete data set of 422 data points for each of the two quasar images, we discarded six outliers.

The number of pairs included in the dispersion estimates depends on $\delta$ . Figure 6 displays selected window functions for D²₃ and D²_4,2. The window function is the number of nearby AB (or BA) pairs in the CLC as a function of time shift. Obviously, larger $\delta$ yields more pairs in the computation, but the overall shape of the window functions remains more or less the same (the curves get smoother as $\delta$ increases). The sampling of the observations may disfavor some time shifts, i.e. the number of pairs of nearby points in the CLC can be very low for certain shifts. Fortunately, there are no major depressions in the curves, so the statistical reliability of the dispersion values should not vary much for the different trial shifts (especially in the interesting range 400-440 days). This is a reassuring and important fact.

We shall plot the dispersion spectra for trial shifts in the range 380-480 days. However, before we present and discuss the behavior of the dispersion curves in this limited range, we show in Fig. 7 a plot of two spectra calculated using the D²_4,2 estimate where the interval goes from 0-600 days. Over the entire range, the dispersion is smaller for the estimate that includes the perturbing polynomial. A higher-order polynomial would account even more for differences in the two quasar signals, and thus decrease the general dispersion even further. One must be wary not to "over-correct'' the B data, though.

We first computed dispersion spectra using various decorrelation lengths, but without any corrections for microlensing (l(t)=l₀). Then, the calculations were repeated, but this time we included polynomials to model long time scale ML variability in the light curves. The results were not very sensitive to the degree (2nd-8th order) of the perturbing polynomial. To limit the number of variables, we fixed the degree to 5th order. Figures 8 and 9 display a selection of spectra ( $\delta = 1$ , 3, 5, 7 days) derived from the two methods.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F8.ps} \end{figure}$	Figure 8: D²_4,2 dispersion spectra with $\delta = 1$ , 3, 5, 7 days. No corrections for additional, possibly microlens-induced variability.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F9.ps} \end{figure}$	Figure 9: D²_4,2 dispersion spectra with $\delta = 1$ , 3, 5, 7 days. A 5th order polynomial was used for modeling of additional variability in the light curves.
Open with DEXTER

The curves are smoother for larger $\delta$ . Note that $\delta = 1$ hardly involves any smoothing. For a given $\delta$ , the two methods (without/with correction for additional variability) yield very similar results. In both cases, the position of the minimum increases slightly for increasing decorrelation length, from 424 to 425.5 days. Secondary minima are almost always found around 410 days, but they are moderated as the smoothing increases.

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{H4015F10.ps} \end{figure}$	Figure 10: Dependence of the $\delta$ parameter for the D²_4,2 statistic. Asterisks: no correction for additional variability. Open squares: 5th order polynomial.
Open with DEXTER

We also computed spectra for $\delta=1-20$ with increments of one, and noted the position of the corresponding minima. Figure 10 plots shift values corresponding to minimum dispersion as a function of the $\delta$ parameter. Here we see more clearly that the minima are shifted towards higher values as pairs with larger time separations are included in the estimates. However, for $\delta > 14{-}15$ , the trend is reversed, but we are probably smoothing too much already. Because the sampling of the observations is rather good (85% of the time separations are less than five days), we can afford to use small decorrelation lengths. This is reassuring, because it reduces the danger of bias from pairs with large time separations. A striking feature is also seen: whether we account for slowly varying microlensing effects or not, the minima are all at 424 days for $\delta = 1$ , 2, 3 and 4.

The same procedure was employed for the D²₃ statistic, which only includes consecutive A and B points whose time separation in the CLC is less than $\delta$ . Hence, no smoothing is performed. These spectra were similar to the curves in Figs. 8 and 9 corresponding to $\delta = 1$ . The results confirmed the trends and features which were highlighted above. Since the number of pairs included in this estimate is lower than the "smoothing'' estimates, it is not as reliable, statistically speaking. The spectra did exhibit more noise, but consistently produced global minima between 423.5 and 424.5 days when $\delta$ was $\leq$ 4 days. As before, larger decorrelation lengths yielded higher time shift values.

To summarize: the different dispersion estimates (D²₃, D²_4,2, both with and without microlensing correction) all produce minima which are concentrated around 424 days for $\delta \leq 4$ . It seems that larger decorrelation lengths introduce bias. Hence, we performed bootstrap runs only for the limited range of $\delta$ -values. As noted earlier, the results were not significantly affected by the nature of the polynomial used to model additional variability. In particular, for the preferred range of $\delta$ (i.e. 1-4 days), the minima in the dispersion spectra were all at 424 days for degrees of order two to eight. We do not want to use very high-order polynomials, as this could suppress some of the intrinsic quasar fluctuations. Hence, it seems justified to fix the degree to five.

Table 3: Time delay results for the D²_4,2 dispersion estimate applied to the complete aperture photometry data set. A 5th order polynomial was used to model additional variability in the quasar light curves. As noted in the text, other polynomials gave very similar results. Estimated errors are 1 $\sigma$ limits.

Table 3 lists the results of the bootstrap procedure. In Fig. 11 we show an example of the distribution of time delays from one of the bootstrap runs. Here, the mode and the median were both 424.5 days, while the mean time delay was 424.7 days. The shapes of the other distributions were similar, and they all had a small skew. We also note that the D²₃ estimate gave similar results with the bootstrap procedure. The estimated time delay using different setups of the Dispersion estimation technique agree well. We thus take the most probable time delay to be the average of the numbers in the table, i.e. 424.7 days, with a mean estimated error of 1.1 days. However, we do not claim that the true error is as low as this.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F11.ps} \end{figure}$	Figure 11: Bootstrap results. The mean time delay from 1000 bootstrap runs is indicated by the large tick mark. Statistic: D²_4,2, $\delta =4$ , l(t) is a 5th order polynomial.
Open with DEXTER

The magnitude difference, $\Delta m_{\rm AB}$ , between the A and B components in the 1992-1997 time span was $\approx$ 0.076 mag (we were not able to compute error bars for this parameter). Hence, B was somewhat brighter than A in the time span which our data covers. This is commonly explained by microlensing in the B component, see e.g. Pelt et al. (1998).

Finally, it is worth noting that all spectra show a local maximum at $\sim$ 418 days. So with this particular data set and the Dispersion estimation method, we can say that a time delay of roughly 418 ( $\pm$ 2) days seems rather unlikely.

$\begin{figure} \par\includegraphics[width=12cm,clip]{H4015F12.ps} \end{figure}$	Figure 12: D²_4,2 dispersion spectra. $\delta =4$ . No corrections for microlens-induced variability. The shapes of the dispersion spectra from the different data sets are very different. The position of some minima correspond to minima in the window function, i.e. the number of pairs in the dispersion estimates reaches a minimum for these time shift values.
Open with DEXTER

5.1.2 Truncated data sets

We now estimate the time delay using only selected data segments. We divide each light curve into four parts, corresponding to the seasons with J.D.-2448000 roughly in the ranges 850-1200 (period 1), 1200-1600 (period 2), 1600-1900 (period 3) and 1900-2300 (period 4), see Fig. 5. Assuming a time delay of around 425 days, it is clear that there is (fortunately) a large overlap between periods 1, 2 and 3 of A and periods 2, 3 and 4 of B, respectively. We thus have the possibility of estimating the time delay between the quasar images from three different data subsets (let us call these S1, S2, S3). The motivation behind this is to see whether the truncated data sets all produce a minimum in the dispersion spectra around 425 days. We shall not perform an exhaustive analysis, though.

Because the number of pairs included in the calculations is much lower for these subsets, the statistical reliability is not as good as in the case where we used the complete data set. We employ the D²_4,2 estimate, $\delta$ in the range 2-5 days, and compute spectra for trial shifts in the range 400-450 days. The effect of correcting for any non-intrinsic quasar fluctuations is also tested. For this we use (only) a 3rd order polynomial (the curves overlap for about 200 days in all three cases, and we do not allow for any extrinsic high-frequency components within this time span).

Figure 12 displays D²_4,2 dispersion spectra computed for the subsets S1, S2 and S3, assuming a constant magnification ratio, l(t)=l₀. Allowing for a time-dependent magnification ratio (to account for possible microlensing-effects) did not significantly change the overall shape of the curves. The minima are sometimes split in two for short decorrelation lengths, so we use $\delta =4$ which reduce the "noise''.

For the first data subset, S1, the deepest minimum is at 430 days. Secondary minima are seen around 424 and 415 days. We checked the window function too, and it contained a prominent minimum for the 430 day shift. It might be that the 430 day candidate is caused by the combination of irregular sampling and an "unfortunate'' time shift. The second subset, S2, yields two minima in the dispersion spectrum, 412 days and 425 days, the latter being the deepest. Here, the window function had no unfavorable time shift. The third plot shows the results using the S3 data, and here the deepest minimum is positioned around 425 days. Secondary minima are found for time shifts of 408 and 434 days. The window function had minima at 408 and 425 days, corresponding exactly to two of the observed local minima in the dispersion spectrum.

Bootstrap runs were performed to get an estimate of the uncertainty. From the distribution of 1000 time delays, the results (mean and standard deviation) were as follows: $430.6 \pm 2.7$ days (S1), $424.9 \pm 3.0$ days (S2), $426.1 \pm 2.3$ (S3). We do not attempt to judge the reliability of the different time delay candidates. It is interesting to see, however, that the local minima may be found in usually one of three regions, i.e. around 410, 425 and 430 days. Also, a local peak in the interval 416-420 days is found in all spectra, thus supporting the statement made in the previous section that the correct time delay is less likely to lie in this range.

For the complete data set, we found in the previous section that $\Delta m_{\rm AB} \approx 0.076$ mag. However, a brief look at the combined light curve (B shifted in time by -425 days and in magnitudes by 0.076 mag) indicated that a single magnitude difference did not optimally align the A and B data. We thus computed $\Delta m_{\rm AB}$ for each of the three data subsets in order to investigate this further. We got $\approx$ 0.089 mag, $\approx$ 0.086 mag and $\approx$ 0.050 mag for subsets S1, S2 and S3, respectively. The fact that the magnitude difference seems to be time-dependent will be discussed in Sect. 6.

5.1.3 Time delay from RES's data

Now that we have performed a time delay determination using our new photometric results, it might be interesting to see whether the "old'' reductions by Schild and collaborators give similar results. We used 537 data points for each quasar image which covered the same period as our data. An extensive analysis is beyond the scope of this paper. We shall thus only comment on the main results from the D²_4,2 dispersion estimate.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F13.ps} \end{figure}$	Figure 13: Dispersion spectra D²_4,2 ( $\delta = 1$ , 3, 5, 7 days) with 5th order polynomial, using RES's data.
Open with DEXTER

Table 4: Time delay results for the D²_4,2 dispersion estimate applied to the RES photometry data. A 5th order polynomial was used to model additional variability in the quasar light curves. Estimated errors are 1 $\sigma$ limits.

The second point above describes a general trend seen in all the different analyzes: the position of the dispersion minimum increases with increasing decorrelation length.

Here we find that the position of the global minimum was 411 days for $\delta \leq 4$ days, and it did not depend on the order of the perturbing polynomial (2nd to 8th order). There are certainly indications of additional variability in the observational data, hence it seems justified to introduce the l(t) polynomial into the dispersion minimization process. The results of the bootstrap procedure are listed in Table 4. The mean time delay is 411.7 days, with a mean estimated error of 1.9 days.

5.2 $\chi ^2$ minimization

The method of Burud et al. (2001) is based on $\chi ^2$ minimization between the data and a numerical model light curve. Microlensing in one or both light curves may be corrected for. We shall in the following carry out a brief, non-exhaustive time delay analysis with this method. Because the procedure is explained in detail in the paper by Burud et al., we only summarize the main features.

The underlying idea is that, in the absence of ML, one can model the two quasar light curves with one model curve, g(t), together with two parameters, $\tau$ and $\Delta m$ , describing the time shift and magnitude offset between images A and B. An arbitrary model curve with equally spaced sampling points is $\chi ^2$ minimized to the two original light curves. The minimization is done only for the observed data points, so only the model curve is interpolated and not the data. Because the data is irregularly sampled and contain noise, a smoothing scheme is necessary. Here, the model light curve, g(t), is smoothed on a time scale, T₁, corresponding to the typical sampling interval of the data. The smoothing term is multiplied by a Lagrange parameter, $\lambda$ , which can be chosen so that the model curve matches the data correctly in a statistical sense for adopted Gaussian statistics (this parameter has no physical meaning). In addition, each data point is given a weight which depends on the relative distances to all other points in the curve. Down-weighting is performed using a Gaussian with ${\it FWHM} = 2\sqrt{2 \ln 2}~T_2$ , where the user may choose the T₂ parameter. The weights are normalized, so that the maximum value of W_i is 1. This would be the case if only one point is within the time interval defined by T₂. According to the authors, a sensible choice of T₂ would be the approximate time scale of the intrinsic quasar fluctuations.

The method was only applied to the complete light curves. We did not attempt to model higher-order ML fluctuations in the light curves, as this is quite an elaborate process. A wide range of parameter setups was tested, and the results proved to be remarkably stable. We present only the main results.

The $\chi ^2$ -values as a function of time shift (in the range 400-450 days) are plotted in Fig. 14. The parameters were as follows: T₁=4 days, T₂=20 days and $\lambda=4000$ . We recognize some of the features from the analysis with the Dispersion estimation technique: The minimum $\chi ^2$ -value occurs for a time shift of 425 days. A secondary minimum is found around 411 days, but the $\chi ^2$ is not as low as the tiny, local minimum at 431-432 days. The overall shape of the $\chi ^2$ distribution remains the same even for large variations in the parameters. The lowest $\chi ^2$ -value is always in the range 424-426 days.

The results form Monte Carlo simulations yielded a time delay of $425.1 \pm 1.3$ days. Also with this method we find that a time delay of roughly 415-420 days seems less likely; the $\chi ^2$ -curve typically has a maximum in this range.

6 Microlensing

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F15.eps} \end{figure}$	Figure 15: Time-shifted B curve values subtracted from A curve values. Dates on the abscissa relate to the (unshifted) A light curve. The A-B differences, $\Delta m_{\rm AB}$ , are clearly time-dependent. Three examples of spline approximations to the ML residuals are also shown. Red, green and blue lines correspond to 6, 9, and 11 nodes, respectively.
Open with DEXTER

We will now briefly investigate the microlensing residual in the quasar light curves. The standard procedure is to shift the light curves in time to correct for the different light travel times, and then subtract them from each other. The last step is not trivial, as the A and B data points are irregularly sampled. This means that when we shift the B data in time by $-\tau$ , the A and B points will generally not overlap. The B data point to be subtracted from a particular A point might be several days away. We have addressed this issue in a simple way.

After having shifted the B curve by -425 days, we check for each data point of image A whether there is a point from the B curve within a certain gap limit of the current A point. If this is the case, then B is subtracted from A, and the result is stored in a "residuals array'' with an averaged time argument. The procedure only makes use of the original, raw data points, and there is one free parameter, namely the gap limit. The exact value of this parameter depends mostly on the spacing of the data, but also on the (assumed) ML time scale. With good sampling the gap limit can be set quite low (a few days) and thus, at least in principle, enable investigation of rapid ML. Lowering the gap limit will obviously decrease the number of points in the residuals array. On the other hand, the A-Bdifferences are then calculated from AB pairs with smaller time separations, which is a good thing if one wants to probe short time scale fluctuations in the residuals.

With a time delay of $\sim$ 425 days, our A and B data overlap for about 3.5 years. The ML investigation will consequently only cover this time span. The large temporal gaps in the residuals data (a result of the lack of observations in the months where 0957+561 was below the horizon) also precludes a continuous "signal'' for the whole period.

6.1 Trends on long time scales

In Fig. 15 we display the A-B residuals ( $\Delta m_{\rm AB}$ ), computed by following the above procedure. We adopted a time delay of 425 days, and the gap limit was set to 2.5 days. There are three seasons which contain an adequate number of points. The results look rather noisy, but there are some significant features. In particular, the third season clearly has a different magnification ratio between A and B, compared to the first two. Moreover, it also varies within the particular season. Variability on shorter time scales may also be seen at certain periods.

The amplitude of the variations in the first two seasons is about 0.05 mag. From the plot we can see that the average magnitude offset for the two first seasons are $\sim$ 0.09 mag, while for the third period the number is roughly 0.05 mag. This is in good agreement with the values we obtained in Sect. 5.1.2.

Here, we only want to get an idea of the general trends in the ML residuals. We thus tried fitting standard cubic splines with different number of nodes into the residuals, see Fig. 15. The (red) curve with 6 nodes "detects'' only changes on long time scales. It is thus a rather conservative guess as to how microlenses in the macro-lensing galaxy affected the light from the quasar images. The general trends are similar to the results of Pelt et al. (1998) - compare with the three last "seasons'' in their Fig. 9. Pelt et al. do not have the points around J.D.-2448000 = 2100 which we do (see Fig. 15). Although very sparse, these data indicate that the curve does not continue to fall off. The other two splines probe finer details in the ML residuals, and are more optimistic approximations. Some of the variability, notably in the gaps, is highly questionable. Still, the $\sim$ 0.05 mag drop in the difference data from the second to the third season (on the order of 300 days) seems significant. We conclude that microlensing variability of approximately 5% amplitude on time scales of less than a year has been significantly observed.

6.2 Short time scale variation

On short time scales (the order of weeks) the residuals are rather noisy. Figure 16 shows the three first seasons of the ML residuals in greater detail. For two of the seasons we also include spline approximations. As can be seen from the middle panel, the data here shows no significant trends, thus no interpolations are attempted. (Compare the scatter here to the the mean formal errors in the aperture photometry of the quasar images, i.e. 17 mmag.) It is not trivial to extract information on the true microlens-induced fluctuations from data such as these, and the interpolated curves are mostly meant to guide the eye. However, on the first plot we can discern a steep negative slope in the residuals around JD-2 448 000 = 900 followed by a positive slope some 50 days later. Optimistically, we can explain this as an "event'' lasting on the order of 70 days, but nothing certain can be said about its amplitude. It could potentially be a strong event, because of the steep gradients. The third period (lower panel) indicates more clearly a U-shaped feature. The amplitude and time scale is hard to assess, because we do not have any points in the "wings''. But is seems that the amplitude is at least 0.05 mag, and the time scale is on the order of 200 days.

7 Summary and discussion

We have presented a re-reduction of the CCD image frames critical to the discussion about time delay and microlensing in the 0957+561 gravitational lens system. Improved computational techniques allow better subtraction of the effects of the lens galaxy, and correction for the aperture crosstalk that arises in aperture photometry of the somewhat overlapping quasar images.

Analysis of the re-reduced photometry for time delay, principally using several variants of the Dispersion technique, gives consistent values around 425 days. The average result from the Dispersion method and the $\chi ^2$ minimization method is 424.9 days, with an estimated mean error of 1.2 days. However, we do not claim that the true error is as small as this. We also note that time delays of roughly 416 to 420 days were never seen in this investigation and are thus less favored by us. This is not in agreement with e.g. Kundic et al. (1997), Pelt et al. (1998) and Colley & Schild (2000).

Analysis of principally the same image frames with fundamentally different reduction and time delay estimation techniques had previously given 404 days (Schild & Thomson 1997, Direct Autocorrelation) and 416.3 days (Pelt et al. 1998, Dispersion estimation procedure), but re-analysis by Oscoz et al. (2001) of the same brightness record gave estimates near 422.6 days. Other smaller data sets for approximately the same observational epochs gave 417 days (Kundic et al. 1997, PRH method, Linear Interpolation) and 425 days (Serra-Ricart et al. 1999, $\delta^2$ method).

In all cases but the first, the quoted errors (typically a day or two) are much smaller than the discrepancies between different data sets or between estimates for the same brightness record. Critical to the discussion is the fact illustrated in Fig. 7 that the FWHM of the dispersion curve has a value of approximately 100 days, and even the local minima, as seen in Figs. 8, 9, 12 and 14, have a FWHM of 10-20 days in spite of the daily data sampling and the available several hundred data points for any test lag (Fig. 6). It is by now evident that something fundamental limits our ability to estimate time delay to the expected limits imposed by the data sampling and the observational errors.

The physical origin of this discrepancy has been attributed by Colley et al. (2003) to the fact that the quasar's luminous structure is time-resolved and microlensed. This combination of microlensing and a time-resolved source might produce multiple time delays whose pattern changes from year to year. Some evidence for this may be seen in Fig. 12, where the relative importance of persistent lags near 410, 425, and 430 days seems to have changed during the observational period. We note, however, that for some subsets these time lags coincide with minima in the window functions.

This puts a new perspective to the understanding of the role of microlensing for a quasar source. Previous discussions of microlensing (see Schmidt & Wambsganss 1998, and references contained therein) have focused upon the role of small accretion discs crossing the network of caustics in the magnification diagram produced by MACHOs in the lens galaxy. If real, microlensing fluctuations on time scales on the order of 70 days (Sect. 6.2) may signal the presence of MACHOs with masses possibly down to planetary masses. The small, 5% amplitude of this short time scale microlensing signal is consistent with previous conclusions that the luminous source may be quite large relative to the Einstein Rings (Refsdal & Stabell 1991, 1993, 1997; Refsdal et al. 2000). Because this is near to the noise level, extremely careful data acquisition and analysis is called for in determining the time delay and microlensing.

We hope in future papers to extend the time delay and microlensing analysis, using an even larger data set. A longer observational base line and maybe more statistical techniques could shed new light on the time delay issue. We also hope that our new reduction scheme (both aperture and PSF photometry, see Ovaldsen 2002) includes some new features which could be of interest to other researchers.

New aperture photometry of QSO 0957+561; application to time delay and microlensing

1 Introduction

2 Data set

3 Aperture photometry scheme

3.1 Field star photometry and calibration

3.2 Photometry of the two quasar components

3.2.1 Galaxy subtraction

3.2.2 Crosstalk correction

4 Photometric results

4.1 Stellar photometry

4.2 Quasar photometry

4.2.1 Binning and censoring

4.2.2 Crosstalk and galaxy contamination

4.2.3 Quasar light curves

5 Time delay

5.1 Dispersion estimates

5.1.1 Complete data set

5.1.2 Truncated data sets

5.1.3 Time delay from RES's data

5.2 $\chi ^2$ minimization

6 Microlensing

6.1 Trends on long time scales

6.2 Short time scale variation

7 Summary and discussion

References

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F14.ps} \end{figure}$	Figure 14: $\chi ^2$ as a function of time shift.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{H4015F16.eps}\par\vspace{... ...s} \vspace{2mm} \includegraphics[width=8.8cm,clip]{H4015F18.eps} \end{figure}$	Figure 16: The ML residuals from the three first seasons of Fig. 15. For two periods we have fitted a cubic spline (upper panel; 9 nodes, lower panel; 4 nodes).
Open with DEXTER

New aperture photometry of QSO 0957+561; application to time delay and microlensing

2 Data set

3 Aperture photometry scheme

3.1 Field star photometry and calibration

3.2 Photometry of the two quasar components

4 Photometric results

5 Time delay

5.1 Dispersion estimates

5.1.2 Truncated data sets

5.2 minimization

6 Microlensing

6.1 Trends on long time scales

6.2 Short time scale variation

7 Summary and discussion

5.2 $\chi ^2$ minimization