SPITZER OBSERVATIONS OF A GRAVITATIONALLY LENSED QUASAR, QSO 2237+0305

QSO 2237+0305 was chosen for this study as it holds several records among gravitationally lensed quasars: it was one of the first four-image lenses discovered (Huchra et al. 1985); its lens galaxy has the lowest redshift, z_l = 0.0395 (Huchra et al. 1985); and it was the first to show gravitational microlensing (Irwin et al. 1989). This last fact is a result of the second: a nearby lens galaxy causes a large velocity of the quasar relative to the magnification patterns created by stars in the lens galaxy projected onto the source plane; this large relative velocity results in a shorter timescale for microlensing. The discovery of microlensing in this system and its short microlensing timescale made it a "rosetta stone" for studies of the size of the quasar emission region: the time-dependent microlensing magnification is sensitive to the size of the source, effectively resolving the quasar on submicroarcsecond scales. Larger sources smooth over the microlensing magnification pattern and thus experience smaller and more gradual variations in magnification (Refsdal & Stabell 1991). In unlensed quasars, only the SED can be compared to models (e.g., Sanders et al. 1989; Blaes et al. 2001; Malkan 1983), while for QSO 2237+0305 the size as a function of wavelength can be compared to models as well, in principle, giving much stronger constraints on the emission mechanism.

Despite this promise, the interpretation of the first optical microlensing events in QSO 2237+0305 was puzzling: one study showed the inferred size of the big blue bump was consistent with the accretion disk model (Wambsganss et al. 1990), while another study showed the size of the emission region was too small (Rauch & Blandford 1991). The latter result led to other newer models which require more theoretical development, e.g., Barvainis (1993) and Czerny et al. (1994). With a much larger data set and more sophisticated analysis of the microlensing light curves, Kochanek (2004) showed that thermal emission from an accretion disk is consistent with the size inferred from microlensing. However, microlensing in a sample of gravitationally lensed quasars has led to a different conclusion: the size of the optical/UV emission region inferred from microlensing is too large compared with the size of quasar accretion disk models inferred from fitting the SEDs (Pooley et al. 2007; Morgan et al. 2009).

This confused state of affairs of microlensing of the big blue bump partly stems from the fact that the absolute size is difficult to constrain as there are degeneracies between the mass of the microlenses and the sky velocity of the quasar relative to the magnification pattern. However, the relative size versus wavelength is much easier to constrain since it is not as subject to these degeneracies (Wambsganss & Paczynski 1991). In particular, the first results for the wavelength-dependent relative size seem to be in good agreement with the accretion disk model for a different lensed quasar (Poindexter et al. 2008), although the absolute size is still discrepant. The relative size of the optical/UV/X-ray emission region for QSO 2237+0305 is well constrained by microlensing (e.g., Wyithe et al. 2000; Kochanek 2004; Anguita et al. 2008), taking advantage of the long timescale data set collected by the Optical Gravitational Lensing Experiment (OGLE; Udalski et al. 2006; Woźniak et al. 2000). An intensive monitoring campaign with the Very Large Telescope (VLT) promises to give a very detailed picture of the relative sizes of the big blue bump and broad-line regions as a function of wavelength (Eigenbrod et al. 2008a, 2008b). In this paper, we will not attempt to resolve the absolute size problem, but rather we will argue that the standard two-component model provides good agreement with the wavelength dependence of the microlensing flux ratios, adding credence to the standard model.

The first microlensing study of the IR bump was carried out with QSO 2237+0305 to distinguish synchrotron and dust emission models for the IR bump (Agol et al. 2000). The synchrotron emission model provides an alternative, albeit less natural, explanation for the IR bump. The synchrotron emission region responsible for the IR bump has to be compact to avoid self-absorption; thus it should show strong variability due to microlensing. On the other hand, the dusty torus model must be extended to avoid sublimation, and thus should vary weakly due to microlensing. Agol et al. (2000) found that the mid-IR flux ratios were consistent with no microlensing (Schmidt et al. 1998); this despite the fact that the optical source was simultaneously undergoing strong microlensing events. These observations ruled out strong microlensing magnification of the mid-IR emission region, which was one of the first clearcut demonstrations that the IR emission region in radio-quiet quasars is due to thermal emission by dust, not synchrotron emission (Wyithe et al. 2002). Here, we extend these results to observations near the 1 μm dip where both the dusty torus and accretion disk contribute to the flux.

