© ESO, 2011

1. Introduction

To determine the expansion rate of the Universe, H₀, accurately, it is important to scale the extragalactic distance ladder and to measure w(z), the redshift evolution of the dark energy equation of state parameter (e.g., Frieman et al. 2008).

Some of the most popular methods in use to measure H₀ have recently been reviewed (Freedman & Madore 2010). The huge observational and theoretical efforts invested in these measurements have led to random errors between 3% and 10%, depending on the methods. However, all these methods rely at some level on each other and are in the end all based on the same local standard candles. In addition, the accuracy required on H₀ to measure w(z) and H(z) is of the order of 1% (Hu 2005; Riess et al. 2009, 2011). It is therefore of interest (i) to explore methods that are fully independent of any standard candle and (ii) to combine the different methods to further reduce the current error bar on H₀.

Strong gravitational lensing of quasars and the so-called “time delay method” in multiply imaged quasars is independent of the traditional standard candles (Refsdal 1964) and is based on well understood physics: general relativity. However, it requires to measure the time delays from long-term photometric monitoring of the many lensed quasars, which has long been a serious observational limitation. While the attempts to measure accurate time delays have been numerous, only few quasars have been measured with an accuracy close to the percent (e.g., Goicoechea 2002; Fassnacht et al. 2002; Vuissoz et al. 2008). In addition, the lens model necessary to convert time delays into H₀ often remains poorly constrained and hampers a breaking of the degeneracies between the properties of the lensing galaxy and H₀. In recent years, successful attempts have been made to constrain the lens model as much as possible to break these degeneracies (Suyu et al. 2009, 2010).

COSMOGRAIL, the COSmological MOnitoring of GRAvItational Lenses, aims both at measuring precise time delays for a large sample of strongly lensed quasars, and at obtaining and using all necessary observations to constrain the lens models. The present paper describes the COSMOGRAIL results for the quadruply imaged quasar HE 0435-1223, using deep VLT spectra and deconvolved HST images.

HE 0435-1223 (α(2000): 04h38min14.9s; ) was discovered by Wisotzki et al. (2000) during the Hamburg/ESO Survey (HES) for bright quasars in the Southern Hemisphere. It was identified two years later as a quadruply imaged quasar by Wisotzki et al. (2002). The redshift of the source is z_s = 1.689 (Wisotzki et al. 2000) and that of the lens is z_l = 0.4546  ±  0.0002 (Morgan et al. 2005). The quasar shows evidence for intrinsic variability, which makes it a good candidate for determining the time delays between the different images. The local environment of the lensing galaxy has been studied in detail by Morgan et al. (2005) using the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) and by Momcheva (2009), who has found that the lensing galaxy lies in a group of at least 11 members. The velocity dispersion of this group is σ ~ 496 km s^-1.

Analytical lens models of HE 0435-1223 are given by Kochanek et al. (2006), who have also measured time delays from two years of optical monitoring:  days,  days, and  days. For a fixed H₀ = 72  ±  7   km   s^-1   Mpc^-1 they found that the lensing galaxy must have a rising rotation curve at the position of the lensed images and a non-constant mass-to-light ratio. Moreover, high dark matter surface densities are required in the lens halo. New monitoring data of Blackburne & Kochanek (2010) analysed using a physically motivated representation of microlensing give time delays compatible with those of Kochanek et al. (2006), although these authors do not provide their measured values.

2. Photometric monitoring

2.1. Optical imaging

HE 0435-1223 was monitored during more than six years, from January 2004 to March 2010, through the R filter, using three different telescopes: the Swiss 1.2 m Euler telescope located on the ESO La Silla site (Chile), the Belgian-Swiss 1.2 m Mercator telescope located at the Roque de Los Muchachos Observatory, La Palma, Canary Island (Spain), and the 1.5 m telescope located at the Maidanak Observatory (Uzbekistan). In addition we also use 136 epochs from the two-year long monitoring of Kochanek et al. (2006), from August 2003 to April 2005, obtained with the ANDICAM camera mounted on the 1.3 m Small and Moderate Aperture Research Telescope System (SMARTS) located at the Cerro Tololo Inter-American Observatory (CTIO) in Chile. A summary of the observations is given in Table 1. Note that we used the published photometry for the SMARTS data, i.e., we did not reprocess the original data frames with our own photometric pipeline.

Fig. 1 Part of the field of view of the 1.2 m Swiss Euler telescope, with HE 0435-1223 visible in the centre. The four PSF stars used for deconvolution purposes and the two reference stars used to carry out the flux calibration are indicated.

