Quasar Microlensing Variability Studies Favor Shallow Accretion Disk Temperature Profiles

Quasars play a critical role in galaxy evolution and growth (Boyle & Terlevich 1998; Di Matteo et al. 2005). Because of the massive size of their central black holes, they may also serve as tests of general relativity (Fabian et al. 2000; Chartas et al. 2017; Dai et al. 2019). Despite the astrophysical importance of quasars and significant effort in the field, we still have an incomplete picture of the physical processes governing quasar emission (e.g., Blaes 2004; Antonucci 2013). No model of quasar accretion disks can fully explain the wide range of observed phenomena, including spectral slope scaling with mass (Davis et al. 2007), the ∼1000 Å spectral cutoff (Shull et al. 2012), and the large optical/UV continuum emission region size (Morgan et al. 2010).

Despite some limitations, the thin-disk model of Shakura & Sunyaev (1973, SS73) and general relativistic extensions (e.g., Novikov & Thorne 1973) have provided a coherent framework within which to interpret observed quasar behavior. In the thin disk model, a hot disk accretes onto the central black hole, with energy transported outward through viscous dissipation. The SS73 thin disk model remains in wide use in the present day mostly due to its relatively simple analytic form.

Subsequent studies have developed progressively more realistic quasar accretion models. Abramowicz et al. (1988) proposed a "slim disk" model in which vertical energy transfer within the disk plays an important role. This changes the disk dynamics and leads to advective flows radially within the disk when accretion rates rise above roughly $\dot{m}\equiv {\dot{M}}_{\mathrm{BH}}/(16{\dot{M}}_{\mathrm{Edd}})=0.1$. Although this model and more recent extensions require numerical solutions, they may better describe the soft X-ray spectra in high-luminosity active galactic nuclei (AGN) and can explain long optical-UV lags (Szuszkiewicz et al. 1996; Sądowski et al. 2011). A key advance came when Balbus & Hawley (1991) demonstrated that a small magnetic field can create turbulence in a disk. These magnetorotational instabilities (MRI) finally gave a physically plausible source of the accretion disk viscosity required in the thin and slim disk models. Since this discovery, magnetic field contributions are included in many accretion disk models (e.g., Agol & Krolik 2000; Penna et al. 2010; Sądowski 2016; Y.-F. Jiang et al. 2019, private communication). State-of-the-art models typically use numerical general relativistic magnetohydrodynamic (GRMHD) simulations to predict accretion disk properties (see Porth et al. 2019). Most of these simulations are, however, presently limited to low accretion rates rather than conditions found in typical quasars.

One of the well-known shortcomings of the thin disk model is its inability to correctly predict the size of accretion disks from their observed luminosities. Assuming blackbody emission within each annulus of the disk, we can estimate the size of an optically thick accretion disk using the observed specific flux, reported as the half light radius r_L,

$$\begin{eqnarray}\begin{array}{rcl}{r}_{L} & = & C(\beta ){r}_{\lambda }=C(\beta ){\left(\displaystyle \frac{{\lambda }_{\mathrm{obs}}^{5}{D}_{L}^{2}{F}_{\lambda ,\mathrm{obs}}}{4\pi {{hc}}^{2}\cos (i){\left(1+z\right)}^{4}}\right)}^{1/2}\\ & & \times \ {\left({\int }_{{u}_{\mathrm{in}}}^{\infty }\displaystyle \frac{u{du}}{\exp \left({u}^{\beta }\right)-1}.\right)}^{-1/2}.\end{array}\end{eqnarray} \tag{ 1 }$$

Here ${F}_{\lambda ,\mathrm{obs}}$ is the observed flux per wavelength in a particular band centered at the observed wavelength λ_obs and D_L is the luminosity distance to the quasar. The size r_λ is the scale radius of the disk at this wavelength, defined such that ${kT}({r}_{{\lambda }_{\mathrm{rest}}})={hc}/{\lambda }_{\mathrm{rest}}$ (see Equation (2) of Morgan et al. 2010 or Equation (5) of Li et al. 2019). The integral is over the dimensionless parameter u = r/r_λ. The factor of cos (i) accounts for disk inclination. For later generalization, we leave this as a function of the temperature profile slope β, where β = 3/4 for the Shakura & Sunyaev (1973) thin disk model. The factor C(β) is given in Li et al. (2019) as the numerical solution to

$$\begin{eqnarray}&&{\int }_{{u}_{\mathrm{in}}}^{C(\beta )}\displaystyle \frac{u\,{du}}{(\exp ({u}^{\beta })-1)}=\displaystyle \frac{1}{2}{\int }_{{u}_{\mathrm{in}}}^{\infty }\displaystyle \frac{u\,{du}}{(\exp ({u}^{\beta })-1)}\end{eqnarray} \tag{ 2 }$$

