A&A 488, 481-490 (2008)
DOI: 10.1051/0004-6361:200809866

COSMOGRAIL: the COSmological MOnitoring of GRAvItational Lenses^,

VII. Time delays and the Hubble constant from WFI J2033-4723

C. Vuissoz¹ - F. Courbin¹ - D. Sluse¹ - G. Meylan¹ - V. Chantry^2, - E. Eulaers² - C. Morgan^3,4 - M. E. Eyler⁴ - C. S. Kochanek³ - J. Coles⁵ - P. Saha⁵ - P. Magain² - E. E. Falco⁶

1 - Laboratoire d'Astrophysique, École Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland
2 - Institut d'Astrophysique et de Géophysique, Université de Liège, Allée du 6 août 17, Sart-Tilman, Bât. B5C, 4000 Liège, Belgium
3 - Department of Astronomy and the Center for Cosmology and Astroparticle Physics, The Ohio State University, Columbus, OH 43210, USA
4 - Department of Physics, United States Naval Academy, 572C Holloway Road, Annapolis MD 21402, USA
5 - Institute of Theoretical Physics, University of Zürich, Winterthurerstrasse 190, 8057 Zürich, Switzerland
6 - Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge MA 02138, USA

Received 28 March 2008 / Accepted 2 July 2008

Abstract
Gravitationally lensed quasars can be used to map the mass distribution in lensing galaxies and to estimate the Hubble constant $H_{\rm0}$ by measuring the time delays between the quasar images. Here we report the measurement of two independent time delays in the quadruply imaged quasar WFI J2033-4723 (z=1.66). Our data consist of R-band images obtained with the Swiss 1.2 m EULER telescope located at La Silla and with the 1.3 m SMARTS telescope located at Cerro Tololo. The light curves have 218 independent epochs spanning 3 full years of monitoring between March 2004 and May 2007, with a mean temporal sampling of one observation every 4th day. We measure the time delays using three different techniques, and we obtain $\Delta t_{B-A} = 35.5$ $\pm$ 1.4 days (3.8%) and $\Delta t_{B-C} = 62.6^{+~4.1}_{-~2.3}~{\rm days} ~\ (^{+~6.5\%}_{-~3.7\%})$, where A is a composite of the close, merging image pair. After correcting for the time delays, we find R-band flux ratios of F_A/F_B = 2.88 $\pm$ 0.04, F_A/F_C = 3.38 $\pm$ 0.06, and F_A1/F_A2 = 1.37 $\pm$ 0.05 with no evidence for microlensing variability over a time scale of three years. However, these flux ratios do not agree with those measured in the quasar emission lines, suggesting that longer term microlensing is present. Our estimate of $H_{\rm0}$ agrees with the concordance value: non-parametric modeling of the lensing galaxy predicts $H_{\rm0}$ = 67 ⁺¹³_-10 km s^-1 Mpc^-1, while the Single Isothermal Sphere model yields $H_{\rm0}$ = 63 ⁺⁷_-3 km s^-1 Mpc^-1 (68% confidence level). More complex lens models using a composite de Vaucouleurs plus NFW galaxy mass profile show twisting of the mass isocontours in the lensing galaxy, as do the non-parametric models. As all models also require a significant external shear, this suggests that the lens is a member of the group of galaxies seen in field of view of WFI J2033-4723.

Key words: gravitational lensing - cosmology: cosmological parameters - quasars: individual: WFI J2033-4723

1 Introduction

When a quasar is gravitationally lensed and we observe multiple images of the source there are light travel time differences between the images. Any intrinsic variation of the quasar is observed at different times in each image with a measurable ``time delay'' between them. Refsdal (1964) first showed that time delays provide a means of determining the Hubble constant $H_{\rm0}$ independent of any local distance calibrator, provided a mass model can be inferred for the lensing galaxy. Conversely, one can also assume $H_{\rm0}$ in order to study the distribution of the total mass in the lensing galaxy.

Figure 1: The 6.3$^\prime$ $\times$ 6.3$^\prime$ field of view around WFI J2033-4723. This image is the central part of a combination of 418 R-band frames obtained with the 1.2 m EULER Telescope with a total exposure time of 48 h and a mean seeing of 1.3 $^{\prime \prime }$. The three stars PSF1-3 used to model the Point Spread Function (PSF) and the 7 reference stars S4-10 used for frame alignment and flux calibration are indicated.

Open with DEXTER

During the past 25 years, time delays have been measured in only 17 systems, at various accuracy levels (see Oguri 2007, for a review). As the error in the time delay propagates directly into $H_{\rm0}$, it is important to make it as small as possible. Unfortunately, most existing time delays have uncertainties of the order of 10% that are comparable to the current uncertainties in $H_{\rm0}$. This uncertainty can be reduced by increasing the sample of lenses with known time delays, and by simultaneously fitting all lenses in the sample with a common value for $H_{\rm0}$ (Saha et al. 2006b; Coles 2008; Saha et al. 2006a).

COSMOGRAIL is an optical monitoring campaign that aims at measuring time delays for a large number of gravitationally lensed quasars to accuracies of a few percent using a network of 1- and 2-m class telescopes. The first result of this campaign was the measurement of the time delay in the doubly imaged quasar SDSS J1650+4251 to an accuracy of 3.8% based on two observing seasons of data (Vuissoz et al. 2007). COSMOGRAIL complements a second monitoring group whose most recent results are a delay for HE 1104-1805 (Poindexter et al. 2007). In this paper we present time-delay measurements for the quadruply imaged quasar WFI J2033-4723 based on merging 3 years of optical monitoring data from the two groups. In a companion effort, Morgan et al. (2008) analyzed the merged data for the two-image lens QJ0158-4325, succeeding in measuring the size of the source accretion disk but failing to measure a time delay due to the high amplitude of the microlensing variability in this system.

