Constraining AGN Torus Sizes with Optical and Mid-Infrared Ensemble Structure Functions

TL;DR

This paper introduces a novel ensemble structure function method to measure the size of AGN dusty tori using optical and MIR variability data, providing a less demanding alternative to reverberation mapping.

Contribution

The paper develops and validates a new SF-based approach to constrain AGN torus sizes from optical and MIR light curves, demonstrating its effectiveness with SDSS quasar data.

Findings

Torus size correlates with luminosity as log R_eff (pc) = 0.51 * log(L_bol/10^46 erg/s) - 0.38.

The MIR torus sizes are approximately 1.4 times larger in W2 band than in W1 band.

The SF method agrees well with dust reverberation mapping results.

Abstract

We propose a new method to constrain the size of the dusty torus in broad-line active galactic nuclei (AGNs) using optical and mid-infrared (MIR) ensemble structure functions (SFs). Because of the geometric dilution of the torus, the MIR response to optical continuum variations has suppressed variability with respect to the optical that depends on the geometry (e.g., size, orientation, opening angle) of the torus. More extended tori have steeper MIR SFs with respect to the optical SFs. We demonstrate the feasibility of this SF approach using simulated AGN light curves and a geometric torus model. While it is difficult to use SFs to constrain the orientation and opening angle due to insensitivity of the SF on these parameters, the size of the torus can be well determined. Applying this method to the ensemble SFs measured for 587 SDSS quasars, we measure a torus $R-L$ relation of ${\rm log}\,R_{\rm eff}\, {(\rm pc)} = 0.51_{-0.04}^{+0.04} \times {\rm log}\,(L_{\rm bol}/10^{46}\,\rm erg\,s^{-1}) -0.38_{-0.01}^{+0.01}$ in the WISE $W1$ band, and sizes $\sim1.4$ times larger in the $W2$ band, which are in good agreement with dust reverberation mapping measurements. Compared with the reverberation mapping technique, the SF method is much less demanding in data quality and can be applied to any optical+MIR light curves for which a lag measurement may not be possible, as long as the variability process and torus structure are stationary. While this SF method does not extract all information contained in the light curves (i.e., the transfer function), it provides an intuitive interpretation for the observed trends of AGN MIR SFs compared with optical SFs.

Tables

Table 1. Table 2: Parameters used to generate simulated optical and MIR light cuves.

parameter	description
${\mathrm{SF}}_{\infty }$	optical ${\mathrm{SF}}_{\infty }$ in Eq. 2
$\tau$	optical damping timescale (day) in Eq. 2
$\beta$	optical slope in Eq. 2; fixed to 1.0
$f_{g,\mathrm{AD}}$	mean flux in the optical; fixed to 1.0
$\theta$	inclination angle; $\theta =0^{\circ }$ means face-on
$R_{\mathrm{in}}$	torus inner radius (pc)
$Y$	ratio of outer to inner radius; fixed to 1.5
$\sigma$	half-opening angle (degree)
$p$	radial powerlaw index; fixed to –1.0
$f_{w,\mathrm{torus}}/f_{g,\mathrm{AD}}$	mean MIR/optical flux ratio

Equations

$$k(\Delta t)={\sigma }^2\exp (-|\Delta t|/\tau ),$$

$$\mathrm{SF}(\Delta t)=\text{SF∞}\sqrt{1-e^{-(|\Delta t|/\tau {)}^{\beta }}},$$

$$\mathrm{SF}(\Delta t)=\sqrt{\frac{1}{N_{\Delta \mathrm{t},\mathrm{pairs}}}\sum _{i=1}^{N_{\mathrm{\Delta t},\mathrm{pairs}}}\Delta f^2},$$

$$\Delta f^2=(f(t)-f(t+\Delta t){)}^2-{\sigma }_f(t{)}^2-{\sigma }_f(t+\Delta t{)}^2.$$

$${\text{SFw}}^2={\text{SFw,torus}}^2+{\text{SFw,AD}}^2+2\rho \times \text{SFw,torus}\text{SFw,AD},$$

$$F_{\mathrm{MIR}}(t)=\text{famp}\int F_{\mathrm{optical}}(t-\tau )\text{ϕ(τ)}d\tau ,$$

$${\tau }_{\mathrm{lag}}=\frac{\int \tau \varphi (\tau )d\tau }{\int \varphi (\tau )d\tau },$$

Full text

black hole physics — galaxies: active — quasars: general — surveys

††software: astropy (Astropy Collaboration et al., 2013), celerite (Foreman-Mackey et al., 2017), emcee (Foreman-Mackey et al., 2013)

1 Introduction

Active galactic nuclei (AGNs) are powered by accreting supermassive black holes (SMBHs) that produce enormous quantities of radiation across the entire electromagnetic wavelength range. It is widely accepted that the inner region of an AGN consists of a corona that produces compact X-ray emissions and an accretion disk that emits thermal blackbody spectrum in the UV/optical bands (Shakura & Sunyaev, 1973; Haardt & Maraschi, 1993), although the exact structure and geometry of the X-ray corona and accretion disk are still debated (e.g., Cackett et al., 2021). The central SMBH is surrounded by a dusty torus that absorbs UV/optical photons and re-emits them in the infrared, producing a MIR bump in their spectral energy distribution (SED; Elvis et al. 1994; Lyu et al. 2017). The inclination angle and covering factor of the dusty torus are key ingredients in the AGN unification model that determine the appearance of AGNs as unobscured/Type 1s and obscured/Type 2s (Antonucci, 1993; Netzer, 2015).

Variability is a ubiquitous feature of AGNs (Ulrich et al., 1997). It offers a unique means to investigate the inner structure of distant AGNs that cannot be spatially resolved, such as measuring the size of the accretion disk (a few light days; e.g., Mudd et al. 2018; Yu et al. 2020) and the dusty torus ($\sim 0.01-1$ pc; e.g., Suganuma et al. 2006; Kishimoto et al. 2007; Koshida et al. 2014; Minezaki et al. 2019; Lyu et al. 2019; Yang et al. 2020) through reverberation mapping. The variability of the continuum emission of AGNs in the optical bands can be well characterized by the damped random walk (DRW) model over months to years timescales, although its physical origin is still unclear (Kelly et al., 2009; MacLeod et al., 2010; Zu et al., 2013; Burke et al., 2021; Stone et al., 2022). The DRW model describes AGN variability as a stochastic process using an exponential covariance function

[TABLE]

where $k(\Delta t)$ measures the correlation between two observations separated by time $\Delta t$, $\sigma$ is the driving rms variability amplitude, and $\tau$ is the characteristic decorrelation (damping) timescale for returning to the mean flux. In general, the rms variability, measured as a function of $\Delta t$ (known as the structure function; SF), can be expressed by a generalized power exponential covariance as