and converts the scale radius to a half light radius. To a good approximation, we can set u_in = 0 for optical and infrared λ_obs. For a standard thin disk this gives C(3/4) = 2.44.

Equation (1) provides estimates of the quasar size that are systematically smaller by a factor of ∼3–4 than the measurements from microlensing (Pooley et al. 2007; Morgan et al. 2010; Blackburne et al. 2011) and reverberation mapping (Cackett et al. 2007; Bentz et al. 2010; Edelson et al. 2015). A number of groups have recently posited a variety of new models, some of which might resolve this discrepancy. Dexter & Agol (2011) demonstrated that an inhomogeneous disk could reproduce the correct disk size, although the level of inhomogeneity required may not be physically plausible. The disk structure may be different at large radii (Jiang et al. 2016), leading to an adjustment in inferred luminosity sizes. A modified disk atmosphere, such as that proposed in Hall et al. (2018), would give rise to a nonblackbody spectrum that could also yield luminosity-based sizes consistent with microlensing sizes. Alternatively, the disk could have mass outflow in the form of disk winds (e.g., Proga 2005; Slone & Netzer 2012; Yuan et al. 2012; Laor & Davis 2014; Tombesi et al. 2015; Li et al. 2019), which could also match measured disk sizes. More complex physics, included in GRMHD simulations, can potentially also reproduce measured continuum emission region sizes.

A more sophisticated picture comes from constraining the shape of the disk continuum emission through modeling or observation. Keeping for now the assumption of blackbody emission within each annulus, and examining regions far from the inner disk edge, the effective surface temperature profile takes on the relatively simple form of

$$\begin{eqnarray}&&T(r)\propto {r}^{-\beta },\end{eqnarray} \tag{ 3 }$$

where T is the disk temperature, r is the radius from the center of the black hole, and β is the temperature profile slope. The thin disk model predicts a temperature profile slope of β = 3/4. The slim disk model of Abramowicz et al. (1988) predicts a shallower temperature profile slope of β ≈ 0.5 for high accretion rates (estimated from Figure 11 in Abramowicz et al. 1988). Other models include a nonzero torque on the inner disk boundary (Agol & Krolik 2000; Penna et al. 2012), and the MRI model of Agol & Krolik (2000) explicitly predicts a steeper value of β = 7/8. Sincell & Krolik (1998) account for the structure of a dissipative accretion disk and found that, although the ultraviolet spectrum is not a local blackbody, the spectral slope can be nearly consistent with the thin disk model. Other numerical models explore the impact of a radially varying disk viscosity (Penna et al. 2013) and MHD spinning black holes (Schnittman et al. 2016) though neither reports the effective temperature profile slope.

Similarly, attempts to empirically measure or estimate accretion disk temperature profiles using observational constraints have led to very inconsistent results. In the rest-frame optical and ultraviolet, Vanden Berk et al. (2001) and Gaskell (2008) estimated a temperature profile slope of β = 0.57 using the spectral continuum slope (Sun & Malkan 1989). A comparably shallow temperature profile slope is inferred in several similar studies (Davis et al. 2007; Bonning et al. 2013; Xie et al. 2016), but the authors emphasize that dust extinction in the host galaxy can significantly affect the observed temperature profile. This host extinction is difficult to correct for and gives rise to high uncertainty in the intrinsic temperature profile slope.

Reverberation mapping (RM) studies have measured the disk temperature profile through continuum interband lags (e.g., Cackett et al. 2007, 2018; Jiang et al. 2017; Edelson et al. 2019). These studies find profile slopes consistent with the thin-disk model, although Edelson et al. (2019) demonstrate that a wide range of profile slopes provide adequate fits to the observed lags. Furthermore, Cackett et al. (2018) find that the broad-line region can contribute to the interband lags, requiring a careful treatment of different AGN components during reverberation mapping analysis. Although the results are consistent with the thin-disk model, reverberation mapping does not presently exclude shallower or steeper temperature profile slopes.

