© ESO, 2013

1. Introduction

In the current cosmological paradigm, only a handful of parameters seem necessary to describe the Universe on the largest scales and its evolution over time. Testing this cosmological model requires a range of experiments, characterized by different sensitivities to these parameters. These experiments, or cosmological probes, are all affected by statistical and systematic errors and none of them on its own can uniquely constrain the cosmological models. This is due to the degeneracies inherent in each specific probe, implying that the probes become truly effective in constraining cosmology only when used in combination.

The latest cosmology results by the Planck consortium beautifully illustrate this (Planck Collaboration 2013). In particular, the constraints obtained by Planck on the Hubble parameter H₀, on the curvature Ω_k, and on the dark energy equation of state parameter w rely mostly on the combination of the baryonic acoustic oscillations measurements (BAO) with the cosmic microwave background (CMB) observations.

Strong gravitational lensing offers a valuable yet inexpensive complement to independently constrain some of the cosmological parameters, through the measurement of the so-called time delays in quasars strongly lensed by a foreground galaxy (Refsdal 1964). The principle of the method is the following. The travel times of photons along the distinct optical paths forming the multiple images are not identical. These travel-time differences, called the time delays, depend on the geometrical differences between the optical paths, which contain the cosmological information, and on the potential well of the lensing galaxy(ies). In practice, time delays can be measured from photometric light curves of the multiple images of lensed quasar: if the quasar shows photometric variations, these are seen in the individual light curves at epochs separated by the time delay.

A precise and accurate measurement of such a time delay, in combination with a well-constrained model for the lensing galaxy, can be used to constrain cosmology in a way which is very complementary to other cosmological probes (see, e.g., Linder 2011). A recent and remarkable implementation of this approach can be found in Suyu et al. (2013a) that uses the time-delay measurements from Tewes et al. (2013b). We note, however, that to obtain a robust cosmological inference from this time-delay technique, particular attention must be paid to any possible lens model degeneracies (Schneider & Sluse2013a,b; Suyu et al. 2013b).

So far, only a few quasar time delays have been measured convincingly, from long and well-sampled light curves. The international COSMOGRAIL¹ (COSmological MOnitoring of GRAvItational Lenses) collaboration is changing this situation by measuring accurate time delays for a large number of gravitationally lensed quasars. The goal of COSMOGRAIL is to reach an accuracy of less than 3%, including systematics, for most of its targets.

In this paper, we present the time-delay measurement for the two-image gravitationally lensed quasar SDSS J1001+5027 (α₂₀₀₀ = 10:01:28.61, δ₂₀₀₀ = +50:27:56.90) at z = 1.838 (Oguri et al. 2005). The image separation of Δθ = 2.86′′ (Oguri et al. 2005) and the high declination of the target make it a relatively easy prey for medium-sized northern telescopes and average seeing conditions. The redshift of the lensing galaxy z_l = 0.415 has been measured spectroscopically (Inada et al. 2012).

Our paper is structured as follows. Section 2 describes our monitoring, the data reduction, and the resulting light curves. In Sect. 3 we present a new time-delay point estimator. We add this technique to a pool of four other existing algorithms, to measure the time delay in Sect. 4. Finally, we summarize our results and conclude in Sect. 5.

2. Observations, data reduction, and light curves

2.1. Observations

We monitored SDSS J1001+5027 in the R band for more than six years, from March 2005 to July 2011, with three different telescopes: the 1.2-m Mercator Telescope located at the Roque de los Muchachos Observatory on La Palma (Spain), the 1.5-m telescope of the Maidanak Observatory in Pamir Alai (Uzbekistan), and the 2-m Himalayan Chandra Telescope (HCT) located at the Indian Astronomical Observatory in Hanle (India). Table 1 details our monitoring observations. In total we obtained photometric measurements for 443 independent epochs, with a mean sampling interval below four days. Each epoch consists of at least three, but mostly four or more, dithered exposures. Figure 1 summarizes the image quality of our data. The COSMOGRAIL collaboration has now ceased the monitoring of this target, to focus on other systems.

Fig. 1 Distribution of the average observed FWHM and elongation ϵ of field stars in the images used to build the light curves of SDSS J1001+5027.

