Gravitationally lensed quasars provide a wealth of resources to observers. Their utility in cosmology was realized quite early on (e.g., Refsdal 1964), and a number of collaborations (e.g., COSMOGRAIL, Courbin et al. (2005) and H0LiCOW 6 ) continue to pursue measurements of lensed quasar time delays to make independent measurements of the Hubble Constant, H 0 (e.g., Kochanek 2002;Vuissoz et al. 2008;Fohlmeister et al. 2013;Bonvin et al. 2017;Suyu et al. 2017) and a range of other useful cosmographic measurements. Quasar microlensing, also predicted quite some time ago (e.g., Chang & Refsdal 1979), provides additional motivation for monitoring lensed quasars because the analysis of microlensing variability (e.g., Kochanek 2004;Morgan et al. 2006Morgan et al. , 2010;;Poindexter et al. 2007;Hainline et al. 2013;MacLeod et al. 2015) and chromatic flux ratio anomalies (e.g., Pooley et al. 2007;Bate et al. 2008Bate et al. , 2011;;Mediavilla et al. 2011;Pooley et al. 2012;Schechter et al. 2014) can be analyzed to probe the central engines of the quasars and the properties of the lens galaxy.
The quadruply lensed quasar WFI J2033-4723 (hereafter WFI2033; 20 h 33 m 42 08, -47°23′43 0 [J2000.0]) was discovered during a wide-field imaging survey for lensed quasars in the southern hemisphere using the MPG/ESO 2.2 m telescope (Morgan et al. 2004). It has a source redshift of z s =1.66, a lens redshift of z l =0.661 (Eigenbrod et al. 2006), and a maximum image separation of 2 5. Vuissoz et al. (2008, hereafter V08) used three seasons (2004)(2005)(2006)(2007) of monitoring data from the Small and Moderate Aperture Research Telescope System (SMARTS) 1.3 m telescope at the Cerro Tololo Inter-American Observatory (CTIO) and the 1.2 m Leonhard EULER Swiss telescope at La Silla, Chile to measure time delays of t | | days between the merged A1+A2=A, B, and C images. V08 found no evidence of variability due to extrinsic factors such as microlensing. Recently, however, Giannini et al. (2017) made a robust detection of microlensing in their four-season monitoring campaign using the 1.54 m Danish telescope at La Silla, a result that we independently corroborated in this investigation. Most recently, Motta et al. (2017) used the single-epoch chromatic microlensing technique to make estimates of the size of the central engine and broad line region (BLR) in WFI2033.
In this paper, we combine 9 new seasons of WFI2033 monitoring data with the 4 seasons of data from V08 to create a 13-season set of light curves. We present our observational data and reduced light curves in Section 2, and we analyze the full combined light curves in Section 3 to confirm the V08 time delays and measure the A1-A2 delay for the first time. In Section 4, we describe our microlensing analysis technique to include the properties of our strong lens models for WFI2033. In Section 5, we present the results of our analysis, and we discuss their implications for accretion disk theory. Throughout our discussion, we assume a flat cosmology with Ω M =0.3, Ω Λ =0.7, and H 70 km s Mpc 0 1 1 = -- (Hinshaw et al. 2009).
We observed WFI2033 in the V-(F555W), I-(F814W), and H-(F160W) bands using the Hubble Space Telescope (HST) as an element of the the CfA-Arizona Space Telescope Survey (CASTLES, 7 Lehár et al. 2000). The V-and I-band images were taken using the Wide-Field Planetary Camera 2 (WFPC2). The H-band images, originally presented in V08 and displayed here in Figure 1, were taken using the Near-Infrared Camera and Multi-Object Spectrograph (NICMOS). We fit the astrometry and photometry of the lens in the HST imagery with the imfitfits (Lehár et al. 2000) routine, using a de Vaucouleurs model for the lens galaxy G1, an exponential disk model for the quasar host galaxy and point sources for the quasar images. Our astrometric and photometric fits, consistent with those made independently by V08, are presented in Table 1.
On the 1.3 m SMARTS telescope, we used the optical channel of the dual-beam ANDICAM instrument (DePoy et al. 2003), which has a plate scale of 0 369 pixel -1 and a 6 5×6 3 field of view. The mean sampling of the SMARTS data is one epoch every eight days, with three 300 s exposures at each epoch using the R-band filter. The R-band filter has an effective wavelength of 658 nm, translating to a rest-frame wavelength of 2473 Å in the UV. The SMARTS data set consists of 117 epochs between 2004 March and 2017 August.
On the 1.2 m EULER telescope, we used the EulerCAM camera, which has a plate scale of 0 2149 pixel -1 and a 15 0×15 0 field of view. The mean sampling of the EULER data is one epoch every five days, with five 360 s exposures at each epoch using the "RG" or "Rouge Genève" filter. The RG filter is a modified broad Gunn R filter, with an effective wavelength of 660 nm translating to a rest-frame wavelength of 2481 Å. The new EULER data set consists of 178 epochs between 2010 October and 2016 December.
The details of our photometric measurement technique are discussed in Kochanek et al. (2006), but we provide a brief summary of that process here. We use five reference stars, located at (-2 4, 61 4), (-16 1, 15 3), (56 0, 70 3), (60 2, -0 9), and (-34 4, 82 5), with respect to image A1. These reference stars are used as the basis for a three-component elliptical Gaussian point-spread function (PSF) model, which we apply to the blended quasar images to obtain the relative brightness of each component at each epoch. When applying the PSF model, we hold the relative positions of the lens components fixed to the astrometry from our HST H-band images. We model the lens galaxy using a nested Gaussian with fixed effective radius and flux to approximate a de Vaucouleurs profile. For the effective radius, we used our measurement from the HST fits r eff =0 83±0 1, and for the flux we use the value that minimizes the total χ 2 in the residuals following galaxy model subtraction when summed over all epochs. We measured a very small color offset of 0.002 mag between the EULER and SMARTS photometry that we applied to the SMARTS data when creating our combined light curves and the data provided in Tables 2 and3. Both the EULER and SMARTS images are characterized by a median stellar FWHM (seeing) of 1 2. Since the merging A1/A2 pair are separated by only 0 72, deconvolving the flux from these two images was challenging. For seeing conditions worse than 1 5 and 1 62 for SMARTS and EULER, respectively, we were unable to reliably resolve any of the quasar's images, so we were forced to discard images taken under these conditions. We also discarded 31 of the 326 total observing epochs from SMARTS and EULER due to bright sky or cloudy observing conditions. In Figure 2, we display our new light curves alongside the published light curves from V08. Since V08 were
Using the polynomial light curve fitting technique of Kochanek et al. (2006), we measured the time delays between
) . The Image A1-C columns give the magnitudes of the quasar images relative to the comparison stars.
(This table is available in machine-readable form.) the combined image A=A1+A2, image B, and image C. In the Kochanek et al. (2006) method, the intrinsic and extrinsic variability in the light curves are fitted by Legendre polynomials, and the polynomial order is chosen using the F-test. In the case of the delay between images A and B and images C and B, we found that a N 10th source = -order polynomial provided a sufficient fit for the source variability and that a N 2nd = m -order polynomial was appropriate for approximating and removing the microlensing variability. In Figure 3, we show the χ 2 statistic for the time-delay fits. ) . The Image A1-C columns give the magnitudes of the quasar images relative to the comparison stars.
(This table is available in machine-readable form.)
images A and C). In Figure 4, we show the light curves shifted by these delay values. These new measurements are fully consistent with those of V08. Using our newly reduced light curves, we also obtain a tentative measurement of the delay between images A1 and A2. With a N 5th source = -and N 2nd = m -order polynomials for the intrinsic and microlensing variability, respectively, we find that image A1 leads image A2 by t t t 3.9
days. We display these results in Figure 5. While V08 were not able to measure the A1-A2 delay, they did constrain the expected range of that delay to -1>Δt A1A2:model >-3 days using a series of lens galaxy mass models. While significantly coarser than our measurement of the B-A and B-C delays, the A1-A2 measurement is consistent with the V08 lens models, although this pair will have the largest fractional uncertainties from microlensing-induced variability (Tie & Kochanek 2018).
In the present paper, we generate a series of lens galaxy models in which the expected A1-A2 delay is -1.6>Δt A1A2(model) > -3.3 days, also consistent with our A1-A2 measurement. With these updated time-delay measurements, we proceed to the primary goal of this investigation, the analysis of extrinsic variability from microlensing in the reduced light curves. A full analysis of the updated delays will be published in V. Bonvin et al. (2018, in preparation).
In essence, our Bayesian Monte Carlo technique for microlensing analysis is an attempt to reproduce the observed microlensing variability using a large set of models for the physical conditions that might have led to this variability (Kochanek 2004). All of this hinges on our ability to accurately model the conditions in the lens galaxy through which the quasar's light must pass. We started by applying the LENSMODEL code of Keeton (2001) to the astrometry from our HST observations to yield a range of models for the stellar and dark matter content in the lens galaxy at the positions of the lensed images. Following V08, we adopted a two-component model for the lens galaxy (G1 in Figure 1). Since this system is now known to exhibit microlensing of both the continuum and the BLR (Sluse et al. 2012;Motta et al. 2017), we did not use the HST flux ratios or those from our ground-based observations as a constraint on the lens galaxy mass models. We required our models to reproduce the astrometry of the lensed images, and we allowed the position, effective radius, ellipticity, and position angle of the lens galaxy to vary within the uncertainties of the photometric model presented in Section 2.1. Consistent with V08, we were unable to model the astrometry of the lensed images unless we included the influence of the neighboring galaxy G2, the east-west shear from which cannot be created by G1 since it has an ellipticity position angle of only 25°east of north. We modeled G2 as a singular isothermal sphere whose properties were also allowed to vary within the uncertainties of our HST photometric and astrometric fits.
Since the dark matter content in the lens galaxy is unknown, we created a series of 10 models for the lens galaxy in which the dark matter fraction varies across an order of magnitude. We began by modeling the lens galaxy using only a de Vaucouleurs profile. In each subsequent model, we decreased the monopole moment of the stellar de Vaucouleurs component by 10% of the constant M/L model mass, and we added an extended, concentric (Navarro et al. 1996, NFW) component to model the dark matter. We parameterized this series using the quantity f M/L , representing the fraction of the lens galaxy mass relative to the constant mass-to-light (M/L) ratio model. From this model sequence, we extract the total convergence κ, the 4.
Using the inverse ray-shooting technique (as described in Kochanek et al. 2006), we generated 40 random realizations of the expected microlensing magnification conditions in the vicinity of each image for each of our 10 macro models Note. Convergence κ, shear γ, and the fraction of the total surface density composed of stars κ * /κ at each image location for the series of mass models. The parameter f M/L =1.0 is the mass of the de Vaucouleurs model for the visible lens galaxy relative to its mass in the absence of dark matter. The χ 2 per degree of freedom for each model is provided in the N parameterized by f M/L . The magnification patterns are 8192×8192 pixels representing a projected source plane scale of 20 1 M e Einstein radii or 8.66×10 17 cm. This implies a pixel scale of 1.06×10 14 cm in the source plane. We assumed an initial stellar mass function (IMF) of dN(M)/dM∝ M -1.3 with a dynamic range of 50, which approximates the microlensing-based Galactic bulge IMF of Gould (2000), although Wyithe et al. (2000) and Congdon et al. (2007) show that microlensing statistics are not especially sensitive to choice of IMF.
Armed with 400 sets of magnification patterns for a range of lens galaxy mass models, we used the Monte Carlo light curve fitting technique of Kochanek (2004) as modified by Morgan et al. (2008) to model unresolved image pairs. Since the light curves are a full 13 seasons in length, we binned them using a window of δt bin =20 days to reduce computation time. The 20 day binning window was sufficiently short to avoid overly smoothing the microlensing variability while adequately reducing the run time for the Monte Carlo routine. The date of each 20 day bin was set as the mean Heliocentric Julian Date (HJD) of the measurements included in that bin. Since the light curves from V08 do not provide individual measurements for the A1 and A2 images, we adjusted the statistical weight of the A=A1+A2 data points to appropriately account for the combined fluxes in these cases. Also, for the combined cases we used either the A1 or A2 magnification pattern (with equal frequency) when creating a simulated light curve for the combined A1+A2 image.
Prior to each Monte Carlo trial, we convolved each set of magnification patterns with a Gaussian surface brightness profile at a range of trial source sizes r M M 14.5 log cm 18.0, 1
where r s ˆis the radius of the accretion disk scaled by M 1 2 * á ñ , the mean mass of a lens galaxy star. Although we used a Gaussian profile, the exact choice of photometric emission model is unimportant, since Mortonson et al. (2005) showed that microlensing statistics are largely a function of the halflight radius of the emitting region, not the exact properties of the emission profile. In a given trial, a convolved pattern is run past a model point source on a random trajectory and at a random transverse speed 10 km s -1 v M M 10 e 1 2 6 * á ñ - ( ˆ) km s -1 . Changes in magnification with time are logged at the epochs of the observed data and a running comparison with the light curves is made. The quality of the fit is tallied in real time using a χ 2 statistic, and, to save computational time, fits with χ 2 /N dof >3.0 are discarded since they will not contribute significantly to the Bayesian integrals in our post-run analysis. During the curve fitting process, we allowed for 0.07 and 0.03 mag of systematic error in the photometry of images A1 and A2 and B and C, respectively. We also allowed for 0.5 mag of uncertainty in the intrinsic flux ratios between the lensed images, since both the continuum and the BLR are affected by microlensing in this lens (Motta et al. 2017) and to allow for the influence of substructure in the lens. We attempted 10 7 fits per set of magnification patterns for a grand total of N trials =4× 10 9 trials, requiring approximately two weeks of run time on the US Naval Academy High Performance Computing (HPC) cluster. In Figure 6, we display two of the best fits from our Monte Carlo analysis to the time-delay corrected difference light curves of WFI2033. Consistent with the findings of Giannini et al. (2017), we easily see ∼0.2 mag of microlensing variability in the difference light curves from image C. The microlensing in images A1, A2, and B is less pronounced, with 0.1 mag of extrinsic variability over the 13 seasons of monitoring.
Using Bayes' theorem, the probability of the parameters given the data D is
where ξ p is the collection of physical variables, ξ t is the collection of trajectory variables, and P(ξ t ) and P(ξ p ) are the prior probabilities for the trajectory and physical variables, respectively. The physical variables are parameters of the local magnification tensors (mean convergence κ and mean shear γ), the local properties of the stars (surface density of stars κ å , mass of the average microlens M * á ñ), the scale radius of the source r s , and the effective velocity of the source v e . The probability distribution for any variable of interest can be obtained by marginalizing over the other variables of the simulation. For example, to find the probability density for the scale radius r s , p p p P r D P D r P P r d , , 3 ) is the probability of fitting the data in a particular trial, p P ( ) sets the priors on the microlensing variables ξ p and ξ t , and P(r s ) is the (uniform) prior on the scale radius. The total probability is then normalized so that P r D dr 1
. We initially do the calculation in "Einstein units," where all lengths depend upon a factor of the unknown mean stellar mass, M M 1 2 * á ñ We construct the prior on the transverse velocity, dP v d v log e e ( ) ( ), using the method described in Kochanek (2004). For the peculiar velocity components of both lens and source, we make redshift-based estimates from the models of Mosquera & Kochanek (2011). We estimate the velocity dispersion of the lens galaxy from its Einstein radius, assuming the galaxy is a singular isothermal sphere with relaxed dynamics, which Treu & Koopmans (2004) and Bolton et al. (2008) show is a good approximation.
We display the probability density for the scale radius in physical units r s in Figure 8. In this plot, we also show a probability density for r s obtained by assuming a uniform prior on the median lens galaxy stellar mass of
A brief inspection of the plot reveals that the results without a prior on the microlens mass M * á ñ are robustly consistent with the results using the uniform mass prior. As a final step we must correct the scale radius for the disk's inclination i by multiplying by i cos 1 2 -(
) , which is necessary because we have assumed a face-on disk in our simulations and microlensing amplitudes depend on the projected area of a source rather than the shape. We adopt the measurement made without the mass prior r i log cm cos 0.5
at λ rest =2481 Å, where i is the inclination angle.
In Figure 9, we plot the size of the accretion disk in WFI2033 on the Accretion Disk Size-Black Hole Mass Relation (Morgan et al. 2010) assuming an inclination angle i=60°, where we have corrected the scale radius r s at the wavelength corresponding to the center of the rest-frame ˆ. The thinner curve indicates the prior probability distribution for the true source velocity v e , which we construct using the method described in Kochanek (2004). Right: the convolution of the prior on v e with the probability density for v e ˆyields the probability density for M * á ñ, which we use to convert the source size measurement into physical units, independent of a prior on M * á ñ. Our estimate for the median lens galaxy stellar mass is M M 0.08 0.05 0.36 ˆ. Right: relative probability density for the thin disk scale radius r s in physical units for WFI2033 at λ rest =2632 Å. The solid line represents the probability density arising from the microlensing simulations using the prior on the true effective velocity v e , while the dotted line shows the result of imposing a uniform prior on the mean microlens mass of
(2017) result is strongly dependent upon priors, especially the assumption of a median microlens mass M M 0.3 * á ñ = ; nevertheless, the independence of the techniques provides robust support for our result. Blackburne et al. (2011) also estimated the scale radius of the accretion disk in WFI2033 using single epoch, multiwavelength photometry, but their results are, by their own admission, anomalous, as they also predict with highest likelihood an accretion disk with an inverted (increasing toward the outer edge) temperature profile. Like Motta et al. (2017), our size measurement is significantly smaller than the Blackburne et al. (2011) estimate.
In Figure 9, we made several updates in addition to the new WFI2033 measurement. Measurements of the accretion disk scale radius using the microlensing variability technique of Q0957+561 (Hainline et al. 2012), SBS0909+532 (Hainline et al. 2013) were added, and updates to the QJ0158-4723 (Morgan et al. 2012), andSDSS0924+0219 (MacLeod et al. 2015) measurements were also included. With changes to 2 out of the 11 existing points and the addition of 3 new measurements, the updated Accretion Disk Size-Black Hole Mass Relation (Morgan et al. 2010) ). There are now 14 systems in which the accretion disk size has been measured using the microlensing variability (e.g., Kochanek 2004) technique. With the exception of SBS0909+532, in which the luminosity-based size estimate is marginally larger than the microlensing-based size measurement, microlensingbased size measurements are consistently larger than the luminosity-based thin disk size estimates by an average of ∼0.55 dex. The SBS0909+532 luminosity-based thin disk size estimate is somewhat suspect, however, since Sluse et al. (2012) and Lehár et al. (2000) found very different photometric fits for the lens galaxy in this system, leading to significant uncertainty in the magnification and, consequently, the intrinsic luminosity. Very recent continuum emission region reverberation mapping studies in local, lower luminosity AGNs (e.g., McHardy et al. 2016;Edelson et al. 2017;Fausnaugh et al. 2018) have revealed similar discrepancies between observed accretion disk size measurements and the predictions of thin disk theory. In Morgan et al. (2010), we proposed that real accretion disks lack the necessary surface brightness to produce their observed luminosity from the smaller area of a simple thin disk model, and we remain confident in that conclusion. We were nevertheless intrigued to find that the slope of r s versus M BH is remarkably consistent with the predictions of thin disk theory (r M s BH 2 3 µ ), so it is the intercept in the accretion disk size black hole mass relation that is inconsistent with thin disk theory, rather than the slope. In Morgan et al. (2010), we suggested that the most promising explanation for the discrepancy is that accretion disks may have shallower temperature slopes than that predicted by thin disk theory T(r)∝r 3/4 , and this hypothesis remains fully viable. We are hopeful that our ongoing lensed quasar monitoring campaign in the infrared (J-, H-, and K-bands), corresponding to optical emission in the rest frame, will allow for measurements of accretion disk temperature profiles in the near future.
The Astrophysical Journal, 869:106 (14pp), 2018 December 20Morgan et al.