Temperature profile measurements also arise naturally from single-epoch microlensing measurements. Chromatic flux ratio anomalies can lead to estimates of the quasar accretion disk size, which scales as r ∝ λ^1/β. The findings from different studies have not, however, converged to a single value. Some studies find evidence for a shallower profile than the thin-disk model (Rojas et al. 2014; Bate et al. 2018) though more studies point to a steeper profile (Blackburne et al. 2011; Muñoz et al. 2011, 2016; Motta et al. 2012, 2017; Jiménez-Vicente et al. 2014; Blackburne et al. 2015). Still other studies based on multiepoch, multifilter curves conclude that a thin-disk model can explain observations (e.g., Eigenbrod et al. 2008; Mosquera et al. 2011). Studies such as that by Bate et al. (2018) have begun to explore the systematic uncertainties that have lead to such a wide spread in measurements. They demonstrated that in systems with low chromatic variation, the single-epoch method will return an overly steep profile, and they propose a means to quantify this effect in future work.

Morgan et al. (2010) developed an alternate, model dependent, means to constrain the temperature profile using the observed difference between the accretion disk size measured by microlensing, r_μ (e.g., Kochanek 2004; Bate et al. 2008) and that inferred from the quasar luminosity, r_L. By modeling the expected ratio r_L/r_μ as a function of temperature profile slope, constraints on β can be visualized. In Morgan et al. (2010) the authors found evidence for a shallower-than-thin-disk temperature profile, but they did not exclude thin disk models. With a larger sample size, we are now able to employ this technique to tightly constrain the slope of the accretion disk temperature profile.

This paper is organized as follows. In Section 2 we present the discrepancy between microlensing and luminosity-based continuum sizes in our sample of quasars with multiepoch microlensing lightcurves. In Section 3 we show how microlensing and luminosity sizes might be brought into agreement with a different temperature profile slope β. We also explore the impact of contaminating light sources. We compare our results to alternate disk models, focusing on two recent analytic models from Hall et al. (2018) and Li et al. (2019). We summarize our conclusions and discuss future observational programs in Section 4. In all calculations we adopt a flat cosmological model with Ω_M = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹ (Hinshaw et al. 2009).

Accretion disk microlensing size measurements are consistent between the single and multiepoch approaches (e.g., Pooley et al. 2007; Morgan et al. 2010; Blackburne et al. 2011). The multiepoch Monte Carlo lightcurve analysis technique of Kochanek (2004) is, however, effectively independent of the most sensitive statistical prior, the mean mass of a microlens star in the lens galaxy. In this technique, the quasar is modeled as a disk of varying size crossing a range of possible stellar magnification fields. By coupling these results with cosmological velocity models, one does not need to assume a mean microlens star mass. We therefore restrict our analysis to only those systems with microlensing sizes measured from multiepoch lightcurves, given in Table 1. These values and uncertainties are found from the median and 16th to 84th percentile range of the reported posterior microlensing size distributions. These 15 systems are shown in Figure 1 with both microlensing and luminosity sizes as a function of black hole mass, where the size offset is readily apparent (Kochanek 2004; Morgan et al. 2008, 2010, 2012, 2018; Dai et al. 2010; Hainline et al. 2012, 2013; MacLeod et al. 2015; Cornachione et al. 2020).

Zoom In Zoom Out Reset image size

Table 1.  Microlensing Size Measurements of the Continuum Size $\mathrm{log}({r}_{\mu })$, Luminosity Continuum Size Estimates for a Standard Thin Disk (β = 3/4), $\mathrm{log}({r}_{L})$, and the Ratio, $\mathrm{log}({r}_{L}/{r}_{\mu })$, for Our Sample of 15 Lensed Quasars

Quasar	$\mathrm{log}({r}_{\mu }/\mathrm{cm})$	$\mathrm{log}({r}_{L}/\mathrm{cm})$	$\mathrm{log}({r}_{L}/{r}_{\mu })$
QJ0158	${15.87}_{-0.15}^{+0.30}$	15.40 ± 0.05	−0.47 ± 0.23
HE0435	${16.07}_{-0.73}^{+0.55}$	15.16 ± 0.10	−0.91 ± 0.65
SDSS0924	${15.80}_{-0.48}^{+0.28}$	15.03 ± 0.10	−0.77 ± 0.39
FBQ0951	${16.32}_{-0.40}^{+0.42}$	15.77 ± 0.04	−0.55 ± 0.42
SDSS1004	${15.37}_{-0.29}^{+0.29}$	15.13 ± 0.18	−0.24 ± 0.34
HE1104	${16.34}_{-0.29}^{+0.22}$	15.79 ± 0.13	−0.55 ± 0.28
PG1115	${16.99}_{-0.38}^{+0.34}$	15.42 ± 0.11	−1.57 ± 0.38
RXJ1131	${15.35}_{-0.22}^{+0.21}$	14.80 ± 0.04	−0.54 ± 0.22
SDSS1138	${15.44}_{-0.60}^{+0.59}$	15.04 ± 0.08	−0.40 ± 0.60
SBS1520	${16.08}_{-0.22}^{+0.23}$	15.56 ± 0.05	−0.52 ± 0.23
Q2237	${16.12}_{-0.30}^{+0.27}$	15.74 ± 0.18	−0.39 ± 0.34
WFI2033	${16.41}_{-0.30}^{+0.30}$	15.31 ± 0.10	−1.10 ± 0.32
SBS0909	${15.66}_{-0.30}^{+0.30}$	15.76 ± 0.04	0.10 ± 0.30
Q0957	${16.57}_{-0.50}^{+0.50}$	15.76 ± 0.03	−0.81 ± 0.50
WFI2026	${16.13}_{-0.34}^{+0.29}$	15.47 ± 0.12	−0.66 ± 0.34

