© ESO, 2011

1. Introduction

The determination of the time delay between different images of gravitationally lensed quasars is an important step in different kinds of studies: in deriving H₀, the expansion rate of the Universe (e.g. Vuissoz et al. 2008), for microlensing studies (e.g. Paraficz et al. 2006a), and for detailed studies of the structure of a lensed quasar (e.g. Goicoechea 2002; Morgan et al. 2008a).

However, previous time delay determinations have been far from homogeneous, not only because they are based on different methods, but also because of their varying levels of reliability. Their accuracy depends, among other factors, on the amplitude and shape of the quasar’s intrinsic variations, the perturbations of the light curves by microlensing effects, the photometric error bars, the typical time sampling of the monitoring, the total time span of the observing campaign, and the intervals between observing seasons. Moreover, the published error bars are generally internal errors only and the way in which they are determined varies from study to study. Unfortunately, once a time delay has been published, the value may be used for years without verification, even when the authors of the original article caution the reader that the result is not very well-constrained (e.g. Jakobsson et al. 2005).

Hence, we are of the opinion that it would be useful to re-evaluate published time delays in a number of systems using the same methods as for all of the delays, as well as for the estimate of the error bars. An idea of the robustness of the time delay values can thus be obtained, not only internally with our methods, but also by comparing our results with published ones. We do not claim that our methods are superior to other ones. However, a critical reanalysis of published results using two fundamentally different approaches allows us to sort the results in terms of the reliability and independence of the method and to determine which lensed systems may be useful for determining H₀. Our main purpose is thus to examine whether the published light curves allow the determination of reliable time delays. If the answer is positive, we attempt to estimate realistic error bars, and correct some of the published values for small systematic errors.

For the determination of H₀ by means of gravitational lensing to be competitive with more classical methods (e.g. Riess et al. 2009), we need to reach at least a comparable accuracy of  ~5% in H₀. As time delays from different lensed systems should be combined to obtain H₀, we can assume that the error in the individual time delays contributes to the statistical error in H₀. Since the time delay uncertainties are only one of several sources of error in H₀ determinations (to be added e.g. to uncertainties in the dark matter distribution in the lens), they should in any case not exceed 5%.

Sections 2 and 3 present the two methods used for time delay determination. These methods are applied to the 11 lensed systems, for which the main results are described in Sect. 4. The lenses of our sample are those for which accurate astrometry has been determined by means of the deconvolution of near-infrared Hubble Space Telescope images (Sluse et al. 2011). A summary of these results is presented in Sect. 5, together with our conclusions.

2. Numerical model fit (NMF)

We revised and improved the method described in Burud et al. (2001). The basic idea can be summarized as follows: for a series of given time delays, the method minimizes the difference between the data and a numerical-model light curve with equally spaced sampling points, while adjusting the two parameters of the difference in magnitude between the light curves and a slope that models slow linear microlensing variations. The model is smoothed by introducing the convolution of the model curve with a Gaussian r(t) of full width at half maximum comparable to the typical sampling of the observations, and this smoothing term is weighted by a Lagrange multiplier λ. The function to be minimized is: $\begin{array}{ccc}\mathit{S}& \mathrm{=}& \end{array}$(1)with $\begin{array}{ccc}{\mathit{\chi }}^{\mathrm{2}}& \mathrm{=}& \sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{A}}}{\left(\frac{{\mathit{d}}_{\mathit{A}}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}\mathrm{-}\mathit{g}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}}{{\mathit{\sigma }}_{\mathit{A}}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathrm{\right)}}\right)}^{\mathrm{2}}\\ & \mathrm{+}& \end{array}$(2)as used in Burud et al. (2001) where d_k(t_i) and σ_k(t_i) are the data for image k(k = A,B) with the associated error bar, g(t_i) the model curve, Δt the time delay, and Δm and α the parameters representing the difference in magnitude and slope between the light curves.

The optimal time delay is the one that minimizes the reduced ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ between the model and the data points. It is important to insist on the difference between χ² as defined in Eq. (2) and on the other hand the reduced ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}_{\mathit{A}}\mathrm{+}{\mathit{N}}_{\mathit{B}}}{\mathit{\chi }}^{\mathrm{2}}\mathit{,}$(3)in which χ² is divided by the number of data points in common¹ between the light curves of the quasar images for a given time delay. Indeed, the longer the time delay one tests, the fewer points these light curves have in common, which tends to reduce the χ² and result in a bias towards longer time delays, hence the use of ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ to avoid this bias.

A second important difference from the original version of the method is technical: for computational reasons, the length of the model curve should be a power of two, which in some cases proves to be too long in comparison to the data, thus falsifying the balance between data and smoothing terms. In the original version the smoothing term was applied to the full length of the model. We adapted the program in such a way that the part of the model that is only needed to complete the length until the next power of two, is no longer taken into account in the minimization process. In this way, the method becomes independent of the number of data points in the light curve, which was not the case in its original form.

