© ESO, 2014

1. Introduction

Follow-up campaigns of gravitationally lensed quasars are unique astrophysical tools because they reveal properties of the accretion flow in distant AGN, the structure of lensing galaxies, and the physics of the Universe as a whole (e.g., Schneider et al. 1992, 2006). To model a lens mass distribution, it is necessary to derive accurate constraints on relative positions of the lensed quasar images and the main lensing galaxy, as well as information on the structure of the galaxy’s light (e.g., Lehár et al. 2000; Chantry et al. 2010; Sluse et al. 2012). Other important constraints are the flux ratios between quasar images. These macrolens (intrinsic) flux ratios are difficult to measure because the quasar light may be affected by differential microlensing magnification and dust extinction within the intervening galaxy (e.g., Wambsganss 1990; Falco et al. 1999). Fortunately, one can use flux ratios of emission line cores to fit macrolens ratios and extinction properties (e.g., Wucknitz et al. 2003; Mediavilla et al. 2009; Motta et al. 2012). The typical Einstein ring of stars (in the source plane) has a size similar to that of the accretion disk, which is the reason why microlenses (stars) strongly affect the X-ray, UV, and optical continuum emission (e.g., Yonehara et al. 1998; Morgan et al. 2012; Mosquera et al. 2013). However, emission lines in the quasar optical spectra would have central features (cores) that are likely produced in regions much more extended than the typical Einstein ring of intervening stars, that is, emission from the narrow-line region (NLR) and the outer parts of the broad-line region (BLR). Thus, emission line cores are presumably unaffected by microlensing.

After macrolens-extinction parameters are inferred from line-core flux ratios, it is possible to pave the way for mass-modeling and continuum-microlensing studies (after correcting continuum flux ratios for macrolens and extinction effects; Wucknitz et al. 2003). Microlensing spectra (microlensing magnification ratios of the continuum at many wavelengths) allow astronomers to discuss the accretion disk structure (e.g., Mediavilla et al. 2011; Motta et al. 2012), while microlensing (extrinsic) variability in broadband light-curves of the quasar images often complicate the measuring of time delays between intrinsic variations. However, intensive microlensing simulations in the time-domain can separate both kinds of variability and yield accurate time delays (and sizes of variable sources; e.g., Hainline et al. 2013). Optical spectroscopy of a lens system also provides the redshift of the lensing galaxy, which is a basic quantity to relate lensing mass/cosmology and time delays (e.g., Eigenbrod et al. 2006, 2007).

A main goal of our Liverpool quasar lens monitoring (LQLM) project is to build a rich time-domain database for several lens systems (Goicoechea & Shalyapin 2012). This robotic optical monitoring is complemented with other relevant observations, such as deep optical imaging and spectroscopy of the quasar fields (e.g., Shalyapin & Goicoechea 2013). Here, we focus on the imaging and spectroscopy of the doubly imaged quasar SDSS J1339+1310 (Inada et al. 2009). The outline of the paper is the following: Sect. 2 presents new astro-photometric and spectrophotometric data obtained with the 2.5 m Nordic Optical Telescope (NOT) and the 10.4 m Gran Telescopio Canarias (GTC). Section 3 is devoted to the analysis and interpretation of data, determining the lens redshift, extracting macrolens-extinction solutions and the microlensing signal, and discussing possible lens models. Our results and conclusions are summarized in Sect. 4. Some concluding remarks are also incorporated into this last section.

We use a standard cosmology (h = 0.71, Ω_M = 0.27 and Ω_Λ = 0.73; e.g., Spergel et al. 2003) throughout this paper.

2. Observations and data reduction

The lens system SDSS J1339+1310 consists of two quasar images, A and B, separated by 1.̋69, and a lensing galaxy (G) placed between the images of the lensed quasar (Inada et al. 2009). Inada et al. (2009) also reported a redshift for the source quasar: z_s = 2.24, which is equal to the average redshift of the two images in the latest Sloan Digital Sky Survey (SDSS) data release DR10¹. This average redshift basically comes from the very prominent Lyα emission, while the C iv (λ1549), C iii] (λ1909) and Mg ii (λ2798) emission lines have z_s(A) = 2.231 and z_s(B) = 2.234. The spectroscopic redshift of the lensing galaxy is unknown. There is only a rough photometric estimate of z_l ~ 0.4 (Inada et al.). Here, we present deep NOT I-band imaging (Sect. 2.1) and GTC optical spectroscopy (Sect. 2.2) that are used to study the physical properties of SDSS J1339+1310 (see Sect. 3).

2.1. Deep I-band imaging at the NOT

Using a 480 s I-band exposure with the 2.2 m UH88 Telescope under about 0.̋8 seeing, GALFIT software (Peng et al. 2002) allowed Inada et al. to separate the two quasar images from the lensing galaxy, and to fit the galaxy light distribution to the convolution of a 2D Sérsic profile and the point spread function (PSF). To improve the astro-photometric data for the lens system, I-band frames of SDSS J1339+1310 were taken on 13 June 2010 at the NOT with ALFOSC. This camera provides a pixel scale of and a field of view of ~. We obtained twelve 300 s exposures, which produced a deep combined frame. The total exposure time in the I band is 1 h = 12 × 300 s, which is ten times longer than the previous one in Inada et al. (correcting for the NOT/UH88 diameter ratio). Our NOT combined (stacked) frame is also characterised by a seeing of only 0.̋64.

Fig. 1 Stacked frame of SDSS J1339+1310 with the NOT in the I band. This corresponds to exposures taken with ALFOSC on 13 June 2010 (total exposure time = 3.6 ks, seeing = 0.̋64). In addition to the original signal (left panel), we also consider a model that includes four additive ingredients: two PSFs (quasar images), the convolution of a 2D de Vaucouleurs profile and the PSF (lensing galaxy), and a constant term (sky background). We display the residual fluxes after subtracting only the quasar images and the background (middle panel), as well as the residual signal after subtracting the global model (right panel). The fit to the model is made with IMFITFITS software (see main text).

The combined frame incorporates the lens system and bright stars, which enabled us to finely sample the stellar PSF and perform an accurate PSF-fitting astro-photometry of SDSS J1339+1310. We used two different software packages, GALFIT and IMFITFITS (McLeod et al. 1998), to fit the light distribution of the lens system. This was modelled as two PSFs (quasar images) plus a standard galaxy profile convolved with the PSF. In the left panel of Fig. 1, we display the region of interest where both quasar images are clearly resolved. The middle panel shows the residual signal after subtracting the two images and a constant sky background from IMFITFITS, and the lensing galaxy is now evident. In the right panel of Fig. 1, we also show the final residual light after subtracting the global model (quasar images + lensing galaxy + sky background) from IMFITFITS. The 2D galaxy profile is assumed to follow a de Vaucouleurs law, the reduced chi-square value is 1.08, and the right panel does not contain substantial residues in its central zone.