Note. All sizes are scaled to 2500 Å in the quasar rest frame and a nominal inclination angle of 60°. Microlensing size errors are found from the 16th and 84th percentile values of the posterior size distributions. Luminosity size uncertainty is calculated from the intrinsic flux uncertainty, which includes both photometric error and magnification uncertainty.

Download table as:  ASCIITypeset image

In lensed quasars we observe individual image fluxes. Using these and magnifications from strong lensing models we calculate a quasar's intrinsic luminosity. To estimate the magnification, most of these studies generated a range of galaxy mass models with varying smooth matter fractions (see, for example, Section 4.1 of Morgan et al. 2018). These sequences give both the typical magnification and the magnification uncertainty range. The estimated lens magnification uncertainty usually dominates the photometric measurement error, leaving the typical intrinsic flux uncertainty on the order of 0.1 dex. Adopting these flux values and uncertainties in each system and, for the moment, assuming a standard thin disk model, we can find the luminosity-based continuum emission region sizes, reported as a half light radius r_L, using Equation (1). These luminosity-based size estimates are given in Table 1, along with the ratio $\mathrm{log}({r}_{L}/{r}_{\mu })$. Since r_L and r_μ measure the same physical quantity, this ratio must be unity for any valid disk model.

Of these 15 quasars, the thin disk luminosity sizes are smaller than the microlensing sizes in all but one system as seen in Figure 2. Only SBS 0909+532 has a luminosity size larger than the microlensing size, but the intrinsic luminosity of this system and luminosity-based size estimate are uncertain because of significant systematic differences between the lens models and magnifications predicted by the astrometric fits in Lehár et al. (2000) and Sluse et al. (2012). Nevertheless, we include this system for completeness. The error-weighted mean offset between the luminosity and microlensing size, $\mathrm{log}({r}_{L}/{r}_{\mu })$, is −0.57 ± 0.08 dex, assuming standard thin disk theory (β = 3/4). We use this ratio as a constraint in subsequent modeling.

Zoom In Zoom Out Reset image size

The ratio $\mathrm{log}({r}_{L}/{r}_{\mu })$ shows no significant dependence on black hole mass. We do find, however, that the measured microlensing sizes scale with black mass as ${r}_{2500\mathring{\rm A} }\propto {M}_{\mathrm{BH}}^{0.66\pm 0.15}$ (Morgan et al. 2018). This scaling is nearly identical to the thin disk prediction of $r\propto {M}_{\mathrm{BH}}^{2/3}$, a curious result given the clear inconsistency between the thin disk luminosity sizes and microlensing measurements. This empirical scaling relation with black hole mass must still be reproduced by an alternate disk model, a point we explore in Section 3.3.

In this work we explore several possible solutions to the size inconsistency described in Section 2. First, we quantify the ratio offset, $\mathrm{log}({r}_{L}/{r}_{\mu })$, as a function of β, adopting the effective surface temperature profile in Equation (3), $T(r)\propto {r}^{-\beta }$. We calculate the dependence of the luminosity continuum sizes on the temperature profile slope. While microlensing size measurements are effectively independent of continuum region shape, we consider the impact on these measurements due to contaminating flux from large scale emission. We also recognize that quasar accretion disks may not be homogenous, optically thick emitters, in which case the generalized temperature profile slope of Equation (3) need not describe the surface temperature. We depart from the simplest assumptions and consider alternate disk models that could reconcile microlensing size measurements with observed fluxes. Finally, we show that the empirical scaling of continuum emission size with black hole mass provides another test for any proposed disk model.

3.1. Size Ratio versus Temperature Profile Slope

Based on the findings of Mortonson et al. (2005), size, not shape, dominates microlensing variability. This means that the microlensing continuum size measurement, r_μ, is effectively independent of β. Although Mortonson et al. (2005) and Congdon et al. (2007) do find a small difference between microlensed images due to shape, any changes in r_μ with β will have a minor influence on the observed size ratio. The dominant β dependence, therefore, comes from the luminosity size estimate, r_L.