This method has not only been implemented for two light curves, but also for deriving time delays from three and four light curves simultaneously. The strength of this simultaneous approach lies not only in the improved constraints on the model, but also in that we assume the coherence between pairs of time delays, differences in magnitudes, and the slope parameter values.

The robustness of the measured time delay is tested in two ways. First of all, we iteratively attempt to find the three parameters of the model light curve: the spacing of the model curve’s sampling points, the range of the smoothing term, and the Lagrange multiplier. The results should be independent of these parameters as long as we remain in a certain range adapted to the data.

In a second step, we wish to test the influence of each individual point of the light curve on the time delay. This is achieved by means of a classical jackknife test: for a light curve consisting of N data points, we recalculate N times the time delay in the light curve of N − 1 data by successively leaving out one data point at a time. Time delays should not change drastically because of the removal of a single point from the light curve. If they do, we know which data point is responsible for the change and we can have a closer look at it.

Errors are calculated by means of Monte Carlo simulations. Normally distributed random errors with the appropriate standard deviation are added to the model light curve and the time delay is redetermined. We note that errors are not added to the data as they already contain the observed error, so adding another would bias the results. The model, to which the measurement errors are added, is assumed to provide a more accurate description of the real light curve of the quasar than the data. This procedure is repeated at least 1000 times, preferably on different combinations of smoothing parameters. The mean value of the time delay distribution that we obtain is considered to be the final time delay and its dispersion represents the 1σ error bar. When we have a markedly asymmetrical distribution, we take its mode as the final time delay and use the 68% confidence intervals to obtain error bars. In this paper, all quoted uncertainties are 1σ error bars except where mentioned explicitly.

The advantages of this method are manifold. First, none of the light curves is taken as a reference curve; they are all treated on an equal basis. Second, a model light curve is obtained for the intrinsic variations in the quasar, which is also the case for the polynomial fit method described by Kochanek et al. (2006), but not for the minimum dispersion method developed by Pelt et al. (1996). This is important when calculating the error bars, thus avoiding adding random errors to the data. Finally, since the model is purely numerical, no assumption is made about the quasar’s intrinsic light curve, except that it is sufficiently smooth, and we only interpolate the model, never the data.

3. Minimum dispersion (MD)

The second method we use is derived from the minimum dispersion method by Pelt et al. (1996) with a number of adjustments as described in Courbin et al. (2011). The main improvements to the original Pelt et al. (1996) method consist in:

• 1.

  No light curve is taken as a reference, they are all treated on anequal basis;

• 2.

  A flexible modelling of microlensing by polynomials up to third order, per light curve or per observing season.

Since no model light curve is constructed, computation time is a lot shorter than for the NMF method. By using two methods based on completely different principles, we are able to check whether the derived time delays are independent of the method, thus testing their robustness (i.e. independence of the particular way in which the data are analysed).

4. Application to 11 lensed quasars

We now present the main results of our time delay analysis for each of the published light curves of 11 gravitationally lensed quasars.

Fig. 1 Light curves at 8.4 GHz and 15 GHz for JVAS B0218+357 after transforming the flux density measurements into magnitudes. The B curve has been shifted by 1 magnitude for clarity.

• JVAS B0218+357: we have used the data set publishedby Cohen et al. (2000),consisting of 51 flux density measurements at 8.4GHz and 15 GHz. The authors obtained a time delay$\mathrm{\Delta }{\mathit{t}}_{\mathit{AB}}\mathrm{=}\mathrm{10.}{\mathrm{1}}_{-1.6}^{\mathrm{+}\mathrm{1.5}}$ days, where A is the leading image, thus confirming independently two values published earlier by Biggs et al. (1999) and Corbett et al. (1996) of Δt_AB = 10.5 ± 0.4 days and Δt_AB = 12 ± 3 days, respectively. After transforming the flux densities onto a logarithmic scale as shown in Fig. 1, we applied the NMF method to the 8.4 GHz and 15 GHz light curves. Using the entire light curve did not give a clear and unique solution. The jackknife test shows that certain data points can change the value of the time delay. After eliminating three of these points, the 9th, 12th and 35th, from the 8.4 GHz light curve, and choosing appropriate smoothing parameters, we obtain a time delay $\mathrm{\Delta }{\mathit{t}}_{\mathit{AB}}\mathrm{=}\mathrm{9.}{\mathrm{8}}_{-0.8}^{\mathrm{+}\mathrm{4.2}}$ days at 68% confidence level. The larger error bars for higher values of the time delay are due to a secondary peak in the histogram (see Fig. 2) around Δt_AB ~ 14 days.

  Fig. 2 Histogram of 1000 runs of the NMF method for the 8.4 GHz data light curve of JVAS B0218+357 leaving out three deviating points. One day gaps in the histogram are artefacts due to the quasi-periodicity of the data.

  Taking into account all points of the 15 GHz light curve provided a comparable value of the time delay of $\mathrm{\Delta }{\mathit{t}}_{\mathit{AB}}\mathrm{=}\mathrm{11.}{\mathrm{1}}_{-1.1}^{\mathrm{+}\mathrm{4.0}}$, even though we noted that the importance of the secondary peak around Δt_AB ~ 14 days was significantly lower after we had eliminated outlying points in both the A and B curve, in the same way as for the 8.4 GHz curve. This suggests that the secondary peak around Δt_AB ~ 14 days is probably caused by artefacts in the data, hence we can confirm with confidence the previously published results: combining the values based on the 8.4 GHz and 15 GHz light curves gives a time delay of $\mathrm{\Delta }{\mathit{t}}_{\mathit{AB}}\mathrm{=}\mathrm{9.}{\mathrm{9}}_{-0.9}^{\mathrm{+}\mathrm{4.0}}$ days.

  The MD method confirms the strong influence of these deviating points: the secondary peak around Δt_AB ~ 14 days even completely disappears when they are removed from both the 8.4 GHz and the 15 GHz light curves. The 8.4 GHz curve gives a time delay of Δt_AB = 12.6 ± 2.9 days, and the 15 GHz data lead to Δt_AB = 11.0 ± 3.5 days, which gives a combined result of Δt_AB = 11.8 ± 2.3 days, all in agreement with the above-mentioned values.

  Even if all of these time delay values for this object are in agreement with each other, the data do not allow a precision of the order of 5% in the delays, which would be necessary for a useful estimate of H₀.