In Section 2, we discuss the observations and data reductions. Although the primary focus of this paper is on probing the relative source size versus wavelength, there are two problems related to the SED that may be addressed with our data as well. (1) How similar is the SED of QSO 2237+0305 to other quasars and Seyfert galaxies? Ground-based observations indicated that QSO 2237+0305 contained hotter dust than any other quasar (Agol et al. 2000), while the observations presented here show that the ground-based observations at one wavelength were in error. The large intrinsic luminosity and high magnification, μ ∼ 16, (Schmidt et al. 1998) make this QSO an excellent candidate for spectroscopy and allow us to compare the spectrum of a high-redshift quasar with nearby Seyfert galaxies. In Section 3.1, we show that QSO 2237+0305 looks very similar to other quasars and Seyfert galaxies. (2) How well does the IR SED match dusty torus models? In Section 3.2, we show that the overall shape agrees qualitatively, but the quantitative agreement is poor.

In Sections 3.3–3.4, we present the microlensing results and interpretation for QSO 2237+0305, demonstrating that two size scales are required to fit the microlensing flux ratios, as predicted by the accretion disk/dusty torus model. In Section 4, we discuss the implications for quasars in general and in Section 5 we summarize.

Basic Calibrated Data (BCD) images were obtained from the Spitzer archive, pipeline version S14.0.0. QSO 2237+0305 was observed with the IRS modules Short-Low (SL2, 5.2–8.7 μm, and SL1, 7.4–14.5 μm), Long-Low (LL2, 14.0–21.3 μm), and Long-High (LH, 18.7–37.2 μm), for a full observed wavelength coverage of 5.2–37.2 μm. Rogue pixels were eliminated using the IRSCLEAN_MASK software package provided by the Spitzer Science Center. We created our own rogue pixel maps (pixels with anomalous behavior) for each spectral order by measuring two quantities from a series of sky images for each order (for the LH data we chose only portions of the image that did not contain the target): (1) the scatter in each pixel with time; (2) the difference between the value of a pixel and the median of a region within a 5×5 pixel region surrounding it. We then flagged pixels which had either excessive scatter or consistently had values much larger than the median-smoothed image, and included these in the rogue pixel map. This procedure resulted in similar maps to those generated automatically by the IRSCLEAN_MASK software, but was better at flagging more rogue pixels so that we did not have to flag any pixels by hand. Rejected pixels were interpolated from the surrounding pixels.

The cleaned spectra were co-added and extracted using the Spectroscopy Modeling Analysis and Reduction Tool (SMART; Higdon et al. 2004). The method of sky subtraction depended on the resolution of the data. For the low-resolution data, sky subtraction was performed by subtracting one nod position from the other before extraction. For high-resolution data, the narrow width of the slit required that separate sky images be subtracted from the QSO images in each nod position. A separate set of images 130 arcseconds away from the QSO was taken for this purpose. This same procedure applied to a standard star (HR 7341) yielded a spectrum which matched between each of the IRS orders/modules, matched the calibration spectrum within 5%, and gave a spectrum which obeyed the Rayleigh–Jeans limit, so we are confident of the relative calibration of our data, but expect that the absolute calibration has an uncertainty of 5%. We rereduced the data with later versions of the pipeline which resulted in fluxes that differ by as much as 20% in the overlapping region between different orders, while the 14.0.0 pipeline did not have this problem.