The luminosity continuum size estimate, however, varies significantly with changes in β. We use luminosity sizes from Equation (1) and numerically calculate $\mathrm{log}({r}_{L}/{r}_{\mu })$ versus β. We plot the resulting curve as a solid line in Figure 3, fixing the inner disk edge to r_in = 0. The yellow circle indicates the constraint from Section 2 of $\mathrm{log}({r}_{L}/{r}_{\mu })=-0.57\pm 0.08$ at β = 3/4, while the dashed horizontal lines indicate the region in which continuum luminosity and microlensing sizes are consistent, $\mathrm{log}({r}_{L}/{r}_{\mu })=0\pm 0.08\,\mathrm{dex}$. Interpreted as constraints on β, this implies 0.37 < β < 0.43 (68% confidence). This value of β ≈ 0.4 is significantly shallower than the thin disk value of β = 0.75.

Zoom In Zoom Out Reset image size

We have so far approximated the inner disk edge at r_in = 0, but a finite inner edge will partially reconcile the size inconsistency. To quantify the impact of the inner edge, we modify the temperature profile of Equation (3) as

$$\begin{eqnarray}&&T(r)\propto {r}^{-\beta }{\left(1-\sqrt{\displaystyle \frac{{r}_{\mathrm{in}}}{r}}\right)}^{(1/4)}.\end{eqnarray} \tag{ 4 }$$

This is holds exactly for the thin disk model; though, it only captures the approximate behavior for general values of β. Based on the disk sizes and black hole masses in Table 1, we choose to fix the inner edge at a typical value of r_in = 0.1r_λ, 2500 _Å. This is roughly at or immediately outside of the innermost stable circular orbit of the central black hole, conservatively assuming a nonspinning (Schwarzchild) black hole. Under this modification, we find the curve indicated by the dotted line in Figure 3. For high values of β the inner edge significantly alters the expected ratio. In the gray region constrained by the data, however, the addition of the inner edge makes a negligible change in our inferred β.

Although microlensing sizes do not change with disk shape, they can be impacted by contaminating flux from large scales. Emission lines in the broad-line region, for example, account for significant quasar flux. Even if a line center does not fall within our observing band, wings of the lines can contribute flux on the order of 30% at λ_rest = 2500 Å (e.g., Maoz et al. 1993; Sluse et al. 2007). This effect has been analyzed in previous multiepoch microlensing studies (e.g., Dai et al. 2010) and is shown to reduce the measured microlensing size r_μ by 20%–50% for 30% contamination. This shows that the measured microlensing size decreases roughly in proportion to the broad-line region contamination, ${r}_{\mu }\propto 1-{F}_{\mathrm{BLR}}/{F}_{\mathrm{tot}}$. Though we expect the contamination in our aggregated sample to be less than 30%, we adopt this value as an approximate upper limit.

Broad-line contamination will also reduce the estimated luminosity size r_L, although by a smaller percentage because our models scale as ${r}_{L}\propto {F}_{\lambda ,\mathrm{obs}}^{1/2}$, following Equation (1). For 30% contamination this decreases r_L by 16%. The scale factor C(β) also changes because the continuum half light radius is no longer formally determined by the accretion disk alone, although this effect is smaller than size decrease due to lower flux. We plot the combined effect of this contamination on $\mathrm{log}({r}_{L}/{r}_{\mu })$ as the dotted–dashed line in Figure 3. Broad-line contamination brings the inferred β closer to the thin disk prediction, giving 0.44 < β < 0.52. While a more significant broad-line contamination would bring the inferred β even closer to the thin disk prediction, a contamination percentage significantly greater than 30% seems physically unlikely.

Light emitted from the accretion disk may also be scattered back into the line of sight at larger scales (e.g., Dai et al. 2010). This can come from scattering off of the outer disk edge or torus or as polar scattering from the surrounding medium (Goosmann & Gaskell 2007). We approximated this scattered light as a constant additive term, which will not be microlensed because of the large scattering scales. This will affect the microlensing size in the same way as broad-line region contamination, decreasing r_μ proportionally to the scattered light percentage. Because scattered light originates from the disk, however, the disk luminosity, and therefore luminosity size, r_L, are effectively unchanged, although we still adjust C(β) as with broad-line contamination.