• SBS 0909+523: we used the data set published by Goicoechea et al. (2008), which contains 78 data points spread over two observing seasons. Their analysis leads to a time delay Δt_BA = 49 ± 6 days where B is the leading image, confirming the previously reported delay $\mathrm{\Delta }{\mathit{t}}_{\mathit{BA}}\mathrm{=}\mathrm{4}{\mathrm{5}}_{-1}^{\mathrm{+}\mathrm{11}}$ of Ullán et al. (2006).

  The NMF method, when applied to the entire light curve, gives a delay Δt_BA ~ 47 days, as displayed in Fig. 3, which is within the error bars of the previously published delay. On closer inspection however, we note that this delay strongly depends on two points that are outside the general trend of the lightcurve for image B and fall right at the end of the time interval covered by the A data points for this time delay value: the 63rd and the 64th data points. Recalculating the delay while omitting these two points gives a different result of Δt_BA ~ 40 days or even lower values, which is not within the published ranges. The same happens if we only take into account the second observing season, which is the longer one: the delay then shortens to Δt_BA ~ 40 days. The parameter modelling slow linear microlensing is also significantly smaller in this case. Visually, both results, with and without the two problematic points, are acceptable. Nevertheless, although both values have proven to be independent of the two smoothing parameters, the NMF method is sensitive to the addition of normally distributed random errors with the appropriate standard deviation at each point of the model curve. This is because the dispersion in the data points is too small compared to the published error bars. That we obtain ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\mathrm{\ll }\mathrm{1}$ also highlights some possible problems in the data reduction or analysis.

  Fig. 3 Light curves of SBS0909+523: A is shifted by a time delay Δt = 47 days and a difference in magnitude Δm =  −0.6656. The slope parameter α = 6.0587 × 10-5 corresponds to a model of slow linear microlensing. The encircled points are the 63rd and 64th observations that were omitted from later tests.

  The MD method gives similar results. When using all data points, two possible delays can be seen, depending on the way microlensing is modelled: Δt_BA ~ 49 and Δt_BA ~ 36. When the two aforementioned data points are left out, we only find Δt_BA ~ 36, independently of how microlensing is handled.

  In all cases, with or without these points, leaving more or less freedom for the microlensing parameters, the large photometric error bars result in very large error bars in the time delay when adding normally distributed random errors to the light curves, so that delays ranging from Δt_BA ~ 27 to Δt_BA ~ 71 are not excluded at a 1 σ level.

  We conclude that this light curve does not allow a reliable determination of the time delay. To determine whether these two data points that do not follow the general trend are due to genuine quasar variations and thus crucial for the time delay determinations or whether in contrast, they are affected by large errors and contaminate the published results, we will need new observations and an independent light curve.

• RX J0911+0551: data for this quadruply lensed quasar were made available by Paraficz et al. (2006b), but had been previously treated and analysed by Burud (2001) and Hjorth et al. (2002), who proposed time delays of Δt_BA = 150 ± 6 days and Δt_BA = 146 ± 4 days respectively, where B is the leading image of the system and A the sum of the close components A1, A2, and A3.

  Fig. 4 Light curves of RXJ0911+0551: A, which is the sum of the close components A1, A2 and A3, has been shifted by one magnitude for clarity.

  Using all data except the first point, which has too strong an influence on our slope parameter because of its isolation as can be seen in Fig. 4, we can at first sight confirm the published delays: the NMF method gives Δt_BA = 150 ± 2.6 days and the MD method results in Δt_BA = 147.4 ± 4.6 days. However, the histogram in Fig. 5 shows a secondary peak at Δt_BA ~ 157 days. Investigating this peak further, we come to the conclusion that some points have a very strong influence on the delay: the first observing season, and especially the first ten points of the light curve, indicate a shorter time delay. According to Burud (2001), these points were added to supplement the regular monitoring data of the Nordic Optical Telescope. However, the first three points of the regular NOT monitoring in the B curve are similarly crucial. Omitting these three points leads to larger error bars of Δt_BA = 151.6 ± 7.0 days using the NMF method. Finally, recalculating the time delay in the regular monitoring data only and without the first three points in the B curve, gives Δt_BA = 159 ± 2.4 days with the NMF method. The MD method results in this case in a histogram with two Gaussian peaks, one around Δt_BA ~ 146 days and one around Δt_BA ~ 157 days, implying a mean time delay of Δt_BA = 151.4 ± 6.7 days. Only a new and independent light curve of similar length could tell us with more confidence which of these values is correct and which is possibly biased (e.g. by microlensing).

  Fig. 5 Sum of three histograms of 1000 runs each for RXJ0911+05, using three different combinations of smoothing parameters.