WFI J2033-4723 (20$^{\rm h}$33$^{\rm m}$42 $.\!\!^{\rm s}$08, -47$^$23$^\prime$43 $.\!\!^{\prime\prime}$0; J2000.0) was discovered by Morgan et al. (2004) and consists of 4 images of a z =1.66 quasar with a maximum separation of 2.5 $^{\prime \prime }$. The lens galaxy was identified by Morgan et al. (2004) and has a spectroscopic redshift of $z_{\rm lens}$ = 0.661 $\pm$ 0.001 (Eigenbrod et al. 2006). The lens appears to be the member of a group, with at least 6 galaxies within 20 $^{\prime \prime }$ of the lens (Morgan et al. 2004), and we will have to account for this environment in any lens model.

We describe the monitoring data and its analysis in Sect. 2. In Sect. 3 we present the near-IR Hubble Space Telescope (HST) observations that we used to obtain accurate differential astrometry of the lens components and surface photometry of the lens galaxy. We estimate the time delays in Sect. 4 and model them using parametric (Sect. 5) and non-parametric (Sect. 6) models for the mass distribution of the lens galaxy. We summarize our results in Sect. 7.

   
2 Photometric monitoring

Our 3-year photometric monitoring of WFI J2033-4723 was carried out from March 2004 to May 2007 with the 1.2 m EULER telescope and the 1.3 m SMARTS telescope located in Chile at La Silla and the Cerro Tololo Interamerican Observatory (CTIO), respectively. WFI J2033-4723 was monitored from both sites for three years during its visibility window from early March to mid-December.

The 1.2 m EULER telescope is equipped with a 2048 $\times$ 2048 CCD camera which has 0.344 $^{\prime \prime }$ pixels and produces an image with an 11$^\prime$ field of view. The mean sampling of the EULER data is one epoch every five days, where each epoch consists of five dithered 360 s images taken with an R-band filter. The worst gaps due to weather or technical problems are 2-3 weeks. The EULER data set consists of 141 epochs of data obtained between May 2004 and May 2007. The image quality varies between 0.9 $^{\prime \prime }$ and 2.0 $^{\prime \prime }$ FWHM over 3 years, with a median of 1.4 $^{\prime \prime }$.

The 1.3 m SMARTS telescope is equipped with the dual-beam ANDICAM (DePoy et al. 2003) camera. Here we use only the optical channel which has 0.369 $^{\prime \prime }$ pixels and a 6.5$^\prime$ $\times$ 6.3$^\prime$ field of view. The mean sampling of the SMARTS data is one epoch every eight days, with three 300 s exposures at each epoch. The SMARTS data set consists of 77 epochs of data obtained between March 2004 and December 2006. The seeing on the images varies between 0.5 $^{\prime \prime }$ and 2.0 $^{\prime \prime }$, with a median of 1.4 $^{\prime \prime }$.

The combined data set consists of 218 observing epochs comprising 956 images covering the common field of view shown in Fig. 1. The average temporal sampling when WFI J2033-4723 was visible is one observation every 4 days over a period of three years, one of the best sets of monitoring data available for a lensed quasar.

Figure 2: Result of the simultaneous deconvolution of the 956 R-band images (EULER+SMARTS) of WFI J2033-4723. The pixel size of this image is half the pixel size of the EULER detector, i.e., 0.172 $^{\prime \prime }$, and the resolution is 2 pixels Full-Width-Half-Maximum, i.e., 0.344 $^{\prime \prime }$. The field of view is 22 $^{\prime \prime }$ on a side. Two galaxies G2 and G3 are seen to the West and North of the main lensing galaxy G1. G3 is part of a group included in the lens modeling, while G1 and G2 are modeled individually (see Sect. 5).

Open with DEXTER

The EULER data are reduced using the automated pipeline described in Vuissoz et al. (2007) and the SMARTS data with the SMARTS pipeline, using standard methods. The reduced frames are then aligned and interpolated to a common reference frame, one of the best-seeing (1 $^{\prime \prime }$) EULER images, taken on the night of 5 April 2006. The 10 stars (PSF1-3 and S4-10) shown in Fig. 1 are used to determine the geometric transformation needed for each EULER and SMARTS image to match the reference frame. The transformation includes image parity, rotation, shifting and rescaling. These 10 stars are also used to determine the photometric zero point of each image relative to the reference image. After interpolation, cosmic rays are removed using the L.A.Cosmic algorithm (van Dokkum 2001). We check that no data pixels are mistakenly removed by this process.

The light curves of the quasars are measured by simultaneously deconvolving all the images using the MCS algorithm (Magain et al. 1998). This method has already been successfully applied to the monitoring data of several lensed quasars (e.g. Burud et al. 2002b; Hjorth et al. 2002; Burud et al. 2002a; Vuissoz et al. 2007). The deconvolved images have a pixel scale of 0.172 $^{\prime \prime }$ (one-half the pixel scale of the EULER data) and are constructed to have a Gaussian PSF with a 2 pixel (0.344 $^{\prime \prime }$) FWHM. The Point Spread Function (PSF) for each of the 956 images is constructed from the three stars PSF1-3 (see Fig. 1). During the deconvolution process, the relative positions of the quasar images are held fixed to the values derived from fitting the HST images in Sect. 3, while their fluxes are allowed to vary from frame to frame. The flux, position and shape of the lensing galaxy are the same for all frames, but the values vary freely as part of the fit.