Fig. 2 Result from the simultaneous deconvolution of the ground-based frames. G is the lensing galaxy and G22 (Morgan et al. 2005) is its closest neighbour on the plane of the sky. The grey scale in the deconvolved image is set to display all light level above 3  ×  σsky. The FWHM resolution of the deconvolved image is 034.

Fig. 3 R-band light curves of the four lensed images of HE 0435-1223 from December 2003 to April 2010. The magnitudes are given in relative units as a function of the Heliocentric Julian Day (HJD), along with their total 1σ error bars. These light curves are available in tabular form at the CDS.

Fig. 4 Left: combination of the four original HST/NIC2 F160W frames of HE 0435-1223. The field of view is 9  ×  9 arcsec. Middle: deconvolved image, where the lensing galaxy is modelled as a de Vaucouleurs profile (see text). The nearest galaxy on the plane of the sky, G22, is also indicated. Right: residual map in units of the noise. The colour scale ranges from  − 4σ (white) to +4σ (black).

2.2. Image processing and deconvolution photometry

The data from the three telescopes used by the COSMOGRAIL collaboration are analysed with the semi-automated reduction pipeline described in Vuissoz et al. (2007). The main challenge was that we had to assemble data from different telescopes: each camera has a different size, resolution and orientation on the plane of sky.

The pre-reduction for each observing epoch consists of flat-fielding using master sky-flats. The Euler image with the best seeing was taken as the reference frame to register all other frames. This reference frame was taken on the night of November 11, 2005 and has a seeing of 082. Two reference stars were then chosen in the field of view (Fig. 1) to compute the geometrical transformations between the images. This transformation involves a spatial scaling and a rotation. The reference stars were also used to compute the relative photometric scaling between the frames taken at different epochs. Eventually, the L. A. Cosmic algorithm (van Dokkum 2001) was applied separately to every frame to remove cosmic rays. All images were checked visually to make sure that no pixel was removed inappropriately, especially in the frames with good seeing.

The photometric measurements were carried out using “deconvolution photometry” with the MCS deconvolution algorithm (Magain et al. 1998). This software has been successfully applied to a variety of astrophysical problems ranging from gravitationally lensed quasars (e.g., Burud et al. 2000, 2002) to the study of quasar host galaxies (e.g., Letawe et al. 2008), or to the search for extrasolar planets using the transit technique (e.g., Gillon et al. 2007a,b). Image deconvolution requires accurate knowledge of the instrumental and atmospheric point spread function (PSF). The latter was computed for each frame from the four stars labelled PSF 1-4 in Fig. 1. These stars are from 0 to 1 mag brighter than the quasar images in HE 0435-1223 and are located within 2′ from the centre of the field, which minimizes PSF distortions.

Because it does not attempt to achieve an infinitely high spatial resolution, the MCS algorithm produces deconvolved images that are always compatible with the sampling theorem. This avoids deconvolution artefacts and allowed us to carry out accurate photometry over the entire field of view. Moreover, the deconvolved image was computed as the sum of extended numerical structures and of analytical point sources whose shape is chosen to be symmetrical Gaussians. In the case of gravitationally lensed quasars, the numerical channel of this decomposition contains the lensing galaxy. The photometry and astrometry of the quasar images were returned as a list of intensities and positions of Gaussian deconvolved profiles. Finally, the deconvolved image can be computed on a grid of pixels of arbitrary size. In the present work, the pixel size in the deconvolved frames is half the pixel size of the Euler data, i.e., 0172. The spatial resolution in the deconvolved frames is two pixels full-width-at-half-maximum (FWHM), i.e., 035.

With the MCS software, dithered images of a given target can be “simultaneously deconvolved” and combined into a single deep and sharp frame that matches the whole dataset at once, given the PSFs and the noise maps of the individual frames. In doing this, the intensities of the point sources are allowed to vary from one frame to the next while the smooth background, which includes the lensing galaxy, is held constant in all frames. The result of the process is shown in Fig. 2, where the point sources are labelled as in Wisotzki et al. (2002). Prior information on the object to be deconvolved can be used to achieve the best possible results. In the case of HE 0435-1223 the relative positions of the point sources are fixed to the HST astrometry obtained in Sect. 3.

Figure 3 shows the deconvolution light curves obtained for each quasar image of HE 0435-1223, where the 1σ error bars account both for the statistical and systematic errors. The statistical part of the error was taken as the dispersion between the photometric points taken during each night. The systematic errors were estimated by carrying out the simultaneous deconvolution of reference stars in the vicinity of HE 0435-1223.

Finally, a small scaling factor was applied to the light curves of all telescopes, including the published light curves of Kochanek et al. (2006), to match the Euler photometry. These shifts are all smaller than 0.03 mag.

3. HST NICMOS2 imaging