Fig. 2 R-band image centered on SDSS J1001+5027. The image is the combination of the 210 best exposures from the Mercator telescope, for a total exposure time of 21 h. We use the stars labeled 1, 2, and 3 to model the PSF and to cross-calibrate the photometry of each exposure. The position of the two lensing galaxies G1 and G2 are indicated in the zoomed image in the upper left. They are most clearly seen in the deconvolved images presented in Fig. 3.

Fig. 3 Two ways of modeling the light distribution for extended objects during the deconvolution process. On the left is shown a single 360-s exposure of SDSS J1001+5027 obtained with the Mercator telescope in typical atmospheric conditions. The other panels show the parametric (top row) and pixelized light models (bottom row) for the lens galaxies as described in the text. The residual image for the single exposure is also shown in each case, as well as the average residuals over the 120 best exposures. The residual maps are normalized by the shot noise amplitude. The dark areas indicate excess flux in the data with respect to the model. Gray scales are linear.

2.2. Deconvolution photometry

Fig. 4 R-band light curves of the quasars images A and B in SDSS J1001+5027 from March 2005 to July 2011. The 1σ photometric error bars are also shown. For display purpose, the curve of quasar image B is shown shifted in time by the measured time delay (see text). The light curves are available in tabular form from the CDS and the COSMOGRAIL website.

The image reduction and photometry closely follows the procedure described in Tewes et al. (2013b). We performed the flat-field correction and bias subtraction for each exposure using custom software pipelines, which address the particularities of the different telescopes and instruments.

Figure 2 shows part of the field around SDSS J1001+5027, obtained by stacking the best monitoring exposures from the Mercator telescope to reach an integrated exposure time of 21 h. The relative flux measurements of the quasar images and reference stars for each individual epoch were obtained through our COSMOGRAIL photometry pipeline, which is based on the simultaneous MCS deconvolution algorithm (Magain et al. 1998). The stars labeled 1, 2, and 3 in Fig. 2 are used to characterize the point spread function (PSF) and relative magnitude zero-point of each exposure.

The two quasar images A and B of SDSS J1001+5027 are separated by 2.86″, which is significantly larger than the typical separation in strongly lensed quasars. In principle, this makes SDSS J1001+5027 a relatively easy target to monitor, as the quasar images are only slightly blended in most of our images. However, image B lies close to the primary lensing galaxy G1. Minimizing the additive contamination by G1 to the flux measurements of B therefore requires a model for the light distribution of G1. In Fig. 3, we show two different ways of modeling these galaxies. Our standard approach, shown in the bottom panels, consists in representing all extended objects, such as the lens galaxies, by a regularized pixel grid. The values of these pixels get iteratively updated during the deconvolution photometry procedure. Because of obvious degeneracies, this approach may fail when a relatively small extended object (lens galaxy) is strongly blended with a bright point source (quasar), leading to unphysical light distributions. To explore the sensitivity of our results to a possible bias of this kind, we have adopted an alternative approach of representing G1 and G2 by two simply parametrized elliptical Sersic profiles, as shown in the top row of Fig. 3. For both cases, the residuals from single exposures are convincingly homogeneous. Only when averaging the residuals of many exposures to decrease the noise can the simply parametrized models be seen to yield a less good overall fit to the data, since they cannot represent additional background sources nor compensate for small systematic errors in the shape of the PSF.

We find that the difference between these approaches in terms of the resulting quasar flux photometry is marginal; it is insignificant regarding the measurement of the time delay. In all the following we will use the quasar photometry obtained using the parametrized model (top row of Fig. 3) which is likely to be closer to reality than our pixelized model in the immediate surroundings of image B.

2.3. Light curves

Following Tewes et al. (2013b), we empirically corrected for small magnitude and flux shifts between the light curve contributions from different telescopes/cameras to obtain minimal dispersion in each of the combined light curves. In the present case we chose the photometry from the Mercator telescope as a reference, and for the data from the Maidanak and HCT telescopes, we optimized a common magnitude shift and individual flux shifts for A and B.