Figure 3: Our R-band light curves obtained for WFI J2033-4723 as well as for the reference star S6 (see Fig. 1). The magnitudes are given in relative units. The filled symbols correspond to the EULER observations while the SMARTS data points are marked by open symbols. Components A1 and A2 were summed into one single component A. The curves have been shifted in magnitude for visibility, as indicated in the figure.

Open with DEXTER

Figure 2 shows an example of a deconvolved image. It is clear that we will have no difficulty estimating the fluxes of components B and C separately. Components A₁ and A₂, however, are separated by only 0.724 $^{\prime \prime }$, which is only twice the resolution of our deconvolved images, and remain partially blended after deconvolution. Since the delay between these images should be very small, we will sum the fluxes of the two images and consider only the light curve of the total flux A=A₁+A₂. The resulting R-band light curves are displayed in Fig. 3.

We also display in Fig. 3 the deconvolved light curve of the isolated star S6, which has roughly the same color as WFI J2033-4723. Each point is the mean of the images taken at a given epoch and the error bar is the 1$\sigma$ standard error of the mean. The light curve is flat, with a standard deviation over the 3 years of monitoring of $\sigma_{{\rm tot}} = 0.010$ mag about its average, which is roughly consistent with the mean error bar of $\sigma_{{\rm mean}} = 0.006$ mag of the individual epochs.

The dispersion of the points in the star's light curve reflects both statistical errors and systematic errors from the photometric calibrations and the construction of the PSF. To the extent that all the light curves suffer from the same systematic errors, we can correct the quasar's light curves by subtracting the residuals of the star's light curve from each of them. We then define the uncertainty in a quasar's light curve as the quadrature sum of the uncertainties in the two light curves. This procedure will increase the photon noise but should minimize the systematic errors.

   
3 HST near-IR imaging

Figure 4: Left: deep NICMOS2 image, taken in the F160W-band. This image is a combination of 4 frames, for a total exposure time of 46 min. North is up and East to the left. The field of view is 4 $^{\prime \prime }$ on a side. Middle: simultaneous deconvolution of the individual NICMOS images (see text), using the MCS deconvolution algorithm. The PSF in this image is an analytical Gaussian with 2 pixels Full-Width-at-Half-Maximum ( FWHM), i.e., the resolution is 0.075 $^{\prime \prime }$. The pixel size is 0.035 $^{\prime \prime }$, i.e., oversampled by a factor of two compared to the original pixel size. Right: residual map of the deconvolution, with the cut levels set to $\pm$$5\sigma$. Only minor structures are seen in the center of the sharpest objects, which is acceptable given the quality of the NICMOS PSF.

Open with DEXTER

We determine the relative positions of the lens components and the light profile for the main lens galaxy G1 and its closest neighbor G2 (see Fig. 2) by analyzing our HST images of the system. These data were obtained in the framework of the CASTLES program (Cfa-Arizona Space Telescope LEns Survey), which provides HST images for all known gravitationally lensed quasars. We deconvolve these images using a modified version of the MCS deconvolution algorithm for images with poor knowledge of the PSF (Magain et al. 2007). We previously used this approach to unveil the faint Einstein ring and the lensing galaxy hidden in the glare of the quasar images of the so-called ``cloverleaf'' HE1413+117 (Chantry & Magain 2007). We analyze the Near Infrared Camera and Multi-Object Spectrometer (NICMOS) F160W (H-band) images obtained with the NIC2 camera. The data consist of four dithered MULTIACCUM images with a total exposure time of 2752 s and a mean pixel scale of 0.07568 $^{\prime \prime }$ (Krist & Hook 2004). We calibrate the images using CALNICA, the HST image reduction pipeline, subtract constants for each quadrant of NIC2 to give each image a mean background of zero, and create a noise map for each frame.

The images are simultaneously deconvolved in a manner similar to that used for the EULER and SMARTS data. The NICMOS frames lack any PSF stars, so we construct the PSF using the quasar images themselves in the iterative method of Chantry & Magain (2007). We first estimate the PSF of each frame using Tiny Tim (Krist & Hook 2004) and then we deconvolve them to have the final Gaussian PSF. During the deconvolution, each image is decomposed into a set of point sources and any extended flux. The latter is then reconvolved to the resolution of the original data and subtracted from the four initial frames, leading to images with far less contamination by extended structures. Four new PSFs are estimated from these new images, and we carry out a new deconvolution. The process is repeated until the residual maps are close to being consistent with the estimated noise (e.g. Courbin et al. 1998). In this case, convergence is reached after three iterations and the final reduced $\chi^{2}$ is 3.59. The final deconvolved image shown in Fig. 4 has half the pixel scale of the initial images and a Gaussian PSF with a FWHM of 0.075 $^{\prime \prime }$.

As part of the MCS deconvolution we also fit analytical models to the main lens galaxy (G1) and its nearby companion G2. The main lens is an early-type galaxy (Eigenbrod et al. 2006), as its companion is likely to be, so we use elliptical de Vaucouleurs profiles for both. The uncertainties are estimated by the scatter of the measurements from a separate set of fits to each independent image. We also estimate that there are systematic errors in the astrometry from the NICMOS pixel scale and focal plane distortions of order 2 milli-arcsec based on our earlier fits to the NICMOS data of H1413+117 (Chantry & Magain 2007). These systematic errors are compatible with the Lehár et al. (2000) comparison of NICMOS and VLBI astrometry for radio lenses.

