SIZE OF THE ACCRETION DISK IN THE GRAVIATIONALLY LENSED QUASAR SDSS J1004+4112 FROM THE STATISTICS OF MICROLENSING MAGNIFICATIONS

The flux variability of the images of a gravitationally lensed quasar is a combination of the intrinsic variability of the source, correlated by a time delay between the different images, and gravitational microlensing that depends on the random distribution of stars in the lens galaxy. This last variability is uncorrelated between images (Chang & Refsdal 1979; see also the review by Wambsganss 2006). The analysis of lensed quasars' light curves has important applications in cosmology (determination of time delays to infer the Hubble constant; Refsdal 1964) and in the study of quasar structure (Chang & Refsdal 1979, 1984; see also Kochanek 2004 and Wambsganss 2006). In this paper we will focus in the last application, using microlensing statistics to determine the quasar accretion disk size (Fohlmeister et al. 2007, 2008; Pooley et al. 2007; Fohlmeister 2008; Mosquera et al. 2009, 2013; Morgan et al. 2010; Blackburne et al. 2011, 2014, 2015; Mosquera & Kochanek 2011; Sluse et al. 2011; Jiménez-Vicente et al. 2012, 2014; Motta et al. 2012; Hainline et al. 2013; Jiménez-Vicente et al. 2015a, 2015b; MacLeod et al. 2015; Mediavilla et al. 2015; Muñoz et al. 2016).

The wide-separation lensed quasar SDSS J1004+4112 is a rare example of a quasar lensed by an intervening galaxy cluster (Wambsganss 2003; Inada et al. 2006). It is lensed into five images (Inada et al. 2005), with a maximum image separation of $14\buildrel{\prime\prime}\over{.} 62$. The quasar has a redshift of z_s = 1.734 and the redshift of the galaxy cluster is ${z}_{l}=0.68$ (Fohlmeister et al. 2008). The lag between images A, B, and C has been experimentally measured by Fohlmeister et al. (2008) who obtained a time delay of 40.6 days between components A and B and of 822 days between A and C (the largest delay measured for a gravitational lens system). It has not yet been possible to measure the time delay of image D. From a model of SDSS J1004+4112 in which the mass distribution of the system is revisited, Oguri (2010) predicts a time delay between A and D of 1218 days

After correcting for the time delays and mean magnitude differences, Fohlmeister et al. (2008) found clear indications of microlensing flux variability in the residuals of the light curves (i.e., the differences between the observed light curves and the modeled intrinsic variability of the quasar). They fitted the residual light curves of image A using a model of microlensing flux variability based on magnification maps to estimate the accretion disk size of the lensed quasar, obtaining a small size below the predictions of thin disk theory.

In this paper we present seven years, corresponding to eight observational seasons, of optical monitoring data for the four brightest images of SDSS J1004+4112 spanning 2505 days from 2003 December to 2010 October, partly coincident with the epochs studied by Fohlmeister et al. (2008), but which significantly extend the coverage by around 1200 additional days. Our aim is to use these data to study the existence of possible microlensing events and to estimate the size of the accretion disk of the lensed quasar. We have been able to carry out a microlensing analysis including image D for the first time thanks to the availability of sufficiently long observations to correct for the time delay shifting the light curves.

Instead of the light curve fitting method (see e.g., Kochanek 2004) we will follow the single epoch method combined with the flux ratios of a sufficiently large source in the quasar as to be insensitive to microlensing to establish the baseline for no microlensing magnification (see, e.g., Mediavilla et al. 2009). Based on single epoch spectroscopy Motta et al. (2012) and Jiménez-Vicente et al. (2014) obtained disk size estimates larger than the results of Fohlmeister et al. (2008). In the present work we will extend the single epoch method to all the epochs in the available light curves, increasing the statistical significance. The basic idea is to compare the histogram of microlensing magnification obtained from the observations corresponding to a time interval with the simulated predictions of microlensing variability for sources of different sizes. This comparison will allow us to evaluate the likelihood of the different values adopted for the size. We will use the optical light curves, with coverage extended by us, to infer microlensing flux variability and the mid-IR data from Ross et al. (2009) to determine the baseline for no microlensing variability.

The paper is organized as follows. In Section 2 we present the data and in Section 3 the light curves of each image. The estimate of the quasar accretion disk size based on the statistics of microlensing magnifications is discussed in Section 4. Finally, the main results are summarized in Section 5.