We used deep near-IR HST images of HE 0435-1223 to derive the best possible relative astrometry between the quasar images and the lensing galaxy and to constrain the light distribution in the lensing galaxy. The data are part of the CASTLES project (Cfa-Arizona Space Telescope LEns Survey) and were acquired in October 2004 (PI: C. S. Kochanek) with the camera 2 of NICMOS, the Near-Infrared Camera and Multi-Object Spectrometer. They consist of four dithered frames taken through the F160W filter (H-band) in the MULTIACCUM mode with 19 samples and calibrated by CALNICA, the HST image reduction pipeline. The total exposure time amounts to approximately 44 min and the pixel scale is 0075652.

The MCS deconvolution algorithm was used to combine the four NIC2 frames into a deep sharp IR image. We followed the iterative technique described in Chantry et al. (2010) and Chantry & Magain (2007), which allowed us to build a PSF in the absence of a stellar image in the field of view. The method can be summarised as follows. First, we estimated the PSF using Tiny Tim software (Krist & Hook 2004) and carried out the simultaneous deconvolution of the four F160W frames using a modified version of the MCS software (Magain et al. 2007). This produces a first approximation of the extended channel of the deconvolved image, i.e., the lensing galaxy and the lensed quasar host galaxy. We reconvolved the latter by the PSF and subtracted it from the original data. A new estimate of the PSF was built on the new image that now contained only the quasar images. The process was repeated until the residual image was satisfactory (for more details see Chantry et al. 2010). Figure 4 shows the result. In this image the pixel size is half that of the original data and the resolution 0075 (FWHM), unveiling an almost full Einstein ring.

In the final deconvolved image, the lensing galaxy was modelled analytically rather than numerically to minimise the number of degrees of freedom. We found that the best-fit profile is an elliptical de Vaucouleurs with the parameters as given in Table 3. The astrometry of the quasar images relative to the lensing galaxy, corrected for the known distortions of the NIC2 camera and for the difference of pixel scale between the x and y directions is summarised in Table 2. Based on our previous work using deconvolution of NICMOS images (Chantry & Magain 2007), we estimate that the total error bars, accounting for residual correction of the distortions amounts to 2 mas. Our results agree well with previous measurements from HST/ACS (Morgan et al. 2005) or HST/NIC2 imaging (Kochanek et al. 2006), also shown for comparison in Table 2.

4. Time delay measurement

4.1. Curve shifting method

Our method to measure the time delays is based on the dispersion technique of Pelt et al. (1996): the light curves are shifted in time and in magnitude to minimise a global dispersion function. In addition, the light curves are distorted on long time scales to account for slow microlensing variations. This was made by adding low-order polynomials to either the full curves or to specific observing seasons.

Pelt et al. (1996) has defined several dispersion statistics between pairs of light curves. We implemented a dispersion estimate similar to (see Eq. (8) of Pelt et al. 1996), which performed a linear interpolation between points of one of the curves over a maximum range of 30 days. In the case of four light curves, we defined a total dispersion that is the sum of the dispersions computed using the 12 possible permutations of two curves among four. Each pair was considered twice so to avoid the arbitrary choice of a reference light curve. The photometric error bars were taken into account to weight the influence of the data points in the dispersion. We then minimised the total dispersion by modifying the time delays and the microlensing polynomials.

4.2. Microlensing and influence on the time delay

Simulated light curves that mimic the observed data were used to estimate the robustness of the method. The error bars on the time delays were calculated using Monte Carlo simulations, i.e., redistributing the magnitudes of the data points according to their photometric error bars. The width of the resulting time delay distributions gives us the 1σ error bars.

Because of microlensing we do not have access to the intrinsic variations of the quasar. We represent microlensing in three of the light curves as a relative variation with respect to the fourth light curve, taken as a reference. We tested each of the four light curves in turn as a reference and kept the one that minimised the residual microlensing variations to be modelled in the 3 others. This is best verified with component B as a reference.

With B as a reference light curve, we note that microlensing in C and D remains smooth and can therefore be modelled with a low-order polynomial drawn over the full length of the monitoring. However, A contains higher frequency variations that need to be accounted for in each season individually, as illustrated in Fig. 5. In doing this, we obtain fairly good fits to the light curves, as shown in the residual signal. To quantify the quality of these residuals, we applied the so-called one-sample runs test of randomness, a statistical test to estimate whether successive realizations of a random variable are independent or not. In practice the test was applied to a sequence of residuals to decide whether a model is a good representation of the data. For most seasons in our curves the number of runs was between 1σ and 3σ lower than the value expected for independent random residuals. Thus, although our microlensing model is not fully representative of the real signal, the deviations from the data points remain small.