For each order of each module, we fitted a Gaussian to the distribution of the difference in flux between the two nods divided by the sum of the squares of the uncertainties; in all cases the standard deviation of this distribution differed from unity indicating that the uncertainties were misestimated. We scaled the SMART uncertainties by the standard deviation of this Gaussian. We then fitted the median-smoothed spectrum from all nods and orders with a fifth-order polynomial, and cleaned the data of points which disagreed by >3σ from this fit, as well as points for which the two nods disagreed by >3σ. This procedure automatically removed data near the edges of each order, which are notoriously unreliable, and also removed other outliers which may be due to improperly cleaned cosmic rays or rogue pixels.

QSO 2237+0305 was observed for 38.4 s in each of the four IRAC wave bands (3.6, 4.5, 5.8, and 8.0 μm) in full array mode with a 16 position spiral dither pattern with 2 s exposures at each position. Post-BCD mosaics were obtained from the Spitzer archive (pipeline ver. 13.0.2), which we used in our analysis. As QSO 2237+0305 is only a few pixels across at the resolution of IRAC (see Figure 1), the four lensed images and lens galaxy are unresolved; however, we were still able to derive the flux ratios of the four images.

Zoom In Zoom Out Reset image size

For comparison with microlensing models, we derived the fluxes of the lensed quasar images from the IRAC data. As the pixel size of the IRAC images, 12, is comparable to the separation of the quasar images, this required a multicomponent model fit. Since the relative positions of the quasar images are known extremely accurately, and since the IRAC point-spread function (PSF) is known fairly precisely, we derived the flux ratios of the four lensed images of the quasar with PSF fitting. The main uncertainty in the fitting is the contribution of the lens galaxy to the flux in the IRAC bands; we addressed this by using the Hubble Space Telescope (HST) H band as a model from the CASTLES survey (Kochanek et al. 2009), assuming no color gradients between the H band and IRAC bands. This is likely a good approximation as extinction and intrinsic colors should vary weakly in the IR since the stellar emission is well into the Rayleigh–Jeans tail.

We created a model composed of (1) four quasar images; (2) the one star 10'' from the center of the lens galaxy; (3) the lens galaxy flux, scaled from a deconvolved HST H-band image; and (4) a uniform sky background. We held the relative positions of the quasar images (and the nearby star) fixed to the values derived from the HST H-band image, while we allowed the absolute position to vary (given the uncertain absolute pointing of Spitzer).

The lens galaxy was isolated in the HST H-band image by masking the quasar images and stars (within a circle 14.25 pixels from the location of each point source), and the masked region for each quasar image was replaced with an elliptical Sersic model fit to the remaining H-band data, while the masked region near the star was filled in with the median flux near its location. As the Spitzer IRAC PSF is derived at five times the pixel resolution (024), we rotated and compressed the HST image to fit the Spitzer images at five times the resolution. We then convolved the HST image with the IRAC PSF for each band, multiplying by a constant factor to scale to each IRAC band, and added to this the five point fluxes by interpolating the Spitzer IRAC PSF to the location of each point source and multiplying by their respective fluxes. Finally, we added in a constant flux to represent the sky.

These model components give a total of nine free parameters to fit: five point sources, the extended galaxy flux scaling factor, the sky flux, and the R.A. and decl. of image A, which was taken as the reference point. We computed the χ² of this model by comparing with the post-BCD mosaic and uncertainties from the Spitzer IRAC pipeline. We optimized the model parameters using the Levenberg–Marquardt method, and then found the uncertainties on each parameter from a Markov chain computed using the method described in Tegmark et al. (2004). The best-fit χ² for the four IRAC channels was (314, 209, 91, 162) for 236 degrees of freedom (dof; nine model parameters to fit the flux of 245 pixels which is the region covered by the HST H-band image). Formally, these fits range from very good to poor, which may indicate that the model is inadequate (e.g., possibly the galaxy has color gradients between 2.2 and 3.6 μm), or the error bars are misestimated. We also computed error bars on the model parameters using the covariance matrix evaluated at the best fit and by finding the region with Δχ² < 1 for each parameter while marginalizing over the other parameters; each of these techniques gave error bars nearly identical to the Markov chain. We converted these values to fluxes in mJy, as well as flux ratios, and report the derived fluxes and errors in Table 2, as well as the flux ratios in Table 3. For the galaxy, we report the entire model flux within 5 pixels (6'') of the center of the galaxy except for the contribution from quasar images. The V-band data is from data taken by the OGLE collaboration 1 day before the Spitzer observations, and the errors reported are relative flux errors, not absolute (Udalski et al. 2006). In addition, we report the continuum spectral slope, α_ν measured at 5400 Å for $f_\nu \propto \nu ^{\alpha _\nu }$, for all four images measured with a VLT observation on 2005 November 11 (Eigenbrod et al. 2008b).