[TABLE]

where a parameter $\beta$ is introduced to describe different shapes of the SF (e.g., Zu et al., 2013; Kozłowski, 2016).

In the optical regime, the SF follows a power-law form at $\Delta t\ll\tau$ with a slope of $\beta\sim 1$ (i.e., a DRW process) and levels out to a constant $\rm SF_{\infty}$ at $\Delta t\gg\tau$. Its amplitude ${\rm SF}_{\infty}$ is found to increase with decreasing luminosity and wavelength, and the decorrelation timescale is correlated with black hole mass over $\sim 6$ orders of magnitude, offering new insights into the origin of AGN variability (e.g., Vanden Berk et al., 2004; Morganson et al., 2014; Kozłowski, 2016; Li et al., 2018; Burke et al., 2021; Stone et al., 2022).

Compared with optical studies, far less is known about the variability properties in the MIR. The MIR SF shape appears similar to that in the optical, but the power-law slope is notably steeper and the variability amplitude smaller. Moreover, its slope appears uncorrelated with AGN luminosity and wavelength, while its rms variability is anti-correlated with luminosity (e.g., Kozłowski et al., 2010, 2016; Wang & Shi, 2020). These results have been qualitatively understood as the smoothing effect, in which the MIR variability is caused by reprocessing of the variable UV/optical emission, whereas the responses are diluted by the extended structure of the torus. However, a detailed study connecting the shape of the MIR SF to the torus structure and geometry is still absent from the literature.

In this paper, we investigate the ensemble variability of quasars in the optical and MIR bands, making use of optical light curves from all available ground-based surveys in the SDSS Stripe 82 region (as compiled in Yang et al., 2020) and MIR light curves from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010; Mainzer et al. 2011). We will measure the ensemble SFs and directly infer the torus geometry from these second-order variability statistics. This paper is structured as follows. The data and sample are described in Section 2. The ensemble SF measurements are presented in Section 3. In Section 4 we illustrate our method through simulated light curves and infer the torus parameters for SDSS Stripe 82 quasars. The conclusions are summarized in Section 5. Throughout this paper, timescales by default are in the restframe of the quasar. For luminosity calculations, we adopt a flat $\Lambda$CDM cosmology with $\Omega_{\Lambda}=0.7$ and $H_{0}=70\,{\rm km\,s^{-1}Mpc^{-1}}$.

2 Sample and Data

To study quasar variability in the optical and MIR bands and measure the corresponding SFs, we utilize the quasar sample and their multiband light curves compiled in Yang et al. (2020) (Y20). The Y20 sample contains 587 spectroscopically confirmed SDSS quasars at $0.2\lesssim z\lesssim 2$ in the Stripe 82 region, which is selected from a parent sample of 7384 quasars by applying several cuts on variability significance and data quality111We have tested that using the entire 7384 quasars does not improve the SF measurements, given that over 80% of the sample have MIR variability amplitudes below the noise level, making the SF measurements at short timescales noisy.. This sample has well measured spectroscopic properties, such as BH mass and bolometric luminosity ($L_{\rm bol}$) from the Shen et al. (2011) catalog. The optical light curves from various ground-based imaging surveys have an observed baseline of $\sim 20$ years. The WISE light curves used in Y20 have a baseline of $\sim 10$ years (2009 Dec – 2011 Feb, 2013 Dec – 2019 Dec) and a cadence of 6 months (see Table 1 in Y20 for the detailed survey information). Y20 performed cross-correlation analysis for each quasar and measured the time lag of the light curves between the $g$ band and the WISE $W1$ band. The lag is then converted to the average size of the torus assuming the time delay is the light-crossing time from the accretion disk to the dusty torus. Remarkably, the size of the torus follows a tight $R-L$ relation over $\sim 4$ orders of magnitude in quasar luminosity when combined with local measurements of low-luminosity AGNs (Kishimoto et al., 2007; Koshida et al., 2014; Lyu et al., 2019; Minezaki et al., 2019; Yang et al., 2020).

In this work, we supplement the Y20 data with the latest NEOWISE 2022 data release (2020 Dec – 2021 Dec) which extends the WISE baseline to $\sim 12$ years ($\sim 19$ observing epochs). Each WISE data point is the coadded (median) photometry of $\sim 8$ observations obtained within a short $\sim 1$-day window over every 6 months. The uncertainty on the coadd photometry ($\sigma_{f}$) is derived by randomly perturbing each individual detection based on its photometric error (assuming a Gaussian-distributed error) within the 1-day window for 1000 times which produces 1000 median values, then use the semi-amplitude of the enclosed 16th and 84th percentiles on the sampled median to represent the $1\sigma$ uncertainty. We measure the optical SF in the bluest $g$ band222Y20 compiled light curves in the $g$, $r$, $i$, $V$ and unfiltered bands from various surveys and converted them into $g$ band to increase the cadence. Here we only use the actual $g$ band data. since it is dominated by the accretion disk continuum. The MIR SF analysis is performed in the most sensitive $W1$ band (3.4\,\textit{\rm\mu m}) that traces the hottest dust associated with the inner edge of the torus, as well as the $W2$ band with a slightly longer wavelength (4.6\,\textit{\rm\mu m}). We will mainly explicate our methods using the $W1$ data, since its light curve is of much higher quality than that in the $W2$ band. While we also present the $W2$ results in the main text, the details of measurements in the $W2$ band are presented in Appendix A. We will directly compare the torus geometry probed by the SFs to that measured from dust reverberation mapping to test the feasibility of our method.

3 Ensemble Structure Function

The structure function measures the rms variability of quasars at different timescales, which is a powerful model-independent tool to empirically quantify AGN variability. Following earlier works (e.g., Kozłowski, 2016), we measure the SF as

[TABLE]

where $\Delta f$ is the magnitude or flux difference between two epochs separated by $\Delta t$ (in the rest-frame) after subtracting the photometric noise $\sigma_{f}$ as

[TABLE]

The error bars on $\rm SF(\Delta t)$ are calculated as the standard error on the mean (i.e., std/$\sqrt{N-1}$) for all the $\Delta f$ pairs within a given $\Delta t$ bin. In addition, we define the observed noise term ($N_{\rm err}$) as the mean of all $\sigma_{f}(t)^{2}+\sigma_{f}(t+\Delta t)^{2}$ then take the square root.