A detailed comparison between previous astro-photometric parameters (see Tables 4 and 5 of Inada et al.) and new results from the NOT observations is shown in Table 1. Inada et al. used a Sérsic profile for the galaxy, while we obtained solutions for two standard galaxy profiles: Sérsic and de Vaucouleurs. In Table 1, the positions of B and G are specified as relative to the image A, and positive values of x and y are westward and northward of image A, respectively. The single-epoch flux ratio B/A = 10^{0.4 × (m_A − m_B)} could not represent the I-band ratio at the same emission time, since each image is most likely variable on a timescale equal to the time delay between A and B. However, considering the central wavelength of the I filter at 7970 Å (rest-frame wavelength λ = 2460 Å) and a time delay below 50 d (rest-frame variability timescale Δτ ≤ 15 d; see our results in Sect. 3.3), we can estimate the expected flux variation of an image (in mag) from the structure functions in Fig. 11 of MacLeod et al. (2012). These give an upper limit of 0.04 mag at λ ~ 2000−3000 Å and Δτ ≤ 15 d, which is taken as a simultaneity error, and thus, added in quadrature to the photometric error in m_A − m_B. In Table 1, the light distribution of the lensing galaxy is described by four structure parameters: r_eff, n, e and θ_e are the effective radius, the concentration index, the ellipticity, and the position angle, respectively.

Table 1 lists some differences between the results in the second and third columns from GALFIT (Sérsic profile). There are significant discrepancies between the UH88 and NOT astrometric solutions, and the new astrometric uncertainties are half the previous ones. In addition, the new Sérsic concentration index is consistent with n = 4, that is, a de Vaucouleurs law, and there is a noticeable change in the orientation of the lensing galaxy (θ_e). The relatively short exposure and the more extended (and elongated) PSF with the UH88 Telescope probably led to an incorrect value for θ_e. The third and fourth columns in Table 1 prove that both profiles (Sérsic and de Vaucouleurs) work in a similar way, and the fifth column includes the solution from IMFITFITS using a de Vaucouleurs law. We adopt this last solution for subsequent analysis. IMFITFITS minimizes the sum of squared differences between observed values and fitted values provided by a model, whereas GALFIT minimizes the sum of squared residuals normalized by the expected errors. Although the outputs of both packages are similar, we prefer to assign the same weight to all squared residuals.

2.2. GTC optical spectroscopy

We also used OSIRIS (on the GTC) to obtain spectra of the sources in SDSS J1339+1310. These observations are part of an ongoing project with the new 10.4 m telescope (Shalyapin & Goicoechea 2013). We took a single 2700 s spectroscopic exposure of SDSS J1339+1310 on 13 April 2013, using the grism R500R (wavelength range and resolving power of 4800−10 000 Å and ~500) and a long slit. The slit width and position angle were 1.̋0 (4 pixel) and −31.̊7, which means that it was oriented along the line joining the B image and the galaxy, but also including an appreciable amount of light from the off-axis image A. In Fig. 2, we show the positions of the three components A, B and G on the slit. The three associated circles have a diameter of 1.̋4, which equals the seeing at 7000 Å. In our observations, there were two main reasons for a tiny differential atmospheric refraction. First, the mean airmass was quite small (1.050). Second, the slit orientation was close to the parallactic angle (see Fig. 2). The vertical shift of components (due to differential atmospheric refraction) throughout the entire spectral range does not exceed 0.̋05, or equivalently, 0.2 pixel.

After a standard data reduction with IRAF², we obtain a wavelength-calibrated spectrum for each slice along the slit. From spectral calibration lamps, it is also possible to infer the spectral resolution for a slit width of 1.̋0 at a central wavelength. We calculated FWHM = 15.6 Å for the Ne i line at 7174 Å (resolving power of 460). We then used additional observations of the standard star GD140 to calibrate the flux scale, and thus, derived wavelength- and flux-calibrated spectra. There are several approaches to extract individual spectra of quasar images and galaxies in lens systems on a long slit (e.g., Fassnacht & Cohen 1998; Ofek et al. 2006; Sluse et al. 2007). In particular, Sluse et al. (2007) extracted spectra of three close quasar images along a slit by fitting three 1D PSF profiles in the spatial direction for each wavelength bin. Our extraction method relies on the scheme of Sluse et al., but it incorporates some procedures to accurately account for off-axis quasar images and galaxies. First, we used the NOT astro-photometric solution in the last column of Table 1 to make an ideal 2D light model with two points (A and B) and a de Vaucouleurs profile (G). In a second step, the ideal model is convolved with a 2D PSF. Tests with the star GD140 favour a Moffat over a Gaussian PSF, therefore we adopted a 2D Moffat PSF. Third, we set the slit mask on the model frame (see the white rectangle in Fig. 2), focused on the 2D light distribution for this mask, and integrated it along the vertical axis. In principle, the 1D model depends on the position of the image B, the two structural parameters (slope and width) of the circular Moffat PSF, and the total fluxes of A, B, and G. However, the slit was centred on B, which means that we deal with an on-axis B image. In addition, Jahnke et al. (2007) showed that the slope and width of the Moffat profile are directly correlated when fitting a given shape, and thus, we set a slope β = 2.5 and allowed the width to vary. The final 1D model contains five free parameters: the horizontal position of B, the width of the PSF, and the fluxes of the galaxy and quasar images.

Fig. 2 Positions of the three components A, B and G of SDSS J1339+1310 on the slit for GTC-OSIRIS. The circles represent seeing disks (diameter = seeing), whereas E, N and Z indicate east, north, and zenith.

For the transmission of the slit, we used an empirical approach. We focused on the through-slit acquisition frame to study the homogeneous sky background in columns far from the region illuminated by the lens system. Instead of a constant brightness along these columns, the vertical brightness distribution is well represented by an exp [− (Y/σ)⁴] law with σ = 2.1 pix (the slit axis corresponds to Y = 0). This describes the slit transmission significantly better that the simple step function. In the top panel of Fig. 3, we show the 1D brightness profile at the central wavelength (~7050 Å) in the spatial region including the lens system for the wavelength- and flux-calibrated spectral frame. The filled circles are the observed values, while the solid lines trace the 1D model solution (black) and its individual contributions A, B and G. This solution is related to the 2D light distribution (before integrating across the slit) that is drawn in the bottom panel of Fig. 3. The 2D distribution has a spatial resolution of 0.1 pix, and it displays the smoothing effect of the PSF and the dip in transmission near the slit edges.

Fig. 3 1D brightness profile (along the slit) at ~7050 Å. Top panel: observational data (filled circles) and the 1D model solution (black line). The contribution of each source is also detailed. Bottom panel: 2D brightness map for the 1D model solution. The two horizontal dashed lines are the slit edges.

