SIZES AND TEMPERATURE PROFILES OF QUASAR ACCRETION DISKS FROM CHROMATIC MICROLENSING^*

Microlensing of the images of gravitationally lensed quasars by stars in the foreground galaxies causes time-variable anomalies in their flux ratios (see the review by Wambsganss 2006). In-depth study of these anomalies offers a unique glimpse into the properties of both the deflector and the quasar itself. Because the amplitude of microlensing anomalies depends on the angular size of the source, it can be used to estimate the size of quasar emission regions at the level of microarcseconds, far beyond current resolution limits (e.g., Rauch & Blandford 1991). The temperature structure of quasar accretion disks on these same scales causes chromatic microlensing, where blue light from the inner regions is more strongly microlensed than red light from farther out (Wambsganss & Paczynski 1991). Thus, multiwavelength observations of quasar microlensing have a unique capability to measure the structure of accretion disks and put direct constraints on accretion models.

Several recent studies have made use of microlensing to explore the size scale of the accretion disk, which emits most of the optical light. Pooley et al. (2007) follow up on earlier results (Blackburne et al. 2006; Pooley et al. 2006) to show that the flux ratio anomalies caused by microlensing are in general stronger in X-rays than at optical wavelengths, and infer that the source of optical light is larger than the X-ray source, and larger on average than expected for a standard thin accretion disk model (Shakura & Sunyaev 1973; Novikov & Thorne 1973). Morgan et al. (2010) estimate disk sizes for several lensed quasars using time-series observations of microlensing variability, rather than single-epoch flux ratio anomalies. They find only marginally better agreement with thin disk models, and note that flux-based size predictions are systematically too small.

Several other studies move beyond single-wavelength size estimates to estimate the dependence of the source size on wavelength and thus constrain the temperature profile of the disk. Bate et al. (2008) and Floyd et al. (2009) place upper limits on the radii of accretion disks using single-epoch broadband flux ratio measurements; they also put weak constraints on the temperature profile slope. Mosquera et al. (2010) make size and slope estimates using narrowband photometry. Other studies combine multiwavelength observations with time-series light curves (Anguita et al. 2008; Poindexter et al. 2008; Mosquera et al. 2009; Poindexter & Kochanek 2010), or even employ spectroscopic monitoring (Eigenbrod et al. 2008). Additionally, efforts have been made (Morgan et al. 2008; Dai et al. 2010) to estimate the size of the source of the non-thermal X-ray flux by combining Chandra measurements with optical light curves. Because of the effort involved in these observations, each of these studies focuses on an individual lensed quasar.

The effects of microlensing also depend on properties of the foreground galaxy—the mass function of the stars, and the surface density of clumpy stellar matter as a fraction of the total (smooth and clumpy) surface density. The amplitude of microlensing perturbations decreases for very small stellar fractions, and also decreases, perhaps counterintuitively, for large stellar fractions (Schechter & Wambsganss 2002). This has led to results that effectively rule out the 100% stars case as well as the 100% smooth matter case, with a best-fit stellar mass fraction of ∼10% (Schechter & Wambsganss 2004; Pooley et al. 2009). The mass function of the microlens stars changes the variability characteristics of microlensing light curves; Poindexter & Kochanek (2010) exploit this to estimate the mean stellar mass in the lensing galaxy of the quadruple lens Q2237+0305.

With the goal of measuring accretion disk sizes and temperature profiles, we present single-epoch observations of the flux ratios of a sample of 12 lensed quasars in eight broadband optical and infrared (IR) filters in the "big blue bump" region of the spectrum. The low cost of our single-epoch observing strategy relative to long monitoring campaigns allows us to cover a factor of six in wavelength, while roughly tripling the number of lensed quasars with observations of this kind. We also present X-ray flux ratios from the Chandra X-ray Observatory for 5 of the quasars. Comparing these flux ratios, as well as X-ray ratios from the literature for six of the seven remaining lenses, to those predicted by smooth lens models, we use a Bayesian analysis method to determine posterior probability distributions for the sizes of the quasar emission regions in each filter. This allows us to estimate the overall size scales of the accretion disks, as well as the scaling of the size with black hole mass and wavelength. As we will see, the relatively large size of our sample is crucial for drawing conclusions about accretion disk structure.

Our sample is entirely made up of quadruply lensed quasars; they have the advantage of having three flux ratios apiece (rather than one for double lenses) and are easier to model. They are as follows:

HE 0230−2130 (Wisotzki et al. 1999),

MG J0414+0534 (Hewitt et al. 1992),

HE 0435−1223 (Wisotzki et al. 2002),

RX J0911+0551 (Bade et al. 1997),

SDSS J0924+0219 (Inada et al. 2003),

HE 1113−0641 (Blackburne et al. 2008),

PG 1115+080 (Weymann et al. 1980),

RX J1131−1231 (Sluse et al. 2003),

SDSS J1138+0314 (Eigenbrod et al. 2006),

SDSS J1330+1810 (Oguri et al. 2008),

WFI J2026−4536 (Morgan et al. 2004), and

WFI J2033−4723 (Morgan et al. 2004).

For brevity, we will usually exclude the declination from their names from now on. Following Pooley et al. (2007), we categorize the quasar images by their macro-magnification and parity: highly magnified minima and saddle points of the light travel time surface are designated HM and HS, while their less-magnified counterparts are LM and LS, respectively.

The variation of accretion disk size with black hole mass and wavelength provides direct constraints on disk models. The standard thin disk model predicts a profile for the effective temperature that falls approximately as the β = 3/4 power of the radius:

$$\begin{equation} T_{\rm eff}(r) = \left(\frac{3 G^2 M_{\rm BH}^2 m_p f_{\rm Edd}}{2 c \sigma _B \sigma _T \eta r^3}\right)^{1/4} g(r_{\rm in}/r)^{1/4}. \end{equation} \tag{ 1 }$$

where M_BH is the black hole mass, m_p is the proton mass, σ_B is the Stefan–Boltzmann constant, and σ_T is the Thomson scattering cross-section. The dimensionless factors f_Edd and η are the ratios of the quasar's bolometric luminosity to its Eddington luminosity and to the accretion rate $\dot{M}c^2$, respectively. The function g is defined as g(x) = 1 − x^1/2, and r_in is the innermost radius of the accretion disk. As we will see, at optical wavelengths the disk size is much larger than r_in, so for simplicity we set g to unity. Likewise, the general relativistic corrections of Novikov & Thorne (1973) are small at the radii of interest. By integrating the specific flux at fixed wavelength over the entire multi-temperature disk, we calculate the half-light radius—that is, the radius at which half of the light at a given wavelength is emitted—predicted by the thin disk temperature profile:

$$\begin{eqnarray} r_{1/2} &= 2.44 \nonumber \left[\frac{45 G^2 M_{\rm BH}^2 m_p f_{\rm Edd} \lambda ^4}{4 \pi ^5 h_P c^3 \sigma _T \eta }\right]^{1/3} \sqrt{\cos {i}} \\ &= 1.68 \times 10^{16} {\rm cm} \left(\frac{M_{\rm BH}}{10^9\, M_\odot }\right)^{2/3} \left(\frac{f_{\rm Edd}}{\eta }\right)^{1/3} \left(\frac{\lambda }{\mu {\rm m}}\right)^{4/3}. \end{eqnarray} \tag{ 2 }$$

In this equation, h_P is Planck's constant, and i is the unknown inclination angle. We include the factor of (cos i)^1/2 to account for the inclination of the disk to our line of sight and set it to a reasonable average value of (1/2)^1/2 in the second line. For a thin disk, the half-light radius scales with the rest wavelength as λ^4/3; for a more general disk temperature profile T_eff ∝ r^−β, the half-light radius will depend on λ^1/β.

In Sections 2 and 3, we describe the X-ray and optical/IR observations and flux ratio measurements. In Section 4, we estimate the systematic uncertainties in the measured flux ratios. We describe our lens models in Section 5, and in Section 6 we describe our Bayesian analysis method. We discuss our results and give concluding remarks in Sections 7 and 8. Throughout this work, we calculate distances and time delays using a geometrically flat universe with Ω_M = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹.