Fig. 5 Light curves obtained with all four telescopes and shifted by time delays of ΔtBA = 8.4 days, ΔtBC = 7.8 days and ΔtBD =  −6.5 days. The relative microlensing representations applied on curves A, C and D are shown as continuous curves with respect to the dashed blue line (see text). A fifth-order polynomial was used over the seven seasons to model microlensing on the quasar images C and D, while seven independent third-order polynomials were used for image A. The lower panels show residuals obtained by subtracting a simultaneous spline fit (grey) from the light curves.

We tested the robustness of our curve-shifting method in several ways. First, we modelled the microlensing variations using polynomial fits of different orders. Second, we fitted these polynomials either across each individual season or across groups of seasons. Finally, we masked the seasons with the worst residual signal (Fig. 5). All these changes had only a negligible impact on the time delay measurements. We note that this is not the case when considering only two or three seasons of data, which shows the importance of a long-term monitoring with good temporal sampling.

4.3. Final results

Our results are summarised in Table 4 and are compared with the previous measurements of Kochanek et al. (2006), who have used pure a polynomial fit to the light curves and two seasons of monitoring. Using the same data but with our modified dispersion technique, we obtained very similar time delays as Kochanek et al. (2006), but larger error bars. We prefer keeping a minimum possible number of degrees of freedom (e.g., in the polynomial order used to represent microlensing), in accordance with the Occam’s razor principle, even to the cost of apparently larger formal error bars.

We also note that Kochanek et al. (2006) give their time delays with respect to A, which, with seven seasons of data, turns out to be the most affected by microlensing. As a consequence the error bars on these time delays are dominated by residual microlensing rather than by statistical errors. The time delays used in the rest of our analysis are therefore measured relative to B.

Finally, we used the mean values of our microlensing corrections to estimate the macrolensing R-band flux ratios between the four quasar images, assuming that no long-term microlensing affects the data. We found m_B − m_A = Δm_BA = 0.62  ±  0.04, Δm_BC = 0.05  ±  0.01, and Δm_BD =  −0.16  ±  0.01, which is well compatible with the ratios measured at seven wavelengths by Mosquera et al. (2011). However, these authors report significant wavelength dependence of the image flux ratios, which led us not to use flux ratios as a constraint in the lens models.

5. Constraining the mass profile of the lensing galaxy

The goal of the present section is to constrain the radial mass profile of the lensing galaxy as much as possible, which is the main source of uncertainty on the determination of H₀ with the time delay method. One way of doing this is to use the information contained in the Einstein ring that is formed by the host galaxy of the lensed quasar (Suyu et al. 2010, 2009; Warren & Dye 2003). This works well when a prominent Einstein ring is visible. Unfortunately, given the depth of the current HST images of HE 0435-1223, the radial extent of the ring is too small to efficiently apply this technique. We propose instead to use information on the dynamics and on the stellar mass of the lens, using deep optical spectroscopy.

5.1. Stellar population and velocity dispersion of the lensing galaxy

A deep VLT spectrum of the lensing galaxy is available from Eigenbrod et al. (2006). While the spectrum was originally used to measure the redshift of the galaxy, it turns out to be deep enough ( ⟨ S/N ⟩   ~  20) to measure the stellar velocity dispersion and the mass-to-light ratio.

We analysed the data using full spectrum fitting with the ULySS package (Koleva et al. 2009). The method consists in fitting spectra against a grid of stellar population models convolved by a line-of-sight velocity distribution. A single minimisation allowed us to determine the population parameters (age and metallicity) and the kinematics (redshift and velocity dispersion). The stellar mass-to-light ratio was derived from the age and metallicity for the considered model. The simultaneous fit of the kinematical and stellar population parameters reduced the degeneracies between age, metallicity and velocity dispersion (Koleva et al. 2008).

We performed the spectral fit using the PEGASE-HR single stellar population models (SSP, Le Borgne et al. 2004), where the observed flux, F_λ, is modelled as follows: (1)The models were built using the Elodie.3.1 (Prugniel & Soubiran 2001; Prugniel et al. 2007) spectral library and the Kroupa (Kroupa 2001) and Salpeter (Salpeter 1955) initial mass functions (IMF). L(v_sys,σ) is a Gaussian function of the systematic velocity, v_sys, and of the velocity dispersion, σ. S(t, [Fe/H] ,λ) is the model for the SSP and depends on age and metallicity. P_n is a polynomial of degree n, which models any residual uncertainty in the flux calibration and extinction correction. In addition, the quasar spectrum might not be perfectly subtracted from the galaxy. To mimic this effect, we included an additive polynomial, Q_m. The order of the additive and multiplicative polynomials is the minimum required to provide an acceptable χ², i.e., in our case (n = 10,p = 1).