The photometric monitoring presented in this paper took place between 2003 December and 2010 October. We monitored SDSS J1004+4112 in the Johnson-Bessell's R-band using the 82 cm telescope (IAC80) at the Instituto de Astrofísica de Canarias Teide Observatory (Tenerife, Canary Islands, Spain). Two different CCDs were used. From 2003 to 2005 a Thomson 1024 × 1024 chip was employed, giving a field of view of about 75, with a pixel size of 043. Since 2005 a new CCD, CAMELOT, was installed. CAMELOT hosts a E2V 2048 × 2048 chip with 0304 pixels, corresponding to a 10.4 × 10.4 arcmin² FOV. A standard R broadband filter was always used for the observations, fairly close to the Landolt R (Landolt 1992). Due to the instrumental change a constant gap of ∼0.2 mag appeared in the light curves that has been removed taking the new instrumental system as a reference. The combined data set consists of 109 epochs (i.e., 109 nights) and the average sampling rate is once every 23rd day. This large average sampling rate arises due to the relatively large seasonal gaps. There are seven seasonal gaps, the largest one stretching over 11 months from 2007 June to 2008 May, i.e., a time period of more than 300 days. The mean observational cadence of the four images is ∼9 days for the first two seasons, ∼10 days for the third and fourth season, ∼6 days for the fifth season and ∼18 days for the last two seasons.

Figure 1 shows an image taken with the IAC80 of the four brightest quasar images A, B, C, and D of SDSS J1004+4112. Photometric data were obtained by applying a completely automatic IRAF task, pho2com, developed by Serra-Ricart et al. (1999). This code yields accurate photometry by simultaneously fitting a stellar two-dimensional profile to each quasi-stellar object (QSO) component by means of DAOPHOT software. To remove inconsistent data due to instrumental problems or other sources of error related to the data reduction we have, in first place, removed the observations in which a sudden change in magnitude simultanously appears in all the images. In a second step, the magnitudes between two consecutive points were compared and if the difference in magnitude was greater than twice the standard deviation, the point was discarded. Only around 10 data points out of 109 had to be rejected and around 100 data points were left for analysis.

Zoom In Zoom Out Reset image size

In Figure 2 we show the resulting light curves of the quasar images A–D split over eight observing seasons. The dashed vertical line shown at 3700 days represents the change of the old to the new CCD and one can clearly see that the error bars of the observations with the new CCD are smaller by almost a factor of two. The images C (yellow light curve) and D (red light curve) are shifted by 0.3 mag and 1.0 mag respectively so that they do not overlap with each other.

Zoom In Zoom Out Reset image size

Due to our larger uncertainties, we have not been able to improve the time delay measurements of Fohlmeister et al. (2008), therefore we use their time delay estimations (for images A, B, and C) to shift the light curves and remove intrinsic variability. For calculating the time delay between image A and B they used the polynomial fitting method and the analysis of their data yields a time delay of ${\rm{\Delta }}{t}_{\mathrm{BA}}=40.6\pm 1.8$ days. They also measured the time delay for the wide separated (1462) image C relative to the close image pair A/B. Using the dispersion spectra method they found ${\rm{\Delta }}{t}_{\mathrm{CA}}=822\pm 7$ days and ${\rm{\Delta }}{t}_{\mathrm{CB}}=780\pm 6$ days. Image D should lag the other three images but to date no features can be seen in the light curve of image D than can be matched to the first season of image A. They derived a lower limit on the time delay between images A and D of ${\rm{\Delta }}{t}_{\mathrm{DA}}\gt 1250$ days (3.4 years) (Fohlmeister et al. 2008; Fohlmeister 2008). We will use the model-predicted AD time delay of 1218 days of Oguri (2010), which is slightly smaller than the lower limit reported by Fohlmeister et al. (2008).

Figure 3 shows a comparison of the IAC data and the data used by Fohlmeister et al. (2008) for the quasar images A and B. The light curves show, in general, a good agreement with those of Fohlmeister et al. (2008) after applying global shifts to match both photometries in the overlapping regions. The exception is a region of image B at JD 3750 where the data of Fohlmeister et al. seem to be slightly below the expected values if one compares with the light curve for image A.

Zoom In Zoom Out Reset image size