In this section, we compare the Spitzer SED of QSO 2237+0305 to Seyfert galaxies and quasars to show that it looks like a typical radio-quiet active galaxy in the IR. Figure 2 shows the Spitzer spectrum of QSO 2237+0305 which has been binned so that each bin has a signal to noise greater than 100. The excellent agreement between the IRS and IRAC results, which had completely independent flux calibration, bolsters our confidence in the accuracy of our reported fluxes.

Zoom In Zoom Out Reset image size

We fitted a power law of the form F_ν ∝ ν^α to the IRS spectrum, and we find α = −0.96 ±.02, giving a spectrum which is nearly flat in νF_ν.

Figure 3 shows the full spectrum of QSO 2237+0305; the isotropic luminosity is defined as (νL_ν)_rest = 4πD²_L(νF_ν)_obs/μ, where rest/obs refer to the rest-frame/observed frequencies, D_L is the luminosity distance of the quasar, and μ is the total magnification of the quasar. Included are our data from both IRS and IRAC, optical and near-IR data points from the CfA-Arizona Space Telescope LEns Survey (Falco et al. 2001), OGLE (Udalski et al. 2006), and Eigenbrod et al. (2008b), X-ray data from Dai et al. (2003), two data points from the Galaxy Evolution Explorer (GALEX) archive (Martin et al. 2005), and a submillimeter data point at 850 μm from Barvainis & Ivison (2002). To compute the total luminosity, we have assumed the cosmological parameters from the WMAP 5 year data set (Dunkley et al. 2009) as well as a total macrolensing magnification of μ = 16 (Schmidt et al. 1998). The optical data have been corrected for extinction in the Milky Way assuming a Galactic reddening of E(B − V) = 0.068 with R_V = 3.1 extinction curve. As the light from the quasar passes through different portions of the bulge of the lens galaxy, we need to make additional extinction corrections for the four lensed images. We have used the flux ratios of the broad lines averaged over time from the data set of Eigenbrod et al. (2008b) to derive the relative extinction of the four quasar lensed images. We find image B is extincted relative to image A by ΔE(B − V) = 0.02 ± 0.05, while images C and D are reddened with respect to images A and B by ΔE(B − V) = 0.10 ± 0.04,  0.18 ± 0.03, respectively. Since images A and B have small (possibly zero) relative extinction, we assume that each of these images has zero extinction in the lens galaxy, and simply correct images C and D. We have not attempted to correct the data for microlensing, so the overall uncertainty is at least 0.2 mag. From the QSO 2237+0305 SED, we find a total luminosity of L_tot = 4.0 × 10⁴⁶ erg s⁻¹.

Zoom In Zoom Out Reset image size

In Figure 3, we plot the SED of QSO 2237+0305, and for reference compare it to the composite radio-quiet quasar SED from Elvis et al. (1994), normalized to 1.3 μm. Although QSO 2237+0305 appears underluminous in the X-ray, UV, mid-IR, and submillimeter relative to the composite, this behavior is well within the range of SEDs in the Elvis sample, and it is likely that the composite is affected by selection biases at these wavelengths where many quasars were not detected. The SED of QSO 2237+0305, shown in Figure 4, looks fairly typical compared to a composite spectrum of Palomar–Green (PG) quasars with weak far-IR emission (Netzer 2007). There are minor differences such as an extra bump near 6–7 μm and a peak associated with the hottest dust at slightly longer wavelengths (∼2.8 ± 0.3 μm), but these differences are well within the range of variation within the PG quasar sample. If we "fitted" the QSO 2237+0305 SED by scaling the Netzer composite spectrum by an arbitrary factor, we find a χ² = 573 for 32 dof, which is formally a very poor fit, but the discrepancy is dominated by the disagreement in the cutoff at short wavelengths and the bump near 7 μm.