• FBQS J0951+2635: we used the data set containing 58 points published by Paraficz et al. (2006b) and presented in Fig. 6.

  Fig. 6 Light curves of FBQ0951+2635. The B curve has been shifted by 0.8 mag for clarity.

  Jakobsson et al. (2005) published a time delay Δt_AB = 16 ± 2 days, a result that is only based on the last 38 points of the light curve when the system had been observed more intensively. They found various possible time delays according to the method, the smoothing, and the data points included, so we performed the same tests.

  We can confirm that the time delay is very sensitive to the choice of smoothing parameters in the NMF method, especially when using the entire light curve, but is still more sensitive to the data points used: leaving out a single point completely changes the time delay. We calculated time delays in the light curve using between 55 and 58 data points and we found delays ranging from Δt_AB = 14.2 ± 4.5 days to Δt_AB = 26.3 ± 4.7 days. Taking into account three possible smoothing combinations and four sets of data (leaving out one more data point in each set) leads to a combined histogram of 12000 Monte Carlo simulations, as shown in Fig. 7. It is clear that a mean value with error bars Δt_AB = 20.1 ± 7.2 days is not of any scientific use: the error bars are too large relative to the time delay. Moreover, the histogram is quite different from a normal distribution. There is no significant concentration of the results, which would allow the determination of a meaningful mode, independently of the chosen binning. One can see that different time delays are possible and can be divided in two groups: shorter values of  ~10.5,  ~15, and  ~18.5 days, and longer values of  ~26.5 and 30.5 days.

  Fig. 7 Sum of twelve histograms of 1000 runs each for FBQ0951+2635, using three different combinations of smoothing parameters for four sets of data consisting of 58, 57, 56, and 55 data points.

  When using only the third observing season, which is more finely sampled, a relatively stable time delay Δt_AB = 18.8 ± 4.5 is measured, but once a single point is left out (for example the 19th point of this third season, a point that deviates from the general trend in spite of a small error bar), the result completely changes towards longer values (Δt_AB = 25.0 ± 4.9 days) and becomes sensitive to smoothing. As the measured time delay should not depend on the presence or absence of a single point, we can only conclude that this light curve, even if it consists of three observing seasons, does not allow a precise determination of this delay.

  The MD method entirely confirms the large uncertainty in this time delay: using all data points we find a time delay of Δt_AB = 21.5 ± 6.8 days, whereas the third season only leads to Δt_AB = 19.6 ± 7.6 days, with peaks in the histogram around Δt_AB ~ 12 and Δt_AB ~ 28 days.

  According to Schechter et al. (1998) and Jakobsson et al. (2005), there are spectroscopic indications of possible microlensing, so this might explain the difficulty in constraining the time delay for this system. Longer and more finely sampled light curves might help us to disentangle both effects. However, at the present stage, we can conclude that this system is probably not suitable for a time delay analysis.

• HE 1104-1805: we used the data published by Poindexter et al. (2007), which combine their own SMARTS R-band data with Wise R-band data from Ofek & Maoz (2003) and OGLE V-band data from Wyrzykowski et al. (2003). The three data sets are shown in Fig. 8. Table 1 lists the four time-delay values published for HE1104-1805, where B is the leading image.

  Fig. 8 Light curves of HE1104-1805, combining the OGLE V-band data, the Wise R-band data and the SMARTS R-band data.

  Table 1

  Published time delays for HE1104-1805.

  We performed tests with both methods on different combinations of the data using all telescopes or only one or two of them. Unfortunately, the results seem to be sensitive to this choice, as they are to the way in which microlensing is treated: both the OGLE and Wise data sets were analysed to find a time delay Δt_BA ~ 157 days, whereas SMARTS data converge to a higher value of Δt_BA ~ 161 days or more as is shown in Fig. 9. In addition, Poindexter et al. (2007)’s smaller value is recovered with the MD method when including OGLE and Wise data but only for some ways of modelling microlensing. We therefore conclude that we can neither make a decisive choice between the published values, nor improve their error bars, which are large enough to overlap.

  Fig. 9 Sum of three histograms for HE1104-1805, using three different combinations of smoothing parameters combining only OGLE and SMARTS data. OGLE data point to a time delay ΔtBA ~ 157 days, whereas SMARTS data converge to a longer value of ΔtBA ~ 161 days.