The relative astrometry and photometry of the lens components and of the lensing galaxies G1 and G2 are given in Table 1. Coordinates are relative to image B, in the same orientation as Fig. 4. The photometry is in the Vega system. For each measurement, we give the $1 \sigma$ internal error bars, to which a systematic error of 2 milli-arcsec should be added. The models for G1 and G2 are presented in Table 2, with the effective semi-major and semi-minor axes of the light distribution a₀ and b₀. Each measurement is accompanied by its $1 \sigma$ error bar.

Table 1: Relative astrometry and photometry for the four components of WFI J2033-4723 and for the lensing galaxies G1 and G2.

Table 2: Shape parameters for the main and secondary lensing galaxies.

   
4 Time delay measurement

We measure the time delays between the blended light curve of A₁/A₂ and images B and C using three different techniques: i the minimum dispersion method of Pelt et al. (1996); ii the polynomial method of Kochanek et al. (2006); and iii the method of Burud et al. (2001). Since WFI J2033-4723 shows well-defined variations (see Fig. 3), it is already clear by visual inspection that $\Delta t_{B-A} \sim 35$ days and $\Delta t_{B-C} \sim 65$ days.

   
4.1 Minimum dispersion method

	Figure 5: Example of a dispersion curve as obtained from the minimum dispersion method, for components B and A. The position of the parabola minimum gives the time delay. Each point is separated by 2 days, i.e. about half the data mean sampling. The time delay indicated here is for only one realization of the boot-strap procedure (see text).
Open with DEXTER

Figure 6: Variation of the dispersion function minimum  $D_{{\rm min}}$ (red solid curves), as a function of the magnitude shift used for the normalization (see Sect. 4.1). Each $D_{{\rm min}}$ corresponds to a different estimate of the time delay, indicated on the blue dotted curves. The final time delay is the one with the lowest  $D_{{\rm min}}$. The top panel is for the B-A time delay, the bottom one for B-C. Time delays indicated here are for only one realization of the boot-strap procedure (see text).

Open with DEXTER

In the minimum dispersion method, time delays are computed for pairs of light curves using a cross-correlation technique that takes into account irregular sampling. The two light curves are first normalized to have zero mean. Then, one of the light curves is used as a reference and the second curve is shifted relative to it by a range of time delays. For each delay, we calculate the mean dispersion between the shifted points and their nearest temporal neighbors in the reference light curve. The best time-delay estimate is the one that minimizes this dispersion function. Since the mean sampling of our curves is one epoch every four days and since there is a limit to the number of time delays that can be tested independently, we test time delays in steps of 2 days. Figure 5 shows an example of a dispersion curve where we have then fit a parabola and set the best time delay to be the one corresponding to the minimum of the parabola.

There is, however, a complication in the step of normalizing the light curves, arising from sampling the light curve of each lensed image over a different time period of the intrinsic source variability (Vuissoz et al. 2007). We solve this problem by computing the dispersions as a function of a small magnitude shift $\Delta m$ in the normalization, measuring both the minimum dispersion  $D_{\rm min}(\Delta m)$ and the best fitting time delay $\Delta t(\Delta m)$ as shown in Fig. 6. Our final value for the delay is the one corresponding to the shift $\Delta m$ that minimizes the overall dispersion.

We then estimate the uncertainties by randomly perturbing the data points, based on a Gaussian distribution with the width of the measurement errors, and computing the dispersions and time delays again. We define the uncertainties by the $1 \sigma$ dispersion in the results for 100 000 trials (Vuissoz et al. 2007). The resulting uncertainty estimates are symmetric about the mean, so our final estimates based on this method are

$$\Delta t_{B-A} = 35.6 \pm 1.3 ~\ {\rm days} ~\ (3.6\%)$$

$$\Delta t_{B-C} = 64.6 \pm 3.4 ~\ {\rm days} ~\ (5.3\%).$$ (1)

While we have not taken microlensing effects into account for this analysis, it should matter little, as the method is not very sensitive (Eigenbrod et al. 2005) to the very low amplitude microlensing variability observed for this system (see Sect. 4.2).

   
4.2 Polynomial fit of the light curves

	Figure 7: Top: best polynomial fit to the light curves, which are shifted vertically for display purpose. Bottom: the residuals of the fit.
Open with DEXTER

In the polynomial method (Kochanek et al. 2006), the intrinsic light curve of the source is modeled as a polynomial that is simultaneously fit to all three light curves. Each quasar image has an additional low order polynomial to model slow, uncorrelated variability due to microlensing. We increase the source polynomial order for each season until the improvement in the $\chi^{2}$ statistics of the fits cease to be significant. This results in using polynomial orders of 11, 10, 17 and 4 for the four seasons of data. The low amplitude microlensing variations are modeled with a simple linear function for the four seasons. Figure 7 shows the best fits to the data and the residuals from the model. The effects of microlensing in this system are very small, with variations of only $\sim$0.01 mag over three years. As with the minimum dispersion method, we estimate the uncertainties by randomly perturbing the light curves 100 000 times and using the standard deviation of the trials as the error estimates to find that

$$\Delta t_{B-A} = 35.0 \pm 1.1 ~\ \hbox{days} ~\ (3.0\%)$$

$$\Delta t_{B-C} = 61.2 \pm 1.5 ~\ \hbox{days} ~\ (2.4\%).$$ (2)

   
4.3 Numerical modeling of the light curves

	Figure 8: Top: the light curves of the quasar images shown along with the best numerical model. Bottom: the residuals of the fit.
Open with DEXTER