Figure 4 shows the combined 6.5-season long light curves, from which we measure a time delay of Δt_AB = −119.3 days (see Sect. 4). In this figure, light curve B has been shifted by this time delay to highlight the correspondence and temporal overlap of the data. We observe strong intrinsic quasar variability, common to images A and B. In the period 2006 to 2007, the variability in image A is as large as 0.25 mag over a single year. In addition to this large scale correspondence, several small and short scale intrinsic variability features are common to both curves, for instance around December 2005 and January 2010. Our data unambiguously reveal, already to the eye, an approximate time delay of Δt_AB ≈ −120 days, with A leading B.

2.4. An apparent mismatch between the light curves of the quasar images

The apparent flux ratio between the quasar images, as inferred from the time-shifted light curves shown in Fig. 4, stays roughly in the range from 0.40 to 0.44 mag over the length of our monitoring. Strong gravitational lens models readily explain different magnifications of the quasar images, yielding stationary flux ratios or magnitude shifts between the light curves. Figure 4 hints, however, at a moderate correlation between some variable flux ratio and the intrinsic quasar variability. In particular, the amplitude of the quasar variability, in units of magnitudes, appears to be smaller in B than in A. Potential reasons for this mismatch include the effects of microlensing by stars of the lens galaxy, or a contamination of the photometry of B by some additive external flux. We find that one has to subtract from curve B about 20% of its median flux to obtain an almost stationary magnitude shift of about 0.66 mag between the light curves. As this contamination would be several times larger than the entire flux of galaxy G1, we conclude that plausible errors of our light models for G1 cannot be responsible for the observed discrepancy between the light curves.

3. A new time-delay estimator

Although an unambiguous approximation of the time delay of SDSS J1001+5027 can be made by eye, accurately measuring its value is not trivial, and is made more difficult by the extrinsic variability between the light curves. Even more obvious features of the data, such as the sampling gaps due to non-visibility periods of the targets, could easily bias the results from a time-delay measurement technique. The impact of these effects on the quality of the time-delay inference clearly differs for each individual quasar lensing system and dataset. To check for potential systematic errors, we feel that a wise approach is to employ several numerical methods based on different fundamental principles.

In the present section we introduce a new time-delay estimation method, based on minimizing residuals of a high-pass filtered difference light curve between the quasar images.

Fig. 5 Difference light curves of SDSS J1001+5027 as obtained by the new difference-smoothing technique introduced in this paper. The curves are shown for the best time-delay estimate found with this technique (top panel, ΔtAB = −118.6 days), and for a wrong time-delay value (bottom panel, ΔtAB = −100.0 days). The difference light curves di are shown as colored points. They are smoothed using a kernel of width s = 100 days to compute the fi (black points). The error bars on the black points show the uncertainty coefficients σfi. The points in the difference light curve di are color-coded according to the absolute factors of their uncertainties σdi by which they deviate from fi. In both panels, light curve A is used as reference, and light curve B is shifted in flux by the same amount.

3.1. The difference-smoothing technique

This technique is a point estimator that determines both an optimal time delay and an optimal shift in flux between two light curves, while also allowing for smooth extrinsic variability. The correction for a flux shift between the light curves explicitly addresses the mismatch described in Sect. 2.4, whatever its physical explanation. This flux shift may be due to a contamination of light curve B by residual light from the lensing galaxy, from the lensed quasar host galaxy, or by microlensing resolving the quasar structure.

We consider two light curves A and B sampled at epochs t_i, where A_i and B_i are the observed magnitudes at epochs t_i, (i = 1,2,3,...,N). We select A as the reference curve. Light curve B is shifted in time with respect to A by some amount τ, and in flux by some amount Δf. Formally, this shifted version B′ of B is given by