We obtained SSP-equivalent ages and metallicities of t ~ 3 Gyr and  [Fe/H]  ~ 0.0 dex respectively. The corresponding rest-frame B-band stellar mass-to-light ratio is using a Kroupa IMF and using a Salpeter IMF. The uncertainties in the age and metallicity were estimated via Monte Carlo simulations and propagated in the error in M_⋆/L_B.

To compute the physical velocity dispersion we subtracted quadratically the instrumental broadening from the measured profile, neglecting the dispersion of the models since they are based on high-resolution templates. The instrumental broadening was measured both from the PSF stars used to carry out the spatial deblending of the spectrum (Eigenbrod et al. 2006) and from the lamp spectra. We obtained the rest-frame physical stellar velocity dispersion of the lensing galaxy: σ_ap = 222  ±  34 km s^-1 in an aperture of 1′′, i.e., 5.7 kpc.

5.2. Numerical integration of the Jeans equations

In this section we model the lensing galaxy using a 3D spherical potential well formed of two components, one for the stellar part of the mass and one for the dark matter halo. We then perform a numerical integration of the Jeans equations in 3D to predict a theoretical velocity dispersion and a total mass for the model. The assumption of spherical symmetry is sufficient for our purpose, as illustrated by the study of the lensed quasar MG 2016+112, where Koopmans & Treu (2002) introduced an anisotropy parameter and showed that it has almost no influence on the inferred mass slope.

The luminous component of the model is a Hernquist profile (Hernquist 1990): (2)where ρ _⋆ (0) is the central density and r _⋆  is a scale radius chosen so that the integrated mass in a cylinder of radius R_eff (effective radius) is equal to half the total stellar mass M _⋆ . The profile has a maximum radius of r _⋆ ,max = 20   r _⋆ .

The dark matter halo is modelled as a generalised Navarro, Frenk & White (NFW) profile (Navarro et al. 1996): (3)where γ_DM is the inner slope of the profile, r_s is the scaling radius and ρ_h(0) is the central mass density. For γ_DM = 1 the model closely follows the standard NFW profile. Its total mass, M_h, is given in the truncation radius r_h,max = 10   r_s.

Following the usual convention, the integrated stellar and dark matter masses are related by the baryonic fraction f_b, (4)Because all integrations in this work were carried out numerically, f_b can easily be computed in any aperture. We chose to compute it in the Einstein radius, which is also where lensing gives the most accurate mass measurement. The velocity dispersion of the stellar component was computed by solving the second moment of the Jeans equation in spherical coordinates (Binney & Merrifield 1998). The velocity dispersion is then (5)where Φ(r) is the total gravitational potential. Equation (5) is solved numerically as follows. First, the density of each mass component (stars and dark halo) is sampled by N particles, using a Monte-Carlo method. Second, the potential is computed using a treecode method (Barnes & Hut 1986). For self-consistency, the density is computed by binning the particles in spherical shells. Then, a velocity is allocated to each particle at a distance r from the galaxy center, following a Gaussian distribution of variance σ _⋆ (r). Finally, the velocity dispersion of the model, σ_Jeans, is computed numerically, by integrating all the particle velocities in an aperture that matches exactly the slit used to carry out the observations.

The advantages of the numerical representation of the lens models are multiple. They allow us (i) to compute velocity dispersions for any combination of density profiles, even non-parametric ones; (ii) to account for the truncation radius of the halo; (iii) to test the dynamical equilibrium of the system by evolving the model with time. Computing the velocity dispersion as well as the total mass in a cylinder along the line of sight of the observer is then straightforward.

In all calcultations we fixed the effective radius to the one observed for the lensing galaxy, i.e., R_eff = 8.44 kpc. The remaining parameter space to explore is composed of the halo scale radius, r_s, the slope of the profile, γ_DM, and the baryonic fraction, f_b, within the Einstein radius.

5.3. Applying the dynamical and stellar population constraints

Using our HST photometry, we measured the total rest-frame B-band galaxy luminosity by converting its total H-band magnitude, m_F160W = 16.20. Using a k-correction of 1.148, and a galactic extinction¹ of E(B − V) = 0.059 (Schlegel et al. 1998), we find L_B = 1.04  ×  10¹¹ L_B, ⊙ . For a given galaxy model we can therefore compute the total stellar mass-to-light ratio, M _⋆ /L_B, as well as the baryonic fraction that we can compare with the observed ones. This requires a choice of IMF. For a Salpeter IMF we measured , while for a Kroupa IMF we measured . In our computation of the baryonic fraction, the total mass is the mass in the Einstein radius, M(<R_E).