Our last approach is based on the method described in Burud et al. (2001), which determines the time delay between a pair of light curves using a gridded numerical model for the source light curve. For a series of time delays, we fit the data with a flux ratio between the two curves, and a linear trend for microlensing on each full light curve. The numerical source model is smoothed on the scale of the temporal sampling, based on a smoothing function weighted by a Lagrange multiplier. The best time delay is the one that minimizes the $\chi^{2}$ between the model and data points.

This method has several advantages: first, none of the data light curves is taken as a reference: they are all treated on an equal basis. Furthermore, as the model is purely numerical, no assumption is made on the shape of the quasar's intrinsic light curve (except that it is sufficiently smooth). Finally, a model light curve is obtained for the intrinsic variations of the quasar, as for the polynomial fit method (see Sect. 4.2).

We have applied this method to the two pairs of light curves of WFI J2033-4723. The resulting fits to the light curves and their residuals are shown in Fig. 8. Using a Monte Carlo method to estimate the uncertainties, we find from 7000 trials (adding normally distributed random errors with the appropriate standard deviation on each data point):

$$\Delta t_{B-A} = 36.0 \pm 1.5 ~\ \hbox{days} ~\ (4.2\%)$$

$$\Delta t_{B-C} = 61.9~ ^{+~6.7}_{-~0.5} ~\ \hbox{days} ~\ (^{+~11\%}_{-~1\%}).$$ (3)

We note a secondary peak in the $\Delta t_{B-C}$ Monte Carlo distribution, around 69 days, in addition to the main peak at 61.9 days. There is, however, no evidence of such a secondary peak in the results of the minimum dispersion method and the polynomial fitting technique. This translates into an asymmetrical error bar on the result obtained with the numerical fitting of the light curves, and is taken into account in our final estimate of the time delay between quasar images B and C.

4.4 Final time delays

	Figure 9: Light curves of the three quasar images, shifted by their respective time delay and flux ratio. The blue circles correspond to image A, the red triangles to B and the green squares to C.
Open with DEXTER

For the final delay estimate we adopt the unweighted mean of the three methods, and we take as uncertainties the quadrature sum of the average statistical error and the dispersion of the results from the individual methods about their mean. Our final estimate of the time delays is

$$\Delta t_{B-A} = 35.5 \pm 1.4 ~\ \hbox{days} ~\ (3.8\%)$$

$$\Delta t_{B-C} = 62.6~ ^{+~4.1}_{-~2.3} ~\ \hbox{days} ~\ (^{+~6.5\%}_{-~3.7\%}).$$ (4)

We cannot measure the time delay between the individual A₁ and A₂ light curves, but values larger than $\Delta t_{A1-A2}=2$ days are incompatible with any of the models we consider in the following section. We can nevertheless estimate the flux ratio between A₁and A₂. After correcting for the time delays, we find that the R-band flux ratios between the images are

$$\frac{F_{A}}{F_{B}} = 2.88 \pm 0.04, ~\ \frac{F_{A}}{F_{C}} = 3.38 \pm 0.06, ~\ \frac{F_{A1}}{F_{A2}}= 1.37 \pm 0.05.$$ (5)

Figure 9 displays the light curves obtained for the three quasar images, after shifting by the time delays and flux ratios. Note that these flux ratios differ from those measured by Morgan et al. (2004) from the MgII broad emission lines, probably due to long-term microlensing (on a longer scale than our monitoring 3-year baseline), as discussed in the next section.

Table 3: Result of the parametric lens modeling.

   
5 Parametric modeling

5.1 Observational constraints

We constrain the mass models of WFI J2033-4723 using the positions of the four lensed images, the position of the lens galaxy G1 and the two delay measurements, for a total of 12 constrains. Except when indicated, we do not use the image flux ratios because they can be affected by extinction (Jean & Surdej 1998; Falco et al. 1999) and milli-lensing by substructures (Mao & Schneider 1998). We can also constrain the structure of G1 given its ellipticity e, position angle $\theta_e$ and effective radius $R_{\rm e}$ to the extent that these properties are correlated with those of its dark matter halo. Although a possible mismatch between the light and mass distributions is not impossible, we adopt the formal errors of 0.002 $^{\prime \prime }$ on the position of the lens G1 (Table 2). This is motivated by the small offset between the centers of mass and light found by Yoo et al. (2006) in a sample of four lensing galaxies.

Finally, WFI J2033-4723 is located in a group of galaxies, labeled G2-G6 in Fig. 10. We include G2 as a singular isothermal sphere (SIS) in all our models since it is close (4 $^{\prime \prime }$) and of similar luminosity to G1. As Morgan et al. (2004), we are unable to successfully model the system without including G2. When enough observational constraints are available we further add galaxies G3-G6 as a SIS mass distribution located at the barycenter  $G_{\rm group}$ of the group. In all models we include an external shear of amplitude $\gamma$ and position angle  $\theta_{\gamma}$ that represents the gravitational perturbations due to mass unaccounted for explicitly. We also experiment with including mass at the position of object X (Fig. 4, 2 $^{\prime \prime }$ North of G1) and find that doing so does not improve the models.

	Figure 10: Environment of WFI J2033-4723 as seen in an HST/ACS F814 (I-band) image. The main lens galaxy G1 and the close companion G2 were included in our analysis of the NICMOS image, and here we have labeled additional group members as G3-G6.
Open with DEXTER