In Figure 2 one cannot see strong differences between the light curves A and B that could be directly related to microlensing. Thus, we can estimate the amplitude of the intrinsic variability of the quasar performing a single polynomial fitting to the B light curve. Despite the similarity between the A and B light curves, weak microlensing could be still present due to a large source size or a location of both images in regions with low microlensing magnification, even strong microlensing magnification but with different sign in A and B. Consequently, although we use the polynomial fitting to make a tentative estimate for the intrinsic variability of the source, some contribution from microlensing variability cannot be completely discarded. Hence, we should compare with the A–B, C–B and D–B modeled magnification pattern/histograms for all the images. We obtain a source variability of ∼0.7 mag for all the data and of 0.5 mag for the first four seasons whereas Fohlmeister et al. (2008) found an intrinsic variability of ∼0.7 mag. Figure 4 shows the simulated quasar variability (black solid line) with the A, B, C, and D light curves overimposed. In the three panels in Figure 5 the residuals of the A, C and D light curves are shown. We have calculated the residuals from: ${\rm{\Delta }}{m}_{X}={m}_{X}-{m}_{\mathrm{Bfit}}-{({m}_{X}-{m}_{B})}_{\mathrm{midIR}}$ where X = A, C, D. The mid-IR data have been taken from Ross et al. (2009). The mid-IR emission is supposed to arise from a large enough region as to be not affected by microlensing and, hence, we use the ${({m}_{x}-{m}_{B})}_{\mathrm{midIR}}$ offsets to set the microlensing baseline (see Mediavilla et al. 2009). In Figure 5 one can see that the residuals of the A light curve are relatively constant. Only in the fourth season, variability induced by microlensing may occur with amplitudes of order 0.15 mag. As opposed to this, microlensing variability is clearly visible in the residuals of the C and D light curves with amplitudes of the order of ∼0.5 mag and ∼0.7 mag. From these curves that represent the differential (with respect to B, the image less prone to microlensing) microlensing of the A, C, and D images, we have obtained the microlensing variability histograms (Figure 6), i.e., the frequencies in which each microlensing amplitude appears in the microlensing variability light curves. We adopt a bin size of 0.2 mag.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Microlensing Estimate of the Quasar Accretion Disk Size

Taking into account that microlensing is sensitive to the size of the source (Morgan et al. 2010; see also the review by Wambsganss 2006) we will use our determinations of microlensing magnification amplitude to estimate the size of the accretion disk in the SDSS J1004+4112 lensed quasar.

We obtained microlensing magnification maps for each image using the inverse polygon mapping method described by Mediavilla et al. (2006, 2011). We used the following parameters obtained from a singular isothermal sphere plus external shear (SIS+${\gamma }_{e}$) model fitted to the images coordinates: convergence $\kappa =0.48$ and shear $\gamma =0.57$ for image A, $\kappa =0.47$ and $\gamma =0.39$ for B, $\kappa =0.38$ and $\gamma =0.33$ for C and $\kappa =0.71$ and $\gamma =0.83$ for image D. We used a surface mass density in stars ${\kappa }_{* }$ of 10% (Mediavilla et al. 2009) and generated 2000 × 2000 pixel² magnification maps with a size of 24 × 24 Einstein radii². The value of the Einstein radius for this system is $2.35\times {10}^{16}\sqrt{M/0.3{M}_{\odot }}$ cm = $9.1\sqrt{M/0.3{M}_{\odot }}$ lt-days at the lens plane. The ratio of the magnification in a pixel to the average magnification of the map gives the microlensing magnification at the pixel and histograms of normalized to the mean maps deliver the relative frequency of microlensing magnification amplitude for a pixel-size source.

To model the structure of the unresolved quasar source we considered a circular Gaussian intensity profile of size r_s, $I(R)\propto \exp (-{R}^{2}/2{r}_{s}^{2})$. It is generally accepted that the specific shape of the radial profile is not important for microlensing flux variability studies because the results are essentially controlled by the half-light radius rather than by the detailed profile (Mortonson et al. 2005). The characteristic size r_s is related to the half-light radius by ${R}_{1/2}=1.18{r}_{s}$. We convolve the magnification maps with Gaussians of 14 different sizes over a logarithmic grid which spans approximately from r_s ∼0.2 to 55 lt-days for a mean stellar mass $\langle M\rangle =0.3{M}_{\odot }$. The source sizes can be scaled to a different mean stellar mass, $\langle M\rangle$, using ${r}_{s}\propto \sqrt{M}$. After convolution we normalized each magnification map by its mean value. The histograms of the normalized map represent the histograms of the expected microlensing variability. Thus, we obtain 14 different microlensing histograms corresponding to different source sizes for each of the images A, C, and D. Finally, convolving the histograms of A, C, and D with the histogram of B we built the microlensing difference histograms A–B, C–B and D–B for different values of r_s to be compared with the experimental histograms obtained from the observed light curves (see Figure 6). The thick solid line in each panel indicates the best fit of the model to the observed frequency distributions of microlensing magnifications.