We compared the model properties to the data by computing a chi-square, . This includes the measured lens velocity dispersion and the total mass in the Einstein radius, M(<R_E) = (3.16  ±  0.31)  ×  10¹¹   M_⊙. We allowed for a 10% error on the total mass in the Einstein radius to account for the weak dependence of this mass upon the choice of a lens model, which led to (6)

Fig. 6 Radial slope, γDM, of the dark matter halo as a function of the baryonic fraction in the Einstein radius, RE. The colour code represents the value of  for all acceptable models, i.e.,  (see text; Eq. (6)). The observational limits drawn on fb using our spectrum of the lensing galaxy are indicated as shaded vertical regions. The blue region corresponds to a Kroupa IMF and the red regions corresponds to a Salpeter IMF. In this figure, all models have rs = 4  ×  RE = 25 kpc.

Fig. 7 Effect of change in rs, the scale radius of the dark matter halo, on the γDM − fb relation. The contours show the region containing 68% of the models. The blue area has rs = 25 kpc (our adopted choice for the final modelling), the green region has rs = 15 kpc, and the red region is for models with small scale radii, rs = 5 kpc.

Fig. 8 Effet of a change in Reff on the γDM − fb relation. The contours show the region containing 68% of the models. The value adopted in this paper is Reff = 8.44 kpc, as measured from HST images. This corresponds to the area in red. The models in the blue region have Reff = 4.22 kpc, and the region in green has Reff = 12.88 kpc.

We used our numerical integration of the Jeans equations to sample the γ_DM − f_b plane, for a given choice of the scale radius, r_s, of the dark matter halo. Figure 6 shows the value of for 960 000 galaxy models with r_s = 4  ×  R_E = 25 kpc. For a given choice of IMF a precise measurement of the baryonic fraction, f_b, can constrain the dark matter slope γ_DM well, at least in the case of a Kroupa IMF.

Figure 6 is drawn for r_s = 4  ×  R_E = 25 kpc, which gives the best  values. However, models with are also obtained with other values of r_s. We therefore test the influence of a change in r_s by building models for a range of scale radii, r_s = 5,15,30 kpc (i.e., r_s = 0.75 × , 2.25 × , 4.50  ×  R_E). Figure 7 summarises the results and shows no influence of r_s on the γ_DM − f_b relation unless r_s becomes unrealistically small.

The Hernquist profile in our model also has a scale parameter, r _⋆ , which is fixed to the observed value of the effective radius, R_eff = 8.44 kpc. Although we have HST images that allow us to measure R_eff and although our error bars on this parameter are small, we note that studies of lensing galaxies by different authors often find different R_eff values from the same data. This is because different methods are used to subtract the quasar images and to mask (or not) the Einstein ring. We therefore computed models for a range of different R_eff and show the influence on the γ_DM − f_b plane in Fig. 8, for R_eff = 4.22,8.44,12.88 kpc. This reflects an (unlikely high) 50% change in R_eff with respect to the actual measurement and translates into a 10 − 20% shift of the γ_DM − f_b relation. This is comparable to the width of the coloured areas in Fig. 8. Because our true measurement error on R_eff is certainly much better than 50%, we conclude that this parameter does not affect the models much, given the accuracy of the observations.

6. Towards H₀ with HE 0435-1223

The gravlens software (Keeton 2004) was used in combination with the dynamical and stellar population constraints to reproduce the lensing configuration of the quasar images and to attempt converting the time delays into H₀. We adopted a total potential well composed of the main lensing galaxy, the nearby galaxy G22 (see Morgan et al. 2005) plus an external shear.

The main galaxy is composed of a projected Hernquist + cuspy halo model identical to those in Eqs. (1) and (2).

For the Hernquist profile, we fixed the ellipticity and PA of the lensing galaxy to the observed ones.

We then estimated how well each of the models described in Sect. 5.2 reproduces the observed image configuration. In doing this, we allowed only the external shear (γ, θ_γ), the Einstein radius of G22 and H₀ to vary. The potential well of G22 was assumed to lie at its observed position and was modelled as a singular isothermal sphere (SIS). We assumed a conservative value of , following the results of Morgan et al. (2005), who have found R_E(G22)  =  018.

Fig. 9 Distribution of our lens models as a function of the dark matter slope and H0. The colour code gives the value of (see text; Eq. (7)). Only the models with are shown. The results in the left panel are for a Salpeter IMF, and for a Kroupa IMF in the right panel. The lensing χ2 itself is not included in the figure. However, it is systematically lower for galaxies with Kroupa IMFs than for galaxies with Salpeter IMFs.

Fig. 10 Same as Fig. 9, but the colour code now gives the slope of the total (dark+luminous) mass profile, γ′. A logarithmic slope of γ′ = 2 corresponds to an isothermal profile, shown in green.