• PG 1115+080: we used the data taken by Schechter et al. (1997). They published a time delay Δt_CB = 23.7 ± 3.4 days between the leading curve C and curve B and estimated the delay between C and the sum of A1 and A2 at Δt_CA ~ 9.4 days. Barkana (1997) used the same data, which are shown in Fig. 10, but a different method to determine the delays. His value of $\mathrm{\Delta }{\mathit{t}}_{\mathit{CB}}\mathrm{=}\mathrm{25.}{\mathrm{0}}_{-3.8}^{\mathrm{+}\mathrm{3.3}}$ is compatible with Schechter’s one. Morgan et al. (2008b) published new optical light curves for this quadruply imaged quasar in order to study microlensing in the system. Unfortunately, these light curves cannot be used to determine a time delay independently because of the clear lack of features in the variability of the quasar and the inconsistency of the individual error bars relative to the dispersion in the data.

  Fig. 10 Light curves of PG 1115+080, where the A curve is the sum of the A1 and A2 component. The A and the B curve have both been shifted by 1.9 and  − 0.3 mag, respectively.

  A first series of tests on Schechter’s data with the NMF method led to time delays of Δt_CA ~ 15 days and Δt_CB ~ 20.8 days with a minor secondary peak around Δt_CB ~ 23.8. We then corrected the published data for the existing photometric correlation between the quasar images and the two stars used as photometric references, as mentioned by Barkana (1997). This caused the shorter time delay to shift either towards Δt_CA ~ 11 days or Δt_CA ~ 16 days, and transformed the longer delay into two nearly equally possible results of Δt_CB ~ 20.8 days or Δt_CB ~ 23.8, as indicated by the two main peaks in the histogram in Fig. 11. Adding observational errors to the model light curves and taking into account four different ways of smoothing results in Δt_CA = 11.7 ± 2.2 (see Fig. 12) and $\mathrm{\Delta }{\mathit{t}}_{\mathit{CB}}\mathrm{=}\mathrm{23.}{\mathrm{8}}_{-3.0}^{\mathrm{+}\mathrm{2.8}}$ (see Fig. 11).

  Fig. 11 Sum of four histograms of 1000 runs each for ΔtCB in PG1115+080, using four different combinations of smoothing parameters.

  Fig. 12 Sum of four histograms of 1000 runs each for ΔtCA in PG1115+080, using four different combinations of smoothing parameters.

  The MD method confirms a delay of Δt_CB ~ 20.0 days, but finds a second solution around Δt_CB ~ 12.0 days, which results in Δt_CB = 17.9 ± 6.9. The value for Δt_CA is also shorter than the one obtained with the other methods namely Δt_CA = 7.6 ± 3.9 days and has larger error bars. Unfortunately, the length and quality of this light curve do not allow one to choose between the possible time delays that differ according to the method used, but are generally lower than published values.

• JVAS B1422+231: for this quadruply lensed quasar, we used the data published by Patnaik & Narasimha (2001), consisting of flux density measurements at two frequencies, 8.4 and 15 GHz. Their results for the time delays were based only on the 15 GHz data without image D, which is too faint. These data are shown in Fig. 13. They obtained Δt_BA = 1.5 ± 1.4 days, Δt_AC = 7.6 ± 2.5, and Δt_BC = 8.2 ± 2.0 days when comparing the curves in pairs.

  Fig. 13 Light curves of JVAS B1422+231 for the A, B, and C components after transforming the flux density measurements at the 15 GHz frequency into magnitudes. The C curve has been shifted by half a magnitude.

  The NMF method allows time delays to be tested for the three light curves simultaneously, thus imposes coherence on the results. Given the error bars in the published results, which are large compared to the time delays, we can confirm the results here, but we emphasize that they include two distinct groups of solutions between which we cannot decide based on the actual light curves: for the shortest delay, we either have Δt_BA ~ 1.0 day or Δt_BA ~ 2.0 days, as shown in Fig. 14. The choice between both solutions is sensitive to the smoothing parameters: the importance of the first group Δt_BA ~ 1.0 day is lower, and even disappears completely, with greater smoothing.

  Fig. 14 Sum of four histograms for ΔtBA of 1000 runs each for JVAS B1422+231, using four different combinations of smoothing parameters.

  For Δt_BC, the situation is similar but even less clear: low smoothing parameters lead to a range of possible solutions between Δt_BC ~ 6 and Δt_BC ~ 10 days, which are all within the error bars of the published results. However, Monte Carlo simulations of reconstructed light curves with higher smoothing parameters give a time delay Δt_BC = 10.8 ± 1.5 days with a secondary peak around Δt_BC ~ 8 days.

  The MD method gives a completely different result: it converges towards time delays that invert the BAC-order into CAB but with error bars large enough not to exclude the BAC order of Δt_BA =  −1.6 ± 2.1 days, Δt_AC =  −0.8 ± 2.9, and Δt_BC =  −2.4 ± 2.7 days.

  New observations are clearly necessary to reduce the uncertainties in both the different solutions and the error bars that are too large in comparison with the time delays to be useful to any further analysis.

