Issue 
A&A
Volume 556, August 2013



Article Number  A22  
Number of page(s)  10  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201220352  
Published online  19 July 2013 
COSMOGRAIL: the COSmological MOnitoring of GRAvItational Lenses ^{⋆,}^{⋆⋆}
XIII. Time delays and 9yr optical monitoring of the lensed quasar RX J1131−1231
^{1} Laboratoire d’astrophysique, École Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland
email: malte.tewes@epfl.ch
^{2} Department of Astronomy, The Ohio State University, 140 West 18th Av., Columbus, OH 43210, USA
^{3} Center for Cosmology and Astroparticle Physics, The Ohio State University, 191 West Woodruff Av., Columbus, OH 43210, USA
^{4} Institut d’Astrophysique et de Géophysique, Université de Liège, Allée du 6 Août, 17, 4000 Sart Tilman, Liège 1, Belgium
^{5} Instituut voor Sterrenkunde, Katholieke Universiteit Leuven, Celestijnenlaan 200B, 3001 Heverlee, Belgium
^{6} ArgelanderInstitut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
^{7} School of Physics and Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK
Received: 8 September 2012
Accepted: 31 May 2013
We present the results from nine years of optically monitoring the gravitationally lensed z_{QSO} = 0.658 quasar RX J1131−1231. The Rband light curves of the four individual images of the quasar were obtained using deconvolution photometry for a total of 707 epochs. Several sharp quasar variability features strongly constrain the time delays between the quasar images. Using three different numerical techniques, we measured these delays for all possible pairs of quasar images while always processing the four light curves simultaneously. For all three methods, the delays between the three close images A, B, and C are compatible with being 0, while we measured the delay of image D to be 91 days, with a fractional uncertainty of 1.5% (1σ), including systematic errors. Our analysis of random and systematic errors accounts in a realistic way for the observed quasar variability, fluctuating microlensing magnification over a broad range of temporal scales, noise properties, and seasonal gaps. Finally, we find that our timedelay measurement methods yield compatible results when applied to subsets of the data.
Key words: gravitational lensing: strong / quasars: individual: RX J11311231 / cosmological parameters
Based on observations made with the 1.2m Swiss Euler telescope (La Silla, Chile), the 1.3m SMARTS telescope (Las Campanas, Chile), and the 1.2m Mercator Telescope. Mercator is operated on the island of La Palma by the Flemish Community, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.
Light curves are available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/556/A22
© ESO, 2013
1. Introduction
Using the time delays between multiple images of gravitationally lensed sources to measure cosmological distances (Refsdal 1964) has several advantages: there is no need for any primary or secondary calibrator, and there are no effects from the intergalactic or interstellar medium. The method, originally proposed for gravitationally lensed supernovae, has so far exclusively been applied to quasars lensed in most cases by individual massive galaxies. Exceptions are SDSS J1004+4112 and J1029+2623, two quasars lensed by galaxy clusters, with long time delays (Fohlmeister et al. 2008, 2013). The quasar lens timedelay method is now recognized as a tool that complements other cosmological probes, in particular for constraining H_{0} as well as the darkenergy equationofstate parameter, w (e.g., Suyu et al. 2012; Linder 2011; Moustakas et al. 2009). In spite of its advantages, the method has long faced two severe limitations to its effectiveness in constraining cosmology.
Overview of our optical monitoring of RX J1131−1231.
First, time delays between the gravitationally lensed images of a quasar are hard to measure. Some claimed timedelay measurements turned out to be erroneous (see, for instance, the controversy around Q0957+561: Vanderriest et al. 1989; Press et al. 1992; Schild & Thomson 1997; Kundic et al. 1997, and references therein). Understandably, early light curves tended to be short and sparse, often too short to clearly demonstrate that microlensing variability was not interfering with their analysis. Microlensing is seen as an uncorrelated extrinsic variability in the quasar images, which results from the timevariable magnification created by stars in the lensing galaxy (e.g., Chang & Refsdal 1979; Schmidt & Wambsganss 2010). In the best cases, light curves spanned a few years (see, e.g., Wyrzykowski et al. 2003; Hjorth et al. 2002; Burud et al. 2002a,b). One consequence is that the numerical methods used to measure time delays from these light curves were exceedingly “optimistic” in their assumptions about extrinsic variability. More recent measurements with better data frequently yielded delays inconsistent with the error estimates of the earlier measurements.
Second, given the measured time delays and lens and quasar image astrometry, there is a famed degeneracy between the timedelay distance, which is a scale parameter inversely proportional to H_{0}, and the spatial distribution of mass responsible for the stronglensing phenomenon. The delays constrain only a combination of this timedelay distance and the surface density of the lens near the images (Kochanek 2002). This can be overcome with independent constraints on the structure of the lens. Suyu et al. (2009, 2010) convincingly showed that it is possible to control the effects of model degeneracies for B1608+656 (Myers et al. 1995, CLASS survey), a quadruply imaged quasar with accurate radio time delays (Fassnacht et al. 2002). To do this, the authors combined (1) detailed HST images of the lensed quasar host galaxy; (2) a velocitydispersion measurement of the lens galaxy; and (3) information about the contribution of intervening galaxies along the line of sight, from galaxy number counts calibrated with numerical simulations.
In parallel to the advances in lens modeling, the observational situation has impressively evolved as well. Two observational groups, the COSmological MOnitoring of GRAvItational Lenses (COSMOGRAIL) and Kochanek et al. (2006), have been intensely monitoring ≈20 lenses for roughly ten years. In 2010, our two groups decided to merge their observational efforts, with the COSMOGRAIL group focusing on the analysis of time delays and the Kochanek et al. (2006) group focusing on the analysis of microlensing. While preliminary results have been published both before and after this merger (Kochanek et al. 2006; Vuissoz et al. 2007, 2008; Morgan et al. 2008a,b), exquisite data spanning almost a decade of continuous observation are now being released, for instance for the quadruply imaged quasar HE 04351223 (Courbin et al. 2011; Blackburne et al. 2011).
In this paper, we present nine years of optical monitoring of the quadruply imaged quasar RX J1131−1231 (Sluse et al. 2003), and measure its time delays with the techniques of Tewes et al. (2013). RX J1131−1231 is one of the most spectacular lenses of our sample. The redshift of the lensing galaxy is z_{lens} = 0.295, while the quasar is at z_{QSO} = 0.658. This low quasar redshift means (1) that the photometric variations are fast, numerous, and strong because it is a lowerluminosity quasar (see, e.g., MacLeod et al. 2010), and (2) that the host galaxy of the lensed quasar is seen as a full Einstein ring with many spatially resolved structures in HST images. Similarly, the lensing galaxy is sufficiently bright to allow a precise measurement of its velocity dispersion and possibly of its velocity dispersion profile. These characteristics facilitate both the timedelay measurement and the lens modeling. The latter, with stateoftheart inferences of cosmological constraints based on our timedelay measurements of RX J1131−1231, are presented in Suyu et al. (2013).
Fig. 1
Distributions of stellar FWHM and elongation ϵ of all images contributing to the light curves, by telescope. 
Fig. 2
Part of the field of view of the Swiss 1.2m Euler telescope around RX J1131−1231. The widefield image is a combination of 600 exposures of 360s each, corresponding to a total exposure time of 2.5 days. The inset shows a single 360s exposure from the EulerCAM instrument under good conditions (FWHM ~ 1′′, ϵ = 0.1, no Moon). All images are taken in the R band. We mainly use the stars labeled 1 to 4 to build a PSF model for each exposure. Stars 1 to 5 are used for the photometric calibration. 
Our observations of RX J1131−1231 and their reduction are described in Sects. 2 and 3, while the light curves are presented in Sect. 4. In Sect. 5, we apply three different curveshifting techniques to the light curves and infer our best measurements of the delays along with realistic random and systematic error bars. Our results are summarized in Sect. 6.
2. Observations
We have been monitoring the quadruply lensed quasar RX J1131−1231 (J2000: 11^{h}31^{m}52^{s}, − 12°31′59″) since December 2003 with three different telescopes in the R band (~600–720 nm). Table 1 summarizes the observational strategy and instrumental characteristics. The light curves presented in this paper cover nine observing seasons (2004–2012, see Fig. 4), with 707 monitoring epochs in total. The average sampling within the seasons is 2.9 days, and the median sampling is 2.0 days. The mean seasonal gap in the combined light curves is 132 days with a standard deviation of two weeks.
The majority of the measurements for this southern target came from the Swiss 1.2m Euler telescope and the Small & Moderate Aperture Research Telescope System (SMARTS) 1.3m telescope, both located in Chile under equally good atmospheric conditions. Figure 1 shows the distributions of the stellar full width at half maximum (FWHM) and the elongation ϵ of each observing epoch, as measured by SExtractor (Bertin & Arnouts 2010) using field stars. The SMARTS 1.3m telescope is guided, in contrast to Euler and Mercator, which are solely tracking. This accounts in part for the broad elongation distribution of the Euler and Mercator images.
The original imaging instrument C2 of the Euler telescope was replaced in September 2010 by EulerCAM, a liquidnitrogen cooled 4k × 4k e2v 23184 CCD yielding a pixel scale of 0215. At the same time the focusing procedure was improved. Figure 1 shows that after two years of observations, images from EulerCAM tend to be statistically sharper than C2 images, and the FWHM is more similar to the SMARTS 1.3m data. All images from the SMARTS 1.3m telescope were obtained through the optical channel of the ANDICAM^{1} camera (DePoy et al. 2003). See Kochanek et al. (2006) for details about the SMARTS data.
Figure 2 shows the central part of a deep combination of Euler monitoring images, totaling 2.5 days of exposure. In the best individual exposures the four quasar images are well resolved; images A and B have the smallest separation of 12.
3. Image reduction
In this section we describe the reduction procedure leading to the light curves of the multiple quasar images of RX J1131−1231. This same procedure will homogeneously be applied to other lenses monitored by COSMOGRAIL and SMARTS, hence we provide a general description of our methods.
All exposures are corrected for bias/readout effects and flat fielded using standard methods. This prereduction is performed individually for each telescope and instrument. For Euler and Mercator we have developed custom python pipelines that enable efficient human inspection and validation of each step. Particular attention is given to a semiautomated selection and combination of skyflats, which is controlled through a simple graphical interface highlighting the temporal evolution of instrumental contaminations such as dust. Using difference images, we check that all sources in the flatfield exposures are adequately masked and that the contamination does not evolve significantly within sets of flatfields that are stacked as masterflats. The period over which we stack flatfields can be from single days to several weeks. Some exposures of the Mercator telescope are prereduced using superflats, combinations of masked science exposures obtained in a single night. We have developed variants of these pipelines adapted to the particularities of other COSMOGRAIL telescopes.
From here on, all images are processed by a single deconvolution photometry software package. The principle ideas for this reduction procedure are the same as employed for previously published COSMOGRAIL light curves (Vuissoz et al. 2007, 2008; Courbin et al. 2011). Over the years we have reworked the steps, and implemented the procedure in the form of a python pipeline linked to a relational database, containing one entry per exposure.
3.1. PSF construction
We first build a conservatively smooth sky model, obtained using SExtractor, for each exposure. This sky model is not critical, as the subsequent photometry procedure will fit for any residual sky level around all sources of interest.
We then align the images, separately for each instrument, and individually estimate the point spread function (PSF) of each exposure. This is done by fitting several field stars using a common model, composed of (i) a simply parametrized profile and (ii) a regularized fine pixel array. The details of this procedure, which is part of a generalpurpose deconvolution package based on ideas reported in Magain et al. (1998), are described in Cantale et al. (in prep.). For most of the exposures of RX J1131−1231, we used the stars labeled 1–4 in Fig. 2 to build the PSF model. The pipeline allows us to easily explore and compare different choices of PSFstars. The final choice of stars is empirically selected to yield the least scatter in the final quasar light curves. For RX J1131−1231 the situation is close to optimal, since (1) the stars 1–4 are bright but still in the linear regime of the CCDs; (2) they surround the lens at modest angular distances; and (3) all stellar companions or background objects can be identified and masked, yielding clean fitting residuals. Cosmicray hits and CCD artifacts are automatically masked using a variant of the L.A. Cosmic algorithm (van Dokkum 2001). We visually supervise the PSF construction and manually adapt the selection of PSFstars for problematic frames. The construction of the PSF models is the predominant technical issue affecting the final quality of the light curves.
3.2. Photometric normalization of the exposures
We next compute a multiplicative normalization coefficient for each exposure, based on the photometry of field stars. These coefficients will constitute the reference for the differential photometry of the quasar images.
In a first step, this computation is made separately for each instrument. We measure the instrumental fluxes N_{⋆ ij} (in photons) of star i in each exposure j, by fitting the PSF_{j} of this exposure to each star, leaving only the fluxes, astrometric shifts, and background levels of each star as free parameters. For each star and exposure, we then compute an individual calibration coefficient (1)where med_{j} denotes the median over all exposures of this instrument. Then, for each exposure j, we obtain the normalization coefficient through (2)The advantages of this simple and fast method is that it does not select a single exposure as the “reference” for all stars. For each instrument, the distribution of these coefficients c_{j} is highly unimodal, because the exposure time is kept constant.
To select the normalization stars we first inspect the normalized stellar light curves for any anomalies or variability on scales of months or years, and, if required, we iteratively repeat the coefficient computation with a refined selection of stars. In terms of the smoothness of the final quasar light curves, the best results are generally obtained by computing these coefficients from a few wellexposed stars with a similar and high signaltonoise ratio, as opposed to using larger numbers of noisier stars. Furthermore, the ideal normalization stars should be of similar color to the quasar. This minimizes potential differential effects due to varying atmospheric conditions given the relatively broad Rband filters used for the monitoring.
In a second step, we rescale the coefficients of each instrument so that a star whose color is closest to that of the quasar receives the same normalized median flux across all instruments. Residual adjustments to this normalization between instruments are performed later, using the quasar light curves in temporal regions where the data from different telescopes overlap.
3.3. Photometry of the quasar images
We obtain deblended light curves of the quasar images by MCS deconvolution photometry, following Magain et al. (1998). This algorithm fits a single model consisting of (1) some number of point sources and (2) a fine pixel array (hereafter the pixel channel) simultaneously to all exposures. For RX J1131−1231, four point sources model the quasar images, while the pixel channel represents both the lensing galaxy and the lensed host galaxy (Einstein ring). This model is scaled by the normalization coefficients from Eq. (2) and convolved by the exposurespecific PSF, before it is fit to the data. By construction, the relative astrometry of the four point sources and the pixel channel are common to all exposures. Only the fluxes of the point sources, the absolute astrometric shift, and a flat residual sky level around the lens are let free to vary between exposures. Iteratively, all these parameters are optimized together with the structure of the regularized pixel channel to minimize a single global χ^{2}. This algorithm was previously successfully applied in the discovery paper of RX J1131−1231 (Sluse et al. 2003), as well as in past COSMOGRAIL publications (e.g., Eigenbrod et al. 2007; Vuissoz et al. 2008; Courbin et al. 2011; Sluse et al. 2012) or for similar monitoring data (Burud et al. 2002a; Hjorth et al. 2002; Jakobsson et al. 2005; Morgan et al. 2012). The light curves of the quasar images are a direct output of this procedure.
Errors in the astrometry or the structure of the pixel channel might degrade or bias the photometry of the quasar images. We observe that the influence of the astrometry on the light curves is weak. Even for bright quasar images, alterations of the relative position of the point sources by as much as 005 do not significantly modify the light curves. Larger errors tend first to smoothly bias the magnitude measurements, before introducing additional scatter. To obtain light curves for gravitational lenses with image separations on the order of the resolution of the best monitoring data, we find no difference between using the tight astrometric constraints from HST images (Morgan et al. 2006; Chantry et al. 2010) or letting the deconvolution algorithm freely optimize the astrometry of the model.
We reach similar conclusions for the impact of the pixel channel and its mandatory regularization. Simple numerical experiments show that to first order, the effect of different plausible solutions is similar to very small additive shifts to the flux of the quasar light curves, with only marginally perceptible effects even for faint quasar images. In practice, we systematically constrain this pixel channel as well as the point source astrometry using only a subset of images with the best resolution.
3.4. Photometric error estimation
In this section we describe how we compute a rather formal best case error estimate for the photometry of each quasar image in each exposure. The normalized flux of a quasar image in a given exposure can be written f_{⋆} = N_{⋆}·c, where N_{⋆} is the measured number of photons, and c is the normalization coefficient from Eq. (2). We assume that the random error on f_{⋆} has two independent sources: (i) the shot noise σ_{N⋆} of the quasar image itself; and (ii) the noise σ_{c} of the normalization coefficient as computed in Sect. 3.2.
We compute σ_{N⋆} following the standard “CCD equation” (see e.g. Howell 2006, Chap. 4.4) (3)where S is the sky level, R is the CCD read noise (both in photons per pixel), and n_{pix} is the number of (equivalent) pixels of the software aperture. The dark current is negligible in our images. The MCS deconvolution of point sources corresponds to PSF fitting, and following Heyer & Biretta (2004, Sect. 6.5.1), we use for n_{pix} the reciprocal of the PSF sharpness, (4)where PSF_{lm} is the fraction of light in the total PSF at pixel lm. In this simple best case computation, we assume that the PSF is perfectly known, and that the sky level is well determined by the fitting procedure. We do not take into account the surface brightness of the lens and host galaxies, which is generally negligible with respect to the sky level. Furthermore, we do not compute noise terms that are due to the blending with other point sources.
As a sanity check, we tested our deconvolution photometry algorithm on simulated point source images created with SkyMaker (Bertin 2009) for various realistic signaltonoise ratios, elongations, and seeing conditions. Using two or more wellexposed stars to build the PSF, the statistical scatter obtained by our photometry pipeline is only marginally larger than the bestcase shot noise computed through the above procedure. Using similar synthetic images, we also verified that deconvolving two partially blended sources of different fluxes does not significantly bias their measured flux difference in any direction.
To estimate the random uncertainty of the normalization coefficient, we use the standard error of the mean (SEM) of the individual coefficients obtained from the n normalization stars in the exposure (see Eq. (1)) (5)Finally, we obtain the error estimated σ_{f⋆ CCD} on the normalized flux through (6)Note that for our COSMOGRAIL and SMARTS monitoring data, the shot noise term from the CCD equation dominates this error budget. The contribution from σ_{c} is typically ≪3% of σ_{f⋆ CCD}. Such small calibration errors are important only when measuring very short delays, where they would introduce positively correlated noise into the curves.
3.5. Combination of points per night
In each monitoring night, we observe the lenses m times (m = 3 or 5) over ≈30 min, yielding flux measurements f_{⋆ k}, with k = 1, .. ,m for each quasar image. To reduce the CPU cost and to reject outliers, we bin these measurements by epochs, separately for each instrument and telescope, before measuring the time delays. We attribute to each epoch the medians of the photometric measurements f_{⋆ k} within that night, and the mean of the Heliocentric Julian Dates (HJD).
This approach also allows us to obtain an empirical photometric error estimate for each quasar image, using a measure of the spread of the image’s f_{⋆ k}. For increased robustness against outliers, we quantify this spread using the median absolute deviation from the median, or MAD (Hoaglin et al. 1983), (7)To estimate the usual σ of an (assumed) Gaussian distribution, the MAD is rescaled as (8)These procedures give us two different error estimates for each epoch and quasar image: (1) the median of the errors estimated individually for each exposure med(σ_{f⋆ CCD}); and (2) the more empirical σ_{f⋆ MAD}. We observe, as expected, that both error estimations are highly correlated, but also that the empirical σ_{f⋆ MAD} is typically twice as high.
The uncertainty we finally assign to each epoch and quasar image uses the higher of med(σ_{f⋆ CCD}) and σ_{f⋆ MAD}. Indeed, a realistic error estimate cannot be lower than either of them. We divide this estimate by the square root of the number of exposures m of the epoch, so that the final error estimate for the median photometric measurement on a given instrument and night is (9)Figure 3 displays the distribution of the resulting error estimates as a function of Rband magnitude and instrument. The data points from all four quasar images of RX J1131−1231 are shown, yielding a broad spread in magnitudes.
Fig. 3
Photometric 1σ error estimates of each observing epoch as a function of approximate Rband magnitude. These errors are estimated following Eq. (9). 
Fig. 4
Optical monitoring of RX J1131−1231, as obtained from deconvolution photometry. From top to bottom are shown the Rband light curves for the quasar images A, B, C, and D along with the 1σ photometric error bars. Colors encode the contributing instruments. Curves B and D have been shifted by +0.3 mag for display purposes. The light curves are available in tabular form from the CDS and the COSMOGRAIL website. 
3.6. Combination of telescopes
The lens galaxy and the lensed images of the quasar host galaxy differ in color from the quasar itself. As a consequence, the different R filters and CCD response functions of the monitoring instruments might “see” these contaminating sources with slightly different amplitudes even with a perfect calibration of the quasar fluxes. This can result in small mismatches between the light curves obtained using different combinations of telescopes and instruments.
We correct our light curves for such small effects (occasionally as high as 10% of the flux of the faintest quasar images) in this final step. We typically select the instrument with the longest or the highestquality curve as a reference. For RX J1131−1231, this is the SMARTS light curve. For each of the other instruments, we optimize additive magnitude and flux shifts (i.e., multiplicative and additive flux corrections) to minimize a dispersion measure between each instrument’s light curve and the reference light curve. We compute this dispersion following the curveshifting technique presented in Tewes et al. (2013), but evaluate it between the light curves of different instruments, instead of different quasar images. Provided that the colors of the quasar images are not differentially reddened by absorption in the lens galaxy, we optimize a single common magnitude shift per instrument, and individual flux shifts for each quasar image and instrument. This is adequate for RX J1131−1231, as previously suggested by Sluse et al. (2007) and Chartas et al. (2009, 2012).
4. Light curves of RX J1131−1231
All these steps were applied to the COSMOGRAIL and SMARTS data for RX J1131−1231. Figure 4 shows the final nineyearlong light curves of the four quasar images. These data are available in machinereadable form at the CDS and on the COSMOGRAIL website^{2}. The light curves are dominated by intrinsic quasar variability, with some features on scales as short as a few weeks. It can be readily seen in the 2008 season for instance, that the delays between A, B, and C must be very short, while D is delayed by slightly less than 100 days. Intriguingly, looking only at the first season of A, B, and C, one might think that A is significantly delayed with respect to B and C. We attribute this discrepancy to microlensing variability, which manifestly changes the magnitude difference between the A and B images from the first to the second season. We discuss this “event” in more detail in Sect. 5.3. Prominent microlensing variability on long timescales is ubiquitous because the flux ratios between the quasar images evolve by as much as a magnitude. These microlensing effects in RX J1131−1231 have been analyzed in Morgan et al. (2010) and Dai et al. (2010).
Lastly, we observe that the photometric error estimates, obtained from Eq. (9), match the observed scatter well in the smooth sections of curves from the individual telescopes. They are certainly not conservatively large, but we stress that the scale of these error estimates has no direct influence on the uncertainties that we compute for the timedelay measurements in the next Section. Our results are robust against a deliberate increase of these error estimates by up to a factor of 5.
5. Timedelay estimation
In this section we infer the time delays of RX J1131−1231 from the light curves shown in Fig. 4, closely following the curveshifting and uncertainty evaluation procedures described in Tewes et al. (2013). We summarize the principal ideas below. A major difficulty and potential source of bias for curveshifting methods is the presence of extrinsic variability in the light curves, on top of the intrinsic quasar variability common to all four images. The main source of the extrinsic signal is variable microlensing magnification due to the motions of the stars in the lens galaxy. As shown by Mosquera & Kochanek (2011), microlensing can affect light curves over a broad range of time scales. For RX J1131−1231, Mosquera & Kochanek (2011) estimated a time scale of ≈11 years for the crossing of a stellar Einstein radius, and ≈3 months for the source radius to cross a caustic. Other potential sources of extrinsic variability, such as variable quasar structure and spurious additive flux contaminating the photometry, are summarized in Tewes et al. (2013). In the present paper we do not aim to separate or analyze the extrinsic variability in physical terms. We simply consider it as a hindrance for the timedelay measurement.
The curveshifting methods of Tewes et al. (2013) try to minimize the bias due to such extrinsic variability in different ways. They all rely on iterative optimizations of time shifts of the four light curves, and yield selfconsistent point estimates of the time delays between all six image pairs :