The X-ray data in this survey are from the Advanced CCD Imaging Spectrometer (ACIS) instrument aboard the Chandra X-ray Observatory. For five of the lenses among our sample (specifically HE 0230, MG J0414, RX J0911, SDSS J0924, and WFI J2033) we adopt the X-ray flux ratios reported by Pooley et al. (2007). For PG 1115 we adopt the X-ray flux ratios of the 2008 January 31 observation (ObsID 7757) from Pooley et al. (2009), as they are more contemporaneous with our optical/IR observations. For the same reason, for RX J1131 we use observations from the Chandra archive. Specifically, we use the 2007 February 13 observation (ObsID 7787) for comparison with optical data and the 2007 April 16 observation (ObsID 7789) for comparison with IR data (see Section 3). The high level of X-ray microlensing variability in these two systems in particular (Chartas et al. 2009; Pooley et al. 2009) underscores the importance of making our multiwavelength observations as contemporaneous as possible (see Section 4.2 for a detailed discussion).

The remaining lenses in our sample, excluding SDSS J1330, which has not been observed in X-rays, were observed using ACIS between 2006 December 17 and 2007 June 28. The data were reduced using the procedure described by Pooley et al. (2009). We produced 0.5–8 keV images of each lens with a resolution of 00246 pixel⁻¹, and fit four Gaussian components and a constant background to each image. The background level was fixed at a level determined from a nearby source-free region, and the relative positions of the four Gaussians were fixed to the Hubble Space Telescope (HST) positions provided by the CASTLES database⁴ (Falco et al. 2001). The width of the Gaussian point-spread function (PSF) was allowed to vary for each lens, but was constrained to be the same for all four quasar components. The fits were performed using Cash (1979) statistics, which are appropriate when small numbers of photons have been collected. In the cases of HE 1113, SDSS J1138, and WFI J2026, the least-squares optimization settled on a wider PSF than for any of the other observations in our previous work. Given the low numbers of total counts in these observations and the small separations of the images, we fixed the width of the PSF in these cases to the average of the other X-ray PSFs. The flux in image A2 of WFI J2026 was very sensitive to the PSF, so we allocated a generous 1.5 mag of uncertainty to its flux ratio, added in quadrature to the other sources of uncertainty (see Section 4.2). The other flux ratios were not sensitive to the PSF size. Table 1 summarizes the Chandra observations and the measured flux ratios.

Table 1. X-ray Fluxes and Flux Ratios

Quasar			Image Flux Ratiosa			LM Unabs. F0.5–8 keVb
ObsIDc	Date	Exposure (s)	HS/LM	HM/LM	LS/LM	(10−14 erg cm−2 s−1)
HE 0435−1223			B/A	C/A	D/A	FA
7761	2006 Dec 17	10 130	0.375$^{_{+0.049}}_{^{-0.045}}$	0.378$^{_{+0.049}}_{^{-0.044}}$	0.363$^{_{+0.048}}_{^{-0.043}}$	17.75$^{_{+4.43}}_{^{-4.40}}$
HE 1113−0641			D/A	B/A	C/A	FA
7760	2007 Jan 28	15 180	0.78$^{_{+0.32}}_{^{-0.32}}$	0.63$^{_{+0.41}}_{^{-0.41}}$	0.20$^{_{+0.16}}_{^{-0.16}}$	7.8$^{_{+2.3}}_{^{-2.3}}$
RX J1131−1231d			A/C	B/C	D/C	FC
7787	2007 Feb 13	5 190	5.717$^{_{+0.417}}_{^{-0.382}}$	3.176$^{_{+0.244}}_{^{-0.224}}$	0.624$^{_{+0.064}}_{^{-0.058}}$	45.34$^{_{+3.99}}_{^{-3.88}}$
7789	2007 Apr 16	5 190	5.580$^{_{+0.391}}_{^{-0.359}}$	2.922$^{_{+0.218}}_{^{-0.200}}$	0.374$^{_{+0.044}}_{^{-0.040}}$	47.44$^{_{+4.14}}_{^{-4.03}}$
SDSS J1138+0314			D/C	A/C	B/C	FC
7759	2007 Feb 13	19 080	1.3$^{_{+0.6}}_{^{-0.5}}$	3.2$^{_{+1.0}}_{^{-1.0}}$	1.0$^{_{+0.4}}_{^{-0.4}}$	1.6$^{_{+0.9}}_{^{-0.9}}$
WFI J2026−4536			A2/B	A1/B	C/B	FB
7758	2007 Jun 28	10 170	2.0$^{_{+1.7}}_{^{-1.0}}$	5.8$^{_{+1.8}}_{^{-1.7}}$	0.40$^{_{+0.06}}_{^{-0.15}}$	5.3$^{_{+1.2}}_{^{-1.2}}$

Notes. We only report new measurements here; see references in the text for other flux ratios. ^aHS: highly magnified saddle point; HM: highly magnified minimum; LS: less magnified saddle point; LM: less magnified minimum. ^bThe unabsorbed flux of the LM image is computed from the best-fit power-law model described in Section 2. ^cThe observation identifier of the Chandra data set. ^dThe data for RX J1131−1231 come from the Chandra archive.

Download table as:  ASCIITypeset image

Between 2007 February and 2008 May, we undertook an optical observing campaign to obtain multi-band, contemporaneous images of our sample of lenses. At IR wavelengths, we used the J, H, and K_s filters with Persson's Auxiliary Nasmyth Infrared Camera (PANIC), on the 6.5 m Baade telescope at Las Campanas Observatory. In the optical, we used the Sloan u'g'r'i'z' filters with the Raymond and Beverly Sackler Magellan Instant Camera (MagIC), which was on the neighboring Clay telescope when we began our campaign, but was moved to Baade while the observations were still underway. This change of location allowed us to use both instruments during the same observing run, and even during the same night. MagIC and PANIC have fields of view 24 and 21 on a side, respectively, large enough to include stars appropriate for use as PSF templates. The instruments' pixel scales are 0069 and 0125 pixels, which more than adequately sample the PSF. The details of the observations are listed in Table 2.

Table 2. Optical Observations

Quasar	Date	Instrument	Filters	Exposures (s)	Seeing
HE 0230	2007 Sep 16	IMACS	u'g'r'i'z'	720; 240; 240; 240; 480	043
	2007 Jul 29	PANIC	JHKs	540; 540; 540	064
MG J0414	2007 Sep 21	IMACS	r'i'z'	360; 240; 480	067
	2007 Sep 22	PANIC	JHKs	810; 720; 720	076
HE 0435	2007 Sep 16	IMACS	u'g'r'i'z'	720; 240; 240; 270; 480	067
	2007 Sep 22	PANIC	JHKs	810; 1080; 1080	106
RX J0911	2007 Feb 13	MagIC	g'r'i'z'	240; 240; 300; 480	070
	2007 Feb 14	MagIC	u'	720	063
	2007 Apr 7	PANIC	JHKs	540; 405; 405	072
SDSS J0924	2007 Feb 13	MagIC	u'g'r'i'z'	720; 240; 240; 180; 480	061
	2007 Apr 7	PANIC	JHKs	540; 540; 540	064
HE 1113	2007 Feb 14	MagIC	g'r'i'z'	240; 180; 120; 240	049
	2007 Apr 7	PANIC	JHKs	450; 450; 450	045
PG 1115	2008 Feb 1	MagIC	g'r'i'z'	120; 120; 150; 240	055
	2008 Feb 2	MagIC	u'	720	087
	2008 Feb 4	PANIC	JHKs	540; 540; 540	056
RX J1131	2007 Feb 13	MagIC	u'g'r'i'z'	480; 120; 120; 120; 240	065
	2007 Apr 6	PANIC	JHKs	486; 405; 405	076
SDSS J1138	2007 Feb 13	MagIC	u'g'r'i'z'	1080; 480; 480; 360; 960	075
	2007 Apr 8	PANIC	JHKs	540; 540; 540	047