Theoretically, the SF can be measured for light curves at arbitrary observing segments assuming that AGN variability is a stationary process, and should be insusceptible to sampling problems. However, in reality, the AGN variability may be non-stationary, or even if it is stationary, the baseline of the survey is not long enough to well constrain the flattening part of the SF at long timescales (e.g., Kozłowski, 2017). The sparse and/or irregular time sampling will also cause spurious features in the SFs owning to the lack of sufficient pairs to robustly measure an average $\Delta f$ (e.g., Emmanoulopoulos et al., 2010). These problems can be partially alleviated by measuring the ensemble SF which is the average of the SFs of many quasars with similar properties (e.g., Vanden Berk et al., 2004).

It is important to emphasize that the accuracy of noise subtraction can significantly affect the ensemble SF measurement. While Equation 4 provides a straightforward approach for subtracting the noise (hereafter the $N_{\rm err}$ method), its accuracy strongly depends on the accuracy of the reported uncertainties of photometric measurements. Inaccurate estimation of magnitude/flux uncertainties could introduce biases in the shape of the SF. In fact, our analysis reveals that the $N_{\rm err}$ method may result in systematic overestimation of the noise term for the $W2$ band, which has lower light curve quality relative to the $W1$ band (see Appendix A for details). Alternatively, a functional fitting method can be used to estimate the noise term. As we will demonstrate in Section 4.1, Equation 2 can adequately describe the MIR SF at timescales of $100\lesssim\Delta t\lesssim 2000$ days. By adding a constant noise term ($C_{\rm noise}$) in quadrature to Equation 2 and fitting the observed ensemble SF (without subtracting the photometric noise beforehand), which flattens at short timescales due to flux errors (i.e., the SF floor; Kozłowski 2016), the noise term and its associated uncertainty can be determined (hereafter the $C_{\rm noise}$ method).

We note that the $C_{\rm noise}$ method is applicable only if the noise term is nearly constant with time, which is the case for the WISE light curves. In contrast, the $g$ band light curves are obtained from multiple surveys conducted during various observing seasons. They contribute differently to the SF at different timescales, leading to a notable time variation in the noise term. Therefore, we always use Equation 4 to subtract the noise in the $g$ band333Our results are insensitive to the shape of the $g$ band SF at $\Delta t\lesssim 10$ days where noise subtraction impacts the SF the most. . The noise terms estimated from the $N_{\rm err}$ and $C_{\rm noise}$ methods show excellent agreement in the $W1$ band. For the sake of simplicity, we present our methods and results using the $W1$ SFs measured from the $N_{\rm err}$ method, which is straightforward to compute and consistent with our treatment for the $g$ band. However, we report the $W2$ results based on the $C_{\rm noise}$ method as default since the $N_{\rm err}$ method likely overestimates the photometric noise. By subtracting the constant SF floor with the $C_{\rm noise}$ method, we remove all contributions from flux uncertainties in the SF measurements. We refer the readers to Appendix A for a detailed comparison of the results obtained from different noise estimation methods and the discussion of the impact of imperfect noise subtraction.

To measure the ensemble SF, we split our quasar sample into two redshift bins, $z<0.77$ and $z>0.77$, based on the median redshift, and further divide each subsample by their median bolometric luminosity ($L_{\rm bol}$). This sample division provides a higher-$L$ subsample and a lower-$L$ subsample in a given redshift bin, where the WISE bands trace the same rest-frame infrared wavelength. As a result, we have four subsamples each containing $\sim 150$ quasars. The median redshift and luminosity for each samples are $z\approx 0.50$, 0.64, 0.93, 1.20 and \textit{{\rm log},L_{\rm bol}}\approx 45.3, 45.7, 45.9 and 46.4 $\rm erg\ s^{-1}$, respectively. These four subsamples will be referred to as S1, S2, S3, and S4, respectively.

In the top panels of Figure 1, we present the ensemble SFs in the $W1$ band measured in magnitude units, as commonly adopted in the literature, to qualitatively show their differences in the optical and MIR bands. The uncertainties on the ensemble SFs are calculated by bootstrapping each individual SF and deriving the semiamplitude of the enclosed 16th and 84th percentiles on the bootstrapped means. The $W1$ SFs increase with $\Delta t$ in a way that resembles the shape of the optical SFs. Additionally, the variability amplitudes decrease with increasing $L_{\rm bol}$ at fixed redshift in both bands, consistent with the scenario that the MIR variability is caused by reprocessing of UV/optical emission. However, it can be clearly seen that the $W1$ SFs are steeper than that in the $g$ band and have smaller amplitudes at all timescales being probed.

A functional fit to the $W1$ SFs using Equation 2 yields a short-timescale slope of $\beta\approx 1.3-2.0$ (Table 1), which is much steeper than that reported in Kozłowski et al. (2016) based on Spitzer light curves at 3.6\,\textit{\rm\mu m} and 4.5\,\textit{\rm\mu m} of $1<z<3$ AGNs ($\gamma\approx 0.45$ where $\beta\sim 2\gamma$). Moreover, the slope appears to increase with $L_{\rm bol}$, while no correlation was found in Kozłowski et al. (2016). This could be due to the insufficient data quality in Kozłowski et al. (2016) since their Spitzer light curves only cover five epochs with large gaps (several years) exist between adjacent observing windows. Analogous findings are observed in the $W2$ band, with the slope and variability amplitude being slightly steeper and smaller, respectively, when compared to the $W1$ band (Appendix A).

The accretion disk (AD) could also contribute to the observed emission at rest-frame \sim 1-2\ \textit{\rm\mu m} and contaminate the torus variability measurements (e.g., Kishimoto et al., 2007, 2008; Lira et al., 2011). This is illustrated in the inset of Figure 1 where we show that the $W1$ to $g$ band flux ratio ($f_{w}/f_{g}$) of our sample (blue line) decreases significantly with $L_{\rm bol}$. Assuming the SED of the AD follows the $F_{\nu}\propto\nu^{1/3}$ power law that extends to \lambda_{\rm rest}\sim 2\,\textit{\rm\mu m} (Kishimoto et al., 2008), we estimate the AD flux in the $W1$ band ($f_{w,\rm AD}$) from the $g$ band flux and show its relative strength to the pure torus flux ($f_{w,\rm torus}$) in the inset of Figure 1 (red line). The non-negligible value444The contamination from the host galaxy has not been taken into account. The actual $f_{w,\rm AD}/f_{w,\rm torus}$ for the brightest quasars at $z>1$ could be even higher as the host galaxies have a stronger contribution at \textit{\lambda_{\rm rest}}\sim 1.5\,\textit{\rm\mu m} compared to \textit{\lambda_{\rm rest}}\sim 2000\,\rm\AA. at the highest luminosities suggests that the AD could contribute significantly to the short-timescale variability in the $W1$ band. Conversely, given the rising shape of the torus SED and the declining trend of the AD SED in the MIR, the AD contamination in the $W2$ band is insignificant.