1.
The freeknot spline technique simultaneously fits onecommon intrinsic spline for the quasar variability, andindependent, smoother extrinsic splines for the microlensing tothe light curves. The curves are shifted in time to optimize this fit.This method is similar to the polynomial method of Kochaneket al. (2006).

2.
The regression difference technique shifts continuous regressions of the curves in time to minimize the variability of the differences between them. This novel method does not involve an explicit model for the microlensing variability.

3.
The dispersionlike technique, inspired by Pelt et al. (1996), shifts the curves to minimize a measure of the dispersion between the overlapping data points. This method has no explicit model for the common intrinsic variability of the quasar, but it includes polynomial models for the extrinsic variability.
Using several independent algorithms allows us to crosscheck our estimates for techniquedependent biases. Indeed, as important as the timedelay point estimates themselves is the reliable estimation of their uncertainties. For this we follow a Monte Carlo approach, by applying the curveshifting techniques to a large number of synthetic light curves with known time delays. These curves are drawn from a light curve model that mimics both the observed intrinsic and extrinsic variability of the real observations, randomizing only the unrecoverable short timescale extrinsic variability. This latter correlated noise locally adapts its amplitude to the scatter of the observed data points. As a result, the synthetic curves are virtually indistinguishable from the real ones (see Figs. 5 and 6 of Tewes et al. 2013, for an illustration).
Fig. 5
Analysis of the intrinsic variance of the timedelay point estimators, found by running the three curveshifting techniques 200 times on the light curves shown in Fig. 4, starting the optimizations from random initial time shifts. The widths of these distributions reflect the failure of the methods to converge to a single optimal solution given the data. We stress that these distributions do not represent probability density functions of the time delays. A small intrinsic variance does not imply that a timedelay estimation is precise or accurate, only that the error surface is relatively smooth. 
Fig. 6
Timedelay measurements along with the 1σ errors for RX J1131−1231 obtained by our three standard techniques from the full nineseasonlong light curves shown in Fig 4. The random and systematic error contributions are given in parentheses for each delay. The error bars represent the random plus systematic errors summed in quadrature. A positive ABdelay Δt_{AB} means that image B leads image A. We consider the measurements from the regression difference technique, which display the lowest bias and variance in our error analysis as the delays that should be used to constrain cosmology and lens models. 
Fig. 7
Results of the Monte Carlo analysis leading to the error estimates for our timedelay measurements shown in Fig. 6. We obtain our uncertainty estimates by applying the curveshifting techniques to 1000 synthetic light curve sets that closely mimic the observed data but have known true time delays. The vertical axes show the delay measurement error, which is compared with the true delays used to generate the synthetic curves (horizontal axes). Separately for each panel, the outcomes are binned according to the true time delays. The bin intervals are shown as light vertical lines. Within each bin, the shaded rods and error bars show the systematic and random errors, respectively, of the delay measurements for each technique. 
5.1. Application to RX J1131−1231
We begin by evaluating the robustness of the timeshift optimization of each method given the data. For this we repeatedly run the point estimators on the observed light curves. Each time, we start the optimizations from random initial conditions, using time shifts uniformly drawn ± 10 days around plausible solutions. The resulting distributions of point estimates are shown in Fig. 5, characterizing what we call the intrinsic variance of the estimates.
These tests are essential to check the reliability of the nonlinear optimization algorithms for a given data set and model parameters. We stress that this procedure should not be interpreted as importance sampling of any posteriors. Generally speaking, overly flexible models dilute the timedelay information and yield higher intrinsic variances. For the freeknot spline technique we chose an average knot step of 20 days for the quasar variability spline, and 150 days for the extrinsic splines. For the dispersionlike method we model the extrinsic variability by independent linear trends on each season. If the light curves sufficiently constrain the time delays, the freeknot spline and regression difference techniques can easily be adjusted to display small intrinsic variance and roughly unimodal distributions, as in Fig. 5. The dispersionlike technique is inherently more sensitive to initial conditions, owing to a much higher roughness of the scalar objective function that it minimizes. As discussed in Tewes et al. (2013), it is important to note that smoothing the objective function does not necessarily lead to more accurate timedelay estimates. For each technique, we use the mean values of the distributions in Fig. 5 as our best timedelay estimates. They correspond to the points in Fig. 6, which summarizes our results for RX J1131−1231.
The remaining part of the analysis is solely about estimating realistic error bars for these point estimates. For this we proceed by blindly running the exact same curveshifting techniques on 1000 sets of fully synthetic light curves with true time delays randomly chosen around the measured ones. Again, we start the methods from random initial shifts. Statistics of the resulting timedelay measurement errors (i.e., measurement − truth) are shown in Fig. 7 as a function of the true delays of the synthetic curves. In this figure, the shaded rods show the mean measurement error, which we call systematic error or bias, while the error bars indicate the standard deviation of the measurement errors, representing the random errors. The total errors for our delay measurements, shown by the error bars in Fig. 6, are computed by adding the maximum bias and the maximum random error of each panel in quadrature. Note that we run the curveshifting techniques only once on each random synthetic curve set, so this error analysis takes into account the intrinsic variance. This is rather conservative, given that we use the mean result of several timeshift optimizations as our best estimates for the observed data. However, the additional contribution to the uncertainty estimates is negligible for methods with low intrinsic variance. Before discussing these results, we show in Fig. 8 the same measurement errors as in Fig. 7, but plotted for each delay against each other to explore potential abnormal correlations.
Fig. 8
Correlations of delay measurement errors for synthetic light curves mimicking RX J1131−1231. The measurement errors are the same as shown in Fig. 7, but this time marginalizing over the true delays. Crosshairs indicate zero error, and the inset shows the scale of each panel. For clarity, only single contours at half of the maximum density are shown for two of the techniques. We observe no correlations (i.e., oblique contours) between the unrelated delay measurements along the short diagonal of this figure. 
5.2. Discussion
In terms of simple lensmodel considerations, RX J1131−1231 is a “longaxis cusp lens”, whose source is located inside a cusp of the tangential caustic curve on the long axis of the potential (Sluse et al. 2003). Images A and D are the saddle points of the arrivaltime surface, while B and C are minima (Blandford & Narayan 1986). Figure 14 of Saha et al. (2006) gives an illustration of this particular lens. Hence, the delays Δt_{AB} and Δt_{AC} are predicted to be small but positive, while Δt_{AD} is negative and large; the possible arrivaltime orders are BCAD or CBAD.
Our measured delays, as shown in Fig. 6, are consistent with these predictions, although the delays between A, B, and C are compatible with being zero given the 1σ errors. Keeton & Moustakas (2009) described and quantified the exciting idea that subhalos of a lens galaxy could be detected through anomalies in observed time delays with respect to predictions from simple lens models. For a cusp lens such as RX J1131−1231, these authors showed that the temporal ordering between the two minima (B and C) could easily be inverted by modest substructure. Our delay measurements do not hint at the presence of such timedelay millilensing in RX J1131−1231 because they are compatible with the smooth lens models of Suyu et al. (2013).
It is reassuring to observe that the three curveshifting methods yield consistent results for all delays despite the very different methodologies, which we see as a success of our error estimation procedure. Keep in mind, however, that the delay measurements are not independent, because they all use the same single set of observed light curves.
We now investigate in more detail the measurementerror analysis shown in Fig. 7, which is based solely on synthetic curves. Overall, we do not observe any strong dependence of a method’s bias or random error on the true time delays used in the simulations. For most image pairs and methods, the random errors are more important than the almost negligible bias. Exceptions are the two results from the dispersionlike technique, which is observed to overpredict Δt_{AB} and underpredict Δt_{BC} by about one day. Remarkably, this same technique also measured the highest (lowest) delay Δt_{AB} (Δt_{BC}) for the observed curves (Fig. 6), when compared to the other methods. We see this as another indication that the bias estimates are reliable. Finally, Fig. 8 shows no sign of abnormal correlations between measurement errors of quasar image pairs, whether they share a common quasar image or not. By construction, our timedelay uncertainty estimates for each image pair marginalize over all other pairs.
Which delay measurement method performs best on RX J1131−1231? Given that the regression difference technique yields both the lowest biases and the smallest random errors, we simply conclude that its measurements are the most informative ones. We will use the delays from the regression difference technique, expressed with respect to quasar image B, to measure the timedelay distance toward RX J1131−1231. Details and results of lens modeling, as applied to RX J1131−1231, are described in detail in Suyu et al. (2013).
5.3. Delay measurement on subsets of the light curves
Fig. 9
Application of the curveshifting techniques to subsections of the full light curves. The square diamonds (top) show measurements using only the first two seasons, the leftward triangles (middle) use the first five seasons, and the rightward triangles (bottom) use the last four seasons. All error bars depict total errors, as in Fig. 6. The blue and greenshaded regions correspond to the interval covered by the total error bars using the full nine seasons, by the spline and regression difference techniques, respectively. 
Our long light curves, featuring several delayconstraining intrinsic variability patterns, suggest that we should independently measure the time delays from subsets of the observing seasons. This analysis represents an invaluable check of the consistency of the timedelay estimation procedure of Tewes et al. (2013), and hence for the results from the full light curves.
We analyze three subsets of the available data: (1) the first two seasons; (2) the first five seasons; and (3) the remaining last four seasons. Case (1) is chosen for comparison with Morgan et al. (2006). We perform the analyses from scratch, without using any knowledge about the delay estimates from the full curves. In particular, we build a new model for the synthetic light curves for each case, independently adjusted so that they best resemble the observations in that time period based on the statistical criteria presented in Tewes et al. (2013, Sect. 7). We do not alter any parameters of the curveshifting techniques.
Figure 9 presents the resulting delay measurements. The data points and error bars depict the individual point estimates and the corresponding 1σ total errors obtained by each technique for the three cases. The shaded regions show the 1σ intervals from the full nine seasons taken from Fig. 6. We observe the following:

1.
Using only the first two seasons, our three methods are systematically biased toward high values of Δt_{AB} and Δt_{AC}. Our interpretation is that some microlensing variability in image A conspiratorially imitates a time shift, particularly around mid 2004. The estimates from the dispersionlike and spline technique roughly reproduce the results obtained from the same two seasons by Morgan et al. (2006) Δt_{AB} = 12.0 ± 1.5, Δt_{AC} = 9.6 ± 1.9 and Δt_{AD} = −87 ± 8, but with significantly wider, yet still too low, error bars. This demonstrates the need of long monitoring programs to measure accurate time delays and minimize the influence of unfortunately placed microlensing events.

2.
For these same two seasons, the regression difference method yields far better estimates for these delays – including adequately sized error bars that encompass the delays estimated from the full light curves. Note that the regression difference technique is the only one that does not involve an explicit model for the microlensing.

3.
For the other two divisions of the data, namely seasons 1−5 and 6−9, which are disjoint and hence independent, our methods yield consistent results. As expected, both subcases show higher error bars than the combined analysis of all seasons.
6. Conclusions
The first part of this paper describes the COSMOGRAIL data reduction procedure, which will be used to reduce all data gathered by our monitoring campaign of gravitationally lensed quasars. In the second part we apply this pipeline to our COSMOGRAIL and SMARTS observations of the quad lens RX J1131−1231, leading to an unprecedented set of nineyearlong light curves of high photometric quality. Several strong and fast intrinsic quasar variability patterns constrain the time delays between the multiple images. Microlensingrelated extrinsic variability is clearly present, as pointed out and analyzed in previous studies (Sluse et al. 2006, 2007; Morgan et al. 2010; Dai et al. 2010; Chartas et al. 2012). However, this distorting signal does not prevent us from measuring accurate time delays, using the three independent algorithms of Tewes et al. (2013).
The best timedelay estimates of RX J1131−1231 are provided by the regression difference technique. It measures the 91day delays between D and the other quasar images to a fractional 1σ uncertainty of 1.5%. This error estimate is obtained by applying the techniques to synthetic curves with known time delays, which contain extrinsic variability features similar to the observed ones. We demonstrate the consistency of our error estimates by independently measuring time delays – including error bars – from subsets of the observed light curves of RX J1131−1231. This experiment also reveals that long multiyear monitoring is essential for reliably measuring time delays, despite progress on the methods.
The results from this paper are used to constrain the timedelay distance toward RX J1131−1231 and deduce stringent implications for cosmology in Suyu et al. (2013). With our timedelay measurement errors of only 1.5%, the accuracy of this cosmological inference is actually limited by the residual uncertainty of (1) the gravitational potential of the lens galaxy and (2) the largescale structures along the line of sight to the quasar (BarKana 1996).
Acknowledgments
We are deeply grateful to the numerous observers who contributed to the data acquisition at the Swiss Euler telescope, the SMARTS 1.3m telescope, and the FlemishBelgian Mercator telescope. We wish to thank Sherry Suyu and Tommaso Treu for precious discussions and comments, and the anonymous referee for her/his useful suggestions. COSMOGRAIL is financially supported by the Swiss National Science Foundation (SNSF). We also acknowledge support from the International Space Science Institute in Bern, where some of this research has been discussed. C.S.K. and A.M.M. are supported by NSF grant AST1009756. E.E. was partially supported by ESA and the Belgian Federal Science Policy (BELSPO) in the framework of the PRODEX Experiment Arrangement C90312. DS acknowledges support from the Deutsche Forschungsgemeinschaft, reference SL172/11.
References
 BarKana, R. 1996, ApJ, 468, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Bertin, E. 2009, Mem. Soc. Astron. It., 80, 422 [NASA ADS] [Google Scholar]
 Bertin, E., & Arnouts, S. 2010, Astrophysics Source Code Library [Google Scholar]
 Blackburne, J. A., Kochanek, C. S., Chen, B., Dai, X., & Chartas, G. 2011 [arXiv:1112.0027] [Google Scholar]
 Blandford, R., & Narayan, R. 1986, ApJ, 310, 568 [NASA ADS] [CrossRef] [Google Scholar]
 Burud, I., Courbin, F., Magain, P., et al. 2002a, A&A, 383, 71 [Google Scholar]
 Burud, I., Hjorth, J., Courbin, F., et al. 2002b, A&A, 391, 481 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chang, K., & Refsdal, S. 1979, Nature, 282, 561 [Google Scholar]
 Chantry, V., Sluse, D., & Magain, P. 2010, A&A, 522, 95 [Google Scholar]
 Chartas, G., Kochanek, C. S., Dai, X., Poindexter, S., & Garmire, G. 2009, ApJ, 693, 174 [NASA ADS] [CrossRef] [Google Scholar]
 Chartas, G., Kochanek, C. S., Dai, X., et al. 2012, ApJ, 757, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Courbin, F., Chantry, V., Revaz, Y., et al. 2011, A&A, 536, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dai, X., Kochanek, C. S., Chartas, G., et al. 2010, ApJ, 709, 278 [NASA ADS] [CrossRef] [Google Scholar]
 DePoy, D. L., Atwood, B., Belville, S. R., et al. 2003, in Proc. SPIE, 827 [Google Scholar]
 Eigenbrod, A., Courbin, F., & Meylan, G. 2007, A&A, 465, 51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fassnacht, C. D., Xanthopoulos, E., Koopmans, L., & Rusin, D. 2002, ApJ, 581, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Fohlmeister, J., Kochanek, C. S., Falco, E. E., Morgan, C. W., & Wambsganss, J. 2008, ApJ, 676, 761 [NASA ADS] [CrossRef] [Google Scholar]
 Fohlmeister, J., Kochanek, C. S., Falco, E. E., et al. 2013, ApJ, 764, 186 [NASA ADS] [CrossRef] [Google Scholar]
 Heyer & Biretta 2004, WFPC2 Instrument Handbook, STScI, Baltimore [Google Scholar]
 Hjorth, J., Burud, I., Jaunsen, A. O., et al. 2002, ApJ, 572, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Hoaglin, D. C., Mosteller, F., & Tukey, J. W. 1983, Understanding Robust and Exploratory Data Analysis (Wiley) [Google Scholar]
 Howell, S. B. 2006, Handbook of CCD Astronomy, 2nd edn. (Cambridge University Press) [Google Scholar]
 Jakobsson, P., Hjorth, J., Burud, I., et al. 2005, A&A, 431, 103 [Google Scholar]
 Keeton, C. R., & Moustakas, L. A. 2009, ApJ, 699, 1720 [NASA ADS] [CrossRef] [Google Scholar]
 Kochanek, C. S. 2002, ApJ, 578, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Kochanek, C. S., Morgan, N. D., Falco, E. E., et al. 2006, ApJ, 640, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Kundic, T., Turner, E. L., Colley, W. N., et al. 1997, ApJ, 482, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Linder, E. V. 2011, Phys. Rev. D, 84, 123529 [NASA ADS] [CrossRef] [Google Scholar]
 MacLeod, C. L., Ivezic, Z., Kochanek, C. S., et al. 2010, ApJ, 721, 1014 [NASA ADS] [CrossRef] [Google Scholar]
 Magain, P., Courbin, F., & Sohy, S. 1998, ApJ, 494, 472 [NASA ADS] [CrossRef] [Google Scholar]
 Morgan, N. D., Kochanek, C. S., Falco, E. E., & Dai, X. 2006 [arXiv:astroph/0605321] [Google Scholar]
 Morgan, C. W., Eyler, M. E., Kochanek, C. S., et al. 2008a, ApJ, 676, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Morgan, C. W., Kochanek, C. S., Dai, X., Morgan, N. D., & Falco, E. E. 2008b, ApJ, 689, 755 [NASA ADS] [CrossRef] [Google Scholar]
 Morgan, C. W., Kochanek, C. S., Morgan, N. D., & Falco, E. E. 2010, ApJ, 712, 1129 [NASA ADS] [CrossRef] [Google Scholar]
 Morgan, C. W., Hainline, L. J., Chen, B., et al. 2012, ApJ, 756, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Mosquera, A. M., & Kochanek, C. S. 2011, ApJ, 738, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Moustakas, L. A., Abazajian, K., Benson, A., et al. 2009, Astro2010: The Astronomy and Astrophysics Decadal Survey Science White Papers [Google Scholar]
 Myers, S. T., Fassnacht, C. D., Djorgovski, S. G., et al. 1995, ApJ, 447, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Pelt, J., Kayser, R., Refsdal, S., & Schramm, T. 1996, A&A, 305, 97 [NASA ADS] [Google Scholar]
 Press, W. H., Rybicki, G. B., & Hewitt, J. N. 1992, ApJ, 385, 404 [NASA ADS] [CrossRef] [Google Scholar]
 Refsdal, S. 1964, MNRAS, 128, 307 [NASA ADS] [CrossRef] [Google Scholar]
 Saha, P., Courbin, F., Sluse, D., Dye, S., & Meylan, G. 2006, A&A, 450, 461 [Google Scholar]
 Schild, R., & Thomson, D. J. 1997, AJ, 113, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, R. W., & Wambsganss, J. 2010, General Relativity and Gravitation, 42, 2127 [NASA ADS] [CrossRef] [Google Scholar]
 Sluse, D., Surdej, J., Claeskens, J.F., et al. 2003, A&A, 406, L43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sluse, D., Claeskens, J.F., Altieri, B., et al. 2006, A&A, 449, 539 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sluse, D., Claeskens, J.F., Hutsemekers, D., & Surdej, J. 2007, A&A, 468, 885 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sluse, D., Chantry, V., Magain, P., Courbin, F., & Meylan, G. 2012, A&A, 538, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suyu, S. H., Marshall, P. J., Blandford, R. D., et al. 2009, ApJ, 691, 277 [Google Scholar]
 Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Treu, T., Blandford, R. D., et al. 2012 [arXiv:1202.4459] [Google Scholar]
 Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Tewes, M., Courbin, F., & Meylan, G. 2013, A&A, 553, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Dokkum, P. G. 2001, PASP, 113, 1420 [NASA ADS] [CrossRef] [Google Scholar]
 Vanderriest, C., Schneider, J., Herpe, G., et al. 1989, A&A, 215, 1 [NASA ADS] [Google Scholar]
 Vuissoz, C., Courbin, F., Sluse, D., et al. 2007, A&A, 464, 845 [Google Scholar]
 Vuissoz, C., Courbin, F., Sluse, D., et al. 2008, A&A, 488, 481 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wyrzykowski, L., Udalski, A., Schechter, P. L., et al. 2003, Acta Astron., 53, 229 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1