SDSS J1330	2008 Feb 1	MagIC	u'g'r'i'z'	360; 120; 120; 150; 240	053
	2008 Feb 3	PANIC	JHKs	540; 540; 540	068
WFI J2026	2008 May 13	MagIC	u'g'r'i'z'	720; 240; 240; 360; 480	054
	2008 May 12	PANIC	JHKs	135; 135; 90	043
WFI J2033	2007 Jun 15	MagIC	u'g'r'i'z'	720; 240; 120; 300; 480	088
	2007 Jul 4	PANIC	JHKs	135; 135; 90	051

Download table as:  ASCIITypeset image

During our 2007 September observing run, MagIC was offline for an upgrade, so we instead used the Inamori Magellan Areal Camera System (IMACS) in its imaging mode. With its f/4 camera, the instrument has a pixel scale of 011 and a 155 field of view, but in order to reduce the readout time we used only a subraster 22 on a side. We used this instrument for our observations of HE 0230, MG J0414, and HE 0435.

The images were corrected for bias and dark current, flattened, and cleaned of cosmic rays using standard techniques. Where multiple exposures were obtained, they were stacked, resulting in a single image per filter per lensed quasar. Postage-stamp i'-band images of the 12 lenses in our sample are provided in Figure 1, together with the residuals from the PSF fitting described in Section 3.

Zoom In Zoom Out Reset image size

Because of the small separations between the quasar images in our sample, we only observed under the best atmospheric seeing conditions, typically below 07. The excellent image quality at the Magellan telescopes was crucial to the success of our survey. The i'- and J-band seeing is reported in Table 2, which outlines the optical/IR observations obtained for this work. The seeing in the other filters was roughly consistent with the rule of thumb: FWHM ∝ λ^−1/5.

Because we are primarily interested in the flux ratios of lensed quasars (indeed, we cannot measure absolute magnifications), we did not, in general, obtain images of standard stars for the purpose of calibrating our photometry. In the case of SDSS J1330, however, we did observe standard stars in the optical bands, so in Section 3 we report calibrated photometry for this lens. We also estimate the photometric zero point in the J and K_s bands, based on images taken the same night which contain several calibrated sources from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006).

Despite the high quality of our data, careful PSF subtraction was necessary to disentangle the compact clumps of point sources and lens galaxies that make up our sample. For each image, we used a nonlinear least-squares fitting algorithm to simultaneously fit the positions and relative fluxes of the four quasar components and the lens galaxy.

For lenses with well-separated components, we performed a simultaneous fit to the images in all eight bands, with the relative positions of the quasar components and the lens galaxy allowed to vary, but constrained to be consistent across all eight filters. The flux ratios of the quasar components were free to vary in each filter independently. The best-fit positions were consistent in every case with the HST positions from the CASTLES survey. In cases where merging pairs or triplets of images would have caused strong correlations between positions and flux ratios, we fixed the relative positions of the quasar components and the lensing galaxy to the CASTLES positions. We modeled the quasar images using nearby relatively bright stars as empirical PSF templates.

In the r' and redder bands, we did not use least-squares minimization to determine the flux ratios of RX J1131. Because of the bright Einstein ring in this system, this strategy would have oversubtracted the quasar components. Instead we fixed all of the quasar fluxes to values that resulted in residuals that looked like an unbroken Einstein ring. Although it was completely ad hoc, we think that this "chi by eye" technique gave flux ratios less affected by systematic errors. Section 4.1 describes our method for estimating the uncertainty of these flux ratios.