• SBS 1520+530: two data sets exist for this doubly lensed quasar: the set made available by Burud et al. (2002b) and the one published by Gaynullina et al. (2005a). Burud et al. (2002c) were the first to publish a time delay for this system Δt_AB = 128 ± 3 days, where A is the leading image, or Δt_AB = 130 ± 3 days when using the iterative version of the method (Burud et al. 2001). Gaynullina et al. (2005b) used an independent data set and found four possible time delays, of which the one with the largest statistical weight, Δt_AB = 130.5 ± 2.9, is perfectly consistent with the previously published time delays.

  Even if the light curve based on Gaynullina et al. (2005a) data contains more than twice as many data points as Burud et al. (2002b)’s older light curve, we decided not to use it because of the lack of overlapping data between the A and the B curves of the quasar after shifting the B curve for the time delay, as can be seen in Fig. 15.

  Fig. 15 Light curves for SBS 1520+530 based on Gaynullina et al. (2005a)’s data after shifting the B curve by ΔtAB = 130.5 days. Hardly any data points overlap between the light curves for images A and B.

  Applying different tests to Burud et al. (2002c)’s light curves, shown in Fig. 16, using the NMF method led to a time delay for which the error bars overlap with the published value: Δt_AB = 126.9 ± 2.3. That the delay is slightly shorter than Burud et al. (2002c)’s value can be explained by our use of the reduced ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ instead of the χ², the latter implying that longer delays are the more likely ones, as explained in Sect. 2. This effect was also noted using the iterative method.

  Fig. 16 Light curves for SBS 1520+530 based on Burud et al. (2002c)’s data.

  Fig. 17 Sum of three histograms of 1000 runs each for SBS1520+530, using three different combinations of smoothing parameters.

  The MD method yields a comparable time delay of Δt_AB = 124.6 ± 3.6 days and confirms the shape of the histogram: the highest peak value at Δt_AB ~ 125 days and a clear secondary peak at Δt_AB ~ 127.5 days. Combining the values from both methods implies that Δt_AB = 125.8 ± 2.1 days.

• CLASS B1600+434: data, as shown in Fig. 18, were made available by Paraficz et al. (2006b) but had been treated and analysed by Burud et al. (2000), who published a final time delay Δt_AB = 51 ± 4 days, where A is the leading image, consistent with the time delay $\mathrm{\Delta }{\mathit{t}}_{\mathit{AB}}\mathrm{=}\mathrm{4}{\mathrm{7}}_{-9}^{\mathrm{+}\mathrm{12}}$ days from Koopmans et al. (2000) based on radio data.

  Fig. 18 Light curves for CLASS B1600+434: 41 data points spread over nearly two years.

  The NMF method leads to a time delay Δt_AB = 46.6 ± 1.1 days, but the histogram in Fig. 19 clearly shows that we cannot use the mean as the final value. The histogram has two distinct values of  ~46 or  ~ 48 days, so we prefer to speak of a delay of either $\mathrm{\Delta }{\mathit{t}}_{\mathit{AB}}\mathrm{=}\mathrm{45.}{\mathrm{6}}_{-0.4}^{\mathrm{+}\mathrm{1.2}}$ days (68% error) or $\mathrm{\Delta }{\mathit{t}}_{\mathit{AB}}\mathrm{=}\mathrm{45.}{\mathrm{6}}_{-0.4}^{\mathrm{+}\mathrm{2.8}}$ days (95% error).

  Fig. 19 Sum of seven histograms of 1000 runs each for B1600+434, using seven different combinations of smoothing parameters.

  These results are in marginal disagreement with the final delay proposed by Burud et al. (2000). However, Burud et al. (2000)’s final result is an average of four time delays, each of them calculated with a different method. Two of these methods also inferred a value of around  ~48 days.

  We could identify at least three explanations of our lower value and smaller error bars in comparison with Burud et al. (2000)’s time delay. The first one is the same again as for SBS 1520+530: our use of the reduced ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ (see formula 3 in Sect. 2) instead of the χ², the latter introducing a bias towards longer delays. The second reason is the technical issue concerning the length of the model curve as explained in Sect. 2, which was found to be crucial for the time delay of this system. Finally, we observed that higher values of the Lagrange multiplier weighting the smoothing term seemed to lead to longer time delay values, which disappeared with lower smoothing. Taking into account these three adjustments, nearly all values around  ~51 days disappear from the histogram.

  This is not the case for the MD method, which explains the slightly longer value of the time delay: Δt_AB = 49.0 ± 1.2 days. Combining these results gives a delay of Δt_AB = 47.8 ± 1.2 days. Even if these error bars imply that the time delay is very tightly constrained, we emphasize that the delay measurement is only based on 41 data points spread over nearly two years, which gives a relatively high weight to every single data point. When adding random errors, neither of the two methods leads to a histogram with a Gaussian shape. A more finely sampled light curve might remedy this situation.