Considering the sparse time sampling of the light curves and the absence of $g$ band observations after 2018, we adopt a statistical approach to remove the AD contamination, instead of subtracting the AD flux from each IR light curves as done in earlier works (e.g., Koshida et al., 2014; Lyu et al., 2019). The idea is to measure the SF in flux density unit (weighted by luminosity distance2 for each object in the bin to account for the redshift effect; see bottom panels in Figure 1) and infer the SF of the disk continuum in the WISE bands (${\rm SF}_{w,\rm AD}$) from that in the $g$ band (${\rm SF}_{g,\rm AD}$), then subtract it from ${\rm SF}_{w}$ to obtain the pure torus component (${\rm SF}_{w,\rm torus}$). Other sources of contamination to the WISE flux, such as a constant host galaxy or nearby objects will also be subtracted out when measuring $\Delta f$ in flux.

The details and feasibility of our subtraction methodology will be further illustrated in Section 4.2 with simulated light curves and SFs. Here we briefly outline our method. Given the weak wavelength dependence of AD variability (e.g., Stone et al., 2022), ${\rm SF}_{w,\rm AD}$ can be estimated as \textit{{\rm SF}_{g,\rm AD}}\times(\nu_{w}/\nu_{g})^{1/3}. If we roughly assume that the flux variations at different timescales follow a Gaussian distribution and associate its standard deviation with the rms variability, then ${\rm SF}_{w,\rm torus}$ can be derived by solving

[TABLE]

where $\rho$ is the covariance (correlation coefficient) between $f_{w,\rm AD}$ and $f_{w,\rm torus}$ as a function of timescale. In practice, this assumption does not hold since AGN variability is red-noise like instead of a Gaussian white noise. However, we will demonstrate in Section 4.2 that it provides a reasonable approximation. The value of $\rho$ depends on the geometry of the torus: an extremely compact (extended) torus with $R\rightarrow 0$ ($R\rightarrow\infty$) will have $\rho\rightarrow 1$ ($\rho\rightarrow 0$). Observationally, it is challenging to accurately constrain $\rho$ with our sparse light curves since the optical and WISE data do not overlap well. Therefore, we adopt $\rho$ measured from simulated light curves which is typically $-0.3\sim 0$ for a broad range of torus geometry at timescales of $\sim 100-2000$ days (see Section 4.2 for details).

The resulting pure-torus SFs in the $W1$ band are shown in Figure 1 (bottom panels). This time the amplitudes of the MIR SFs are greater than those in the optical as the absolute flux in the $W1$ band is typically a factor of few times larger than that in the $g$ band. The AD contamination is small overall, except for the most luminous S4 sample, where its variability in the $W1$ band is comparable to the torus at $\Delta t\sim 100$ days. The functional fitting results with Equation 2 to the $g$ band ($\beta$ fixed to 1.0, i.e., assuming a DRW process) and pure-torus SFs measured in flux unit are summarized in Table 1. The $g$ band SFs can be well described by a DRW model with $\tau\sim 700$ days, consistent with previous dedicated studies on the AGN optical variability (e.g., Stone et al., 2022). The $W1$ slope becomes slightly steeper compared to that measured in magnitude units because of the subtraction of the AD contamination, while its correlation with luminosity remains.

The steeper MIR SF demonstrates that the variability in the MIR is suppressed compared to the optical variability, and the suppression is most prominent at short timescales. This could be attributed to the smoothing effect, where the optical variations of the AD have been smeared out by the extended structure of the torus, due to that different parts of the torus respond asynchronously to the incident UV/optical continuum. The optical and MIR SFs thus contain critical information on the geometry of the torus (e.g., inclination, size, opening angle). For example, consider a ring-like torus with a face-on inclination. In this configuration, the observed MIR light curve would simply be a time-delayed version of the optical light curve, and the shapes of their SFs should be identical. Therefore, the torus probed by the WISE bands must have some vertical or radial extension and/or the inclination angle is non-zero. Higher-luminosity quasars may have a more extended dust distribution in order to produce a steeper $\beta$. In the next section, we will generate simulated light curves using a toy torus model to test this hypothesis.

4 Constraining Torus Geometry

4.1 Idealized Light Curves

We first generate idealized light curves (regularly sampled with a 5-day cadence and a baseline of $\sim 10^{4}$ days in both bands) to illustrate the shape of the MIR SF under different torus geometries. The simulated light curves and the corresponding SFs in the following sections are presented in flux units (arbitrary normalization). The model parameters are summarized in Table 2. Below we describe the models and parameters in detail.

The optical light curve is generated from a DRW model using celerite (Foreman-Mackey et al., 2017) assuming a mean flux of 1.0, a $\sim 10\%$ fractional variability amplitude (\textit{\rm SF_{\infty}}=0.1) and a decorrelation timescale $\tau=700$ days, as motivated by observations (Table 1). The MIR light curve is the convolution of the optical light curve with the torus transfer function $\phi(\tau)$ as

[TABLE]

where $\phi(\tau)$ quantifies the delayed response (at time lag $\tau$) of distant dusty clouds to the continuum emissions from the accretion disk, and $f_{\rm amp}$ is the ratio between the variable part of the optical and IR emissions which is related with the dust amount and reprocessing efficiency. The transfer function depends on the geometry, composition (grain species), and physics (e.g., anisotropic illumination, optical thickness, cloud occultation) of the torus (e.g., Kawaguchi & Mori, 2011; Almeyda et al., 2020), while a simplified top-hat transfer function is often employed in dust lag measurements (e.g., Lyu et al., 2019; Yang et al., 2020).