In order to study the likelihood of the different r_s we compare the microlensing histograms inferred from the model for different values of r_s with the histograms of the data using two different statistics:

• i.  
  histogram product, defined as

  $$\begin{eqnarray}&&{P}_{X}({r}_{s})=\displaystyle \sum _{i=1}^{{n}_{\mathrm{bin}}}\,{h}_{X-B}^{i}{\hat{h}}_{X-B}^{i}({r}_{s}),\end{eqnarray} \tag{ 1 }$$

  where ${h}_{X-B}^{i}$ and ${\hat{h}}_{X-B}^{i}({r}_{s})$ are the observed and modeled histograms and n_bin is the number of bins
• ii.  
  Pearson's ${\chi }^{2}$ test,

  $$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{{n}_{\mathrm{bin}}}\,\displaystyle \frac{{({h}_{X-B}^{i}-{\hat{h}}_{X-B}^{i}({r}_{s}))}^{2}}{{h}_{X-B}^{i}}\end{eqnarray} \tag{ 2 }$$

  with

  $$\begin{eqnarray}&&{P}_{X}({r}_{s})\propto {e}^{-\tfrac{1}{2}{\chi }^{2}}.\end{eqnarray} \tag{ 3 }$$

After multiplying the probability distributions corresponding to A, C and D we finally obtain the probability density function (PDF) of the source size,

$$\begin{eqnarray}&&P({r}_{s})={P}_{{\rm{A}}}({r}_{s})\cdot {P}_{{\rm{C}}}({r}_{s})\cdot {P}_{{\rm{D}}}({r}_{s}).\end{eqnarray} \tag{ 4 }$$

In the case of the histogram product method we have convolved the histograms with normal distributions with these dispersions, not finding any significant difference.

To account for the rms errors in the data in the case of the ${\chi }^{2}$-method we applied a Monte Carlo method, producing random realizations of the histograms of the data using normal distributions of mean the value of the observed histogram and with the following dispersions: ${\sigma }_{{\rm{A}}}=0.08$, ${\sigma }_{{\rm{C}}}=0.14$ and ${\sigma }_{{\rm{D}}}=0.22$ estimated from the data (see Figure 5). In the case of the histogram product-method we convolved the histograms with normal distributions with these dispersions, not finding any significant difference.

Figure 7 shows the resulting normalized probability distribution for both methods. For the product we obtain a disk size of $\langle {r}_{s}\rangle ={7.4}_{-4.7}^{+15.7}\sqrt{M/0.3{M}_{\odot }}$. The expected value using ${\chi }^{2}$ is $\langle {r}_{s}\rangle ={4.4}_{-1.0}^{+1.0}\sqrt{M/0.3{M}_{\odot }}$ lt-days (for 68% confidence estimates). The ${\chi }^{2}$-method seems to constrain better the size although the errors in the histograms may have been underestimated. This is supported by the rather high values of ${\chi }^{2}$. To reach values of ${\chi }^{2}$ close to 1 we need to multiply the typical deviation in Equation (2), $\sqrt{{h}_{X-B}^{i}}$, by a factor of 3. From the corresponding PDF (see Figure 7) we obtain a similar value $\langle {r}_{s}\rangle ={3.6}_{-1.9}^{+2.7}\sqrt{M/0.3{M}_{\odot }}$ but with uncertainties more compatible with the results obtained using the product method.

Zoom In Zoom Out Reset image size