Zoom In Zoom Out Reset image size

The IR spectrum of QSO 2237+0305 also looks very similar to low-redshift Seyfert galaxies taken from a sample of 23 galaxies (V. Gorjian et al. 2009, in preparation). The Spitzer spectrum of QSO 2237+0305 is plotted with the two most similar Seyfert spectra in Figure 5, scaled to match the flux of QSO 2237+0305. The qualitative shape of the SEDs matches well from 4 to 10 μm, although Mrk 509 shows stronger emission features, presumably due to silicates.

Zoom In Zoom Out Reset image size

The similarity of the IR SED of QSO 2237+0305 to other Seyferts and quasars indicates that our microlensing studies of this object will broadly apply to radio-quiet active galaxies.

To improve our physical understanding of the emission from QSO 2237+0305,  we fitted the SED with the models of Fritz et al. (2006). The models utilize the Mathis–Rumpl–Nordsieck (MRN) dust size distribution (Mathis et al. 1977) with scattering and absorption opacities from Laor & Draine (1993). The model fixes the geometry as a torus with an opening angle that is independent of radius and an inner radius of the torus,

$$\begin{equation} R_{\rm min} = 1.3\,{\rm pc}\,L_{46}^{1/2} T_{1500}^{-2.8}, \end{equation} \tag{ 1 }$$

which is set by a dust sublimation temperature of 1500 K, where L₄₆ is the active galactic nucleus (AGN) luminosity in units of 10⁴⁶ erg s⁻¹. The dust density in the model is described by ρ(r, θ) ∝ r^βe^{−γ|cos(θ)|}, with a dust-free cone within polar angle θ < θ_cone. The grid of models covers a range of parameters for the dust with: (1) the ratio of the outer to inner radius, 30 < R_max/R_min < 300; (2) the variation of dust density with radius, −1 < β < 1/2; (3) the equatorial optical depth at 9.7 μm, 0.1 < τ_9.7 < 10; (4) inclination angles, i, ranging from 0.01 (edge-on) to 89° (face-on), and 11°,  21°,  31°,  41°,  51°,  61°,  71°,  81° in between; (5) dust-free cone with size 20° < θ_cone < 60°; and (6) an angular cutoff in density with 0 < γ < 6.

We have scaled the IR portion of the models by an arbitrary constant to optimize the match with our Spitzer SED. The best-fit model is shown in Figure 6, which has parameters R_max/R_min = 100, β = −1/2, τ_9.7 = 10, i = 71°, θ_cone = 20°, and γ = 3. The best fit implies an AGN luminosity of 4.4 × 10⁴⁶ erg s⁻¹, remarkably close to the luminosity found from the QSO 2237+0305 SED (L = 4.0 × 10⁴⁶ erg s⁻¹), despite the fact that we have only fitted the Fritz models to the IR data. This may imply that the dust acts as a fairly good calorimeter of the total AGN flux. The qualitative fit to the data is fair: the model shows a peak near 2.5 μm in the rest frame, and gradual decline toward longer wavelengths, and an emission feature near 10 μm. However, quantitatively the fit is horrible: χ² = 3045 for 94 dof. This is primarily due to the fact that the short-wavelength peak is more prominent in the model than in the observations and the silicate absorption and emission features in the model are offset in wavelength of the observed features. It is possible that optimizing the parameters will improve the fit as the grid is quite coarse and some of the best-fit parameters are at the extreme values of the grid, such as θ_cone. Also, the viewing angle is 71° which is only 1° within the opening angle of the cone; however, viewing angles of 81° and 89° are very similar in shape, but only slightly poorer fits, plotted as dotted lines in Figure 6. Consequently, we do not believe that the fitted parameters are unique or even correct; indeed the simple geometry chosen by Fritz et al. may be wrong. The main point is that a dusty torus model can produce a fair qualitative fit to the SED of QSO 2237+0305; further development of theoretical models will be required to obtain a better quantitative fit.