In this work, we consider a simple geometric model and neglect detailed radiative transfer given that the light curve quality is insufficient to constrain the complicated physics that could be highly degenerate. Our torus model contains five parameters: the inclination angle ($\theta=0^{\circ}$ means a face-on view), the inner radius $R_{\rm in}$ (in units of pc), the outer to inner radius ratio $Y$, the torus half-opening angle $\sigma$, and the radial power-law index $p$. Specifically, since we are measuring the ensemble SF in a narrow wavelength range (\textit{\lambda_{\rm rest}}\sim 1.4-2.8\,\textit{\rm\mu m} in $W1$ and \textit{\lambda_{\rm rest}}\sim 2.0-3.8\,\textit{\rm\mu m} in $W2$), the clouds are expected to be the hottest dust concentrate near the sublimation radius ($T\sim 1500$ K, $R\lesssim 1$ pc), as opposed to the polar dust component ($T\sim 100-200$ K, $R\sim 100-1000$ pc) that dominates the emissions at \lambda\sim 20\,\textit{\rm\mu m} and contribute little to the $W1$ and $W2$ bands (e.g., Hönig et al., 2013; Lyu & Rieke, 2018; Yamada et al., 2023). Therefore, we assume that the torus has a thin-wall geometry by fixing $Y$ to a small value of 1.5555Assuming $T_{\rm in}=1500$ K (i.e., the sublimation temperature for Silicate) and a $T\propto R^{-1}$ radial temperature profile (Kishimoto et al., 2011), Y\sim T_{\rm in}/T_{2.8\,\textit{\rm\mu m}}\sim 1500 K/1000 K according to Wien’s displacement law., and neglect the polar component in our model. In the radial direction, the clouds are distributed following a power-law profile with its index being fixed to –1.0666We find that as long as $Y$ is small (e.g., $\lesssim 2.0$), the choice of $p$ has a negligible impact on the shape of the transfer function. (e.g., Kishimoto et al., 2011; Lyu et al., 2019). In the vertical direction, the clouds have a sharp edge at [$-\sigma$, $\sigma$] and within $\sigma$ they follow a Gaussian distribution with the standard deviation being set to $\sigma$. As a result, most clouds are clustered on the equatorial plane. In addition, we define the characteristic (response-weighted) time delay, comparable to the time lag measured from dust reverberation mapping, as

[TABLE]

and convert it to the effective size of the torus ($R_{\rm eff}$) by multiplying the speed of light. A schematic view of our torus model is shown in Figure 2.

We trace the time delay of each cloud (Figure 2) relative to the observer’s line of sight and measure the transfer function for different torus geometries. In Figure 3 we show four examples of the transfer functions and the resulting MIR light curves assuming the average MIR to optical flux ratio ($f_{w,\rm torus}/f_{g,\rm AD}$) to be 1.0. For a compact torus with \textit{R_{\rm in}}=0.2 pc, the MIR light curve follows the general variability pattern in the optical with a clear time delay; while for an extended torus with \textit{R_{\rm in}}=0.9 pc, the variability in the MIR is significantly diluted. Increasing $\sigma$ and $\theta$ can further smear out the variability features.

In Figure 4 we plot the SFs of the simulated light curves which are derived from 100 ensembles to mitigate random fluctuations. In the first panel we demonstrate how the SF changes with $R_{\rm in}$ under fixed $\theta$ ($30^{\circ}$) and $\sigma$ ($30^{\circ}$). At \textit{R_{\rm in}}=0.01 pc, the MIR SF follows that in the optical at long timescales. However, its slope becomes steeper and resembles a power-law shape below a timescale of $\sim 100$ days, meaning that the rapid flux variations are suppressed. As $R_{\rm in}$ increases, the dilution effect further weakens the short-timescale variability and starts to impact the longer-timescale regime. When focusing on $100\lesssim\Delta t\lesssim 2000$ days (i.e., the timescale probed by observations), the slope $\beta$ becomes steeper as $R_{\rm in}$ increases. This trend resembles the correlation between $\beta$ and $L_{\rm bol}$ we found in Section 3, indicating that the size of the torus may increase with $L_{\rm bol}$. In the second and third panels we fix $R_{\rm in}$ and vary $\theta$ and $\sigma$. Increasing $\theta$ and $\sigma$ has a similar effect on the shape of the SF as increasing $R_{\rm in}$, meaning that there are inherent degeneracies among different torus parameters.

It is worth noting that the MIR SFs do not strictly follow Equation 2. This empirical functional form fails to decribe the shape of the model SF at the shortest timescales, especially for compact torii where the suppression of short timescale variability is stronger than that predicted by Equation 2. This is illustrated in the bottom panels of Figure 4 where we show that the normalized residuals of fitting the model SFs with Equation 2 can be as large as $\sim 50\%$ at $\Delta t\lesssim 100$ days. However, it can still provide a good fit to the data at $\Delta t\gtrsim 100$ days (i.e., the timescales probed by our WISE light curves), thus offering a convenient means to quantify the observed long-timescale variability.

Overall, the shape of the simulated SFs are in good agreement with those of SDSS quasars, lending support that the observed MIR variability is originated from the delayed and smoothed response to the UV/optical variations. Furthermore, the MIR SFs of SDSS quasars at $\Delta t\sim 100-600$ days follow a steep power-law shape, distinct from the model SFs with \textit{R_{\rm in}}<0.1 pc. This result disfavors a compact torus geometry for these quasars.

4.2 Effect of Accretion Disk Contamination

As noted in Section 3, the variability of disk continuum in the MIR can affect the SF measurement for the most luminous quasars. To demonstrate this effect, we generate a simulated $W1$ light curve with \textit{R_{\rm in}}=0.6, $\theta=30^{\circ}$, $\sigma=30^{\circ}$, $f_{w,\rm torus}/f_{g,\rm AD}=3.4$, and add disk contamination represented by $f_{w,\rm AD}=f_{g,\rm AD}\times(\nu_{w}/\nu_{g})^{1/3}$. The parameters are selected to be representative of luminous quasars (\textit{{\rm log},L_{\rm bol}}\sim 46.4) given the f_{w,\rm torus}/f_{g,\rm AD}-\textit{L_{\rm bol}} relation shown in Figure 1 and the $R\propto L^{0.5}$ relation found in the literature (e.g., Koshida et al., 2014; Lyu et al., 2019). The resulting SFs with and without AD contamination are displayed in Figure 5. A noticeable leveling off towards $\Delta t\lesssim 300$ days is evident when the disk contamination is included in the MIR SF measurement; there is also an excess relative to the true SF at $\Delta t\gtrsim 1000$ days. The flattening at short timescales can also be seen in the observed SF, demonstrating that its shape is a superposition of the torus and AD component. Nevertheless, the variability signals at the timescales probed by our light curves are dominated by the torus emission.

To remove the disk contamination using the method outlined in Section 3, we first measure $\rho$ as a function of $\Delta t$ from our idealized light curves (shown as a black curve in the inset of Figure 5). The pure-torus component, ${\rm SF}_{w,\rm torus}$, is then calculated using Equation 5. The AD-subtracted SF is in good agreement with the true model (Figure 5), demonstrating the effectiveness of this method. The slight overestimate ($\sim 5\%$) of the true SF at the longest timescales may be attributed to the significant deviation of AGN variability from a Gaussian process as assumed in Section 3.