Distributions of stellar FWHM and elongation ϵ of all images contributing to the light curves, by telescope. 

In the text 
Fig. 2
Part of the field of view of the Swiss 1.2m Euler telescope around RX J1131−1231. The widefield image is a combination of 600 exposures of 360s each, corresponding to a total exposure time of 2.5 days. The inset shows a single 360s exposure from the EulerCAM instrument under good conditions (FWHM ~ 1′′, ϵ = 0.1, no Moon). All images are taken in the R band. We mainly use the stars labeled 1 to 4 to build a PSF model for each exposure. Stars 1 to 5 are used for the photometric calibration. 

In the text 
Fig. 3
Photometric 1σ error estimates of each observing epoch as a function of approximate Rband magnitude. These errors are estimated following Eq. (9). 

In the text 
Fig. 4
Optical monitoring of RX J1131−1231, as obtained from deconvolution photometry. From top to bottom are shown the Rband light curves for the quasar images A, B, C, and D along with the 1σ photometric error bars. Colors encode the contributing instruments. Curves B and D have been shifted by +0.3 mag for display purposes. The light curves are available in tabular form from the CDS and the COSMOGRAIL website. 

In the text 
Fig. 5
Analysis of the intrinsic variance of the timedelay point estimators, found by running the three curveshifting techniques 200 times on the light curves shown in Fig. 4, starting the optimizations from random initial time shifts. The widths of these distributions reflect the failure of the methods to converge to a single optimal solution given the data. We stress that these distributions do not represent probability density functions of the time delays. A small intrinsic variance does not imply that a timedelay estimation is precise or accurate, only that the error surface is relatively smooth. 