• CLASS B1608+656: light curves for this quadruply lensed system were first analysed by Fassnacht et al. (1999) and subsequently improved in Fassnacht et al. (2002) by adding more data. Using three observing seasons, they published time delays of $\mathrm{\Delta }{\mathit{t}}_{\mathit{BA}}\mathrm{=}\mathrm{31.}{\mathrm{5}}_{-1}^{\mathrm{+}\mathrm{2}}$, Δt_BC = 36.0 ± 1.5, and $\mathrm{\Delta }{\mathit{t}}_{\mathit{BD}}\mathrm{=}\mathrm{77.}{\mathrm{0}}_{-1}^{\mathrm{+}\mathrm{2}}$ days. Their analysis is based on a simultaneous fit to data from the three seasons but treats the curves only in pairs.

  Fig. 20 Light curves for CLASS B1608+656. The B curve has been shifted by half a magnitude for clarity.

  Fig. 21 Histograms for the three time delays in B1608+656, using two different combinations of smoothing parameters for all seasons and three out of the four curves simultaneously.

  We performed several tests on these light curves, which are shown in Fig. 20: first taking into account only the first and the third season separately (the second season not presenting useful structure), then all data for the three seasons simultaneously, using the three and four curve version of our method as described in Sect. 2. This enables us to impose coherence between the pairs of time delays, which was not done by Fassnacht et al. (2002). The results, as illustrated in Fig. 21, confirm the previously published values, within the error bars, of Δt_BA = 30.2 ± 0.9, Δt_BC = 36.2 ± 1.1, and Δt_BD = 76.9 ± 2.3 days. One particularity deserves more attention: the value for Δt_BA changes slightly according to the seasons and the curves considered simultaneously. When leaving out the second season data (featureless), Δt_BA systematically converges towards Δt_BA = 33.5 ± 1.5 days, as shown in Fig. 22. This is consistent with Fassnacht et al. (2002), who already mentioned a time delay of Δt_AC ~ 2.5 days. Even if this slight difference is probably due to microlensing and should be investigated in more detail, we chose to retain the final value, which is the one based on the use of all data.

  Fig. 22 Sum of four histograms of 1000 runs each for ΔtBA in B1608+656, using two different combinations of smoothing parameters on the first and the third season.

  The MD method entirely confirms these results of Δt_BA = 32.9 ± 2.9, Δt_BC = 35.2 ± 2.5, and Δt_BD = 78.0 ± 3.7 days with another indication for Δt_AC ~ 2.5 days. Combining both methods results in the time delays of Δt_BA = 31.6 ± 1.5, Δt_BC = 35.7 ± 1.4, and Δt_BD = 77.5 ± 2.2 days.

• HE 2149-2745: we reanalysed the data set made available by Burud et al. (2002b). These data consist of two light curves, one in the V-band and one in the I-band, as shown in Fig. 23. Burud et al. (2002a) published a time delay Δt_AB = 103 ± 12 days where A is the leading image. This delay is based on the V-band data, but, according to Burud et al. (2002a), agrees with the I-band data.

  Fig. 23 Light curves of HE 2149-2745 in the V-band and the I-band. The B curve has been shifted by one magnitude for clarity.

  Our tests, based on the light curves as such, and using both methods, clearly reveal two possible delays: one around  ~70−85 days and one around  ~100−110 days. Unfortunately, the light curves of images A and B show little structure and hardly overlap, except for some points in the second season, when shifting them for a delay of more than 100 days, especially in the I-band, which makes it very difficult to choose between the two possibilities. Moreover, once we add random errors to the model light curve and perform Monte Carlo simulations, we only obtain a forest of small peaks, spread over the entire tested range of 50−140 days, instead of a Gaussian distribution around one or two central peaks, demonstrating that these results are highly unstable. Leaving out two outlying data points in the B curve only slightly improves the situation: within the forest of peaks in Fig. 24, those in the range 75 − 85 seem to be slightly more important than those over 100 days. Nevertheless, we cannot derive a reliable time delay from these data sets for this system.

  Fig. 24 Sum of six histograms of 1000 iterations each for HE 2149-2745 leaving out two data points, using six different combinations of smoothing parameters.

5. Conclusions

We have presented an improved numerical method, the numerical model fit method, to calculate time delays in doubly or quadruply lensed quasar systems, and applied it to 11 systems for which time delays had been published previously. This allowed the validity of these time delay values to be evaluated in a coherent way. The use of a minimum dispersion method allowed us to check the independence of the results from the method. The results are summarized in Table 2.

We caution that some published time delay values should be interpreted with care: even if we have been able to confirm some values (time delays for JVAS B0218+357, HE 1104-1805, CLASS B1600+434, and for the three delays in the quadruply lensed quasar CLASS B1608+656) and give an improved value for one system, SBS 1520+530, many of the published time delays considered in our analysis have proven to be be unreliable for various reasons: the analysis is either too dependent on some data points, leads to multiple solutions, is sensitive to the addition of random errors, or is incoherent between the two different methods used.