The scattered light percentage can be constrained through polarization measurements. For example, Hutsemékers et al. (2015) explore a scattered light percentage between 17% and 44% in the lensed quasar H1413+117; though, they do not favor any particular value. The findings of Goosmann & Gaskell (2007) suggest that quasars with low typical inclination angles will likely have less scattered light than Seyferts. Based on these ranges, we adopt a typical scattered light percentage of 30% and show the effect of scattering on the temperature profile as the dashed line in Figure 3. Scattering has the largest impact on the inferred β, allowing values as high as 0.48 < β < 0.56, which gives β < 0.75 at ∼4σ significance. While the true scattered light percentage could exceed 30%, there is at present no strong evidence of high scattering fractions. In such a case, scattered light would need to compose 61% of the total flux to bring our temperature slope into agreement with the thin disk prediction at 2σ.

Accounting for this range of systematic shifts, we find an improved constraint on quasar accretion disk temperature profiles 0.37 < β < 0.56. We plot the temperature profile the central value β = 0.47 in Figure 4, including the steeper thin disk temperature profile for comparison. Even if we allow for the case of both 30% scattered light and 30% broad-line region contamination we still find β < 0.75 at 3σ confidence.

Zoom In Zoom Out Reset image size

This range of temperature profile slopes could be produced, for example, by the phenomenological disk wind model of Li et al. (2019). In this model mass outflows from an underlying thin disk increase the sizes measured by microlensing or reverberation mapping techniques. This model is completely analogous to ours, in which a disk wind parameter s gives an effective temperature profile identical to Equation (3) with β = (3 − s)/4. Our estimate of β ≈ 0.47 translates to a disk wind parameter of s ≈ 1.1, consistent with their findings. While the model of Li et al. (2019) does not arise directly from accretion disk physics, a similar temperature profile slope may be produced in other disk wind models (e.g., Slone & Netzer 2012; Yuan et al. 2012; Laor & Davis 2014; Tombesi et al. 2015). In particular, Figure 5 of Laor & Davis (2014) demonstrates that their disk wind model flattens the effective temperature profile.

3.2. Size Ratio Variation with Disk Atmospheres

Our analysis so far assumes that accretion disks are homogenous with optically thick blackbody emission varying with radius with a power-law form for the surface temperature. However, none of these assumptions need necessarily hold. For example, Hall et al. (2014) showed that shock heating could alter a disk's temperature profile and inferred size. In another work Jiang et al. (2016) found that an iron opacity bump changes disk structure at radii comparable to our measured half light radii. This could manifest in a broken power-law temperature profile. Hall et al. (2018) demonstrated that with optically thin scattering atmospheres, disks can have a flatter observed spectral energy distribution than that of the fiducial optically thick disk. Depending on the characteristics of the disk atmosphere, these nonblackbody spectra can produce luminosity size estimates consistent with microlensing size measurements. These models actually produce a steeper surface temperature profile slope than the thin disk, but the spectral slope of the emergent flux is broadened by the nonblackbody atmosphere.

For a comparison to the nonblackbody atmosphere model of Hall et al. (2018), we must numerically calculate r_L from the definition of the half-light radius

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\nu ,1/2}}{{L}_{\nu ,\mathrm{rest}}}=\displaystyle \frac{{\int }_{{r}_{\mathrm{in}}}^{{r}_{L}}2\pi r\,{F}_{\nu }(r){dr}}{{\int }_{{r}_{\mathrm{in}}}^{\infty }2\pi r\,{F}_{\nu }(r){dr}}=\displaystyle \frac{1}{2}.\end{eqnarray} \tag{ 5 }$$

Here all frequencies are in the rest frame of the quasar. This equation retains the assumption of constant flux in each annular ring, but allows a generalized flux spectrum. For F_ν we use the constant atmosphere flux model, Equation (1) of Hall et al. (2018), where the surface temperature T_s(r) is found according to Hall et al. (2018) Section 2.4. We show this numerical T_s(r) for a typical atmospheric density of $\mathrm{log}({\rho }_{\mathrm{atm}}/\rho )=-4$ in Figure 4. As discussed in Hall et al. (2018), the surface temperature approaches the thin disk temperature from above at high radii. At small radii, however, deviations from blackbody emission lead to T_s that is much higher than the thin disk.

We enforce that the model reproduces the observed flux, ${F}_{\nu ,\mathrm{obs}}$, by

$$\begin{eqnarray}&&{L}_{\nu ,\mathrm{rest}}=4\pi {D}_{L}^{2}{F}_{\nu ,\mathrm{obs}}/(\chi (i)(1+z))=2{\int }_{{r}_{\mathrm{in}}}^{\infty }2\pi r\,{F}_{\nu ,\mathrm{rest}}(r){dr}.\end{eqnarray} \tag{ 6 }$$