For any estimate of the time delay τ and of the flux shift Δf, we form a difference light curve, with points d_i at epochs t_i, (3)where the weights w_ij are given by (4)The parameter δ is the decorrelation length, as in Pelt et al. (1996), and σ_Bj denotes the photometric error of the magnitude B_j. This decorrelation length should typically be of the order of the sampling period, small enough to not smooth out any intrinsic quasar variability features from the light curve B. The uncertainties on each d_i are then calculated as (5)where w_ij are given by Eq. (4). To summarize, at this point we have a discrete difference light curve, sampled at the epochs of curve A, built by subtracting from light curve A a smoothed and shifted version of B. We now smooth this difference curve d_i, again using a Gaussian kernel, to obtain a model f_i for the differential extrinsic variability (6)where the weights ν_ij are given by (7)The smoothing time scale s is a second free parameter of this method. Its value must be chosen to be significantly larger than δ. For each f_i, we compute an uncertainty coefficient (8)The idea of the present method is now to optimize the time-delay estimate τ and flux shift Δf to minimize residuals between the difference curve d_i and the much smoother f_i. Any incorrect value for τ introduces relatively fast structures that originate from the quasar variability into d_i, and these structures will not be well represented by f_i. Figure 5 illustrates this phenomenon in the case of SDSS J1001+5027 by showing d_i and f_i for an optimal and an arbitrarily chosen wrong time-delay estimate. In both panels of Fig. 5, the largest deviations between d_i and f_i are due to poorly constrained points with very high σ_di, and are therefore not significant. However, for the incorrect time-delay estimate, a larger number of well-constrained points of d_i significantly deviate from f_i (yellow and red points). To quantify this match between d_i and f_i we define a cost function in the form of a normalized χ², (9)and minimize this using a global optimization.

In the above description, light curves A and B are not interchangeable, thus introducing an asymmetry into the time-delay measurement process. To avoid this arbitrary choice of the reference curve, we systematically perform all computations for both permutations of A and B, and minimize the sum of the two resulting values of .

3.2. The uncertainty estimation procedure

As a point estimator, the technique described above does not provide information on the uncertainty of its result. We stress that simple statistical techniques such as variants of bootstrapping or resampling cannot be used to quantify the uncertainty of such highly non-linear time-delay estimators (Tewes et al. 2013a). These approaches are not able to discredit “lethargic” estimators, which favor a particular solution (or a small set of solutions) while being relatively insensitive to the actual shape of the light curves. Furthermore, they are not sensitive to plain systematic biases of the techniques.

Consequently, to quantify the random and systematic errors of this estimator, for each dataset to be analyzed and as a function of its free parameters, we follow the Monte Carlo analysis described in Tewes et al. (2013a). It consists in applying the point estimator to a large number of fully synthetic light curves, which closely mimic the properties of the observed data, but have known true time delays. It is particularly important that these synthetic curves cover a range of true time delays around a plausible solution, instead of all having the same true time delay. Only this feature enables the method to adequately penalize estimators with lethargic tendencies.

3.3. Application to SDSS J1001+5027

The decorrelation length δ and the width of the smoothing kernel s are the two free parameters of the described technique. In this work, we choose δ to be equal to the mean sampling of the light curves (δ = 5.2 days) and s = 100 days, yielding a point estimate of Δt_AB = −118.6 days for the time delay. The corresponding d_i and f_i difference light curves are shown in the top panel of Fig. 5. Results of the uncertainty analysis will be presented in the next section, together with the performance of other point estimators.

We have explored a range of alternative values for the free parameters (s = 50,100,150,200 and δ = 2.6,5.2,10.4 days), and find that neither the time-delay point estimate from the observed data, nor the error analysis is significantly affected. The time-delay estimates resulting from these experiments stay within 1.2 days around the reference value obtained for δ = 5.2 and s = 100 days. Regarding the uncertainty analysis, we observe that increasing the smoothing length scale s beyond 100 days decreases the random error, but at the cost of an increasing bias, which is not surprising.

4. Time-delay measurement of SDSS J1001+5027

Fig. 6 Error analysis of the four time-delay measurement techniques, based on delay estimations on 1000 synthetic curves that mimic our SDSS J1001+5027 data. The horizontal axis corresponds to the value of the true time delay used in these synthetic light curves. The gray vertical lines delimit bins of true time delay. In each of these bins, the colored rods and 1σ error bars show the systematic biases and random errors, respectively, as committed by the different techniques.