We consider a sequence of standard mass models for G1, including a singular isothermal ellipsoid (SIE), a de Vaucouleurs (dVC) model and a NFW model (Navarro et al. 1997), and we fit the data using LENSMODEL (v1.99g) (Keeton 2001). The results are summarized in Table 3, where Cols. $\char93 1$ and $\char93 2$ describe the model family, and $\char93 3$ the mass parameter (either the Einstein radius b in arcseconds or the mean mass surface density  $\kappa_{\rm s}$, as defined in Keeton 2001). Column $\char93 4$ is for the ellipticity e and PA $\theta_e$ of the lens G1. Note that the measured PA of G1 is $\theta = 27.8^{\circ}$ (Table 2). Column $\char93 5$ gives the external shear amplitude $\gamma$ and PA  $\theta_{\gamma}$, $\char93 6$ the number of degrees of freedom for each model, and $\char93 7$ the resulting reduced $\chi^{2}$. Column $\char93 8$ finally shows the best estimate for h= $H_{\rm0}$/100. A minus sign in this column means that time delays are not used as constraints. All angles are given positive East of North, and values given in parentheses are fitted to the observations. All models assume $\Delta t_{A1-A2}=2$ days and include the companion galaxy G2, with a resulting mass $0.1~ m_{\rm G1} < m_{\rm G2} < m_{\rm G1}$.

5.2 SIE models

Our first model consists of an SIE for G1, an SIS for G2 and an external shear. When we fit only the image positions but include a prior on the position angle $\theta_e$ (from Table 1) we do not find a good fit unless the constraint on the position of the lensing galaxy is relaxed. The prior on the position angle is justified by statistical studies finding correlations between the position angles but not the axis ratios of the visible and total mass distributions (Keeton et al. 1997; Ferreras et al. 2008). With the inclusion of the time delays we have enough constraints to add the group halo to the model. With the position angle of G1 constrained, we obtain poor fits to the data with reduced $\chi^2=3.6$ for $N_{\rm d.o.f.}=2$. When we leave the position angle free, we find good fits but the model PA is $55^\circ$ from the observed. These models have Hubble parameters of $h\simeq 0.63^{+0.07}_{-0.03}$ with the spread dominated by the degeneracies between the ellipticity and the shear.

5.3 De Vaucouleurs models

Next we consider a constant mass-to-light ratio model of the lens galaxy based on a De Vaucouleurs model. The position angle and the effective radius $R_{\rm e}=0.608$ $^{\prime \prime }$ (the geometric mean of the semi-axes in Table 1, corresponding to 4 kpc for h=0.72) are constrained by the values measured for the galaxy. This model does not fit well the lens configuration ( $\chi^2\sim 34.4$), mainly due to the small uncertainty on the lens galaxy position. When we include the time delays we find a good fit ( $\chi^2\simeq 0.01$) as long as we allow G1 to be misaligned with respect to the observed galaxy. As expected from the reduced surface density compared to the SIE model (Kochanek 2002), we find a much higher value for the Hubble parameter, h=0.92.

5.4 NFW models

We use an NFW model with a fixed break radius of $r_{\rm s}=10$ $^{\prime \prime }$ (40 kpc), where the break radius is related to the virial radius through the concentration $c=R_{\rm vir}/r_{\rm s}$. Since $r_{\rm s}$ lies well outside the Einstein radius of the lens, its particular value is not important. This model is not realistic by itself because the shallow $\rho \propto 1/r$ central density cusp of the model will lead to a visible central image. We again find that we can fit the astrometry well even when the position angle of the lens is constrained, but we cannot do so after including the time delays unless we allow the model of G1 to be misaligned relative to the light. In any case, this model leads to a fifth image about 3 mag fainter than A that should be visible on our NICMOS data. This NFW model has a higher surface density near the Einstein ring than the SIE model, so we find a lower value for the Hubble parameter of $h \simeq 0.29$. Using an unphysically small break radius of $r_{\rm s}=1$ $^{\prime \prime }$ raises the density and hence the Hubble parameter to $h\simeq 0.63$.

5.5 De Vaucouleurs plus NFW models

As our final, and most realistic, parametrized model we combine a dVC model constrained by the visible galaxy with an NFW model for the halo. The two components are first constrained to have the same position and position angle, the parameters of the dVC model are constrained by the observations, and the NFW model has a fixed $r_{\rm s}=10$ $^{\prime \prime }$ break radius. This model leads to a poor fit, with $\chi^2=6.33$ for $N_{\rm d.o.f.}=1$. When we free the PA of the NFW model, we find an acceptable fit for $N_{\rm d.o.f.}=0$, but the misalignment of the NFW model relative to the optical is  $40^{\circ}$.

We can further constrain the model by including the three MgII emission line flux ratios from Morgan et al. (2004). We use the line flux ratios instead of those obtained from the light curves, because they should be insensitive to microlensing. With these three additional constraints we still find that the $89.8^{\circ}$ position angle of the NFW model is strongly misaligned from the $26.4^{\circ}$ position angle of the dVC model, indicating a twisting of the mass isocontours. The model has a reduced $\chi^{2}$ of 3.2 for $N_{\rm d.o.f.}=3$. The model flux ratios are significantly different from the constraints. We find F_A/F_B = 2.81, F_A/F_C = 5.01, and F_A1F_A2 = 1.26 while Morgan et al. (2004) report F_A/F_B = 2.55 $\pm$ 0.60, F_A/F_C = 2.02 $\pm$ 0.35, and F_A1/F_A2 = 1.61 $\pm$ 0.35. The match of the flux ratios is better if we do not include the constraints from the time delays. In all cases, F_B/F_C is the most ``anomalous'' flux ratio, as also found by Morgan et al. (2004). While the differences between the line and continuum flux ratios suggests the presence of long-term microlensing, we see no evidence for the time variability in the flux ratios expected from microlensing. We also note that the broad line flux ratios vary with wavelength (Morgan et al. 2004), which suggests that dust extinction may as well be affecting the flux ratios.

   
6 Non-parametric modeling