When applying this method to observed quasars, one can roughly estimate $R_{\rm in}$ from the $R\propto L^{0.5}$ relation and derive a dedicated $\rho-\Delta t$ relation. Figure 5 shows examples with \textit{R_{\rm in}}=[0.01\sim 1.0,dR=0.1] in the inset. It can be seen that the $\rho-\Delta t$ relation is insensitive to $R_{\rm in}$ (as well as $\theta$ and $\sigma$, which are not shown) at $\Delta t\sim 100-1000$ days for \textit{R_{\rm in}}\gtrsim 0.4. This critical size corresponds to \textit{{\rm log},L_{\rm bol}}\gtrsim 46.0, consistent with our high luminosity samples. On the other hand, a compact torus with \textit{R_{\rm in}}\lesssim 0.4 is expected to be associated with low luminosity AGNs, and thus, subtracting the AD contamination or not will have a minimal impact on the SF.

4.3 Realistic Light Curves

The idealized simulations presented in Section 4.1 have demonstrated that the geometry of the torus is imprinted in the shape of the MIR SF, and its difference compared with the optical SF. In principle, this could be used to constrain torus parameters by fitting the measured ensemble SFs to models. However, in practice, the sparse time sampling of the light curve, photometric noise, and limited number of ensembles make it impossible to measure the SF with the precision like the models shown in Figure 4. Additionally, the observed ensemble SF is the average of individual SFs associated with a broad range of torus parameters. It is unclear whether the average SF also corresponds to the average of the underlying torus parameter distributions. Here we simulate realistic light curves and test to what extend can the current data be used to constrain the geometry of the torus.

We first summarize the observational constraints on the distribution of torus parameters that will be used to generate our simulated populations. The effective size of the torus at \textit{\lambda_{\rm rest}}\sim 2\,\textit{\rm\mu m} is $\sim 0.1-1.0$ pc for luminous quasars. At a given luminosity, it follows a log-normal distribution with an intrinsic scatter of only $\sim 0.1$ dex (Koshida et al., 2014; Yang et al., 2020). The opening angle $\sigma$ can be derived from f_{\rm IR}\equiv L_{\rm IR}/\textit{L_{\rm bol}} once corrected for the anisotropic effect as detailed in Stalevski et al. (2016), and is typically $30^{\circ}-50^{\circ}$ for Type 1 quasars (e.g., Ezhikode et al., 2017; Ichikawa et al., 2019). The inclination angle should be randomly distributed within [$0^{\circ}$, $90^{\circ}-\sigma$] according to the unification model, and its average value is expected to be $\sim 25^{\circ}-45^{\circ}$. This simple argument is consistent with inclinations constraint from SED fitting, kinematic modeling and polarization observations (e.g., Fischer et al., 2013; Marin, 2014; Zhuang et al., 2018).

Based on these results, we consider a hypothetical population of Type 1 quasars that varies as a DRW process in the optical band. The decorrelation timescale $\tau$ (in rest-frame days) is assumed to follow a Gaussian distribution $\mathcal{N}(600\sim 700,\ 50)$, and \textit{\rm SF_{\infty}}\sim\mathcal{N}(0.20\sim 0.25,\ 0.02). For a mean optical flux of 1.0 (which will be used in our simulations), this $\rm SF_{\infty}$ corresponds to a fractional variability amplitude of $\sim 16\%$, similar to the observed value for our sample. The torus opening angle, inclination angle, and logarithm of $R_{\rm in}$ are modeled as Gaussian-distributed with $\sigma\sim\mathcal{N}(30\sim 50,\ 10)$, $\theta\sim\mathcal{N}(25\sim 45,\ 10)$, and {\rm log}\textit{R_{\rm in}}\sim\mathcal{N}(-1.0\sim-0.1,\ 0.1), respectively.

In our simulation, we first generate $\overline{\tau}$, \overline{\textit{\rm SF_{\infty}}}, $\overline{\sigma}$, $\overline{\theta}$, and \overline{{\rm log}\,\textit{R_{\rm in}}} from a uniform distribution $U(600,700)$, $U(0.20,0.25)$, $U(30,50)$, $U(25,45)$, and $U(-1.0,-0.1)$, respectively, which we treat as the average properties of the hypothetical quasar populations. Then for each set of average parameters, we create a simulated quasar sample and assign each quasar an optical variability pattern and torus geometry by randomly sampling the aforementioned Gaussian distributions for 150 times. The sample size is chosen to match our S1–S4 samples. Next, for each simulated quasar, we generate light curves with realistic photometric uncertainties (0.02 for $g$ band and 0.04 for $W1$ band) and cadence (by using the observing epochs of randomly selected SDSS quasars). An “observed” ensemble SF ($\rm SF_{obs}$) similar to those shown in Figure 1 is then derived for the 150 simulated light curves in both bands.

The next step is to generate a model SF library in the MIR. We first fit the $g$ band $\rm SF_{obs}$ with Equation 2 and obtain the decorrelation timescale $\tau_{\rm fit}$, which is typically $500-800$ days. The MIR model SFs used in the fitting are then generated from a grid of parameters spanning $\tau=$ the nearest integral multiple of 100 to $\tau_{\rm fit}$ (in the optical), \textit{\rm SF_{\infty}}=0.1 (optical), $\theta=[25\sim 45,\ d\theta=1]$, \textit{R_{\rm in}}=[0.01\sim 1.0,\ dR=0.01], and $\sigma=[30\sim 50,\ d\sigma=2]$ following the method described in Section 4.1. A parameter $f_{\rm norm}$ is introduced in the fitting to account for the relative normalizations. By searching for model SF (i.e., example curves displayed in the top panels of Figure 4) that best matches the measured $\rm SF_{obs}$  (the actual SF data points, not the analytic approximation obtained using Equation 2) in the MIR through MCMC sampling with emcee (Foreman-Mackey et al., 2013), we can test the recovery of the input torus parameters (i.e., $\overline{\sigma}$, $\overline{\theta}$, and \overline{{\rm log}\,\textit{R_{\rm in}}}) by comparing them with that of the best-fit model SF.

Figure 6 shows the comparison between the input and fitted torus parameters. The best-fit parameter and its $1\,\sigma$ uncertainty are determined from the 16-50-84th percentiles of the posterior distribution (assuming a flat prior). The inner radius of the torus can be robustly constrained, given that it has the strongest impact on shaping the SF. However, the recovery of $R_{\rm in}$ becomes less accurate for extended torus. This may be due to a reduction in its variability amplitude at short timescales, which makes it challenging to separate the variability signal from photometric noise. The posterior distribution of the inclination angle and opening angle spans the entire allowed fitting range, and the fitted values exhibit no correlation with their input counterparts. This failure is driven by the combined effect of the stochastic nature of SF measurement, imperfect noise subtraction, and inherent parameter degeneracies. Nevertheless, the effective size of the torus is insensitive to the poorly constrained $\theta$ and $\sigma$, meaning that although different torus parameters are degenerate, the average response-weighted size can be well determined.