In this work, we use five different methods to measure the time delay of SDSS J1001+5027 from the data shown in Fig. 4. All these methods have been developed to address light curves affected by extrinsic variability, resulting from microlensing or flux contamination. Three of the techniques, namely the dispersion-like technique, the regression difference technique, and the free-knot spline technique are described in length in Tewes et al. (2013a) and were used to measure the time delays in the four-image quasar RX J1131−123 (Tewes et al. 2013b).

In the the previous section, we presented our fourth method, the difference-smoothing technique. These first four methods are point estimators: they provide best estimates, without information on the uncertainty of their results. We proceed by quantifying the accuracy and precision of these estimators by applying them to a set of 1000 fully synthetic light curves, produced and adjusted following Tewes et al. (2013a). These simulations include the intrinsic variations of the quasar source, mimicking the observed variability of SDSS J1001+5027, as well as extrinsic variability on a range of time scales from a few days to several years. They share the same sampling and scatter properties as the real observations.

Figure 6 shows the results of this analysis, depicting the delay measurement error as a function of the true delay used to generate the synthetic light curves. As always, this analysis naturally takes into account the intrinsic variances of the techniques, that are due to the limited ability of the employed global optimizers to find the absolute minima of the cost functions.

As can be seen in Fig. 6, the dispersion-like technique is strongly biased for this particular dataset. This could be a consequence of the simplistic polynomial correction for extrinsic variability linked to this technique. For the other techniques, the bias remains smaller than the random error, and no strong dependence on the true time delay is detected.

The final systematic error bar for each of these four techniques is taken as the worst measured systematic error on the simulated light curves (biggest colored rod in Fig. 6). The final random error is taken as the largest random error across the range of tested time delays. Finally, the total error bar for each technique is obtained by summing the systematic and random components in quadrature.

Fig. 7 Time-delay measurements of SDSS J1001+5027, following five different methods. The total error bar shown here includes systematic and random errors.

In the writing process of this paper, Hojjati et al. (2013) proposed a new independent method to measure time delays that is also able to address extrinsic variability. Their method is based on Gaussian process modeling, and does not rely on point estimation. It provides its own standalone estimate of the total uncertainty. We have provided these authors with the COSMOGRAIL data of SDSS J1001+5027, without letting them know our measured values. They find Δt_AB = −117.8 ± 3.2 days.

We include this measurement by Hojjati et al. (2013) as a fifth measurement in our result summary, presented in Table 2 and in a more graphical form in Fig. 7. Not only do their time-delay values agree with our four estimates, but also their error bars agree well with ours, in spite of the totally different way of estimating them.

We have five time-delay estimates from five very different methods, and all these estimates are compatible with each other. We now need to combine these results. In doing this, we exclude the delay from the dispersion-like technique that, as we show, is dominated by systematic errors. While the estimates from the four remaining techniques are obtained with very different methods, they are still not independent, as they all make use of the same data. We therefore simply average the four time-delay measurements to obtain our combined estimate, and we use the average of the total uncertainties as the corresponding uncertainty. This leads to Δt_AB = −119.3 ± 3.3 days, shown in black in Fig. 7.

5. Conclusion

In this paper, we present the full COSMOGRAIL light curves for the two images of the gravitationally lensed quasar SDSS J1001+5027. The final data, all taken in the R band, totalize 443 observing epochs, with a mean temporal sampling of 3.8 days, from the end of 2004 to mid-2011. The COSMOGRAIL monitoring campaign for SDSS J1001+5027 is no longer in progress. It involved three different telescopes with diameters from 1.2 m to 2 m, hence illustrating the effectiveness of small telescopes in conducting long-term projects with potentially high impact on cosmology.

We analyzed our light curves with five different numerical techniques, including the three methods described in Tewes et al. (2013a). In addition, we introduced and described a new additional method, based on representing the extrinsic variability by a smoothed version of the difference light curve between the quasar images. Finally, we also presented results obtained via the technique of Hojjati et al. (2013), based on modeling of the quasar and microlensing variations using Gaussian processes. The technique was blindly applied to the data by the authors of Hojjati et al. (2013), without any prior knowledge of the results obtained with the other four methods.