Our result for each method expressed in terms of the half-light radius at ${\lambda }_{\mathrm{rest}}=2407\,\mathring{\rm A}$ (${R}_{1/2}={8.7}_{-5.5}^{+18.5}\sqrt{M/0.3{M}_{\odot }}$ and ${R}_{1/2}={4.2}_{-2.2}^{+3.2}\sqrt{M/0.3{M}_{\odot }}$ lt-days) is significantly greater than the value obtained by Fohlmeister et al. (2008) fitting the A and B light curves (${R}_{1/2}=2.44\ {R}_{\lambda }={0.6}_{-0.3}^{+0.5}$ lt-days), but in good agreement with the value inferred by these authors from thin disk theory using an experimental estimate of the black hole mass (Fohlmeister et al. 2008). Our estimate for the size is also in good agreement with the results of Motta et al. (2012) and Jiménez-Vicente et al. (2014) for this system, and with the average determinations obtained for a sample of lensed quasars by Jiménez-Vicente et al. (2012, 2014, 2015a, 2015b) when a fraction of mass in stars of 10% is considered.

As we have used an estimated (not measured) time delay for image D, and this image presents the largest microlensing signal in this system, it is important to check that this assumption is not biasing our results. In order to check the dependence of the results with the adopted time delay for image D, we considered ±100 days shifts of the time delays not finding any significant difference in the size calculations. On the other hand, we also repeated the computations not considering image D finding slightly larger but comparable results for the size: $\langle {r}_{s}\rangle ={9.3}_{-5.2}^{+13.8}\sqrt{M/0.3{M}_{\odot }}$ for the product-method and $\langle {r}_{s}\rangle ={5.4}_{-2.7}^{+4.3}\sqrt{M/0.3{M}_{\odot }}$ for ${\chi }^{2}$.

We also checked the impact of interchanging the two images less affected by microlensing, A and B, by performing the polynomial fitting to the A light curve and considering the pairs, B–A, C–A and D–A in the calculations. The results are almost identical for the ${\chi }^{2}$-based method, $\langle {r}_{s}\rangle ={4.6}_{-1.9}^{+5.2}\sqrt{M/0.3{M}_{\odot }}$, and compatible to within errors for the product-method, $\langle {r}_{s}\rangle ={4.6}_{-2.9}^{+5.1}\sqrt{M/0.3{M}_{\odot }}$.

We have presented eight seasons of monitoring data for the four brightest images of the five-image gravitational lens system SDSS J1004+4112 which significantly extend the time coverage of previous works. Taking as reference image B, which is less affected by microlensing, and using the experimental time delays inferred by Fohlmeister et al. (2008) for components A and C and the theoretical time delay modeled by Oguri (2010) for image D, we have removed the intrinsic variability from the light curves in the overlapping region. Using the mid-IR flux ratios between images determined by Ross et al. (2009) as a baseline for no microlensing magnification, we finally obtained the microlensing light curves, A–B, C–B, and D–B. We detected microlensing variability up to 0.7 mag in images C and D. The light curve of A seems to be less affected by microlensing (changes of order 0.15 mag in the fourth season of our data) in good agreement with the results of Fohlmeister et al. (2008).

We have used the statistics of microlensing magnifications along the available seasons to infer probabilistic distributions for the size using two different methods. Using the product of the observed and modeled microlensing histograms we have obtained a half-light radius of ${R}_{1/2}={8.7}_{-5.5}^{+18.5}\sqrt{M/0.3{M}_{\odot }}$ lt-days. A consistent but more restrictive result of ${R}_{1/2}={4.2}_{-2.2}^{+3.2}\sqrt{M/0.3{M}_{\odot }}$ lt-days is obtained using a ${\chi }^{2}$ criterion to compare the observed and modeled histograms. Both results are only marginally consistent with the previous estimate by Fohlmeister et al. (2008) that obtained a smaller size, but in good agreement with the predictions of thin disk theory and with the measurements by Motta et al. (2012) and Jiménez-Vicente et al. (2014).

We thank the anonymous referee for a thorough revision and valuable suggestions. We are especially grateful to the Instituto de Astrofsica de Canarias astronomers for their observations of the data appearing in this paper. The 0.82 m IAC80 Telescope is operated on the island of Tenerife by the Instituto de Astrofsica de Canarias in the Spanish Observatorio del Teide. JJV is supported by the Spanish Ministerio de Economía y Competitividad and the Fondo Europeo de Desarrollo Regional (FEDER) through grant AYA2014-53506-P and by the Junta de Andalucía through project FQM-108. EM and JAM were supported by the Spanish MINECO with the grants AYA2013-47744-C3-1-P and AYA2013-47744-C3-3-P. JAM was also supported by the Generalitat Valenciana with the project PROMETEOII/2014/060.