Before presenting our final spectra, we discuss the influence of possible systematic effects on them. For example, the slit could have a slightly different position angle from what we used, or the relative positions of A and G (with respect to B) could be different from their best-fit values (see the last column of Table 1). The uncertainty of 0.̊1 in the slit position angle is consistent with small displacements of A and G (≤0.01 pix) in the coordinate system with the origin at image B and the horizontal (X) axis coinciding with that of the slit (see Fig. 2). Additional displacements may be associated with the astrometric uncertainties in Table 1, and both systematic effects have a weak influence on the flux estimates for A, B, and G. We also considered that the vertical position of B could be slightly displaced from the slit axis. To check deviations of fluxes with respect to their solutions for an on-axis B image (our final model with Y_B = 0), firstly, we analysed the through-slit acquisition frame and placed a constraint on the positive vertical offset of B. Taking this constraint Y_B< 0.1 pix into account, the galaxy and B image spectra practically do not change, and the deviations of flux for the A image do not exceed 2%. Chromatic vertical offsets in the interval ±0.1 pix are also expected as a result of differential atmospheric refraction (see above), and thus, changes in the flux of A are below a few percent. A maximum systematic error of a few percent is propagated onto B/A.

Fig. 4 GTC-OSIRIS spectra of SDSS J1339+1310. These are derived from a 2.7 ks exposure on 13 April 2013. The red, blue, and green lines trace the fluxes of A, B, and G (see Fig. 2) as a function of observed wavelength. Data from frames in the r (circles) and I (square) bands are also shown for comparison purposes. The dashed lines describe the responses of the two filters, i.e., the LT r-Sloan and the NOT #12 passbands. Vertical dotted lines are associated with emission lines at zs = 2.24, while grey highlighted regions indicate atmospheric artefacts (see main text).

In a first iteration, we fitted the five free parameters of the 1D model to the 1D data for each ~5 Å wide wavelength bin. The two position-structure parameters were then fitted to smooth polynomial functions of the observed wavelength. In a second iteration, the position of B and the width of the PSF for each wavelength bin were estimated through the corresponding polynomial function, leaving only the three fluxes as free parameters. The final observer-frame spectra of A, B and G are plotted in Fig. 4. These spectra are freely available at the GLENDAMA website³. In Fig. 4, the circles and the square show fluxes from our photometric monitoring programme with the 2.0 m Liverpool Telescope (LT) in the r band (Goicoechea & Shalyapin 2012) and the NOT stacked frame in the left panel of Fig. 1. We note that the r-band fluxes correspond to LT frames taken on the same night that we performed the spectroscopic observations with the GTC. The NOT frames were obtained about three years before the GTC exposure, and hence we only include the I-band flux of the non-variable source. The spectral trend of G agrees reasonably well with independent two-colour photometry. We also compared the direct r-band flux ratio B/A and B/A from properly averaged spectrophotometric data (taking into account the blue and red lines, as well as the sensitivity function of the LT r-Sloan filter), and obtained a small deviation of ~3%.

The most prominent features of the quasar spectra in Fig. 4 are the C iv (λ1549), C iii] (λ1909) and Mg ii (λ2798) emission lines, which are observed in both images in the wavelength range 4850−9250 Å. Using z_s = 2.24 (see above), we also draw three vertical dotted lines as reference wavelengths. Each of the six emission lines allows us to estimate the source redshift. Thus, we obtain average values z_s(A) = 2.230 and z_s(B) = 2.231. Both results agree well with each other and with the z_s values for the C iv, C iii] and Mg ii lines in the SDSS DR10 spectra (see the first paragraph in this section).

3. Properties of the lens system from the new observations

3.1. Spectroscopic redshift of the main lensing galaxy

In this subsection, we used the spectrum of G in Fig. 4 to estimate z_l. The software RVSAO/XCSAO (Kurtz & Mink 1998) is an IRAF package that allowed us to obtain the redshift of the main lensing galaxy using the cross-correlation method. Going into details about the procedure, we compared the observed spectrum at λ ≥ 6100 Å with the template spectrum for elliptical galaxies (Kinney et al. 1996). To perform this comparison, we masked some wavelength intervals related to sky-absorption bands and emission lines. Because the RVSAO/XCSAO package only works well for radial velocities below 100 000 km s^-1 (z< 0.4), the observed spectrum was initially blue-shifted by V₀ = 131 409 km s^-1 (rest-frame of a galaxy at z₀ = 0.6), and then cross-correlated with the template. Instead of a relatively high velocity to be subtracted from V₀ (if z_l ~ 0.4), the cross-correlation indicated the existence of a small additional radial motion of 1363 km s^-1. Hence, we infer z_l = 0.609 with formal error <0.001 (84 km s^-1) and r_cc = 3.81, where r_cc measures the ratio between the cross-correlation peak and the signal expected from noise and template mismatch (Tonry & Davis 1979). This redshift determination seems accurate and reliable (r_cc> 3).

Fig. 5 Spectrum of the main lensing galaxy in SDSS J1339+1310. The observed spectrum (see Fig. 4) is smoothed by a 5-point filter and represented by a green line, while the template spectrum for an elliptical galaxy is properly red-shifted (zl = 0.609) and drawn with a magenta line. Grey highlighted regions and vertical dashed lines are associated with sky absorption/emission and galaxy absorption lines, respectively (see main text).

In Fig. 5 we show a smoothed galaxy spectrum (five-point smoothing) together with the wavelength-shifted template (using z_l = 0.609). The two spectra in the 6100−9250 Å region agree well, which we used to match both trends (see above). Although the cross-correlation in the previous paragraph avoids the noisy spectral region at λ< 6100 Å, we also checked the output from a comparison over a wider wavelength range (5100−9250 Å). We obtained the same best value for z_l, but the new uncertainty is larger. The noisy behaviour at observed wavelengths λ< 6100 Å is due to two reasons: in this spectral range, the galaxy is fainter and the random residues (after subtracting the bright quasar images) are higher. To account for systematics (the choice of the spectral region for cross-correlation and others), we finally considered a velocity error of about twice the formal one (using the 6100−9250 Å region), that is, 150 km s^-1. This leads to z_l = 0.609 ± 0.001. In Fig. 5, we also display several absorption lines in the spectrum of G (see the vertical dashed lines). Main features in the spectrum of an elliptical galaxy are the Ca ii HK doublet and break, the Hβ line, and the Mg i b triplet.

3.2. Separating dust extinction and microlensing through spectroscopic data

3.2.1. Dust extinction and macrolens flux ratio