It is noteworthy that the choice of the radius ratio $Y$ also has an impact on the derived $R_{\rm in}$ and $R_{\rm eff}\equiv c\tau_{\rm lag}$. If the dust emitting over a narrow range of wavelengths has a more compact (e.g., $Y=1.1$) or extended (e.g., $Y=2.0$) spatial distribution, whereas $Y$ is fixed at 1.5 in our fitting, we find that the size tends to be slightly underestimated for the former and overestimated for the latter by a factor of $\sim 1.15$. Nonetheless, this small systematic difference has no significant impact on our main conclusions.

4.4 Torus Geometry for SDSS quasars

We now apply the fitting methodology to the observed ensemble SFs of SDSS quasars to constrain their torus geometry. In the fitting, $R_{\rm in}$, $\sigma$ and $\theta$ are allowed to vary within the ranges of [0.01, 1.5], [30, 50], and [25, 45], as motivated by the observational constraints on the average torus parameters in Section 4.3. Note that slightly changing (e.g., by 10∘) the allowed fitting range for $\sigma$ and $\theta$ has limited impact on the size measurements. The best-fit SFs and torus parameters in the $W1$ band are shown in Figure 1. Our simplified torus model provides a successful fit to the data, which further support the notion that the MIR variability of quasars is driven by UV/optical variability and diluted by the extended torus geometry. Given that the inclination and opening angle cannot be constrained with the current dataset (Section 4.3), we focus on the fitting result of the torus size in the following discussion.

In Figure 7 (red points), we present the correlation between $R_{\rm eff}$ and $L_{\rm bol}$ for the S1–S4 sample measured from our ensemble SF analysis in the $W1$ band. To more effectively sample the $R-L$ relation, we further divide our quasars into nine luminosity-based subsamples (0–20th, 10–30th, …, and 80–100th percentile of the luminosity distribution), and show them as black points. The overlapping luminosity range is chosen to guarantee sufficient numbers of objects are included in each bin to facilitate a robust measurement of the ensemble average. For comparison, previous results from dust reverberation mapping studies are also shown in Figure 7 (Koshida et al., 2014; Lyu et al., 2019; Yang et al., 2020). In particular, Lyu et al. (2019) measured time lags in the WISE $W1$ band for 67 PG quasars (44.6\lesssim\textit{{\rm log},L_{\rm bol}}\lesssim 46.6) at $\overline{z}\sim 0.15$, corresponding to \textit{\lambda_{\rm rest}}\sim 2.9\,\textit{\rm\mu m}. The sample of Koshida et al. (2014) contains local low-luminosity AGNs (42.9\lesssim\textit{{\rm log},L_{\rm bol}}\lesssim 45.6) with time lag measured in the $K$ band (\lambda_{\rm rest}\sim 2.1\textit{\rm\mu m}), which is similar to our less luminous quasars (\textit{{\rm log},L_{\rm bol}}<46.0,\ \lambda_{\rm rest}\sim 1.8-2.3\,\textit{\rm\mu m}) but slightly redder than the most luminous ones (\textit{\lambda_{\rm rest}}\sim 1.5\,\textit{\rm\mu m}, corresponding to the $H$ band). However, simultaneous multiband monitoring of local AGNs has revealed that the time delays between the optical (e.g., $V$ band) and $JHK$ bands are similar (e.g., $\tau_{V,H}/\tau_{V,K}$ has a mean of 0.9 and a scatter of 0.13 for our compilation of $\sim 20$ local AGNs (or the same AGN in different observing runs) from Lira et al. 2015; Pozo Nuñez et al. 2015; Oknyansky et al. 2015; Mandal et al. 2018; Landt et al. 2019; Sobrino Figaredo et al. 2020; Mandal et al. 2021), suggesting similar torus sizes in the rest-frame NIR bands. Therefore, our result in the $W1$ band can be directly compared to that of Koshida et al. (2014).

Intriguingly, the torus sizes derived from our ensemble SF method exhibit a strong correlation with luminosity. The best-fit relation for the nine black points (neglecting correlated errors) is {\rm log}\,\textit{R_{\rm eff}}\,(\rm pc)=0.51_{-0.04}^{+0.04}\times{\rm log}\,(\textit{L_{\rm bol}}/10^{46}\,\textit{\rm erg\ s^{-1}})-0.38_{-0.01}^{+0.01}. This trend is not a result of AD subtraction which steepens the SF for the most luminous quasars: a direct fit of ${\rm SF}_{w}$ at $\Delta t>200$ days, where disk variability has a lesser impact, still yields consistent result (slope $\sim 0.43_{-0.04}^{+0.04}$). We also apply the same fitting methodology to the $W2$ SFs (see Appendix A for details) to measure the torus sizes at a slightly longer wavelength. As shown in Figure 7 (magenta points), the size in the $W2$ band is also strongly correlated with AGN luminosity, which can be well described by the $R-L$ relation in the $W1$ band but with an elevated normalization. This suggests that the $W2$ band traces dust located further away from the central SMBH. The size ratio between the two bands ($\sim 1.4$) derived from our SF method is consistent with that measured from dust reverberation mapping for the same sample ($\sim 1.3$; Yang et al. in prep.). The slope of $\sim 0.5$ is consistent with the results in the low luminosity regime (e.g., Koshida et al., 2014; Lyu et al., 2019) and theoretical expectation of the inner edge of the torus being determined by dust sublimation and radiation equilibrium (Barvainis, 1987; Kishimoto et al., 2007). Our result supports that the circumnuclear torus structures in luminous quasars and their Seyfert counterparts are shaped by the same physical process.

However, as luminosity increases, we observe a systematic difference (increase) in sizes in comparison to the torus lags measured in Y20. As discussed in Y20, the $R-L$ relation in their study may be biased low towards the high-luminosity end. This could be due to the difficulty of reliably constraining the increased time delay within a limited baseline, as well as the contamination from the accretion disk which is not subtracted in Y20. While our SF method does not suffer from this limitation, there is a substantial scatter between the input and fitted $R_{\rm eff}$ for such an extended torus (Figure 6). Therefore, we cannot rule out the possibility that our size measurements at the highest luminosities may be affected by small number statistics. To firmly measure the $R-L$ relation via this SF method, a larger sample is needed to provide more independent ensemble SF measurements at the luminous end (e.g., a sample of $\sim 1200$ quasars can double the number of independent ensemble SF measurements).