Aside from the dispersion-like technique, dominated by systematic errors, we find that the four other methods yield similar time-delay values and similar random and systematic error bars. Our final estimate of the time delay is taken as the mean of these four best results, together with the mean of their uncertainties: Δt_AB = −119.3 ± 3.3 days, with image A leading image B. This is a relative uncertainty of 2.8%, including systematic errors.

The present time-delay measurement can be used in combination with lens models to constrain cosmological parameters, in particular the Hubble parameter H₀ and the curvature Ω_k (e.g., Suyu et al. 2013a). The accuracy reached on cosmology with SDSS J1001+5027 alone or in combination with other lenses, will rely on the availability of follow-up observations to measure: (1) the lens velocity dispersion; (2) the mass contribution of intervening objects along the line of sight; and (3) the detailed structure of the lensed host galaxy of the quasar. This translates in practice into one single night of an 8 m-class telescope, plus about four orbits of the Hubble Space Telescope.

Acknowledgments

We thank the numerous observers who contributed to the data from the Mercator and Maidanak telescopes, and we are grateful for the support provided by the staff at the Indian Astronomical Observatory, Hanle and CREST, Hoskote. We also thank A. Hojjati, A. Kim, and E. Linder for running their curve shifting algorithm on our data. S. Rathna Kumar and C. S. Stalin acknowledge support from the Indo-Swiss Personnel Exchange Programme INT/SWISS/ISJRP/PEP/P-01/2012. COSMOGRAIL is financially supported by the Swiss National Science Foundation (SNSF). We thank the referee Peter Schneider for his timely comments that helped to improve this paper.

References

All Tables

All Figures

	Fig. 1 Distribution of the average observed FWHM and elongation ϵ of field stars in the images used to build the light curves of SDSS J1001+5027.
In the text

Fig. 2 R-band image centered on SDSS J1001+5027. The image is the combination of the 210 best exposures from the Mercator telescope, for a total exposure time of 21 h. We use the stars labeled 1, 2, and 3 to model the PSF and to cross-calibrate the photometry of each exposure. The position of the two lensing galaxies G1 and G2 are indicated in the zoomed image in the upper left. They are most clearly seen in the deconvolved images presented in Fig. 3.

In the text

Fig. 3 Two ways of modeling the light distribution for extended objects during the deconvolution process. On the left is shown a single 360-s exposure of SDSS J1001+5027 obtained with the Mercator telescope in typical atmospheric conditions. The other panels show the parametric (top row) and pixelized light models (bottom row) for the lens galaxies as described in the text. The residual image for the single exposure is also shown in each case, as well as the average residuals over the 120 best exposures. The residual maps are normalized by the shot noise amplitude. The dark areas indicate excess flux in the data with respect to the model. Gray scales are linear.

In the text

	Fig. 4 R-band light curves of the quasars images A and B in SDSS J1001+5027 from March 2005 to July 2011. The 1σ photometric error bars are also shown. For display purpose, the curve of quasar image B is shown shifted in time by the measured time delay (see text). The light curves are available in tabular form from the CDS and the COSMOGRAIL website.
In the text

Fig. 5 Difference light curves of SDSS J1001+5027 as obtained by the new difference-smoothing technique introduced in this paper. The curves are shown for the best time-delay estimate found with this technique (top panel, ΔtAB = −118.6 days), and for a wrong time-delay value (bottom panel, ΔtAB = −100.0 days). The difference light curves di are shown as colored points. They are smoothed using a kernel of width s = 100 days to compute the fi (black points). The error bars on the black points show the uncertainty coefficients σfi. The points in the difference light curve di are color-coded according to the absolute factors of their uncertainties σdi by which they deviate from fi. In both panels, light curve A is used as reference, and light curve B is shifted in flux by the same amount.

In the text

Fig. 6 Error analysis of the four time-delay measurement techniques, based on delay estimations on 1000 synthetic curves that mimic our SDSS J1001+5027 data. The horizontal axis corresponds to the value of the true time delay used in these synthetic light curves. The gray vertical lines delimit bins of true time delay. In each of these bins, the colored rods and 1σ error bars show the systematic biases and random errors, respectively, as committed by the different techniques.

In the text

	Fig. 7 Time-delay measurements of SDSS J1001+5027, following five different methods. The total error bar shown here includes systematic and random errors.
In the text