Here z is the quasar redshift, D_L is the luminosity distance, and i is the viewing inclination. The factor χ(i) corrects for anisotropy in the emitted flux, which can be due to the geometry of the emitting region, limb darkening, or general relativistic warping near the inner disk edge (Krolik 1999; Misra & Sriram 2003). For our model we assume a flat thin disk and neglect relativistic effects, which gives $\chi (i)=2\cos (i)$ (Misra & Sriram 2003). We fix the inclination angle to 60° for consistency with reported microlensing sizes. As in Morgan et al. (2010) we adopt a nominal Eddington fraction of f_edd = L/L_E = 1/3 and an efficiency of η = 0.06 for each model.

From the observed fluxes in the 15 quasars, we calculate the estimated luminosity size r_L for the thin disk, our β model (which also holds for the Li et al. 2019 model), and the constant atmosphere disk model in Hall et al. (2018). Here we use the microlensing sizes r_μ (without considering any large-scale flux contamination). For the Hall et al. (2018) model the outer disk atmospheric density at each radius is assumed to be a constant ratio of the central density. We select two scale factors, $\mathrm{log}({\rho }_{\mathrm{atm}}/\rho )=-3$ and $\mathrm{log}({\rho }_{\mathrm{atm}}/\rho )=-4$, which are favored by Hall et al. (2018). We show the weighted mean ratio, expressed as $\mathrm{log}({r}_{L}/{r}_{\mu })$, in Table 2, where we see that our β ≈ 0.4 disk and the constant atmosphere model with $\mathrm{log}({\rho }_{\mathrm{atm}}/\rho )=-4$ yield consistent luminosity sizes.

Table 2.  The Weighted Mean Ratio of Continuum Sizes Estimated from Luminosity to Those Measured by Microlensing, $\mathrm{log}({r}_{L}/{r}_{\mu })$

Model	$\mathrm{log}({r}_{L}/{r}_{\mu })$
Thin Disk	−0.57 ± 0.08
β = 0.4 Disk	−0.02 ± 0.08
Constant Atmosphere Disk $\mathrm{log}({\rho }_{\mathrm{atm}}/\rho )=-3$	−0.29 ± 0.08
Constant Atmosphere Disk $\mathrm{log}({\rho }_{\mathrm{atm}}/\rho )=-4$	−0.10 ± 0.08

Note. We show the change for the thin disk model, a β disk model at β = 0.4 (which corresponds to the disk wind model of Li et al. 2019 with s = 1.4), and the Hall et al. (2018) constant atmosphere model under two different densities.

Download table as:  ASCIITypeset image

3.3. Continuum Size Variation with Black Hole Mass

As another means to distinguish models, we compare predictions for the continuum emission size as a function of black hole mass. Here we use the constant atmosphere model with the range of atmospheric densities $-1\leqslant \mathrm{log}({\rho }_{\mathrm{atm}}/\rho )\leqslant -5$. For the disk wind model, we explore the range of disk wind parameters 0.3 ≤ s ≤ 1.1. Figure 5 illustrates that under these ranges of atmospheric densities and disk wind parameters, both the models of Hall et al. (2018) and Li et al. (2019) can reproduce the empirical relationship.

Zoom In Zoom Out Reset image size

These two models do, however, predict different slopes of r_{2500 Å} versus M_BH. The constant atmosphere model of Hall et al. (2018), although not strictly a power law, gives an approximate slope of $0.74\lt \mathrm{log}({r}_{2500\mathring{\rm A} })/\mathrm{log}({M}_{\mathrm{BH}})\lt 0.88$. The model of Li et al. (2019) predicts a shallower slope of $0.47\lt \mathrm{log}({r}_{2500\mathring{\rm A} })/\mathrm{log}({M}_{\mathrm{BH}})\lt 0.63$. While both are consistent with the measured scaling $\mathrm{log}({r}_{2500\mathring{\rm A} })/\mathrm{log}({M}_{\mathrm{BH}})\,=0.66\pm 0.15$, higher precision measurements may soon enable a distinction between these models.

We have found that the inconsistency between microlensing size measurements and luminosity size estimates of quasar UV continuum emission regions can be reconciled if the quasar accretion disk temperature profile is shallower than that of the thin disk. The size difference cannot be completely explained by contamination from other emission sources nor by the impact of the inner disk edge. We conclude that for a simple optically thick continuum emission region, the effective temperature profile must be much shallower than a standard thin disk to reproduce observed microlensing sizes.