Given the accuracy that is needed for time delays to be useful to further studies, we note that it will be necessary to perform long-term monitoring programs on dedicated telescopes to obtain high-quality light curves of lensed quasars, not only for new systems but also for the majority of the lenses in this sample, for which the time delay has been considered to be known. The COSMOGRAIL collaboration, which has been observing over 20 lensed systems for several years now, will soon be improving the time delay values for some of these systems for which the accuracy is unsatisfactory.

Acknowledgments

We would like to thank Sandrine Sohy for her very precious help in programming. This work is supported by ESA and the Belgian Federal Science Policy (BELSPO) in the framework of the PRODEX Experiment Arrangement C-90312.

References

All Tables

All Figures

	Fig. 1 Light curves at 8.4 GHz and 15 GHz for JVAS B0218+357 after transforming the flux density measurements into magnitudes. The B curve has been shifted by 1 magnitude for clarity.
In the text

	Fig. 2 Histogram of 1000 runs of the NMF method for the 8.4 GHz data light curve of JVAS B0218+357 leaving out three deviating points. One day gaps in the histogram are artefacts due to the quasi-periodicity of the data.
In the text

	Fig. 3 Light curves of SBS0909+523: A is shifted by a time delay Δt = 47 days and a difference in magnitude Δm =  −0.6656. The slope parameter α = 6.0587 × 10-5 corresponds to a model of slow linear microlensing. The encircled points are the 63rd and 64th observations that were omitted from later tests.
In the text

	Fig. 4 Light curves of RXJ0911+0551: A, which is the sum of the close components A1, A2 and A3, has been shifted by one magnitude for clarity.
In the text

	Fig. 5 Sum of three histograms of 1000 runs each for RXJ0911+05, using three different combinations of smoothing parameters.
In the text

	Fig. 6 Light curves of FBQ0951+2635. The B curve has been shifted by 0.8 mag for clarity.
In the text

	Fig. 7 Sum of twelve histograms of 1000 runs each for FBQ0951+2635, using three different combinations of smoothing parameters for four sets of data consisting of 58, 57, 56, and 55 data points.
In the text

	Fig. 8 Light curves of HE1104-1805, combining the OGLE V-band data, the Wise R-band data and the SMARTS R-band data.
In the text

	Fig. 9 Sum of three histograms for HE1104-1805, using three different combinations of smoothing parameters combining only OGLE and SMARTS data. OGLE data point to a time delay ΔtBA ~ 157 days, whereas SMARTS data converge to a longer value of ΔtBA ~ 161 days.
In the text

	Fig. 10 Light curves of PG 1115+080, where the A curve is the sum of the A1 and A2 component. The A and the B curve have both been shifted by 1.9 and  − 0.3 mag, respectively.
In the text

	Fig. 11 Sum of four histograms of 1000 runs each for ΔtCB in PG1115+080, using four different combinations of smoothing parameters.
In the text

	Fig. 12 Sum of four histograms of 1000 runs each for ΔtCA in PG1115+080, using four different combinations of smoothing parameters.
In the text

	Fig. 13 Light curves of JVAS B1422+231 for the A, B, and C components after transforming the flux density measurements at the 15 GHz frequency into magnitudes. The C curve has been shifted by half a magnitude.
In the text

	Fig. 14 Sum of four histograms for ΔtBA of 1000 runs each for JVAS B1422+231, using four different combinations of smoothing parameters.
In the text

	Fig. 15 Light curves for SBS 1520+530 based on Gaynullina et al. (2005a)’s data after shifting the B curve by ΔtAB = 130.5 days. Hardly any data points overlap between the light curves for images A and B.
In the text

	Fig. 16 Light curves for SBS 1520+530 based on Burud et al. (2002c)’s data.
In the text

	Fig. 17 Sum of three histograms of 1000 runs each for SBS1520+530, using three different combinations of smoothing parameters.
In the text

	Fig. 18 Light curves for CLASS B1600+434: 41 data points spread over nearly two years.
In the text

	Fig. 19 Sum of seven histograms of 1000 runs each for B1600+434, using seven different combinations of smoothing parameters.
In the text

	Fig. 20 Light curves for CLASS B1608+656. The B curve has been shifted by half a magnitude for clarity.
In the text

	Fig. 21 Histograms for the three time delays in B1608+656, using two different combinations of smoothing parameters for all seasons and three out of the four curves simultaneously.
In the text

	Fig. 22 Sum of four histograms of 1000 runs each for ΔtBA in B1608+656, using two different combinations of smoothing parameters on the first and the third season.
In the text

	Fig. 23 Light curves of HE 2149-2745 in the V-band and the I-band. The B curve has been shifted by one magnitude for clarity.
In the text

	Fig. 24 Sum of six histograms of 1000 iterations each for HE 2149-2745 leaving out two data points, using six different combinations of smoothing parameters.
In the text