We use the non-parametric PixeLens (Saha & Williams 2004) method as our second probe of the mass distribution. This approach has the advantage that it makes fewer assumptions about the shape of the G1 than the ellipsoidal parametric models. The models include priors on the steepness of the radial mass profile, imposes smoothness criteria on the profile and we restrict to models symmetric about the lens center. We include two point masses to represent G2 and the group. We run 1000 trial models at a resolution of $\sim$ $0.23 \hbox{$^{\prime\prime}$ }$/pixel which are constrained to fit the image positions and the time delays exactly. For each model we vary the Einstein radii of G2 and the group over the ranges 0.03 $^{\prime \prime }$ < $R_{\rm E}({\rm G2})$ < 3 $^{\prime \prime }$ and 0.3 $^{\prime \prime }$ < $R_{\rm E}({\rm group})$ < 5 $^{\prime \prime }$, respectively. Apart from the inclusion of these additional point masses, the method and priors are as explained in detail in Coles (2008). A test, where the technique is used to infer $H_{\rm0}$ from a N-body and hydro simulated lens, and an additional discussion of the priors are described in Read et al. (2007). Figure 11 shows the resulting probability distribution for $H_{\rm0}$ from the 1000 models. The median value of the distribution is

$$\hbox{$H_{\rm0}$ } = 67.1~^{+13.0}_{-9.9} ~\ \hbox{km~s$^{-1}$ ~Mpc$^{-1}$ }$$ (6)

where the error bars are at 68% confidence. As already illustrated by our parametric modeling, the predicted $H_{\rm0}$ value depends on the density gradient of the models. Figure 12 shows the mean surface density contours of the models, and we see a twisting of the contours away from that of the visible galaxy in the outer regions.

	Figure 11: Distribution of $H_{\rm0}$ from 1000 non-parametric models. The bottom-axis gives the Hubble time, the top-axis $H_{\rm0}$ in km s-1 Mpc-1.
Open with DEXTER

Figure 12: Mean mass distribution from 1000 pixellated lens models. The red dots are the quasar images and the blue dot the source position. The thick solid line indicates the observed major axis of the lensing galaxy. Each tick-mark measures 1 $^{\prime \prime }$. The third contour from the outside traces the critical mass density $\Sigma _{\rm crit} = 1.19$ $\times$ $10^{11}~M_{\odot}$ arcsec-2, and the others are drawn logarithmically from the critical one (each contour is 2.5 larger/smaller than the previous one). North is to the top and East to the left.

Open with DEXTER

   
7 Conclusions

By combining data from COSMOGRAIL and the SMARTS 1.3 m telescope we measure two independent time delays in the quadruply imaged quasar WFI J2033-4723 (Morgan et al. 2004). The fractional uncertainties of $\sim$$4\%$ are among the best obtained so far from an optical monitoring. We obtain the light curves of the quasar images using the MCS deconvolution photometry algorithm (Magain et al. 1998) and then measure time delays using three different approaches with a final estimate of

$$\Delta t_{B-A} = 35.5 \pm 1.4 ~\ \hbox{days} ~\ (3.8\%)$$

$$\Delta t_{B-C} = 62.6~ ^{+~4.1}_{-~2.3} ~\ \hbox{days} ~\ (^{+~6.5\%}_{-~3.7\%})$$ (7)

where A is the mean light curve of the blended of quasar images A₁ and A₂. We find little evidence of microlensing in this system, which makes WFI J2033-4723 a very good system for measuring clean time delays.

The parametric models are consistent with concordance value of $H_{\rm0}$ when the lens galaxy has an isothermal mass profile out to the radius of the Einstein ring. As expected, the models allow higher (lower) values as the mass distribution is more centrally concentrated (extended) using de Vaucouleurs (NFW) mass distribution. The non-parametric models predict $H_{\rm0}$ = 67.1  ^+13.0_-9.9 km s^-1 Mpc^-1.

The addition of the time delays as a constraint on the lens models does not alter the mismatch between the observed and predicted image flux ratios. The largest flux ratio anomaly is the 45% difference between the MgII flux ratios found for images B/C. Morgan et al. (2004) also noted that the F_B/F_C flux ratio varies with wavelength, suggesting the presence of chromatic microlensing. The lack of significant variability in the flux ratio over our three year monitoring period suggests either that the effective source velocities in this lens are very low or that the affected images lie in one of the broad demagnified valleys typical of microlensing magnification patterns for low stellar surface densities.

Several galaxies close to the line of sight have a significant impact on the mass modeling. We generally model the potential as the sum of a main lensing galaxy G1, a companion galaxy G2 ($\sim$4 $^{\prime \prime }$ West of G1), and a nearby group ($\sim$9 $^{\prime \prime }$ North of G1). Both the parametric and non-parametric models suggest that the isodensity contours of G1 itself must be twisted, with some evidence that the outer regions are aligned with the tidal field of the group rather than the luminous galaxy. This could indicate that G1 is a satellite rather than the central galaxy of the group (e.g. Kuhlen et al. 2007). The twisting seems to be required even though the angular structure of the potential can be adjusted through the companion galaxy G2, an external tidal shear, and in some cases a group halo. Clarifying this issue requires more constraints such as detailed imaging of the Einstein ring image of the quasar host, measuring the redshifts of the nearby galaxies, and measuring the velocity dispersion of G1.

Acknowledgements