We show in Fig. 9 the value of H₀ for each lens model as a function of its dark matter slope, γ_DM. The colour code in the figure corresponds to the value of , where the baryonic fraction, f_b ± σ(f_b), in the Einstein radius is now included in the calculation: (7)Including f_b is justified by Fig. 6, showing that different values of f_b select different lens models. We display our results for the two most common IMFs in use in stellar population modelling, the Salpeter and the Kroupa IMFs.

The points define a surface with a clear valley which minimum indicates the best dark matter slopes for each IMF. These are γ_DM(Sal) ~ 1.15 and γ_DM(Kro) ~ 1.54 for the Salpeter and the Kroupa IMFs, respectively. Each model shown in Fig. 9 is also required to display a good lensing chi-square, , after fitting of the quasar image positions with gravlens. The values of are systematically lower for the lensing galaxies with Kroupa IMFs. The present lensing and dynamical work therefore favours a lensing galaxy with a Kroupa IMF, as also found by Cappellari et al. (2006) from 3D spectroscopy of early-type galaxies and by Ferreras et al. (2008) using the SLACS sample of strong lenses (Bolton et al. 2008).

The points in Fig. 9 all have , which prevents us, for now, from giving a value for H₀ given the observational uncertainty on f_b and σ_ap. The measurement errors on both parameters will, however, easily improve with deeper and higher resolution spectroscopy of the lens.

With the current observational constraints we rely on previous work done on the total mass slope, γ′, of lensing galaxies. Koopmans et al. (2009) have measured the probability distribution function of γ′ in the SLACS sample of strong lenses. They have found , with an intrinsic spread of σ(γ′) = 0.20, also confirmed in a more recent study by Auger et al. (2010). The models with the best  values in Fig. 9 correspond to a total slope of γ′ = 2.1  ±  0.1 independent of the choice of an IMF. This is well within the rms limits of Koopmans et al. (2009) and leads to 57  <  H₀  <  71 km s^-1 Mpc^-1 with all models equiprobable within this range.

In our analysis we modelled the environment of the lensing galaxy as a SIS that represents galaxy G22, plus an external shear with a PA of the order of  − 15°, and amplitudes in the range 0.06 < γ < 0.08. However, we did not explicitly account for the fact that HE 0435-1223 lies within a group of galaxies (Momcheva 2009). The unkown convergence, κ, associated with the group leads to an overestimate of H₀ by a factor (1 − κ)^-1, meaning that the range of acceptable H₀ values would decrease further if we have underestimated the convergence caused by the group. Current imaging data suggest, however, that the line of sight is in fact sightly under-dense compared with other lensed systems (Fassnacht et al. 2011). However, the missing convergence is about κ ~ 0.01 − 0.02, leading to at most a 2% change (upwards) in H₀. The above statement should, however, be considered with care because the effect of a group is poorly approximated by a simple convergence term. Explicit modelling of the group halo is likely needed to properly account for the modifications induced on the main lens potential. Deep X-ray and/or optical integral field spectroscopy may turn out to be very useful as well in determining the centroid and mass of the group that contains the galaxy lensing HE 0435-1223.

7. Conclusion

We presented seven years of optical monitoring for the four lensed quasar images of HE 0435-1223. We found that the time delays are better expressed with respect to component B, which is the least affected by stellar microlensing in the lensing galaxy. The formal error bars on the time delays are between 5% and 10% depending on the component, which is remarkable given the very short time delays involved. In addition, the delays are robust against different tests performed on the data, including removal of subsets of data and monte carlo simulations. These tests are possible only with very long light curves, as provided by COSMOGRAIL. Most past lens monitorings have two to three seasons and much coarser temporal sampling. The delays are also independent of the way the microlensing variations are modelled. Given the short length of the delays, additionally improving the error bars will only be possible by increasing the temporal sampling of the curves, i.e., by merging all existing data on HE 0435-1223 taken by different groups over the years.

We introduced a method to convert the time delays into H₀ purely based on external constraints on the radial mass slope of the lens. These constraints come from deep optical spectroscopy of the lens, which allowed us to measure its velocity dispersion and its stellar mass-to-light ratio. The present paper describes our approach, but the current observations of the lensing galaxy so far lead to a broad range of values for H₀.

Our methodology complement that presented by Suyu et al. (2009, 2010) well. Combining the two in the future will require follow-up observations such as (i) deep HST imaging to map the Einstein ring with high signal-to-noise; (ii) deep high resolution spectroscopy of the lens over a broad spectral range to narrow down the uncertainty on its velocity dispersion and stellar mass; (iii) integral field spectroscopy to measure all redshifts within 30 − 60′′ around the lens; and (iv) X-ray imaging to pinpoint massive groups along the line of sight.