5 Concluding Remarks

We have measured the optical and MIR ensemble SFs in the $g$ band and WISE $W1$ and $W2$ bands for 587 quasars at $0.2\lesssim z\lesssim 2$ in the SDSS Stripe 82 region. The SFs in the MIR exhibits smaller variability amplitudes and steeper slopes than those in the optical. In addition, the MIR variability amplitude shows a negative correlation with luminosity while the slope shows a positive trend. We built a geometric torus model and demonstrated that the shape of the MIR SF and its correlation with luminosity can be well explained by the geometric dilution effect, in which the MIR variability of quasars is originated from the delayed and smoothed response of the torus to the incident UV/optical variations. We further generated simulated light curves with realistic photometric noise and cadence and demonstrated that the ensemble SFs in the MIR can be used to constrain the size of the torus. The derived size of SDSS quasars shows a remarkable $R\propto L^{0.5}$ relation that is consistent with dust reverberation mapping measurements and theoretical expectations.

Our SF method provides an important independent means to constrain the torus $R-L$ relation, complementary to the dust reverberation mapping technique. It can be applied to light curves with any observation period, regardless of whether the optical data lead or overlap with the MIR data or not, and is less susceptible to the sampling issue, as long as the ensemble AGN variability is roughly stationary. The only requirement is to have sufficient flux pairs and precise photometric measurements to robustly measure the variability amplitude at a given timescale, which could be achieved by denser sampling, a longer baseline, or larger ensemble samples. In the next decade, the synergy between the $\sim 10$ yr Vera C. Rubin Observatory and the $\sim 5$ yr Nancy Grace Roman Space Telescope will further extend the baseline and provide optical and NIR light curves for AGNs up to $\sim 5$ mag fainter with much higher photometric precision than the current ground-based facilities and WISE. This will not only allow more robust SF measurements for large samples of luminous quasars, but also extend the current analysis to much lower luminosities, for which the variability is intrinsically stronger (e.g., Vanden Berk et al., 2004) and the SFs are expected to be distinctly different (even at $\Delta t>100$ days) from the luminous populations given their compact torus sizes.

Finally, while the SF method is straightforward to implement and much less demanding of light curve data quality than is dust reverberation mapping, we emphasize that more geometric information describing the torus can be extracted from the full shape of the transfer function. For example, as illustrated in Figure 4, although the SF is insensitive to the torus orientation given that it only measures the rms variability, the actual orientation can be constrained from the shape of the transfer function. Therefore, dust reverberation mapping using high-quality optical+IR light curves to measure the shape of the transfer function would be far more powerful to infer the full geometry of the torus. In future work, we plan to investigate the recoverability of the transfer function using optical and IR light curves, and the corresponding torus geometry with the model presented here.

We publicly release the code used to generate the torus models, transfer functions and SFs at this github link: https://github.com/bwv1194/geometric-torus-variability.

We thank the anonymous referee for a careful reading of our manuscript and for constructive comments which helped improve the quality of the paper. This work is supported by NASA grant 80NSSC21K1566. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/.

Appendix A Structure Function in the $W2$ band and the impact of imperfect noise subtraction

This appendix presents the details of SF measurements in the $W2$ band compared with the $W1$ results shown in the main text. Additionally, we discuss the impact of imperfect noise subtraction on the SF measurements and the derived torus sizes. Figure 8 displays the ensemble SFs (in magnitude units) in the $W2$ band measured from the $N_{\rm err}$ method (i.e., noise term estimated from the reported photometric uncertainties) and compare them with those in the $W1$ band. As expected, the $W2$ SFs are steeper and have smaller variability amplitude, given that the $W2$ band traces a colder and more extended torus component. Additionally, the $W2$ slope exhibits a similar positive correlation with AGN luminosity. However, the $W2$ SF appears to be excessively steep, even becoming negative at $\Delta t\lesssim 300$ days. Applying the fitting methodology in Section 4.4 to the $W2$ SFs (measured in flux units and excluding negative points) results in a strong $R-L$ relation that is similar to the $W1$ band (blue points in Figure 9). But the derived $W2$ size is $\sim 2$ times larger than the $W1$ size, which significantly exceeds that measured from dust reverberation mapping for the same sample ($\sim 1.3-1.4$; Yang et al. in prep.).

The steep/negative SF in the $W2$ band is likely caused by its elevated flux uncertainty ($\sim 1.4$ times larger than the $W1$ band). Ideally, this elevated photometric uncertainty will cause an increase in both the observed rms variability and $N_{\rm err}$, while the subtraction of the latter from the former will still return an unbiased estimate of the intrinsic variability, which should be strictly greater than 0. However, our analysis reveals that the noise fraction (i.e., $N_{\rm err}$/observed rms) in the $W2$ band can be even greater than 1.0 (left panel in Figure 10), indicating that the photometric uncertainties are overestimated.

The comparison between $N_{\rm err}$ to the noise term derived from the $C_{\rm noise}$ method based on functional fitting (Section 3) in the $W2$ band is shown in Figure 10 (right panel). In order to obtain a more accurate estimation of $C_{\rm noise}$, we imposed a requirement in the fitting that the $W2$ slope $\beta$ should be steeper than the $W1$ slope, which is more accurately constrained. The fitted $C_{\rm noise}$ is $\sim 1.07$ times smaller than $N_{\rm err}$. As a result, the $W2$ SFs measured by subtracting $C_{\rm noise}$ (red points in Figure 8) become only slightly steeper than the $W1$ SFs, and negative values are largely avoided. Consequently, the $W2$ torus sizes are reduced by a factor of $\sim 1.4$ (magenta points in Figure 9). This reduction brings the $W2/W1$ torus size ratio ($\sim 1.4$) into agreement with that determined from reverberation mapping. Our analysis highlights the crucial role of accurately estimating the noise term for the measurement of SF and the determination of torus size. Even a subtle overestimation of the photometric noise (in this case, only $\sim 7\%$ overestimation) in the $W2$ band can substantially affect the inferred SFs, as the noise contribution dominates the observed rms variability at these shortest timescales (Figure 10). The empirical $C_{\rm noise}$ method properly removes the constant SF floor at short timescales due to photometric errors, and does not depend on the accuracy of photometric errors in $W2$.

In contrast, the two noise estimation methods yield consistent results in the $W1$ band (Figure 10), leading to almost identical SFs (dark and light blue points in Figure 8) and $R-L$ relations (red and black lines in Figure 9). This agreement suggests that the reported photometric uncertainties in the $W1$ band are much more accurate than those in $W2$.