Zoom In Zoom Out Reset image size

Figure 7 shows the ratios of the quasar image fluxes as a function of wavelength from Table 2. The V-band data are the data obtained from the OGLE data archive (Udalski et al. 2006) taken at a time closest to our IRAC observations: 2005 November 16 00:48 UT (HJD). The 10 μm points are from Agol et al. (2000), and thus are not simultaneous to our Spitzer observations; however, there appears to be little variability at this wavelength. Also plotted are the model flux ratios from Trott & Webster (2002) which is the most complete model of the lens galaxy of QSO 2237+0305 constructed to date, including the bar and spiral arms. The agreement with the 10 μm flux ratios is expected since these were used as a constraint on the model; however, the model gives very similar flux ratios as an earlier model by Schmidt et al. (1998) which was constructed before the 10 μm observations.

Zoom In Zoom Out Reset image size

The general trend is obvious in Figure 7: for all three pairs of images there is a strong microlensing anomaly in the optical which gradually disappears toward longer wavelengths. This is precisely the behavior expected from the standard model of quasars: the source should be larger at longer wavelengths and thus less affected by microlensing.

Ideally, we would like to use the microlensing flux ratios to constrain the size of the quasar as a function of wavelength; however, only weak constraints can be derived without a detection of time variability due to microlensing (Wyithe et al. 2002). Instead, we use a semiempirical model for the IR SED to predict what the flux ratios should be as a function of wavelength and compare these predictions with the observed flux ratios to confirm the plausibility of this model.

We constructed a semiempirical model for the flux ratios to compare to the observed IRAC flux ratios as follows. We assume that the spectrum consists of a power-law component due to an accretion disk (we are modeling the region from 0.4 to 4.0 μm which is well longward of the peak of the disk spectrum and is only a decade in frequency, so a power law should be an adequate approximation of a disk spectrum) and a single-temperature thermal dust emission component, representing the inner edge of the dusty torus. We fitted the SED from 0.4 to 4.0 μm in the rest frame with these two components, determining their relative strength at each wavelength. The best fit is shown in Figure 8; the model provides a good fit to the four IRAC data points. We have not attempted to correct for microlensing, nor possible time variability as the SED data are not simultaneous. However, this will likely have a small effect on the SED as summing over all four images reduces the impact of microlensing and in the IR quasars are weakly variable.

Zoom In Zoom Out Reset image size

With these two fits, we determined the minimum possible source sizes to reproduce the observed flux with thermal emission as follows. For the power-law component, we assumed a disk geometry with a temperature that is a power law in radius, r, finding T ∝ r^−0.66, and found that the half-light radius should scale with wavelength as

$$\begin{equation} r_{1/2} = 4 \times 10^{16}\,{\rm cm}\,\lambda ^{1.5}, \end{equation} \tag{ 2 }$$

where λ is measured in microns in the rest frame of the quasar. At this radius, the standard disk model is well outside the inner edge and thus is expected to have a temperature dependence of T ∝ r^−3/4, which is close to the measured dependence. We assumed that the dust component either has an emissivity described by optically thin interstellar medium (ISM) dust with the model of Draine (2003) or emits as an optically thick blackbody (BB). These two extremes were chosen to bracket the range of possible behaviors for the hottest dust at the inner edge of the torus (we did not use the best-fit Fritz model due to the different peak wavelength). We found best-fit temperatures of the dust of 1168 K (ISM) or 1264 K (BB). This requires a minimum distance from the quasar of 3.83 pc (ISM) or 0.76 pc (BB) to maintain temperature equilibrium, assuming that the quasar emits isotropically. The total luminosity in this component is 6 × 10⁴⁵ erg s⁻¹ (ISM/BB) implying an emission region of at least 1.3 pc (ISM) or 2.2 pc (BB) in size, consistent with the radiation equilibrium argument within a factor of 3.