The GTC spectrum of each quasar image contains three prominent emission lines at different wavelengths: we detect the C iv, C iii] and Mg ii pairs of lines presented in Sect. 2.2, and each of these three emissions is most likely generated in spatially different regions (e.g., Marziani et al. 2010; Sluse et al. 2011; Motta et al. 2012). According to Motta et al. (2012), the broad wings of a given emission line could be produced in a compact region, whereas its core is probably dominated by photons arising from the NLR and the outer parts of the BLR. Thus, large sources would yield the line cores, which would be unaffected by microlensing. We estimated flux ratios for emission line cores following the approach of Motta et al.⁴. To illustrate our procedure, we focused on the C iii] lines in the optical spectra of A and B. To obtain each emission line profile, we subtracted a quadratic interpolation of the continuum, using two continuum zones close to the line: a zone to its left (5800−5950 Å) and a zone to its right (6350−6500 Å). In Fig. 6 we show the C iii] profiles for A and B. For the A image, the signal strength is enough to fairly appreciate the expected C iii] asymmetry, with a blue-wing excess due to Al iii (λ1857) and Si iii] (λ1892) emissions (e.g., Sluse et al. 2011). The line flux of A is significantly higher than the line flux of B, but the flux ratio B/A is higher than 1 in both continuum zones around the line profiles (see Fig. 4). After continuum removal, we integrated the profiles over a 30 Å interval centred on the peak of A (strongest signal), and obtained the flux ratio for the C iii] core: (B/A)_core = 0.27.

Fig. 6 C iii] emission line profiles. For each quasar image, we take the continuum in the intervals 5800−5950 and 6350−6500 Å, fit it to a quadratic law, and then remove the fitted quadratic trend. The solid line traces the line profile for the image A. The dashed line describes the line profile for image B.

A key issue is to estimate uncertainties of line-core flux ratios. For the C iii] emission, we considered all sources of error: systematics related to the spectral extraction and the analysis of the C iii] lines, as well as the signal noise. The spectral extraction only causes an error smaller than a few percent (see Sect. 2.2), which is negligible compared with the overall contribution of other sources of error (~10%). For the C iii]-line analysis, we tried several different combinations of continuum intervals, continuum shapes (linear or quadratic), and integration methods (direct sum or Gaussian fit), and compared the results. The C iv and Mg ii emission lines were observed close to the spectral edges, and we simply took linear laws to link left-hand and right-hand continuum regions. However, we checked possible systematics associated with the continuum zones and integration methods we used to derive flux ratios, therefore our results for the C iv and Mg ii cores are also reliable. The 1σ measurements of the line-core flux ratios are provided in Table 2. In this table, we also add the magnitude differences Δm_core = (m_A − m_B)_core = 2.5log (B/A)_core and their corresponding errors.

Even considering that the emission line cores are not affected by microlensing effects and do not vary on timescales below 50 d (see time-delay predictions in Sect. 3.3), we may observe a chromatic behaviour of the line-core flux ratios and magnitude differences (see Table 2). This is very likely caused by differential dust extinction in the lensing galaxy at z_l = 0.609 (e.g., Motta et al. 2012, and references therein). The light rays associated with images A and B are characterised by two different impact parameters to the intervening galaxy (see Table 1 and Fig. 2), and thus, conceivably cross two galactic regions with different dust properties. The general formalism for differential extinction in a pair of lensed images was introduced by several authors (e.g., Falco et al. 1999; Wucknitz et al. 2003; Elíasdóttir et al. 2006), and this is easily applicable to analyse our line-core magnitude differences in Table 2. We assumed a relationship $\mathrm{\Delta }\mathit{m}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\mathrm{\Delta }{\mathit{m}}_{\mathrm{AB}}\mathrm{+}\mathrm{\Delta }{\mathit{A}}_{\mathrm{AB}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathit{,}$(1)where Δm_AB = 2.5log (M_BA) is a logarithmic version of the (achromatic) macrolens flux ratio M_BA = M_B/M_A, and ΔA_AB(λ) = A_A(λ) − A_B(λ) is the amount of differential extinction at the wavelength λ in the observer’s frame. We also used two dust extinction laws. First, A(λ) = E(B − V) [a(x)R(V) + b(x)], where E(B − V) = A(B) − A(V) is the colour excess, R(V) = A(V) /E(B − V) is the ratio of total-to-selective extinction in the V band, and a(x) and b(x) are known functions of x = (1 + z_l) /λ in μm^-1. This is the standard extinction in the Milky Way (Cardelli et al. 1989). The second variant was a power law A(λ) = A(V) [xλ_V] ^α, with λ_V = 0.55 μm and α = 1 − 2. This can account for extinction in galaxies with types of dust other than those in our own Galaxy (Elíasdóttir et al. 2006).

Fig. 7 Extinction curve for SDSS J1339+1310. Filled circles with error bars represent the line-core magnitude differences in Table 2. The red lines describe the best fits, and the light grey areas are the 1σ bands. Top panel: galactic extinction. Bottom panel: quadratic extinction (see main text).

For the Galactic extinction, Eq. (1) can be rewritten as $\mathrm{\Delta }\mathit{m}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathit{p}}_{\mathrm{0}}\mathrm{+}{\mathit{p}}_{\mathit{V}}\mathrm{\left[}\mathit{a}\mathrm{\right(}\mathit{x}\mathrm{)}\mathrm{-}\mathit{b}\mathrm{(}\mathit{x}\mathrm{\left)}\mathrm{\right]}\mathrm{+}{\mathit{p}}_{\mathit{B}}\mathit{b}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}$(2)where p₀ = Δm_AB, while p_V = A_A(V) − A_B(V) and p_B = A_A(B) − A_B(B) are the differential visual and blue extinctions. It is straightforward to derive Eq. (2) from Eqs. (2)−(4) of Wucknitz et al. (2003). For the power-law extinction, Eq. (1) has a simpler form $\mathrm{\Delta }\mathit{m}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathit{p}}_{\mathrm{0}}\mathrm{+}{\mathit{p}}_{\mathit{V}}\mathrm{[}\mathit{x}{\mathit{\lambda }}_{\mathit{V}}{\mathrm{]}}^{\mathit{\alpha }}\mathit{.}$(3)Equations (2) and (3) are valid for a general case in which dust along each sight line has different properties (e.g., McGough et al. 2005). We have three free parameters using the Galactic extinction, and two free parameters with the power-law extinction for a given value of α.

The red line in the top panel of Fig. 7 displays the best fit to a Milky Way-like extinction: p₀ = −0.41, p_V = −1.14 and p_B = −1.08 mag. From this best fit with χ² ~ 0, we do not infer an unique extinction scenario, but many possible scenarios that are indistinguishable from each other (e.g., Wucknitz et al. 2003; McGough et al. 2005; Elíasdóttir et al. 2006). For example, two physical solutions are (1) E_A(B − V) = 1.00, R_A(V) = 2.00, E_B(B − V) = 0.94, R_B(V) = 3.34, and (2) E_A(B − V) = 0.60, R_A(V) = 3.00, E_B(B − V) = 0.54, R_B(V) = 5.44. Provided that p₀ = −0.41 mag, these two solutions perfectly fit the three line-core magnitude differences (filled circles in Fig. 7). We also followed a standard approach to generate synthetic sets of magnitude differences and determine errors (e.g., Press et al. 2007). The observational noise was assumed to be normally distributed. Therefore, one obtains a synthetic set by adding a random quantity to each observed difference. This random quantity is a realization of a normal distribution around zero, with a standard deviation equal to the observed error. We produced 1000 synthetic sets and obtained 1000 trios (p₀, p_V, p_B) and the corresponding 68% (1σ) confidence regions. The 68% confidence region for the three parameters jointly was inferred from the trios giving rise to fits to the real data with Δχ² ≤ 3.53. These solutions generate the 1σ band in the top panel of Fig. 7.