In the text 
Fig. 6
Timedelay measurements along with the 1σ errors for RX J1131−1231 obtained by our three standard techniques from the full nineseasonlong light curves shown in Fig 4. The random and systematic error contributions are given in parentheses for each delay. The error bars represent the random plus systematic errors summed in quadrature. A positive ABdelay Δt_{AB} means that image B leads image A. We consider the measurements from the regression difference technique, which display the lowest bias and variance in our error analysis as the delays that should be used to constrain cosmology and lens models. 

In the text 
Fig. 7
Results of the Monte Carlo analysis leading to the error estimates for our timedelay measurements shown in Fig. 6. We obtain our uncertainty estimates by applying the curveshifting techniques to 1000 synthetic light curve sets that closely mimic the observed data but have known true time delays. The vertical axes show the delay measurement error, which is compared with the true delays used to generate the synthetic curves (horizontal axes). Separately for each panel, the outcomes are binned according to the true time delays. The bin intervals are shown as light vertical lines. Within each bin, the shaded rods and error bars show the systematic and random errors, respectively, of the delay measurements for each technique. 

In the text 
Fig. 8
Correlations of delay measurement errors for synthetic light curves mimicking RX J1131−1231. The measurement errors are the same as shown in Fig. 7, but this time marginalizing over the true delays. Crosshairs indicate zero error, and the inset shows the scale of each panel. For clarity, only single contours at half of the maximum density are shown for two of the techniques. We observe no correlations (i.e., oblique contours) between the unrelated delay measurements along the short diagonal of this figure. 

In the text 
Fig. 9
Application of the curveshifting techniques to subsections of the full light curves. The square diamonds (top) show measurements using only the first two seasons, the leftward triangles (middle) use the first five seasons, and the rightward triangles (bottom) use the last four seasons. All error bars depict total errors, as in Fig. 6. The blue and greenshaded regions correspond to the interval covered by the total error bars using the full nine seasons, by the spline and regression difference techniques, respectively. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.