If the average mass of microlenses is 0.1 M_☉, then the Einstein radius is r_E = 5.77 × 10¹⁶ cm projected to the source plane. At 1 μm, the power-law emission (if thermal) is comparable in size to an Einstein radius, while the dust emission component is about 200 times larger than the Einstein radius. Thus, the power-law component will be affected by microlensing, while the dust component should be nearly unaffected by microlensing, with variations of less than 1% (Refsdal & Stabell 1991).

To model the microlensing of the power-law component, we extrapolate the wavelength dependence of the optical flux ratios measured with the VLT (Eigenbrod et al. 2008b) to the IR. We correct the optical flux ratios for differential extinction using the E(B − V) values derived from broad emission lines, as described in Section 3.1. Then, the expected flux ratio as a function of wavelength, r_BA(λ), of images A and B is given by

$$\begin{equation} r_{\rm BA}= {\mu ^{\rm B}_{P} F_{P} + \mu ^{\rm B}_{D} F_{D} \over \mu ^{\rm A}_{P} F_{P} + \mu ^{\rm A}_{D} F_{D}} e^{-\tau _{\rm B}+\tau _{\rm A}}, \end{equation} \tag{ 3 }$$

where μ^A,B_P,D(λ) are the magnifications of the power-law (P) and dust (D) components as a function of wavelength, F_P,D(λ) are the intrinsic (unmagnified) fluxes of these two components, and τ_A,B(λ) are the optical depths through the lens galaxy for each component (we have dropped λ in this equation for simplicity).

Now, as argued above, the dust component should be large enough to be unaffected by microlensing, so r_BA,D = μ^B_D/μ^A_D can be derived from a model for microlensing. We utilize the model of Schmidt et al. (1998), improved upon by Trott & Webster (2002), with more recent modifications based on kinematic data (Trott et al. 2009). The relative strengths of the power-law and dust components we take from the model of the SED (Figure 8), f_DP = F_D(λ)/F_P(λ); as mentioned above; as the SED is constructed from the sum of all four images it should be affected only weakly by microlensing. Finally, the wavelength dependence of the microlensing magnification we take from the extinction-corrected optical flux ratios measured with the VLT, r_BA,P = μ^B_P/μ^A_P. For the extinction correction, we assume a Milky Way extinction curve with R = 3.1. Dividing through the numerator and denominator by μ^A_PF_P, we can rewrite the above equation as

$$\begin{equation} r_{\rm BA} = { r_{\rm BA,{\it P}} + \big(\mu ^{\rm A}_{D}\big/\mu ^{\rm A}_{P}\big) r_{\rm BA,{\it P}} f_{DP} \over 1 + \big(\mu ^{\rm A}_{D}\big/\mu ^{\rm A}_{P}\big) f_{DP}} e^{-\tau _{\rm B} + \tau _{\rm A}}. \end{equation} \tag{ 4 }$$

The remaining unknown in this equation is μ^A_D/μ^A_P(λ) since the microlensing magnification of the power-law component for each image is unconstrained by our data. Fortunately, the left-hand side of this equation is weakly dependent on this ratio. We use the size distribution derived from the power-law emission spectrum (Equation (2)) to compute the probability distribution as a function of wavelength from microlensing simulations of each image using the macrolensing parameters from the model of Schmidt et al. (1998). We utilized the code of Wambsganss (1999) to run the simulations, creating 10 simulations for each image of a size 20r_E × 20r_E, and then convolving the magnification pattern with a Gaussian source with the same half-light radius as derived in Equation (2).