Regarding the power law model, α = 1 leads to a poor fit with χ² = 1.55, where the number of degrees of feedom (d.o.f.) is 1. However, for α = 2, the best fit is acceptable: χ² = 0.93, d.o.f. = 1 (see the bottom panel of Fig. 7). We then used the synthetic sets of magnitude differences to obtain 1000 pairs (p₀, p_V). The 68% confidence region for the two parameters jointly is defined by the pairs producing Δχ² ≤ 2.30 fits to the real data, and the 1σ band is shown in the bottom panel of Fig. 7. The best-fit parameters and their 1σ uncertainties from our Monte Carlo simulations are summarised in Table 3.

For the Galactic extinction, the best-fit value of the macrolens flux ratio is 0.69, but allowed ratios range from 0.23 to ~2. We remark that the macrolens magnitude difference Δm_AB is equal to Δm(x = 0), and the 1σ band is particularly broad at x ≤ 1.5 μm^-1 (see the top panel of Fig. 7). Thus, near-infrared (NIR) spectroscopy is required to better constrain the model at x ≤ 1.5 μm^-1, and find a narrower interval of M_B/M_A values. For the quadratic extinction, we obtain a narrow 1σ interval for M_BA (<20% uncertainty in Table 3), favouring the lowest values from Galactic extinction and values close to the observed ratio for the Mg ii line core.

3.2.2. Microlensing factory

Wucknitz et al. (2003) presented a novel and comprehensive analysis of spectroscopic data of the double QSO HE 0512−3329. Assuming that the emission line magnitude differences are only due to an achromatic macrolens effect and some amounts of differential extinction, Wucknitz et al. fitted the involved parameters, properly corrected the continuum magnitude differences for differential macrolens magnification and dust extinction, and found evidence for chromatic microlensing of the continuum emission. This idea of extracting macrolens/extinction parameters from data of large sources (BLR/NLR), and to use these parameters to correct data of compact continuum sources, is being actively exploited (e.g., Mediavilla et al. 2005, 2009, 2011; Jiménez-Vicente et al. 2012; Motta et al. 2012). To our knowledge, only one anomalous behaviour has been reported (Goicoechea et al. 2005): spectroscopic data of Q0957+561 show evidence in favour of the existence of compact dusty clouds in the lensing galaxy, which affect continuum sources and some substructures of the BLR/NLR. However, the global emission-line region of the double quasar does not experience differential extinction, since the emission line flux ratios basically coincide with the radio flux ratio.

Here, we adopted the perspective of Wucknitz et al., and used the extinction curve fit results to discuss the microlensing magnification ratio of the continuum. First of all, the continuum flux ratio is given by $\mathrm{(}\mathit{B}\mathit{/}\mathit{A}{\mathrm{)}}_{\mathrm{cont}}\mathrm{=}{\mathit{M}}_{\mathrm{BA}}{\mathit{\mu }}_{\mathrm{BA}}{\mathit{ϵ}}_{\mathrm{BA}}\mathit{,}$(4)where μ_BA and ϵ_BA = 10^0.4ΔA_AB are the microlensing magnification and dust extinction ratios. From Eq. (4), we inferred the microlensing signal as μ_BA = (B/A)_cont/ [M_BAϵ_BA]. Because both (Galactic and quadratic) extinction models can account for the line-core flux ratios, we considered their best fits and 1σ bands in Fig. 7. These varied solutions allowed us to verify the reliability of the microlensing signal reconstruction. The continuum flux ratio is depicted in the top panel of Fig. 8. The original signal (B/A)_cont was then corrected with the solutions for Galactic (middle panel) and quadratic (bottom panel) extinctions. The microlensing signal μ_BA does not depend strongly on the extinction model and macrolens/extinction parameters. This is directly related to the wavelength regime we are probing, which coincides with the spectral range for the observed emission lines and where the macrolens/extinction solutions are similar (see the 1σ bands at 1.8 ≤ x ≤ 3.2 in Fig. 7). In addition, the microlensing spectra suggest that B is amplified (relative to A) by a factor of about 3−5, with larger amplifications at shorter wavelengths. In Eq. (4), we implicitly assume that the source spectrum does not change significantly on timescales <50 d (see the upper limit on the time delay predictions in Sect. 3.3; e.g., Wucknitz et al. 2003). Our final results support this initial hypothesis, because intrinsic flux variations by a factor of ~4 (~1.5 mag) over tens of days are not observed in light curves of normal quasars.

Fig. 8 Continuum flux ratio and microlensing signal. Top panel: continuum flux ratio. Middle panel: microlensing magnification ratio using the best fit (red circles) and the 1σ band (light grey region) for Galactic extinction. Bottom panel: microlensing magnification ratio using the best fit (red circles) and the 1σ band (light grey region) for quadratic extinction. In the middle and bottom panels, the microlensing signal from the 1σ bands is linearly interpolated in the spectral regions where the C iv, C iii] and Mg ii emission lines appear. We also show the zero-level for microlensing (dashed lines; μBA = 1) and fits to a simple microlensing model (solid lines; see main text).