We would like to thank all the observers at the EULER telescope for the acquisition of these monitoring data. We also thank Profs. G. Burki and M. Mayor for their help in the flexible scheduling of the observing runs, improving the necessary regular temporal sampling. We are very grateful to Prof. A. Blecha and the technical team of the Geneva Observatory for their efficient help with the electronics of the CCD camera. Many thanks also to Sandrine Sohy for her commitment in the programming part of the deconvolution techniques. Finally, we would like to thank Chuck Keeton for his advice on the use of LENSMODEL. This work is supported by ESA and the Belgian Federal Science Policy Office under contract PRODEX 90195. C.S.K. is funded by National Science Foundation grant AST-0708082. COSMOGRAIL is financially supported by the Swiss National Science Foundation (SNSF).

References

• Burud, I., Magain, P., Sohy, S., & Hjorth, J. 2001, A&A, 380, 805 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Burud, I., Courbin, F., Magain, P., et al. 2002a, A&A, 383, 71 [NASA ADS] [CrossRef] [EDP Sciences]
• Burud, I., Hjorth, J., Courbin, F., et al. 2002b, A&A, 391, 481 [NASA ADS] [CrossRef] [EDP Sciences]
• Chantry, V., & Magain, P. 2007, A&A, 470, 467 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Coles, J. 2008, ApJ, 679, 17 [NASA ADS] [CrossRef]
• Courbin, F., Lidman, C., Frye, B. L., et al. 1998, ApJ, 499, L119 [NASA ADS] [CrossRef] (In the text)
• DePoy, D. L., Atwood, B., Belville, S. R., et al. 2003, in Presented at the Society of Photo-Optical Instrumentation Engineers (SPIE) Conference, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye, & A. F. M. Moorwood, Proc. SPIE, 4841, 827 (In the text)
• Eigenbrod, A., Courbin, F., Vuissoz, C., et al. 2005, A&A, 436, 25 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Eigenbrod, A., Courbin, F., Meylan, G., Vuissoz, C., & Magain, P. 2006, A&A, 451, 759 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Falco, E. E., Impey, C. D., Kochanek, C. S., et al. 1999, ApJ, 523, 617 [NASA ADS] [CrossRef]
• Ferreras, I., Saha, P., & Burles, S. 2008, MNRAS, 383, 857 [NASA ADS] [CrossRef]
• Hjorth, J., Burud, I., Jaunsen, A. O., et al. 2002, ApJ, 572, L11 [NASA ADS] [CrossRef]
• Jean, C., & Surdej, J. 1998, A&A, 339, 729 [NASA ADS]
• Keeton, C. R. 2001, ArXiv Astrophysics e-prints (In the text)
• Keeton, C. R., Kochanek, C. S., & Seljak, U. 1997, ApJ, 482, 604 [NASA ADS] [CrossRef]
• Kochanek, C. S. 2002, ApJ, 578, 25 [NASA ADS] [CrossRef] (In the text)
• Kochanek, C. S., Morgan, N. D., Falco, E. E., et al. 2006, ApJ, 640, 47 [NASA ADS] [CrossRef] (In the text)
• Krist, J., & Hook, R. 2004, The Tiny Tim's user Guide v.6.3 (In the text)
• Kuhlen, M., Diemand, J., & Madau, P. 2007, ApJ, 671, 1135 [NASA ADS] [CrossRef] (In the text)
• Lehár, J., Falco, E. E., Kochanek, C. S., et al. 2000, ApJ, 536, 584 [NASA ADS] [CrossRef] (In the text)
• Magain, P., Courbin, F., & Sohy, S. 1998, ApJ, 494, 472 [NASA ADS] [CrossRef] (In the text)
• Magain, P., Courbin, F., Gillon, M., et al. 2007, A&A, 461, 373 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Mao, S., & Schneider, P. 1998, MNRAS, 295, 587 [NASA ADS] [CrossRef] (In the text)
• Morgan, C. W., Eyler, M. E., Kochanek, C. S., et al. 2008, ApJ, 676, 80 [NASA ADS] [CrossRef] (In the text)
• Morgan, N. D., Caldwell, J. A. R., Schechter, P. L., et al. 2004, AJ, 127, 2617 [NASA ADS] [CrossRef] (In the text)
• Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] (In the text)
• Oguri, M. 2007, ApJ, 660, 1 [NASA ADS] [CrossRef] (In the text)
• Pelt, J., Kayser, R., Refsdal, S., & Schramm, T. 1996, A&A, 305, 97 [NASA ADS] (In the text)
• Poindexter, S., Morgan, N., Kochanek, C. S., & Falco, E. E. 2007, ApJ, 660, 146 [NASA ADS] [CrossRef] (In the text)
• Read, J. I., Saha, P., & Macciò, A. V. 2007, ApJ, 667, 645 [NASA ADS] [CrossRef] (In the text)
• Refsdal, S. 1964, MNRAS, 128, 307 [NASA ADS] (In the text)
• Saha, P., & Williams, L. L. R. 2004, AJ, 127, 2604 [NASA ADS] [CrossRef] (In the text)
• Saha, P., Coles, J., Macciò, A. V., & Williams, L. L. R. 2006a, ApJ, 650, L17 [NASA ADS] [CrossRef]
• Saha, P., Courbin, F., Sluse, D., Dye, S., & Meylan, G. 2006b, A&A, 450, 461 [NASA ADS] [CrossRef] [EDP Sciences]
• van Dokkum, P. G. 2001, PASP, 113, 1420 [NASA ADS] [CrossRef] (In the text)
• Vuissoz, C., Courbin, F., Sluse, D., et al. 2007, A&A, 464, 845 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Yoo, J., Kochanek, C. S., Falco, E. E., & McLeod, B. A. 2006, ApJ, 642, 22 [NASA ADS] [CrossRef] (In the text)