In all cases we modeled the lensing galaxy as a two-dimensional pseudo-Gaussian function (Schechter et al. 1993), with a full-width at half-maximum as a fixed parameter. The widths were chosen using a trial-and-error technique, examining the residuals by eye. (Because of the frequent presence of a faint Einstein ring due to the quasar's host galaxy, least-squares minimization often overestimated the width.) Because the lens galaxy was nearly always faint compared to the quasar images, the goodness of the fit was insensitive to this approximation. The galaxies were mostly round, but for the few exceptions we fixed the axis ratio and position angle by hand, as we did the width.

The fitted fluxes of the lens galaxies are uncalibrated, and therefore of little interest; we do not report them. SDSS J1330 is the exception; because it was discovered relatively recently, and because observing conditions were photometric on the night it was observed, we observed standard stars in the optical bands, for the purpose of calibrating its fluxes. The calibrated magnitudes of this lens galaxy are reported in Table 4. The optical calibration relies on aperture photometry of the standard stars SA 107-351 (for u') and G163-50 (for g'r'i'z'), and places the fluxes on the Sloan u'g'r'i'z' photometric system, which closely approximates the monochromatic AB system (Smith et al. 2002). In the H and K_s bands, we made use of a field containing several 2MASS sources for photometric calibration; we performed aperture photometry on them and compared their measured magnitudes to those reported in the NASA/IPAC Infrared Science Archive⁵ to obtain photometric zero points. This places our calibrated J and K_s fluxes on the 2MASS photometric system (Cohen et al. 2003). In the optical bands, used the atmospheric extinction coefficients of Smith et al. (2007) to correct for differences in airmass between SDSS J1330 and the standard stars; these corrections were all less than 0.04 mag.

The lack of HST positions for SDSS J1330, combined with the relative brightness and elliptical shape of the lensing galaxy, made its fit a special challenge. In order to get the shape of the galaxy right, we adopted an iterative approach. First, we subtracted a rough model of the lens galaxy, fixing its parameters to reasonable guesses. We then fit the positions and fluxes of the quasar components, subtracting them from the image. We added back our rough galaxy model to the residual image and fit a single pseudo-Gaussian to the resulting galaxy image, with all parameters allowed to vary. The best-fit FWHM of the galaxy was 10 along the major axis; the axis ratio and position angle were 0.67 and 245 east of north. Finally, we fit the original image again, fixing the galaxy shape parameters to these best-fit values, in order to find the quasar positions and fluxes. The results are consistent with those of Oguri et al. (2008), but have higher precision; they are reported in Tables 3 (astrometry) and 4 (photometry).

Table 3. Relative Astrometry for SDSS J1330+1810

Component	x	y
A	−1.26 ± 0.01	−1.17 ± 0.01
B	−0.84 ± 0.01	−1.19 ± 0.01
C	≡0	≡0
D	−1.49 ± 0.01	$\phantom{-}0.44 \pm 0.01$
G	−1.04 ± 0.02	−0.20 ± 0.02

Notes. All positions are in arcseconds. Positive x and y point west and north, respectively.

Download table as:  ASCIITypeset image

For each lens, the transformation between pixel coordinates and sky coordinates was determined by registering our i'- and J-band images to the USNO-B astrometric catalog. On average, there were ∼10–15 objects for each fit. The pixel scale and rotation angle thus determined for i' and J were assumed to apply to the remaining MagIC/IMACS data and PANIC data for that lens, respectively. We tested this assumption for a few cases by fitting images in other bands to the catalog and found it to be justified.

Relative photometry for all the lenses in our sample except SDSS J1330 is listed in Table 5. We also list the calculated rms deviation of the magnitudes in each filter from those predicted by the lens models. This number is an indication of how anomalous the flux ratios are. For comparison with the optical flux ratios, we also list the X-ray ratios, in magnitudes.

Aside from statistical measurement uncertainties in the flux ratios, systematic errors arise from several sources. At optical and IR wavelengths, quasar emission lines, quasar variability, or contamination from a lens galaxy or an Einstein ring may contribute. In X-rays, quasar variability is a factor, and delays between X-ray and optical/IR observations combine with microlensing variability to contribute additional systematic uncertainty.

4.1. Optical Uncertainties

4.2. X-ray Uncertainties

Like the optical flux ratios, the X-ray ratios had measurement uncertainties due to measurement noise and the blending of close pairs; both errors were generally larger because the X-ray observations have fewer photons and a broader PSF. There were also contributions from intrinsic quasar variability, as in the optical case. The X-ray ratios had no systematic errors due to emission lines, but they did have errors due to microlensing variability because they were not generally measured contemporaneously with the other wavelengths.

Crucial to our analysis is the assumption that the arrangement of the source and the microlenses is the same for all wavelengths. But when observations are not contemporaneous, the source and the microlenses have the opportunity to reconfigure themselves. To estimate the magnitude of this effect, we again used a structure function; this time it was not empirical, but was derived from microlensing magnification patterns (see our description of these patterns in Section 6.1). For each image, we chose the appropriate pattern (with the stellar mass fraction k_* set to 0.1; see Section 6.1) and ran 1000 tracks across it, in random directions, measuring the variance in the logarithm of the magnification as a function of distance moved.

The conversion to a structure function (with a time delay on the abscissa instead of a distance) required an estimate of the transverse speed of the source relative to the lens. We added four velocity components in quadrature: the velocity dispersion of the stars in the lens galaxy (as estimated by our lens model), the tangential component of the velocity of the Sun relative to the rest frame of the cosmic microwave background (CMB), and the peculiar velocities of the quasar and the lens galaxy (as estimated using Equation (14.10) of Peebles 1980):

$$\begin{eqnarray} v_{\perp }^2 = \nonumber &2\left(\frac{\sigma _L}{(1 + z_L)} \frac{D_{\rm OS}}{D_{\rm OL}}\right)^2 +\left(v_{\rm CMB} \sin {\alpha } \frac{D_{\rm LS}}{D_{\rm OL}}\right)^2 \\ &+2\left(\frac{\sigma _{\rm pec}}{(1 + z_S)^{3/2}} \frac{f(z_S)}{f(0)}\right)^2 +2\left(\frac{\sigma _{\rm pec}}{(1 + z_L)^{3/2}} \frac{f(z_L)}{f(0)}\frac{D_{\rm OS}}{D_{\rm OL}}\right)^2.\nonumber\\ \end{eqnarray} \tag{ 3 }$$

Here z_S and z_L are the redshifts of the source and the lens, respectively, and D_OL, D_OS, and D_LS are angular diameter distances from observer to lens, observer to source, and lens to source. We used these distances to project all velocities to the source plane. The angle α is measured between the Sun's velocity with respect to the CMB rest frame and the line of sight to the lens, and we set v_CMB to 370 km s⁻¹ (Lineweaver et al. 1996). We followed Kochanek (2004) in setting the present-day velocity dispersion of galaxies σ_pec to 235 km s⁻¹. The function f(z) is the cosmological growth factor; we approximate it as f ∝ Ω_M(z)^0.6. The stellar velocity dispersion σ_L of the lens we estimate from the monopole component b of our lens model (see Section 5) using the relation (e.g., Kochanek 2006)

$$\begin{equation} b = 4\pi \frac{\sigma _L^2}{c^2} \frac{D_{\rm LS}}{D_{\rm OS}}. \end{equation} \tag{ 4 }$$

All velocities are corrected for cosmological time dilation. The factors of two convert one-dimensional velocity dispersions to two dimensions. Multiplying the estimated source speed v_⊥ by the delay between our optical/IR observations and their X-ray counterparts, we determined the distance traveled by the source along the caustic pattern, and from the structure function read off the predicted error in our X-ray flux ratio.

For each X-ray flux ratio, we added in quadrature the uncertainty contributions from measurement errors (including extra uncertainty for image A2 of WFI J2026 because of its degeneracy with the PSF width), intrinsic quasar variability, and microlensing variability. These uncertainties are given in Table 5. In general, the largest contributions were systematic uncertainties from microlensing variability, though the statistical measurement errors were also substantial in some cases. As in the optical case, the uncertainty due to quasar variability was relatively insignificant.

We used the lensmodel program of Keeton (2001) to create parametric models of each lens. The models were constrained by the positions of the four lensed images, and that of the lensing galaxy, a total of 10 constraints. We did not use fluxes for constraints, since most of our lenses suffer from flux ratio anomalies. Nor did we use time delays. The positions for the quasar images and lens galaxies came from the CASTLES survey.

Our default model consisted of a singular isothermal sphere for the lensing galaxy, with a quadrupole component of the potential provided by a constant external shear. With the position of the lens fixed, this model had five free parameters: the monopole strength of the lens, the magnitude and direction of the shear, and the two-dimensional position of the source.

In some cases, such as HE 0435, this simple model fits the image positions well. But in several cases the χ² goodness of fit was poor enough to warrant further complexity in the model. In these cases, we made changes to the model motivated by the appearance of the lens galaxy.

The lens galaxy of HE 0230 has a prominent companion, located close to image D. We modeled this companion as a second isothermal sphere, fixing its position to its HST measured value but allowing its mass to vary. Despite only adding one free parameter to the model, this addition improved the fit considerably. We followed a similar strategy for MG J0414, RX J0911, and WFI J2033, each of which displays a faint smudge in its HST image which is arguably a satellite to the lens galaxy. Adding secondary lenses at the positions of the smudges improved the goodness of fit to acceptable levels.

The lens galaxy in PG 1115 does not have a nearby companion, but is a member of a small galaxy group centered to the southwest of the lens. Following the example of Schechter et al. (1997), we explicitly modeled the group as a second isothermal sphere for this lens. We parameterized its position using polar coordinates, and allowed its mass and distance from the main galaxy to vary while fixing its position angle to that of the brightest galaxy in the group. We did not include an external shear in this fit.

Finally, the lens galaxy in SDSS J1330 displays significant ellipticity. In this case, we used an isothermal ellipsoid instead of including an external shear. We allowed the ellipticity and position angle to vary, along with the galaxy position, for a total of seven free parameters. Since there are no HST data for this lens, we used our measured image positions for constraints (see Section 3 and Table 3). SDSS J1330 was the only lens system where an isothermal ellipsoid made for a better fit than an isothermal sphere with external shear.

The salient features of our best-fit models are listed in Table 6. We used the components of the magnification tensor predicted by these models (namely, the local convergence κ and shear γ) to generate the microlensing magnification patterns described in Section 6.1. We also used the predicted scalar magnifications to estimate the microlensing-free flux ratios, and assumed that any discrepancies between the predicted and observed flux ratios was due to microlensing. We believe that this is a good assumption for our sample of lenses because there is not strong evidence for other potential explanations. For example, differential extinction is unlikely to be a factor in these early-type lens galaxies (Kochanek & Dalal 2004); this is supported by the general lack of strong hydrogen absorption in X-rays (Pooley et al. 2007). Likewise, the chromatic variations we see in the flux ratios argues against millilensing by dark matter substructure. In the case of MG J0414, however, mid-IR flux ratios measured by Minezaki et al. (2009) disagree with the predicted ratios. Since the dusty torus region where this light originates is too large to be affected by microlensing, it is likely that millilensing is affecting this lens. In this case, we adopted the reported mid-IR flux ratios as the "baseline" rather than the model ratios, attributing any optical departures from them to microlensing.

Table 6. Lens Model Parameters

Quasar	Primary Lens			Secondary Lens			Magnificationa
	b	γ	ϕγb	b2	x2c	y2c	HM	HS	LM	LS
HE 0230	087	0.11	−600	033	-0.283	+0.974	+9.42	−9.65	+4.95	−1.35
MG J0414	114	0.11	+746	012	-0.385	+1.457	+22.9	−24.2	+6.23	−3.11
HE 0435	120	0.078	−138	...	...	...	+7.49	−7.90	+7.14	−4.73
RX J0911	097	0.27	+72	024	-0.754	+0.665	+11.0	−5.96	+1.97	−4.99
SDSS J0924	087	0.063	+848	...	...	...	+14.9	−13.0	+6.62	−6.55
HE 1113	033	0.040	+377	...	...	...	+15.8	−16.7	+12.6	−9.59
PG 1115	103	...	...	257	−10.866	−5.300	+19.7	−18.9	+5.09	−3.37
RX J1131	186	0.16	−736	...	...	...	+13.2	−22.7	+12.6	−1.05
SDSS J1138	067	0.10	+326	...	...	...	+7.17	−6.68	+5.17	−3.64
SDSS J1330	094	0.16d	−322d	...	...	...	+27.1	−27.2	+8.41	−5.50
WFI J2026	066	0.11	−900	...	...	...	+13.7	−11.5	+3.78	−4.01
WFI J2033	107	0.11	+360	025	+0.229	+2.020	+6.04	−3.80	+3.88	−2.46

Notes. ^aHS: highly magnified saddle point; HM: highly magnified minimum; LS: less-magnified saddle point; LM: less-magnified minimum. ^bPosition angle of external shear γ. All angles are measured in degrees east of north. ^cFixed position of secondary galaxy, relative to main lensing galaxy, in arcseconds. Allowed to vary radially in the case of PG 1115; see the text. The positive directions of x and y are west and north, respectively. ^dEllipticity and position angle of the lens galaxy for the SIE model of SDSS J1330.

Download table as:  ASCIITypeset image

The predicted flux ratios are model dependent to some extent; this systematic error is difficult to quantify, but ought to be small compared to the other uncertainties. Dalal & Kochanek (2002) suggest error bars of 10%; Keeton et al. (2003) call this "quite conservative." We estimated the model uncertainties in the flux ratios as 0.05 mag for well-separated quasar images, or 0.03 mag for members of a close pair or triplet, since these configurations have model-independent asymptotic flux ratios (Keeton et al. 2003, 2005, and references therein). We added this uncertainty in quadrature to those listed in Table 5.

The effect of a finite source size is to smooth out the temporal variations caused by microlensing. The movement of a large source across a network of microlensing magnification caustics is equivalent to a point source moving across a version of the caustic map that has been smoothed by convolution with the surface brightness profile of the source. Since microlensing variability is stronger in X-rays than at optical wavelengths, it is thought that the X-rays come from a very compact region, whereas the optical source is much larger (Pooley et al. 2007, 2009; Morgan et al. 2008; Dai et al. 2010). For single-epoch observations like ours, microlensing magnifications are more naturally treated as flux ratio anomalies, where the flux ratios of the quasar images fail to match those predicted by lens models, than as variability. We expect the flux ratio anomalies to be stronger in X-rays than at optical wavelengths, and stronger in blue light than in red light, most of the time. However, different regions of the magnification maps have different responses to a change in the source size. Indeed, in some regions (e.g., at the center of an astroid caustic), convolution by a large source profile will increase the anomaly. Poindexter et al. (2008) saw an example of this in the doubly lensed quasar HE 1104−1805.

Single-epoch measurements are much less expensive than long-term monitoring campaigns, especially with multiple wavelengths. But they have the disadvantage that they are unable to constrain the local characteristics of the microlensing magnification map, leading to a degeneracy between the natural range of possible magnifications and source size effects. A strong anomaly is evidence for a relatively small source, but a weak or absent anomaly could be due to a large source or to the chance location of a small source on a relatively calm portion of the magnification map. PG 1115 is a real-life example of this: for several years its optical A2/A1 ratio maintained a low-level anomaly that was revealed to be due to a large source size when the same ratio was found to be extremely low in X-rays (Pooley et al. 2006). Our analysis method follows Pooley et al. (2007) in using X-ray flux ratios to break this degeneracy, but in a more quantitative way.

A similar uncertainty exists when using size estimates at a range of wavelengths to estimate the wavelength dependence of the size of quasar accretion disks. Because of the unknown response of the microlensing magnification to changes in source size, a single-epoch observation of a single lensed quasar cannot tell us much about this dependence. We stress, however, that it does allow for estimates with meaningful error bars. The way to decrease this uncertainty is to sample a greater fraction of the caustic map, either via time-series observations of a single lens or by observing a sample of lenses. We have adopted the latter approach.

We have developed a quantitative Bayesian method for estimating the angular size of quasar accretion disks at a range of optical and IR wavelengths using our measured optical/IR and X-ray flux ratios. This approach relates the posterior probability distribution for the source size to the likelihood of observing the reported X-ray flux ratios, assumed to originate in a very compact source, simultaneously with the observed optical/IR ratios from a larger source (see Equation (5)). The bulk of this section is devoted to the calculation of these likelihood distributions.

6.1. Microlensing Magnification Patterns

6.2. Size Estimation Using Optical Ratios

We used a simple Bayesian method to calculate the posterior probability distribution for the half-light radius of the source r_1/2 at each wavelength. Assuming a uniform prior for the source size, the probability is

$$\begin{equation} p(r_{1/2}|\lbrace m_i\rbrace) \propto \mathcal {L}(\lbrace m_i\rbrace |r_{1/2}), \end{equation} \tag{ 5 }$$

where {m_i} is the set of quasar image magnitudes, i runs from 1 to N_im = 4, and $\mathcal {L}$ is the likelihood of having observed the magnitudes {m_i} given a source half-light radius. Each magnitude m_i may be expressed as the sum of the unmagnified quasar magnitude m_S, macro- and micro-magnification terms μ_i and Δμ_i, and measurement noise n_i:

$$\begin{equation} m_i = m_S + \mu _i + \Delta \mu _i + n_i. \end{equation} \tag{ 6 }$$

The source magnitude m_S is not observable, and because we did not flux calibrate our measurements, it is in uncalibrated units. Since we ultimately marginalized over m_S, these limitations are unimportant.

We computed the right-hand side of Equation (5) by comparing the measured image fluxes m_i to histograms of the magnification patterns described in Section 6.1. Suitably normalized, each histogram was an estimate of the distribution of microlensing magnifications:

$$\begin{equation} p_i(\Delta \mu _i | r_{1/2}, k_*), \end{equation} \tag{ 7 }$$

where the pattern sourcing the histogram was created using the stellar mass fraction k_* and was convolved with a source with radius r_1/2.

Figure 2 shows an example of four such magnification histograms, one for each image in SDSS J1330, for a particular source size and stellar mass fraction. The abscissa of this plot denotes the microlensing deviation Δμ in magnitudes. Atop the histograms have been plotted vertical lines, indicating the observed deviation from the model in the i' band. In the case of no microlensing, the lines would be co-located at zero. The presence of four vertical lines misleadingly implies more information than is available; in fact, we have not measured the four image magnifications, but only their three ratios. Thus, the average position of the four lines is unknown. In Figure 2, we have set it to zero, but in reality the four lines are free to slide from side to side in formation. Each horizontal position corresponds to a distinct value of the quasar's unknown magnitude m_S.

Zoom In Zoom Out Reset image size

We shifted each histogram horizontally to obtain p_i(Δμ_i − m_i + μ_i) = p_i(−m_S − n_i). These shifts brought the four vertical lines in Figure 2 together. Because the shifts depended on the measured fluxes, we marginalized over the measurement errors:

$$\begin{equation} p_i(-m_S) = \frac{1}{\sqrt{2\pi }\sigma _i} \int dn_i \exp {\left(\frac{-n_i^2}{2\sigma _i^2}\right)} p_i(-m_S - n_i). \end{equation} \tag{ 8 }$$

This amounted to a convolution of the shifted histograms with the normal distributions implied by the measurement uncertainties σ_i. These shifted and convolved histograms were probability density functions of the quasar magnitude m_S. We multiplied them together and integrated their product to find the likelihood of measuring all the observed flux ratios, independent of the quasar flux:

$$\begin{equation} \mathcal {L}(\lbrace m_i\rbrace |r_{1/2},k_*) = \int dm_S \phi (m_S) \left[\prod _{i=1}^{N_{\rm im}} p_i(-m_S)\right], \end{equation} \tag{ 9 }$$

where N_im = 4 and ϕ(m_S) is the prior on the luminosity of the quasar. For this luminosity function we used a power law of slope −2.95, as estimated for bright quasars by Richards et al. (2006). Our use of a scale-free luminosity function obviated any difficulties due to the uncalibrated quasar flux.

Finally, we marginalized over the stellar mass fraction k_* to find the right-hand side of Equation (5):

$$\begin{equation} \mathcal {L}(\lbrace m_i\rbrace |r_{1/2}) = \int dk_* p(k_*) \mathcal {L}(\lbrace m_i\rbrace |r_{1/2},k_*). \end{equation} \tag{ 10 }$$

For our prior p(k_*) on the stellar mass fraction we interpolated the estimated probability distribution for its value given in Figure 6 of Pooley et al. (2009). It should be noted that although we have written the marginalizations in Equations (8)–(10) as integrals for clarity, they were actually calculated as discrete sums.

We calculated this relative likelihood for 31 half-light radii, spaced logarithmically between 0.01r_Ein and 10r_Ein, where r_Ein is the Einstein radius of a solar-mass microlens. We let the posterior probability of the source size equal this likelihood, implicitly adopting a logarithmic prior on the source size (between these limits; zero elsewhere). For completeness, we also used a linear prior; this simply involved multiplying the likelihood distribution by a power law of unit slope.

Figure 3 shows the posterior distribution resulting from the use of this technique on the i'-band flux ratios of SDSS J1330. This figure demonstrates the ability of this method to rule out large source sizes (e.g., greater than an Einstein radius). The decrease in probability to zero occurs when the growing source size causes two of the shifted histograms to become narrow enough that they no longer overlap at all; at this point the likelihood is zero. Measurement errors widen the histograms via the convolution in Equation (8), causing the probability distribution to widen as well. There is no corresponding decrease in probability for small source radii. This is to be expected, since the optical flux ratio anomalies simply indicate that the source is undergoing microlensing. To set a lower limit on the size, we must discern why the observed anomalies are not stronger than they are—is the caustic map being smoothed out by a large source, or is the source simply lying in a non-perturbed region? Though some size discrimination is possible by comparing the posterior distributions derived from blue filters versus red filters, a more effective method to rule out very small sizes is to compare the optical flux ratios to those in X-rays.

Zoom In Zoom Out Reset image size

6.3. Size Estimation Using Optical and X-ray Ratios

SDSS J1330, the lens we chose for the examples in the previous section, is unique among our sample in that it has never been observed in X-rays. But for the remaining lenses in our sample there is at least one measurement of the X-ray flux ratios. Since the microlensing variability is in general stronger in X-rays than at optical wavelengths, the X-rays must come from a quite compact region (Pooley et al. 2007, 2009; Dai et al. 2010). Any decrease in the anomalies in the optical and IR must be due to the finite size of the source at these wavelengths. Comparing the optical flux ratios to X-ray ratios gives us the ability to set a lower limit on the size of the optical-emitting accretion disk.

We again wished to find the value of the right-hand side of Equation (5), but this time, we wanted the likelihood of observing both the X-ray and optical flux ratios simultaneously, given a source size. Because the X-rays are thought to originate from a very compact region, we assumed that the X-ray ratios were drawn from the original (unconvolved) magnification pattern. Then, for each image, we constructed a family of two-dimensional histograms:

$$\begin{equation} p_i\big(\Delta \mu _i^X, \Delta \mu _i^O | r_{1/2}, k_*\big), \end{equation} \tag{ 11 }$$

where Δμ^O_i and Δμ^X_i indicate values taken from the convolved (optical) and unconvolved (X-ray) versions of the same pattern. These histograms, after appropriate scaling, are joint probability distributions for the optical and X-ray microlensing magnifications; they are the two-dimensional generalization of the histogram in Equation (7). Figure 4 shows two-dimensional histograms corresponding to the images of PG 1115, with the X-ray magnifications drawn from the original pattern and the optical magnifications drawn from the same pattern smoothed by a source of varying half-light radii. Histograms for all four images are shown for a single source size (0.01r_Ein), and the histogram for the highly magnified saddle point image A2 is shown for a variety of source sizes. The maps used for these histograms all had a stellar mass fraction k_* = 0.1.

Zoom In Zoom Out Reset image size

We shifted the histograms along the optical and X-ray axes in a manner analogous to the one-dimensional case, using the measured flux ratios in X-rays and in the chosen optical/IR filter. For very small source sizes, the two patterns from which each histogram was constructed were nearly identical, so the histograms showed a significant correlation between Δμ^O and Δμ^X. Thus, for a very small source size to be assigned any significant probability, the optical and X-ray shifts had to be the same. As the optical source size increased, the histograms' extent in the Δμ^O direction decreased, as did the correlation between the magnification in the two bands, until for very large source sizes the histograms were nearly aligned with the Δμ^X axis. In this case, the optical shifts had to be nearly zero for the probability to be appreciable (see Figure 4).

In the one-dimensional case, we propagated the uncertainties in the measured flux ratios by convolving our shifted histograms with Gaussians representing the uncertainties. In this two-dimensional context, convolution (e.g., with a two-dimensional Gaussian) gave us insufficient information. Since the source size likelihood in every filter depended on the same X-ray flux ratios, the X-ray uncertainties led to correlated errors in the size estimates at different wavelengths. So we instead performed a Monte Carlo integration in order to propagate the full covariance matrix through the analysis described in the next paragraph. We describe the error propagation in the paragraph that follows.

We multiplied the shifted histograms and marginalized over the unknown X-ray and optical source fluxes to find the likelihood of simultaneously measuring the observed X-ray and optical fluxes:

$$\begin{eqnarray} \mathcal {L}\big(\big\lbrace m_i^X,m_i^O\big\rbrace |r_{1/2},k_*\big) &=& \int dm_S^X dm_S^O \phi \big(m_S^X,m_S^O\big)\nonumber\\ &&\times \left[\prod _{i=1}^{N_{\rm im}} p_i\big({-}m_S^X,-m_S^O\big)\right]. \end{eqnarray} \tag{ 12 }$$

In this case, the weighting function ϕ(m^X_S, m^O_S) incorporated a prior both on the luminosity function of quasars and on the correlation between their X-ray and optical fluxes. The latter is often parameterized using α_OX ≡ 0.3838log(L_2 keV/L_2500 Å). Gibson et al. (2008) find that α_OX = −0.217log(L_2500 Å) plus a constant offset, with a scatter of 0.1. This implies that m^X_S = 0.4346m^O_S plus a constant offset, with a scatter of 0.65 mag. So we constructed a band centered on a line with a slope of 0.4346 in the (m^X_S, m^O_S) plane, and with a Gaussian cross-section. In order to be conservative, we doubled the observed scatter and set the width of the Gaussian to correspond to σ = 1.3 magnitudes. We adjusted the constant offset so that the line passed through the mean of the product histogram. We multiplied this function by a power law of slope −2.95 (as in the one-dimensional case) in the direction along the line of correlation to represent the luminosity function. We note that the slope of the luminosity function of X-ray selected bright quasars is very similar at −2.8 (Silverman et al. 2008). As before, we repeated this technique using magnification patterns representing a range of stellar mass fractions and marginalized over this parameter using a weighted average of the resulting likelihood distributions (see Equation (10)).

The Monte Carlo integration which allowed us to propagate the measurement uncertainties into the likelihood distributions consisted of repeating the process of shifting, multiplying, and integrating under the histograms for 1000 sets of flux ratios drawn from Gaussian distributions implied by the measurement uncertainties. This gave us 1000 distributions for the disk size in every filter. From these we took the most likely values (the modes) and constructed the covariance matrix C_ij = Cov(r_1/2,i, r_1/2,j), where i and j refer to optical/IR filters. The diagonal elements of the matrix are simply the variances of the half-light radii due to measurement uncertainties. To this matrix we added a second, diagonal matrix $\mathcal {S}_{ij}$, which expressed the "intrinsic" width of the first of our 1000 likelihood distributions, the one derived from the reported flux ratios. This intrinsic width arose not from measurement uncertainties, but simply from the widths of the magnification histograms. The diagonal elements of $\mathcal {S}$ were the squares of the halved 68% confidence intervals of these distributions, and the remaining elements were zero. We describe in the next section how we used the inverse of the sum of these covariance matrices

$$\begin{equation} \mathcal {W} = (C + \mathcal {S})^{-1} \end{equation} \tag{ 13 }$$

as a weight matrix in a chi-square fit.

The fact that the optical flux ratios demonstrate microlensing magnifications allows us to place upper limits on the size of the optical disk, since a very large disk will be unaffected by microlensing. Likewise, the fact that the X-ray flux ratios are different from (typically, larger than) the optical ratios allows us to place lower limits. For small optical sizes, the differing flux ratios cause the highly correlated histograms to be shifted in a direction different from that of the correlation, lowering the likelihood (see Figure 4). The posterior probability distribution for the i'-band radius of the accretion disk of PG 1115, determined in this manner, is shown in Figure 5, and shows a central peak, with a decrease in probability to the left due to the interaction of the X-ray and optical flux ratios, and a decrease to the right due primarily to the optical flux ratios alone.

Zoom In Zoom Out Reset image size

We tested the accuracy of our analysis method by applying it to artificial flux ratios generated using the predictions of the thin disk model. We calculated the half-light radii that this model predicts for PG 1115 in each filter, then convolved the magnification patterns appropriate for the four images of this lens (with stellar mass fraction k_* = 0.1) with source profiles of these radii. Then we drew microlensing magnifications from 30 random locations in these patterns. We drew X-ray ratios from the same locations in the corresponding non-smoothed patterns. We ran our analysis on these 30 sets of simulated flux ratios, with uncertainties equal to those assigned to the real PG 1115 data, to check if we recovered the input temperature profile. We fit a line to the recovered sizes to determine the logarithmic slope and normalization of the wavelength-dependent half-light radius, using the procedure described in Section 7. On average, the best-fit overall normalization log(r_1/2/cm) was 16.1, with a scatter of 0.2 dex, compared to the input value of 15.7. The average output slope was 1.2, with a scatter of 0.1, versus the input value of 1.33. The general agreement of the recovered parameter values with the input values indicates that the analysis method produces sensible results. The small offset between the input and output normalizations indicates that our analysis may overestimate the sizes slightly.

The result of the analysis described in Section 6 was a posterior probability distribution for the half-light radius of the quasar accretion disk in each of our 12 lensed quasars, in each filter. The medians and modes of these distributions are listed, along with 1σ confidence intervals, in Table 7, and plotted in Figures 6–9 as a function of rest wavelength. The error bars are taken from the diagonal elements of the covariance matrix $\mathcal {W}^{-1}$ defined in Equation (13). The choice of prior for the half-light radii has a small effect on the measured values. In general, the linear prior favors values of r_1/2 ∼ 0.2 dex larger than the logarithmic prior does. This is to be expected, since the linear prior has more probability density at larger sizes. The slope with wavelength is generally unaffected by the choice of prior. Since a logarithmic prior is more appropriate for a scale-free quantity such as the half-light radius, we adopt these values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The half-light radii are measured in terms of the Einstein radius of a solar-mass microlens r_Ein. We converted these angular sizes to physical radii in the source plane using the angular diameter distances D_OS, D_OL, and D_LS (observer-to-source, observer-to-lens, and lens-to-source, respectively). The physical source radius depends weakly on the redshifts of the source and lens, and the chosen cosmology. It also depends on our choice of the stellar mass function, going as the square root of the mean mass (see Section 6.1). We show physical distances on the left axis of Figures 6–9, with the corresponding fraction of an Einstein radius on the right axis.

In order to compare the measured source sizes to the thin disk predictions, we must estimate the Eddington fraction, accretion efficiency, and black hole masses of our sample of quasars. Following Pooley et al. (2007), we set f_Edd = 0.25 (Kollmeier et al. 2006) and η = 0.15 (Yu & Tremaine 2002). We turned to the literature to find black hole mass estimates for most of our sample. For MG J0414, HE 0435, RX J0911, SDSS J0924, PG 1115, and RX J1131 we adopted the virial mass estimates of Peng et al. (2006), which are based on the widths of the broad emission lines C iv, Mg ii, and Hβ. For SDSS J1138 we adopted that of Morgan et al. (2010), calculated using the same method. We adopted the masses estimated by Pooley et al. (2007) for HE 0230 and WFI J2033 using their optical and X-ray luminosities. For the remaining quasars (HE 1113, SDSS J1330, and WFI J2026) we used the luminosity method of Pooley et al. (2007), inferring a bolometric luminosity from the optical flux. The optical estimates used the relation

$$\begin{equation} L_{\rm bol} = 9[\lambda f_\lambda ]_{5100\,{\rm \AA}} 4 \pi d_l^2 \end{equation} \tag{ 14 }$$

(Kaspi et al. 2000), where d_l is the luminosity distance to the quasar, and the rest-frame 5100 Å flux is extrapolated from the magnification-corrected HST NICMOS F160W flux of the LM image using an assumed power-law spectrum f_λ ∼ λ^−1.7(Kollmeier et al. 2006). We used our assumed Eddington ratio f_Edd = 0.25 to calculate black hole masses. The luminosity and virial mass estimates for our sample of quasars are listed in Table 8. We adopted the virial estimates where available, and the luminosity masses elsewhere. We note that in cases where both are available, the luminosity-based estimates are systematically smaller than the virial estimates by ∼0.5 dex; this could lead to underpredictions of the accretion disk size by factors up to ∼2 for the cases where we use the luminosity estimate in the absence of a virial estimate. This systematic bias is a little surprising, since the bolometric luminosity technique is calibrated using the virial method (Kollmeier et al. 2006). It may be partially due to a bias in our magnification correction, but it is difficult to imagine that we have systematically overestimated magnifications by a factor of three.

Using these black hole mass estimates, we plot the prediction of the thin disk model for each quasar in Figures 6–9. Comparing the microlensing sizes to these curves, it is immediately apparent that the measured sizes are systematically larger than predicted by the thin disk model and vary more weakly with wavelength than predicted.

For each lens we fit a power-law model Y(λ) to the medians of the distributions:

$$\begin{equation} Y(\lambda) = Y_c +\nu \log \left(\frac{\lambda }{\lambda _c}\right), \end{equation} \tag{ 15 }$$

where the central wavelength λ_c is the geometric mean of the rest wavelengths and Y_c is the logarithm of the half-light radius of the source at that wavelength, in cm. Since the errors on the medians are correlated, we minimized the statistic

$$\begin{equation} \chi ^2 = \left(\mathbf {y} - Y(\boldsymbol{\lambda })\right)^{\rm T} \cdot \mathcal {W} \cdot \left(\mathbf {y} - Y(\boldsymbol{\lambda })\right), \end{equation} \tag{ 16 }$$

where y is the vector of median sizes measured at wavelengths $\boldsymbol{\lambda }$, and $\mathcal {W}$ is the weight matrix defined in Equation (13). The parameters Y_c and ν of the best-fit power laws are listed in Table 9.

In nearly every case, the measured disk size at the central wavelength λ_c is larger than the thin disk model predicts, by factors ranging from ∼2 to more than 30. The logarithm of this factor is listed for each quasar in Table 9. Comparing the logarithmic offset to our errors in Y_c, we have ruled out the thin disk normalization by at least 3σ in all cases except MG J0414, and as many as 10σ in some cases. This result is consistent with that of Pooley et al. (2007), though our analysis method is more quantitative. The average offset between Y_c and the thin disk prediction at the same wavelength is 0.89 dex, or a factor of 7.5. A chi-square test comparing the measured sizes to predicted sizes yields χ²/N_dof = 46.

In Figure 10, we plot the half-light disk radii of our sample (excluding SDSS J1330) at fixed rest wavelength versus the black hole mass. The disk radii are calculated using the best-fit model (see Equation (15)) at λ = 2500 Å. This wavelength is chosen to match the one used by Morgan et al. (2010), who extrapolate from single-wavelength microlensing observations to create a similar plot. The errors on the half-light radii are the uncertainties on Y_c (see Table 9), and uncertainties of 0.4 dex have been assigned to the black hole masses (but not plotted). We multiplied the mass uncertainties by the thin disk model slope of 2/3 and added them in quadrature to the uncertainties on the half-light radii. With these uncertainties, the best-fit slope for r_1/2 versus M_BH is 0.27 ± 0.17. The χ² for this fit is 17.3 for 9 degrees of freedom. Thus, we find a slope that is about half the expected r_1/2 ∝ M^2/3_BH dependence (see Equation (2)), but with a substantial uncertainty. Perhaps more importantly, we find a significant offset from the thin disk model in the overall normalization, nearly an order of magnitude. The only quasar that we find to be consistent with the prediction is MG J0414.

Zoom In Zoom Out Reset image size

For comparison, we have plotted the microlensing size estimates of Morgan et al. (2010) in Figure 10, including their wavelength corrections but undoing their inclination corrections to match our convention and multiplying by a factor of 2.44 to convert to half-light radii (see Equation (2)). They find a logarithmic offset between their microlensing sizes and the thin disk prediction in the same sense as ours, and about half as large. Though we do not have a good explanation for the offset between our microlensing size measurements and the prediction, the smaller offset between our sizes and those of Morgan et al. (2010) may arise from the difference in analysis techniques; analyzing this possibility is the subject of future work.

Increasing the Eddington fraction of the quasars, or decreasing their accretion efficiencies, would help to reconcile our microlensing sizes with the predictions of the thin disk model, though the adjustments would need to be drastic, because of the weak dependence of the size on these values. Increasing the black hole mass estimates systematically by factors of ∼20 or decreasing the mean stellar mass in the lensing galaxies by factors of ∼50 would also bring about agreement, but systematic offsets at these levels seem unlikely.

With one exception, every quasar in our sample displays a power-law slope ν flatter than the expected 4/3 by 1.2 to 2.6σ or more (see Table 9). Although these do not individually rule out the thin disk slope at high significance, the 10 to 1 preponderance of quasars with too-shallow slopes makes a convincing case against it. Treating the measurement for each quasar as an independent constraint yields an average value ν = 0.17 ± 0.15. Comparing the eleven measured slopes to 4/3 yields a χ²/N_dof = 6.1. Since a multi-temperature blackbody accretion disk with a power-law temperature profile has T_eff ∝ r^−1/ν, applying our results to this model implies that the effective temperature is a steeply falling function of radius. It is worth noting that nearly all the measured slopes are consistent with ν = 0, i.e., a source with a wavelength-independent size, though the existence of such an object seems highly unlikely, from both theoretical (e.g., energy conservation) and observational (e.g., quasar variability) standpoints.

MG J0414, the only quasar that matches the thin disk size prediction (and the only radio-loud quasar in our sample), is also the only quasar with a slope ν consistent with 4/3. Although MG J0414 was the only case where we used empirical mid-IR flux ratios in place of model ratios (as we suspected the presence of millilensing), it is unlikely that its unique behavior is due to this choice, since mid-IR radios for other lensed quasars (including PG 1115) match the models very well (Chiba et al. 2005; Minezaki et al. 2009).

Two lenses, HE 0230 and WFI J2033, have microlensing size estimates that decrease with wavelength. This result does not seem to arise from a gross failure of our analysis method, since the flux ratios themselves become more anomalous with wavelength, which in most cases implies a smaller source size (see Table 5). It cannot be due to microlensing variability, since the optical/IR measurements are coeval. But it should not be concluded that the temperature profile of these quasars is in fact inverted, as the errors do not rule out positive slopes. It is possible that the sources lie on an unusual region of the caustic pattern where the smoothing caused by a larger source causes more anomalous flux ratios. These lenses underscore the importance of having a large sample of objects when using single-epoch measurements.

We have obtained flux ratios in eight optical and IR filters spanning a factor of six in wavelength, for a sample of 12 strongly lensed quasars. Comparing these ratios to each other and to Chandra X-ray flux ratios, both from the literature and newly reported here, we see chromatic variations, which we attribute to the microlensing of a multicolor accretion disk. The standard thin accretion disk model predicts a temperature profile that falls as the β = 3/4 power of radius, implying that the half-light radius of the disk is proportional to λ^4/3, with an overall normalization that depends on the black hole mass and (weakly) on the Eddington fraction and accretion efficiency. This chromatic dependence combines with the dependence of microlensing magnification on source size to produce higher-amplitude microlensing variations at blue wavelengths than at red wavelengths. For single-epoch observations like ours, this usually means that the observed flux ratios at blue wavelengths are more anomalous than at red wavelengths, compared to the ratios predicted by smooth lens models.

In addition to formal statistical uncertainties, we have estimated the systematic errors in the flux ratios due to spectral contamination by broad emission lines, quasar variability combined with lensing time delays, and microlensing variability combined with delays between measurement epochs. We have described a Bayesian method for determining the probability distribution of the half-light radius of the source in each filter. For quasars without X-ray measurements, we are able to place upper limits on the half-light radii, by virtue of the departures of the flux ratios from the model predictions. When X-ray ratios are available, we assume that they originate from a very compact region, and we are able to place both upper and lower limits on the half-light radius, reasoning that the differences between the X-ray and optical ratios arise because of the appreciable optical radius.

The resulting size estimates are larger overall than is predicted by the standard thin disk model by nearly an order of magnitude, on average. The scaling of the half-light radius with mass is consistent with the expected slope, but the scaling with wavelength is shallower than expected. Although the large error bars on the wavelength slopes mean that a single slope measurement would have little weight, the fact that all but one of the measured slopes are too shallow commands attention. Since the scaling of radius with wavelength is the reciprocal of the scaling of temperature with radius, our results indicate that if we assume a blackbody accretion disk with a power-law temperature profile the temperature slope is in general steeper than the standard thin disk model predicts. However, other explanations for the relative insensitivity of radius to wavelength are also possible, though none immediately presents itself.

There is growing evidence from microlensing studies that quasar accretion disks are larger at optical wavelengths than simple accretion models predict. Our result is a confirmation of earlier, more qualitative work by Pooley et al. (2007). It also corroborates the results of Morgan et al. (2010), who use several lensed quasars to study the mass dependence of the accretion disk radius, and find it difficult to reconcile their radii with thin disk theory. Our results for the observer-frame r'-band size of MG J0414 and u'-band size of SDSS J0924 are consistent with the upper limits of Bate et al. (2008) and Floyd et al. (2009), respectively. Likewise, our estimates of the radius of RX J1131 at around 4000 Å in the rest frame, and of HE 0435 at around 6500 Å in the observed frame, are consistent with those of Dai et al. (2010) and Mosquera et al. (2010), respectively.

The picture is less clear regarding the temperature profile slope β. Using time-domain measurements at several wavelengths, Poindexter et al. (2008) find β ∼ 0.6 for the quasar HE 1104−1805, consistent with the thin disk slope of 3/4. They point out that a flatter temperature slope (i.e., β < 3/4) would help to explain the faintness of lensed quasars relative to the expected flux from a Shakura–Sunyaev disk of the size implied by microlensing. However, our result is just the opposite: the low average value of ν implies that the temperature slope β is greater than 3/4. Other single-lens microlensing studies have addressed this question as well. Bate et al. (2008) and Floyd et al. (2009) use single-epoch broadband photometry similar to ours to constrain the slope ν for the quasars MG J0414 and SDSS J0924, respectively; their limits are compatible both with our results and with the thin disk model. Mosquera et al. (2010) use narrowband imaging of HE 0435, finding a slope ν = 1.3 ± 0.3, consistent with our slope within the errors. In addition, several studies have focused on Q2237+0305, because of the short time scale and large amplitude of its microlensing variability, and the large amount of monitoring data that are available. Anguita et al. (2008) find a ratio of r'- to g'-band source sizes consistent with the thin disk model. Eigenbrod et al. (2008) use spectrophotometric monitoring to constrain ν = 1.2 ± 0.3, and Mosquera et al. (2009) use narrowband photometry, finding solutions compatible with ν = 4/3, as well as with higher values. Our result, working with a larger sample size, is the first to indicate that ν is significantly less than 4/3, so further study is clearly needed. We note that the error bars on ν are uniformly large for individual lenses. Increasing the sample size, as we have done, is a good way to improve this; additional strategies include using spectroscopic methods to transform line emission from a nuisance into a tool, improving the coordination between X-ray and optical observing campaigns, and moving to multiwavelength, long temporal baseline monitoring.

The authors thank Joachim Wambsganss for the use of his microlensing simulation code. J.A.B. and P.L.S. acknowledge the support of US NSF grant AST-0607601.