To interpret the microlensing spectra in Fig. 8, we considered a face-on disk accounting for the continuum sources at different wavelengths. The multiwavelength surface brightness of the disk is I(λ,r) = I₀(λ)f [r/R(λ)], with f being an arbitrary decreasing function of r/R(λ) and R(λ) the typical radius of the source observed at λ. For simplicity, we also assumed that microlensing only affects the B image, so the disk is located in a certain region of the microlensing magnification pattern for B. Although chromatic microlensing is presumably associated with high-magnification regions (e.g., Wambsganss & Paczyński 1991), realistic magnification maps show a wide variety of this type of structures (e.g., Wambsganss 1990). Detailed microlensing simulations are beyond the scope of this paper, and therefore we concentrated on a simple, analytical model for the high-magnification region responsible for the observed chromaticity. Taking a Cartesian coordinate frame in which the origin is placed at the center of the disk, the local microlensing magnification for B was expanded as μ_B(x,y) = ∑ _i,jμ_ijxⁱy^j (i,j = 0, 1, 2...). The intensity distribution I(λ,r), r = (x² + y²)^1/2, was then convolved with μ_B(x,y) to obtain (e.g., Jaroszyński et al. 1992; Shalyapin et al. 2002) ${\mathit{\mu }}_{\mathrm{BA}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{\mathit{\mu }}_{\mathrm{00}}\mathrm{+}\mathrm{(}{\mathit{\mu }}_{\mathrm{20}}\mathrm{+}{\mathit{\mu }}_{\mathrm{02}}\mathrm{)}\mathit{K}{\mathit{R}}^{\mathrm{2}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{+}\mathit{O}\mathrm{\left[}{\mathit{R}}^{\mathit{n}}\mathrm{\right(}\mathit{\lambda }\mathrm{\left)}\mathrm{\right]}\mathit{.}$(5)In Eq. (5), K is a constant, n ≥ 4, and the chromaticity is related to the quadratic and higher order terms in the expansion of μ_B (e.g., Shalyapin 2002). We detect moderate amplification contrasts δμ_BA/μ_BA(λ_V) ≤ 1/3 (δμ_BA = | μ_BA − μ_BA(λ_V) |), therefore we only considered the two first terms in the right-hand side of Eq. (5). If R²(λ) ∝ λ^p, p> 0, our final microlensing model is ${\mathit{\mu }}_{\mathrm{BA}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\mathit{C}\mathrm{+}\mathit{Q}\mathrm{(}\mathit{\lambda }\mathit{/}{\mathit{\lambda }}_{\mathit{V}}{\mathrm{)}}^{\mathit{p}}\mathit{.}$(6)This model was fitted to the two typical microlensing spectra (using the best macrolens/extinction solutions; see the middle and bottom panels of Fig. 8), leading to best-fit values of p in the range 2.3−2.4. It is very remarkable that p ~ 2 − 3 is close to the power-law index for a standard accretion disk (p = 8/3; e.g., Shalyapin et al. 2002).

Very recently, Guerras et al. (2013) have claimed that a significant part of the UV Fe ii and Fe iii emissions may originate in the accretion disk. Adopting the choice of Guerras et al. of wavelength regions (see Table 2 in their paper), we also determined the flux ratio for the two pseudo-continuum spectral features, Fe(1) and Fe(2), associated with iron emission. These Fe(1) and Fe(2) features are dominated by Fe iii and Fe ii, respectively, and are diluted by the true continuum between the C iii] and Mg ii lines (Vestergaard & Wilkes 2001). We find (B/A)_Fe(1) ~ 0.8 and (B/A)_Fe(2) ~ 0.4 in a spectral region where differential extinction is almost achromatic (see the two last rows in Table 2). Hence, the Fe iii emission region seems to be more compact than the region responsible for Fe ii photons, since the (B/A)_Fe(1) value is relatively close to (B/A)_cont ~ 1, but the flux ratio for Fe(2) is close to the ratios for the C iii] and Mg ii emission lines.

3.3. Lens models

Inada et al. (2009) modelled SDSS J1339+1310 using UH88 astro-photometric constraints. They took the relative astrometry and B/A in the second column of Table 1, and a singular isothermal ellipsoid (SIE) mass model (e.g., Lehár et al. 2000). The number of model parameters was the same as the number of observational constraints (d.o.f. = 0), and the LENSMODEL software (Keeton 2001) led to a perfect fit with χ² ~ 0. Table 6 of Inada et al. lists the best-fit values of the mass scale b, the ellipticity e and its position angle θ_e, as well as their 1σ uncertainties. Adopting a redshift of the lens z_l = 0.4 (see the beginning of Sect. 2), the discovery paper also reported the predicted time delay between quasar images (Δt). However, apart from the relative astrometry, the results of Inada et al. are unreliable because the authors used biased estimators of M_BA and z_l. The I-band flux ratio (B/A) of ~0.85 was not corrected for extinction and microlensing, which play a role in this system (see Sect. 3.2). In addition, z_l ~ 0.4 is only a crude photometric estimation of the lens redshift (see the true value in Sect. 3.1).

The new NOT relative astrometry in the fifth column of Table 1 and the new constraints on the macrolens flux ratio in the second column of Table 3 allow us to reliably explore the SIE mass model. We considered that the coordinates of G are Gaussianly distributed around and , with . These are the most uncertain astrometric constraints. The coordinates of B were set to and , neglecting uncertainties of ~0.1%. We also took into account the two Gaussian distributions of Δm_AB that produce the 1σ confidence intervals for M_BA in Table 3: (1) ⟨ Δm_AB ⟩ = −0.41, σ(Δm_AB) = 1.17 (Galactic extinction), and (2) ⟨ Δm_AB ⟩ = −1.86, σ(Δm_AB) = 0.19 (quadratic extinction). Both scenarios (1) and (2) are assumed to be equally likely. From all these Gaussian distributions, we generated several tens of trios (x_G, y_G, M_BA), and used the LENSMODEL package to find the SIE solution for each data set (d.o.f. = 0). The intervals for b, e and θ_e are determined by the broad range of values of M_BA, while astrometry variations play a secondary role. Thus, the 1σ uncertainties in M_BA (see Table 3) were propagated onto the distributions of fitted parameters, or in other words, we focused on the 68% of solutions associated with the 1σ intervals for M_BA.

Fig. 9 SIE mass model of SDSS J1339+1310. The green squares describe the solutions using 0.15 ≤ MBA ≤ 0.21, the red circles describe the solutions using 0.23 ≤ MBA ≤ 2.01 (see main text for details). Top panel: mass scale vs. mass ellipticity. Bottom panel: mass ellipticity vs. position angle. In the bottom panel, the 1σ confidence intervals for the ellipticity and the orientation of the galaxy light are also shown as a black cross.

Two 2D distributions of SIE parameters are shown in Fig. 9. In the top panel, we present the boomerang-shaped relationship between the mass scale and the ellipticity. The bottom panel displays the U-shaped relationship between the ellipticity and its position angle. The green squares and red circles correspond to solutions coming from 0.15 ≤ M_BA ≤ 0.21 (1σ interval using quadratic extinction) and 0.23 ≤ M_BA ≤ 2.01 (1σ interval using Galactic extinction), respectively. The distribution of values of M_BA is very broad, which leads to double-valued functional forms. For example, e ~ 0.2 is related to two values of b and θ_e. The higher b and the negative θ_e come from M_BA ~ 0.2−0.3, while M_BA> 1 yields the lower b and the positive θ_e. For the 1σ intervals of individual mass parameters, we obtain 0.812 ≤ b (′′) ≤ 0.893, 0.123 ≤ e ≤ 0.267 and −40.1 ≤ θ_e (°) ≤ 19.7. Using the spectroscopic redshift of G (z_l = 0.609), we also infer 37.8 ≤ Δt (d) ≤ 44.7. Despite the weak global constraint on M_BA, the uncertainties in the mass scale and the expected delay are below 10%.

Koopmans et al. (2006) reported that an SIE model works well for a sample of Sloan Lens ACS (SLACS) systems, finding a very close alignment between stellar and total mass. However, other recent studies of massive early-type galaxies acting as strong gravitational lenses did not agree with these SLACS results, suggesting the existence of external shear as a source of apparent misalignment (e.g., Gavazzi et al. 2012)⁵. The bottom panel of Fig. 9 shows the shape and orientation of the G galaxy’s light (black cross), indicating an offset between light (stars) and mass position angles. Although this misalignment may be true, it likely means that SDSS J1339+1310 is affected by external shear arising from the environment of G (nearby and/or along the line of sight galaxies). Assuming that the mass of the lensing galaxy is aligned with its light distribution, we also used the best-fit astrometric parameters (neglecting their uncertainties ≤ 2%; see above), the global 1σ range for M_BA (see above), and the 1σ bars of e and θ_e in the bottom panel of Fig. 9 to obtain some insight in an SIE + γ lens model (e.g., Lehár et al. 2000). We focused on 12 representative values of M_BA and the two extreme values of both structure parameters (e, θ_e), constructed all possible trios (M_BA, e, θ_e), and then fitted the mass scale of G and the external shear (strength γ and direction θ_γ) to each data set. For the SIE fits, three properties of the ellipsoid (b, e, θ_e) were allowed to be free. Now (SIE + γ fits) we set (e, θ_e) and allowed the new parameters (γ, θ_γ) to vary. Therefore, the degrees of freedom do not change (d.o.f. = 0).

Fig. 10 SIE + γ mass model of SDSS J1339+1310. The green squares and red circles have the same meaning as those in Fig. 9. Top panel: mass scale vs. shear strength. Bottom panel: shear strength vs. shear direction.

Our SIE + γ results are shown in Fig. 10: mass scale vs. shear strength (top panel) and shear strength vs. shear direction (bottom panel). The green squares are associated with M_BA = 0.15, 0.18, 0.21 (quadratic extinction), and the red circles correspond to M_BA = 0.23, 0.50, 0.69, 0.90, 1.00, 1.12, 1.50, 1.80, 2.01 (Galactic extinction). Figure 10 shows that the M_BA value has a strong influence on the derived external shear. For example, the lowest flux ratios yield γ> 0.1 and θ_γ around −45°, while we obtain γ ~ 0.01 − 0.07 and θ_γ> − 25° using M_BA ≥ 1.50 as a constraint. The ellipticity and its position angle also play a relevant role in the model fitting. In addition, it is interesting to note that 0.800 ≤ b (′′) ≤ 0.946, 0.011 ≤ γ ≤ 0.169 and −54.3 ≤ θ_γ (°) ≤ 20.3. From z_l = 0.609, we obtain 38.3 ≤ Δt (d) ≤ 47.6.

A discussion on the origin of the external shear requires accurate data of the galaxy field around SDSS J1339+1310 (see recent papers by Anguita et al. 2009; MacLeod et al. 2009; Goicoechea & Shalyapin 2010; Wong et al. 2011; Koptelova et al. 2014). There are no SDSS galaxies within a 20′′ radius circle around the lens system (≤130 kpc at z ~ 0.6). Moreover, we do not detect extended sources inside this circle in our deep I-band exposure with the NOT. These NOT/SDSS null detections could mean the absence of massive galaxies close to the main lens. If companion galaxies are not responsible for the external shear in our SIE + γ mass model, γ may arise from other nearby objects. We find 60 SDSS galaxies at z ~ 0.6 inside a 3′ radius circle, that is, at distances ≤1.2 Mpc from SDSS J1339+1310, with r-band magnitudes in the range 20.5 to 25.5. Here, the expression “galaxies at z ~ 0.6” refers to extended sources with redshift error bars including the value 0.6. The number of galaxies in each octant of the sky shows a strong anisotropy in the galaxy distribution, with an overdensity of about 3 in the −45 <θ< 0° octant. This octant contains 17 galaxies, while the other seven have a population of 6 ± 2. Interestingly, the anisotropy in the SDSS galaxy map agrees well with the central part of the distribution of values of θ_γ in the bottom panel of Fig. 10.

4. Conclusions

We are conducting an optical follow-up programme for several bright lensed quasars that are visible from the northern hemisphere³. This is done to perform comprehensive studies on the distant quasars and the galaxies that lens them, using deep imaging and spectroscopy in quasar fields, multiband light curves of quasar images, and polarimetric data. In this paper, we concentrated on the poorly known system SDSS J1339+1310 (Inada et al. 2009). We first presented accurate astro-photometric and spectrophotometric measurements with the Nordic Optical Telescope and the Gran Telescopio Canarias and then discussed the physical properties emerging from the new data. Our results and conclusions are the following:

• 1.

  The new astrometric solutions (relative to the image A) forquasar image B and the main lensing galaxy G are different fromand more precise than the previous solutions by Inadaet al. (2009). These were derivedwith 1 mas (B) and 8−10 mas (G) uncertainties. In addition, we obtained a position angle for G that differs by ~50° from the previous value.

• 2.

  We determined the previously unknown spectroscopic redshift of G: z_l = 0.609 ± 0.001, which differs from the rough photometric estimate z_l ~ 0.4.

• 3.

  We estimated the flux ratio B/A for the cores of the C iv, C iii] and Mg ii emission lines, and the observed chromaticity in (B/A)_core was plausibly associated with differential dust extinction in the lensing galaxy (e.g., Motta et al. 2012). This allowed us to derive constraints on the macrolens flux ratio M_BA and the extinction parameters.