To predict the flux ratios at IRAC wavelengths, we have run 10⁴ Monte Carlo simulations sampling the relative extinction, the optical flux ratios, the galaxy macro-lensing model flux ratios, and the relative microlensing magnification for the two images within the uncertainties of each parameter. We have then computed the median and 68% confidence limits at each wavelength from these simulations, which is plotted in Figure 9, for the ratio of images B to A, as well as several other image pairs. We find good agreement between the model predictions and the observed flux ratios. For all but three flux ratios, the data are within the 68% confidence limits of the model flux ratios. Thus, the wavelength dependence of the flux ratios is consistent with the accretion disk/dusty torus model.

Zoom In Zoom Out Reset image size

Agol et al. (2000) observed of QSO 2237+0305 with the Long-Wavelength Spectrometer on the Keck I telescope at 8.9 and 11.7 μm; the flux ratios at these two wavelengths were identical within the errors, and both agreed with the macrolensing flux ratios as predicted by the best lens models. The flux at 11.7 μm reported in that paper agrees well with our IRS spectrum; however, the 8.9 μm Keck data point is higher than the IRS data by about 40% indicating that the Keck data had an incorrect calibration. We have not tracked down the source of this discrepancy, but it cannot be due to microlensing which would have caused a difference in the flux ratios at 8.9 and 11.7 μm. We have more confidence in the calibration of the Spitzer spectrum; consequently the Keck flux at 8.9 μm was likely in error. The qualitative agreement between the Spitzer SED and the Netzer composite indicates that QSO 2237+0305 behaves as a typical quasar in the near-IR.

As mentioned above, the QSO SED, Figure 6, is qualitatively well fitted with an AGN torus model from Fritz et al. (2006); however, the publicly available models have a fixed temperature for the inner edge of the torus at 1500 K which is somewhat higher than the temperature of the inner edge we have estimated (1164–1250 K). In addition, the silicate spectral features of the model are a poor fit to the SED, as is commonly seen in comparing dusty torus models to AGN SEDs (e.g., Nenkova et al. 2002). We have found that some Type II models (obscured quasar) fit the silicate feature well, having an offset silicate feature due to radiation transfer effects, while these models do not fit shorter wavelengths which are obscured. Thus, it may be possible that the dust composition and/or torus opening angle changes with radius causing the silicate feature to appear more like that of Type II quasars, while allowing the inner edge and quasar to be visible so that shorter wavelengths look more like a Type I quasar. Another possibility is that the silicate properties are modified near quasars, either due to changes in the grain size distribution or due to grain porosity, causing the silicate spectral features to be shifted (e.g., Voshchinnikov & Henning 2008). Both possibilities need to be explored in future models of dust grain opacity, as well as computing better physical models for the dusty torus, such as Krolik (2007).

The agreement between the measured and predicted flux ratios is quite good, close to 1σ for all data points except one (Figure 9). The uncertainties in the flux ratio predictions are highly correlated between all wavelengths since microlensing and extinction have a monotonic variation with wavelength, so the case of the ratio of images C to B (for example) is still highly probable. It is clear from Figure 9 that an extrapolation of the optical power law (which can be seen shortward of 1 μm) does not do a good job of predicting the IR flux ratios, while including the unmicrolensed IR bump due to dust emission brings the flux ratios back into agreement with the data (within one standard deviation). If the power-law component had a cutoff around 1 μm so that the IR data were solely due to the extended dust emission, then the flux ratios would change abruptly to the macrolensed values. This is not the case for, e.g., the ratio of image D to image A, indicating that both the power-law and dust emission components are required to fit the IRAC data. This adds evidence to the case for the presence of an accretion disk emission component under the IR bump.

By extrapolating the wavelength dependence of the microlensing flux ratios from the optical, we mostly avoided needing the size of the accretion disk in units of the Einstein radius. The only place the size of the accretion disk enters our analysis is in computing the flux ratios of the disk and torus components for each image in Equation (4) (μ^A_D/μ^A_P, and the same ratio for images B–D). However, this factor cancels out when either the disk or torus components dominate the flux, so our results are weakly sensitive to our choice of source size (Equation (2)) and Einstein radius (M = 0.1 M_☉).