If all these observations are taken for HE 0435-1223 as well as for a few lenses with well measured time delays, the cost in terms of follow-up observations will still remain very modest compared with other existing methods to measure H₀, such as Cepheids and supernovae.

Acknowledgments

We are grateful to all the observers who contributed to the data acquisition at the Euler and Mercator telescopes as well as at Maidanak Observatory. COSMOGRAIL is financially supported by the Swiss National Science Foundation (SNSF). This work is also supported by the Belgian Federal Science Policy (BELSPO) in the framework of the PRODEX Experiment Arrangement C-90312. V.C. thanks the Belgian National Fund for Scientific Research (FNRS). D.S. acknowledges a fellowship from the Alexander von Humboldt Foundation. M.K. has been supported by the Programa Nacional de Astronomía y Astrofísica of the Spanish Ministry of Science and Innovation under grant AYA2007-67752-C03-01 and DO02-85/2008 from Bulgarian Scientific Research Fund.

References

All Tables

All Figures

	Fig. 1 Part of the field of view of the 1.2 m Swiss Euler telescope, with HE 0435-1223 visible in the centre. The four PSF stars used for deconvolution purposes and the two reference stars used to carry out the flux calibration are indicated.
In the text

	Fig. 2 Result from the simultaneous deconvolution of the ground-based frames. G is the lensing galaxy and G22 (Morgan et al. 2005) is its closest neighbour on the plane of the sky. The grey scale in the deconvolved image is set to display all light level above 3  ×  σsky. The FWHM resolution of the deconvolved image is 034.
In the text

	Fig. 3 R-band light curves of the four lensed images of HE 0435-1223 from December 2003 to April 2010. The magnitudes are given in relative units as a function of the Heliocentric Julian Day (HJD), along with their total 1σ error bars. These light curves are available in tabular form at the CDS.
In the text

	Fig. 4 Left: combination of the four original HST/NIC2 F160W frames of HE 0435-1223. The field of view is 9  ×  9 arcsec. Middle: deconvolved image, where the lensing galaxy is modelled as a de Vaucouleurs profile (see text). The nearest galaxy on the plane of the sky, G22, is also indicated. Right: residual map in units of the noise. The colour scale ranges from  − 4σ (white) to +4σ (black).
In the text

Fig. 5 Light curves obtained with all four telescopes and shifted by time delays of ΔtBA = 8.4 days, ΔtBC = 7.8 days and ΔtBD =  −6.5 days. The relative microlensing representations applied on curves A, C and D are shown as continuous curves with respect to the dashed blue line (see text). A fifth-order polynomial was used over the seven seasons to model microlensing on the quasar images C and D, while seven independent third-order polynomials were used for image A. The lower panels show residuals obtained by subtracting a simultaneous spline fit (grey) from the light curves.

In the text

Fig. 6 Radial slope, γDM, of the dark matter halo as a function of the baryonic fraction in the Einstein radius, RE. The colour code represents the value of  for all acceptable models, i.e.,  (see text; Eq. (6)). The observational limits drawn on fb using our spectrum of the lensing galaxy are indicated as shaded vertical regions. The blue region corresponds to a Kroupa IMF and the red regions corresponds to a Salpeter IMF. In this figure, all models have rs = 4  ×  RE = 25 kpc.

In the text

	Fig. 7 Effect of change in rs, the scale radius of the dark matter halo, on the γDM − fb relation. The contours show the region containing 68% of the models. The blue area has rs = 25 kpc (our adopted choice for the final modelling), the green region has rs = 15 kpc, and the red region is for models with small scale radii, rs = 5 kpc.
In the text

	Fig. 8 Effet of a change in Reff on the γDM − fb relation. The contours show the region containing 68% of the models. The value adopted in this paper is Reff = 8.44 kpc, as measured from HST images. This corresponds to the area in red. The models in the blue region have Reff = 4.22 kpc, and the region in green has Reff = 12.88 kpc.
In the text

Fig. 9 Distribution of our lens models as a function of the dark matter slope and H0. The colour code gives the value of (see text; Eq. (7)). Only the models with are shown. The results in the left panel are for a Salpeter IMF, and for a Kroupa IMF in the right panel. The lensing χ2 itself is not included in the figure. However, it is systematically lower for galaxies with Kroupa IMFs than for galaxies with Salpeter IMFs.

In the text

	Fig. 10 Same as Fig. 9, but the colour code now gives the slope of the total (dark+luminous) mass profile, γ′. A logarithmic slope of γ′ = 2 corresponds to an isothermal profile, shown in green.
In the text