• 4.

  After correcting the continuum flux ratio for macrolens and extinction effects, we found evidence for chromatic microlensing of the continuum emission. This microlensing signal indicates that B is amplified (relative to A) by a factor of about 3−5. To interpret the observed microlensing spectrum, which shows larger amplifications at shorter wavelengths, we considered a face-on disk that contains the continuum sources at different wavelengths. We also assumed that microlensing only affecs the B image, and took a simple model for the involved high-magnification region. Using the best macrolens/extinction solutions to correct the continuum flux ratio, the derived microlensing signal is roughly consistent with the existence of a microlensed, standard accretion disk. This needs to be proved with detailed microlensing simulations, since microlensing spectra of some lens systems suggested the presence of non-standard accretion disks (Motta et al. 2012, and references therein).

• 5.

  We also checked whether weak spectral features associated with UV Fe ii and Fe iii emissions (Vestergaard & Wilkes 2001) originated in compact (continuum) regions or not. While the Fe iii emission region seems to be relatively compact, the Fe ii emission is likely produced in a more extended region. This disagrees with a statistical analysis of other lens systems (Guerras et al. 2013).

• 6.

  Using new observational constraints such as relative astrometry, M_BA, and the luminous structure of G, we discussed two different lens models. We started with an SIE mass model for G, and found possible evidence of external shear. We then studied an SIE + γ model, that consisted of the main lensing galaxy (assuming that the mass of G is aligned with its light distribution; SIE component) plus an external shear from nearby objects (γ component). Our NOT stacked frame and the SDSS database did not show any extended object in the vicinity of G, but there is evidence of a galaxy overdensity at z_l ~ 0.6 in a sky octant associated with the derived shear angles. In spite of the large uncertainty in M_BA, the SIE mass scale is reasonably well constrained: b = 0.̋873 ± 0.̋073 (1σ interval). We also used the new lens redshift to estimate the time delay between A and B, which is expected to be Δt = 43 ± 5 d (1σ; A is leading).