Our findings are in excellent agreement with those from spectral indexing, under models with little to no dust extinction, which measure β ≈ 0.57 (Davis et al. 2007; Gaskell 2008). Reverberation mapping (e.g., Edelson et al. 2017, 2019) generally measures slopes in agreement with the thin disk model, although Jiang et al. (2017) found evidence for a shallower temperature profile. Furthermore, most RM findings do not exclude the temperature profile we found here.

The shallower temperature profile we find here is, however, in tension with many chromatic microlensing studies. Most chromatic microlensing measurements give a steeper-than-thin disk profile of β ≈ 1 (Muñoz et al. 2011, 2016; Motta et al. 2012; Blackburne et al. 2015) while some return a still steeper estimate of β ≈ 1.25 or higher (Blackburne et al. 2011; Jiménez-Vicente et al. 2014; Motta et al. 2017). These measurements are inconsistent with our estimate of β < 0.56. Bate et al. (2018) recently ran simulations that found that low chromatic variability in a lensed quasar system can lead to systematically steeper temperature profile slopes. Indeed Rojas et al. (2014) and Bate et al. (2018) analyzed several systems with high chromaticity and measured temperature profile slopes as low as β = 0.43. However, the analysis of other high chromaticity lenses return values of β ≈ 0.9, still steeper than the thin disk prediction (Muñoz et al. 2011; Motta et al. 2012). Bate et al. (2018) point out that their simulations may not be completely general and it remains unclear whether a better understanding of chromaticity alone will be sufficient to reconcile chromatic microlensing findings with ours.

This inconsistency may also highlight the limitations of our own model dependent estimate of β. If the underlying disk structure deviates strongly from the standard thin disk model, for example, with large disk inhomogeneities (Dexter & Agol 2011), then our analysis may need significant modification. We also have not accounted for the width of our filters or potential geometric effects from the nonzero thickness of the disk; though, we expect these to have only a modest impact our inferred β. Furthermore, we have assumed that the size offset comes from a change in temperature profile slope, though a change in temperature profile normalization could also reconcile the size inconsistency. A physical mechanism for such a normalization offset is, however, not immediately apparent and is difficult to produce under the standard thin disk model. This would require, for example, a significant systematic offset in black hole masses.

From Section 3, it is clear that multiple models are consistent with our measurements, even including the constraint from black hole masses. Our value of β ≈ 0.47 is also consistent with the temperature profile slope predicted by the slim disk models of Abramowicz et al. (1988) and Szuszkiewicz et al. (1996). Takeo et al. (2020) presented a model of a rapidly accreting black hole which also predicts a scaling of β = 0.5 (see Equation (26)). Furthermore Dexter & Agol (2011) showed that disk inhomogeneities could also flatten the spectrum and explain the observed size difference, though the authors caution that this model is not necessarily physical. Iron opacity could also lead to a broken power-law temperature profile slope (Jiang et al. 2016); though, we presently lack sufficient observational constraints from microlensing alone to test this model. Comparison of these models to other quasar observations (e.g., Davis et al. 2007; Shull et al. 2012) will be necessary to fully confirm or refute any proposed models.

We have not considered here any of the recent numerical GRMHD models. We expect that these simulations will generally predict disks with nonblackbody spectra and inhomogeneities, likely including a departure from our assumption of constant annular flux. Further publications of GRMHD disk azimuthally averaged surface temperature predictions, or better still, half-light continuum size predictions, along with size scaling with black hole mass, would enable ready comparison with our measurements. Y.-F. Jiang et al. (2019, private communication), showed that their model readily produces values of β ≈ 0.34 above ∼10 gravitational radii, suggesting that this simulation qualitatively agrees with our findings, but further work is needed to produce a robust comparison.

Our results as given here are still model dependent, hinging on annular blackbody flux and a well-behaved disk surface temperature profile slope. Our next step is to conduct multiwavelength analyses of lensed quasars using the multiepoch variability technique, expanding to H-band lightcurves, which will measure the quasar continuum region size in the rest-frame optical. While observationally expensive, this is the most bias-free method available to concretely measure disk temperature profiles in individual systems and will be independent of any accretion disk model. Fortunately the Large Synoptic Survey Telescope (LSST; LSST Science Collaboration et al. 2009) will provide many lightcurves suitable for such an analysis.

We thank Jason Dexter for helpful discussions on quasar accretion disk modeling. We also thank Yan-Fei Jiang for generating the temperature profile in his recent simulation, for providing estimates on the predicted temperature profile slope, and for helpful comments on this work. We are also grateful to Patrick Hall for providing valuable feedback. We thank the anonymous referee for helpful comments that improved this manuscript.

This material is based upon work supported by the National Science Foundation under grant AST-1614018 to M.A.C. and C.W.M.