Future observations will help to better understand the physics of SDSS J1339+1310. For example, we only used the red grism of GTC-OSIRIS in a 2700 s spectroscopic exposure, therefore new longer exposures with both the red and blue grisms will allow us to detect the Lyα+N v and Si iv+O iv] emission lines and evaluate the velocity dispersion in G. The additional emission lines will lead to new line-core flux ratios, which must provide valuable constraints on M_BA and the extinction parameters. This is an ongoing project with the GTC that will be completed in next years. NIR spectroscopy probably also plays a critical role in better constraining the macrolens flux ratio/dust extinction law. Moreover, at NIR wavelengths, we expect to detect narrow emission features arising from the NLR (e.g., O [iii]), which are unambiguously free of microlensing effects. Adaptive optics observations may also improve the information about the structure of the lens system, and deep NIR images are valuable diagnostic tools. In particular, a new deep image in the H band could reveal details about the brightness profile of G and the galaxy field around SDSS J1339+1310 (e.g., Koptelova et al. 2014, and references therein). These H-band data may be used to look for galaxies in the vicinity of G, adequately study the galaxy clustering, and estimate its associated convergence/shear.

We note that our initial intention was to determine the time delay of the system. However, preliminary r-band light curves in 2009, 2012 and 2013 with the Liverpool Telescope show that A and B have distinct behaviours, with sharp events in B that are never seen in A. This does not permit us to measure the time delay in a direct way, therefore we will try to separate intrinsic and extrinsic variability through intensive microlensing simulations (e.g., Hainline et al. 2013) and present the corresponding results in a subsequent paper. Microlensing simulations must explain results in the joint time-wavelength domain, because we catch a strong decline in the flux of B around the date of our spectroscopic observations (0.5 mag in less than 100 d), just when chromatic microlensing is detected. This microlensing analysis offers a unique opportunity of constraining the r-band size and wavelength-dependent structure of the quasar accretion disk.

Acknowledgments

We wish to thank the anonymous referee for her/his constructive comments, which have significantly improved the original manuscript. Based on observations made with the Gran Telescopio Canarias (GTC), installed at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. This article is also based on observations made with the Nordic Optical Telescope (NOT) and the Liverpool Telescope (LT), operated on the island of La Palma by the Nordic Optical Telescope Scientific Association and the Liverpool John Moores University (with financial support from the UK Science and Technology Facilities Council), respectively, in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. This research has been supported by the Spanish Department of Science and Innovation grant AYA2010-21741-C03-03 (Gravitational LENses and DArk MAtter − GLENDAMA project), and the University of Cantabria.

References

All Tables

All Figures

Fig. 1 Stacked frame of SDSS J1339+1310 with the NOT in the I band. This corresponds to exposures taken with ALFOSC on 13 June 2010 (total exposure time = 3.6 ks, seeing = 0.̋64). In addition to the original signal (left panel), we also consider a model that includes four additive ingredients: two PSFs (quasar images), the convolution of a 2D de Vaucouleurs profile and the PSF (lensing galaxy), and a constant term (sky background). We display the residual fluxes after subtracting only the quasar images and the background (middle panel), as well as the residual signal after subtracting the global model (right panel). The fit to the model is made with IMFITFITS software (see main text).

In the text

	Fig. 2 Positions of the three components A, B and G of SDSS J1339+1310 on the slit for GTC-OSIRIS. The circles represent seeing disks (diameter = seeing), whereas E, N and Z indicate east, north, and zenith.
In the text

	Fig. 3 1D brightness profile (along the slit) at ~7050 Å. Top panel: observational data (filled circles) and the 1D model solution (black line). The contribution of each source is also detailed. Bottom panel: 2D brightness map for the 1D model solution. The two horizontal dashed lines are the slit edges.
In the text

Fig. 4 GTC-OSIRIS spectra of SDSS J1339+1310. These are derived from a 2.7 ks exposure on 13 April 2013. The red, blue, and green lines trace the fluxes of A, B, and G (see Fig. 2) as a function of observed wavelength. Data from frames in the r (circles) and I (square) bands are also shown for comparison purposes. The dashed lines describe the responses of the two filters, i.e., the LT r-Sloan and the NOT #12 passbands. Vertical dotted lines are associated with emission lines at zs = 2.24, while grey highlighted regions indicate atmospheric artefacts (see main text).

In the text

Fig. 5 Spectrum of the main lensing galaxy in SDSS J1339+1310. The observed spectrum (see Fig. 4) is smoothed by a 5-point filter and represented by a green line, while the template spectrum for an elliptical galaxy is properly red-shifted (zl = 0.609) and drawn with a magenta line. Grey highlighted regions and vertical dashed lines are associated with sky absorption/emission and galaxy absorption lines, respectively (see main text).

In the text

	Fig. 6 C iii] emission line profiles. For each quasar image, we take the continuum in the intervals 5800−5950 and 6350−6500 Å, fit it to a quadratic law, and then remove the fitted quadratic trend. The solid line traces the line profile for the image A. The dashed line describes the line profile for image B.
In the text

	Fig. 7 Extinction curve for SDSS J1339+1310. Filled circles with error bars represent the line-core magnitude differences in Table 2. The red lines describe the best fits, and the light grey areas are the 1σ bands. Top panel: galactic extinction. Bottom panel: quadratic extinction (see main text).
In the text

Fig. 8 Continuum flux ratio and microlensing signal. Top panel: continuum flux ratio. Middle panel: microlensing magnification ratio using the best fit (red circles) and the 1σ band (light grey region) for Galactic extinction. Bottom panel: microlensing magnification ratio using the best fit (red circles) and the 1σ band (light grey region) for quadratic extinction. In the middle and bottom panels, the microlensing signal from the 1σ bands is linearly interpolated in the spectral regions where the C iv, C iii] and Mg ii emission lines appear. We also show the zero-level for microlensing (dashed lines; μBA = 1) and fits to a simple microlensing model (solid lines; see main text).

In the text

Fig. 9 SIE mass model of SDSS J1339+1310. The green squares describe the solutions using 0.15 ≤ MBA ≤ 0.21, the red circles describe the solutions using 0.23 ≤ MBA ≤ 2.01 (see main text for details). Top panel: mass scale vs. mass ellipticity. Bottom panel: mass ellipticity vs. position angle. In the bottom panel, the 1σ confidence intervals for the ellipticity and the orientation of the galaxy light are also shown as a black cross.

In the text

	Fig. 10 SIE + γ mass model of SDSS J1339+1310. The green squares and red circles have the same meaning as those in Fig. 9. Top panel: mass scale vs. shear strength. Bottom panel: shear strength vs. shear direction.
In the text