Revisiting the dust torus size – luminosity relation based on a uniform reverberation mapping analysis

We investigate the torus size – luminosity relation of Type 1 AGNs based on the reverberation-mapping analysis using the light curves of the optical continuum and the IR continuum obtained with the W1 and W2-bands of the Wide-field Infrared Survey Explorer (WISE) survey. The final sample consists of 446 and 416 AGNs, respectively, for W1 and W2-band light curves, covering a large dynamic range of bolometric luminosity from 10 43 . 4 to 10 47 . 6 erg s − 1 , which show reliable lag measurements based on our quality assessment analysis. After correcting for the accretion disk contamination in the observed IR flux, we constrain the torus size ( R dust ) and AGN bolometric luminosity ( L bol ) relationship with the best-fit slope of 0.39 (0.33) for the W1-(W2) band, which is shallower than expected from the dust radiation equilibrium model. By combining the previous K-band lag measurements, we find that the measured torus size depends on the observed wavelength of the dust radiation, as R dust,K : R dust,W 1 : R dust,W 2 = 1.0:1.5:1.8 ( R dust ∝ λ 0 . 80 ) at L bol = 10 46 erg s − 1 , confirming a stratified structure of the torus, where wavelength-dependent emissions originate from distinct regions of the torus. By investigating the deviation from the best-fit torus size – luminosity relation, we find a moderate correlation between the offset from the R dust – L bol relation and Eddington ratio. This suggests a possible influence of the Eddington ratio on the observed flattening of the R dust – L bol relationship.


INTRODUCTION
Active galactic nuclei (AGNs) are amongst the most luminous objects in the Universe, powered by the accretion of matter onto a supermassive black hole (SMBH) and emit radiation across the entire electromagnetic spectrum (Rees 1984).The unification model of AGN posits the presence of a dust torus surrounding the broad emission line region (BLR).Based on different obscuration by the torus, AGNs are classified into two main categories, namely Type 1 and Type 2, characterized by the presence or absence of broad emission lines in their optical/infrared spectrum, respectively (Seyfert 1943;Urry & Padovani 1995).
In AGNs with UV luminosities ranging from 10 42 to 10 44 erg s −1 , the BLR typically resides at a distance of ∼ 0.003 to 0.02 pc from the accretion disk (AD), while the surrounding dust torus can extend from 0.01 to 0.1 pc.Consequently, the central regions of AGNs are compact and pose challenges for direct imaging techniques.Nevertheless, the advancements in near-infrared (NIR) and MIR interferometry have enabled the measurement of torus sizes in more than 40 sources, primarily focusing on nearby and bright AGNs (Kishimoto et al. 2009(Kishimoto et al. , 2011a,b;,b;Burtscher et al. 2013;Gravity Collaboration et al. 2020).
An alternative way to resolve the central part of AGNs is the technique of reverberation mapping (RM; Blandford & McKee 1982;Peterson 1993).The broad-band spectral energy distribution (SED) of an AGN presents an excess emission in the IR, which is attributed to the thermal emission from the dusty torus (Salpeter 1964;Sanders et al. 1989;Elvis et al. 1994).By measuring the delayed response (τ ) of the reprocessed IR continuum relative to the ionizing UV/optical continuum from the AD, the torus size can be obtained as R torus < c×τ , where c is the speed of light.
To measure the inner extent of the torus, the dust RM (DRM) between the optical and the NIR continuum in the K-band has been applied to ∼ 40 AGNs over a relatively narrow range of luminosity, i.e., V-band luminosity L V < 10 45 erg s −1 (Suganuma et al. 2006;Koshida et al. 2014;Mandal et al. 2018;Minezaki et al. 2019;Mandal et al. 2021).Koshida et al. (2014) found that the dust torus radius (R dust ) correlated with the optical luminosity (L) as R dust ∝ L α with the exponent α = 0.5, which is consistent with the expectation when the dust particles are in radiative equilibrium with the dust sublimation temperature (T sub = 1500K; Barvainis 1987;Kishimoto et al. 2007).However, Minezaki et al. (2019) reported a somewhat shallower slope with α = 0.424, proposing that the anisotropic illumination of the AD leads to a smaller dust-sublimation radius in the equatorial plane compared to the polar directions.They also suggested that the dust torus does not respond instantaneously to the changes in the luminosity of AD, leading to larger torus sizes being more frequently observed in less luminous AGNs, which show relatively small timescales of flux variation as the AD of low-luminosity AGNs has dimmed (Koshida et al. 2009).Consequently, the average inner radius of the dust torus tends to be larger for less luminous AGNs.In contrast, the torus radius from the K-band interferometric measurements is roughy a factor of two larger than the dust reverberation radius measured in the same band.This difference could be attributed to the distinction between the flux-weighted and response-weighted radii of the innermost dust torus obtained from interferometry and DRM, respectively (Kishimoto et al. 2011b;Koshida et al. 2014).These results indicate that the torus size and luminosity relation is complex, requiring more detailed investigations.
By expanding the K-band studies, Lyu et al. (2019) and Yang et al. (2020) performed RM studies by utilizing various optical surveys and IR imaging data from the all-sky Wide-field Infrared Survey Explorer (WISE) (Wright et al. 2010;Mainzer et al. 2014).They focused on high-luminosity AGNs with bolometric luminosities (L bol ) exceeding 10 45 erg s −1 , effectively excluding the low-luminosity counterparts.Using the WISE W1 and W2-band lag analysis, Lyu et al. (2019) derived reliable torus sizes for 67 among 87 Palomar-Green (PG) quasars at redshifts z ≤ 0.5, reporting that the torus size correlates with AGN bolometric luminosity with a slope of 0.47 (0.45) for W1-(W2) band lags.Using a large sample of 587 AGNs within the SDSS stripe 82 region, Yang et al. (2020) performed a DRM study, although they did not report an independent torus size-luminosity relationship.More recently, Chen et al. (2023) conducted DRM analysis for 78 AGNs, for which Hβ lag measurements are available, reporting a much shallower slope, α = 0.37 between the torus size and bolometric luminosity.They also uncovered a correlation between the deviations from the best-fit torus size-luminosity relation and the accretion rate (Du et al. 2014), concluding AGNs with higher mass accretion tend to have a more shortened torus size because of the self-shadowing effect of the slim disk (Kawaguchi & Mori 2011;Wang et al. 2014a).
While the thermal radiation from the dust torus dominates the IR flux of Type 1 AGNs, it includes emission from the AD.Consequently, the presence of AD contamination in the observed IR light curves leads to derived dust lags appearing shorter than the actual lag (Koshida et al. 2014;Mandal et al. 2021), because the continuum emission from the AD in the IR is nearly synchronous with the optical emission compared to the lag of the torus IR emission (Tomita et al. 2006;Kishimoto et al. 2008;Lira et al. 2011).Therefore, it is necessary to correct for AD contamination in the observed IR emission in order to reliably measure the dust lag.Notably, while Lyu et al. (2019) performed AD contamination correction in the observed IR fluxes, Chen et al. (2023) and Yang et al. (2020) did not consider this correction, potentially leading to an underestimation of their obtained torus sizes.
In a different approach, Li & Shen (2023) employed an ensemble of structure functions (SFs) based on optical and IR light curves in W1-band of 587 AGNs in the SDSS Stripe 82 region.They determined the torus size-luminosity relation with an exponent of α = 0.51 by binning the luminosity range into nine sub-samples and calculating the average effective torus size from the ensemble SF analysis for each sub-sample.Unlike the RM technique, their method relies on fitting parameters, such as the inclination angle (θ inc ), inner radius (R in ), outer-to-inner radius ratio (Y ), torus half-opening angle (θ opn ), and radial power-law index (p).
The previous studies of DRM showed a notable diversity in analysis method and the slope of the correlation between torus size and luminosity varies among these studies (See Lyu et al. 2019;Yang et al. 2020;Li & Shen 2023;Chen et al. 2023).Thus, it is necessary to further investigate the torus size -luminosity relation using a uniform and consistent analysis method over a large dynamic range of AGN luminosity.Such investigations would also contribute to establishing AGNs as standard candles for constraining cosmological parameters (Yoshii et al. 2014;Hönig et al. 2017).2) name of the filter, (3) observing time period of the light curves, (4) cadence of the light curves, (5) number of data points present in the light curves, (6) sky coverage of the corresponding survey, and (7) limiting magnitude in mag unit.
In this paper, we revisit the torus size-luminosity (R dust -L AGN ) relation with a uniform analysis of the torus lag measurements, using of a large sample of AGNs.First, our sample covers a large dynamic range in luminosities (10 43.4 to 10 47.6 erg s −1 ) with a redshift ranging from 0.02 to 2. Second, we employ a uniform method for cross correlation analysis and perform quality assessment, ensuring consistent and reliable lag measurements.Third, we investigate the effect of the AD contamination on the observed IR fluxes in W1 and W2bands, and correct for the contamination in order to determine the R dust -L AGN relationship.The presentation of the paper is as follows.In Section 2, we describe the sample and data acquired from different surveys and literature.We discuss the construction of the optical and IR light curves in Section 3. In Section 4, we describe the time series analysis.The results are presented in Section 5. We discuss the wavelength dependent lag and scatter in our R dust -L AGN relation in Section 6 and a summary in Section 7. .

SAMPLE AND DATA
In order to investigate the relationship between the torus size and AGN luminosity over a wide range of luminosities, we conducted a selection process of Type 1 AGNs with available optical and IR light curve data.We selected targets from two distinct categories based on their redshifts and luminosities.

Sample 1 (z < 0.35)
To explore the low luminosity end of the R dust -L AGN relationship, we selected targets from the catalog of Liu et al. (2019).This catalog comprises a uniformly defined sample of 14,584 Type 1 AGNs compiled from SDSS DR7 at z < 0.35 and covering a range of AGN bolometric luminosities between 10 41.5 and 10 46.6 erg s −1 .Since the AGNs in this catalog are selected based on the robustness of the broad components of the Balmer lines instead of imposing a lower limit on the FWHM of emission line, the catalog contains many low BH mass (M BH ) AGNs with lower luminosity than other existing Type 1 AGN catalogs.Thus, we can explore the low luminosity end of the R dust -L AGN relation using this catalog.The authors carefully removed the stellar continuum by fitting their SDSS spectra.We used other key parameters (i.e., optical luminosity at 5100 Å, Hβ, [O iii] and Fe ii luminosities, FWHM of Hβ, Eddington ratio etc.) from the catalog.

Sample 2 (z > 0.3)
To cover a moderate to high luminosity range in the R dust -L AGN relationship, we used the SDSS Stripe 82 (S82) region (MacLeod et al. 2012), where a total of approximately 9258 spectroscopically confirmed broadline AGNs are available.These AGNs span a wide range of redshifts and luminosities, useful for investigating luminous AGN population in the R dust -L AGN relation.Yang et al. (2020) performed a DRM study of S82 AGNs using ∼ 20 yr ground-based optical light curves and 10 yr IR light curves in the W1-band from the WISE.They constructed optical light curves using data from the SDSS, Pan-STARRS (PS1; Chambers et al. 2016), Dark Energy Survey (DES; Abbott et al. 2018), Catalina Real-Time Transient Survey (CRTS;Drake et al. 2009), Palomar Transient Factory (PTF;Law et al. 2009), Zwicky Transient Facility (ZTF;Bellm et al. 2019) and All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014).They were able to find lag between the optical and W1-bands for 587 AGNs with redshift in the range 0.3 < z < 2, and luminosity in the range 10 44.5 to 10 47.6 erg s −1 .They employed JAVELIN (Zu et al. 2011) to search for the lags between the optical and W1-band light curves.Thus, Sample 2 comprises a total of 587 AGNs from S82 region.We used the same object ID number as assigned by Yang et al. (2020) for Sample 2.

Light curves data
We gathered optical and IR data from both groundbased and space-based telescopes, as well as from the relevant literature for our target AGNs.

Optical data from ground based surveys
For Sample 1, we compiled the available optical photometric data from different ground-based transient surveys.Initially, we conducted a search for data from the ASAS-SN survey, as it overlaps with the major part of the IR light curves from WISE and NEOWISE.Since 2012, ASAS-SN monitored nearly the entire visible sky to a V-band depth of 17 mag with a cadence of 2-3 days using twenty telescopes each of 14 cm diameter.Therefore, we selected AGNs with g < 17.5 mag from our sample of 14,584 Type 1 AGNs, which resulted in 732 targets.We searched light curve data from the publicly available ASAS-SN light curve servers, i.e., ASAS-SN Photometry Database 1 (Shappee et al. 2014;Jayasinghe et al. 2019).Among the 732 targets, we found ASAS-SN light curve data for 142 AGNs.Thus in Sample 1, we arrived at an initial number of 142 sources.
We also collected V-band photometric light curve data for the 142 targets from the CRTS.This survey uses the 0.7 m Catalina Schmidt Telescope in Tucson, Arizona, the 1.5 m telescope on Mt.Lemmon, Arizona and the 0.5 m Siding Springs Survey (SSS) telescope in Australia to search for optical transients with timescales of minutes to yr in the northern sky.CRTS uses SExtractor (Bertin & Arnouts 1996) to process the images following standard techniques.The limiting V-band magnitude for the 1.5 m and 0.7 m are 21.5 mag and 19.5 mag, respectively.The photometry for sources brighter than 13 V-mag is problematic because of saturation and non-linearity response of the CCD.
In addition, we obtained photometric light curve data from the PTF Data Release 3 (DR3) when available.This survey used the 48 inch Samuel Oschin Telescope at the Palomar Observatory in northern San Diego County, California.The photometry was performed in two broadband filters, SDSS g ′ and Mould R. The survey provides a depth of m g ∼ 21.3 mag and m R ∼ 20.6 mag with an exposure of 60 s.The PTF data were used to fill the gap between the CRTS and ASAS-SN light curves.
Finally, we gathered light curve data from the ZTF Data Release 14 (DR14).ZTF is a wide-field timedomain survey, that uses the 1.2 m Schmidt telescope at the Palomar Observatory with 47 deg2 field of view.It provides photometric light curve data with a cadence of 3-4 days in the gri bands since 2018 March (MJD ∼ 58194).It scans the entire northern visible sky to median depths of g ∼ 20.8 mag and r ∼ 20.6 mag (AB, 5σ) with a 30 s exposure (Masci et al. 2019).ZTF uses pointsource-function-fit photometry to build the light curves.We selected ZTF data with the condition catflags=0 to avoid the data contaminated by clouds, moonlight or bad seeing.
In summary, we collected optical photometric light curve data for 142 Type 1 AGNs in Sample 1 with a baseline from 2005 to 2021.This sample consists of targets with bolometric luminosity within the range of 10 43.4 to 10 46.5 erg s −1 at z < 0.3.
We utilized the publicly available light curve data in the optical band from Yang et al. (2020) for 587 Type 1 AGNs in Sample 2 for our own analysis.

IR data from WISE
WISE surveys the full sky at wavelengths centered at 3.4 (W1), 4.6 (W2), 12 (W3) and 22 µm (W4) using a 40 cm diameter infrared telescope since 2010.We collected profile-fit photometric light curve data of our targets in the W1 and W2-bands from the Near-Earth Object Wide-field Infrared Survey Explorer Reactivation Database (NEOWISE-R), WISE All-Sky and WISE Post-Cryo Database from the NASA/IPAC In-frared Science Archive 2 .The retrieved WISE data covers a time period from 2010 to 2020.There is a gap of 3 yr between 2011 and 2014.We only selected goodquality photometric measurements following the suggestions given at https://wise2.ipac.caltech.edu/docs/release/neowise/expsup/sec2 3.html.
For Sample 2, we acquired W1-band light curve data from Yang et al. (2020), while the light curve data in the W2-band was retrieved from the WISE database.
The survey details utilized in this study are summarized in Table 1, while an overview of the initially selected sample properties can be found in Table 2. Figure 1 presents a comparison between the redshift and luminosity distributions of our initially selected sample and those reported in the literature for DRM.Sample 1 covers a range of luminosity with redshifts below 0.33, while Sample 2 mainly represents high-luminosity and high-redshift AGNs.Therefore, our total sample encompasses the range of luminosity and redshift that has been used to measure the torus size in AGNs to date.

Photometric Calibration of the optical light curves
The optical light curves originate from different surveys, each observed with distinct telescopes and filter properties.For instance, the CRTS observations are unfiltered.We converted the CRTS magnitudes to V-band fluxes following Hovatta et al. (2014).While ASAS-SN provides light curve data in V-band, PTF have data either or both in SDSS-g and Mould R-bands.The ASAS-SN and PTF magnitudes were converted to fluxes using the conversion factors given in Bessell et al. (1998) and Fukugita et al. (1996) in their respective filters.We also utilized ZTF-g band data to construct our optical light curves and used ZTF-r when ZTF-g is not available.The ZTF magnitudes were converted into fluxes using the zero-points provided by the Spanish Virtual Observatory (SVO) filter profile service3 .We took the average of the fluxes if more than one observation were taken in a single night.
Various optical monitoring surveys utilize different filters and adopt various aperture schemes in photometry, generating systematic offset among the continuum light curves.Thus, it is necessary to inter-calibrate heterogeneous data sets and include systematic errors.We intercalibrated the light curves from different surveys with various optical bands using the python-based software  PyCALI4 (Li et al. 2014).PyCALI uses Bayesian statistics to determine the best scaling factors for inter-calibration by modeling each light curve using damped random walk model (DRW; Kelly et al. 2009;MacLeod et al. 2010) for describing AGN variability.In this process, we used ASAS-SN light curve as the reference to scale all the other light curves, because ASAS-SN covers a significant portion of the overlapping region between the optical and WISE light curves.In this inter-calibration we ignored the small optical inter-band continuum lags (accretion disk RM; Jiang et al. 2017;Mudd et al. 2018;Guo et al. 2022), which are expected to be a few days (< 10 days; Jha et al. 2022), negligible compared to IR time scale.In our subsequent analysis, we adopted the total error computed by PyCALI, which accounts for measurement and systematic errors.For Sample 2, we used DES g-band light curves that were finally calibrated using various optical data sources from Yang et al. (2020).

WISE Light Curves
The WISE magnitudes are obtained using profilefitting photometry and circular aperture photometry5 .WISE Vega magnitudes were then converted to AB magnitudes using m AB = m V ega + ∆m, where ∆m = 2.699, 3.339, 5.174, and 6.620 for the W1, W2, W3, and W4bands, respectively (Jarrett et al. 2011), and finally into flux densities following Jarrett et al. (2011).WISE observed each object for 15-17 epochs, with 10-20 frames within a 36 hour window in each epoch during the period 2010-2020.Thus, we calculated the mean value of the IR fluxes within a 6-month window for both W1 and W2-band light curves.
The light curves were not corrected for the host galaxy contribution considering a constant host star-light contamination, which makes a negligible effect on the lag measurements.We obtained the errors in the flux measurements through the standard propagation of errors.In Figure 2, we show light curves for 6 targets in the optical and IR bands.The displayed light curves reveal a correlated flux variations between the optical and IR-(W1, W2) bands for the selected AGNs.A few more light curves for objects with various lags and luminosities are shown in Figures A1 and A2 for Sample 1 and Sample 2, respectively.The light curve data are presented in Table 3.The final calibrated optical and WISE light curves, with a typical cadence of 5 days in the optical and 180 days in the WISE light curves, span over 16-20, and 10 yr, respectively.In this section we describe our lag measurements based on these light curves, quality assessment, and correction for z and AD contamination in WISE light curves.

Lag measurements between the optical and IR light curves
All the AGNs selected for this study exhibit synchronized flux variations over time in both the optical and IR W1 and W2-bands, with the optical observations consistently leading the IR measurements.To determine the time lag between the optical and IR continua, we employed the Interpolated Cross-Correlation Function (ICCF; Gaskell & Sparke 1986;Gaskell & Peterson 1987) as the primary method.It determines time lags by linearly interpolating two light curves and calculating the cross-correlation coefficient (r).We evaluated CCF by setting the lag search range to ∼ 70 % of the baseline defined by the optical light curve, which is sufficient for robust time lag measurements (Grier et al. 2017(Grier et al. , 2019;;Yang et al. 2020).We applied the ICCF method to measure time lags in 142 AGNs from Sample 1 and 587 AGNs from Sample 2, respectively.We first located the most significant peak within the search range and the lag was estimated from the centroid (τ cent ) of the CCF around the most significant peak, defined as (Peterson et al. 1998) where, τ cent was obtained by considering all the points within 80 % of the maximum of the CCF (r max ).The uncertainties in the measured lags were estimated using a model-independent Monte Carlo simulation based on flux randomization (FR) and random subset selection (RSS) as described in Peterson et al. (1998); Wandel et al. (1999);and Peterson et al. (2004).This process was repeated for 5000 iterations, and the centroid of the lag was determined in each iteration.The crosscorrelation centroid distribution (CCCD) was then constructed.As the CCCD was not Gaussian, the lag was evaluated from the median of the distribution and the lower and upper uncertainties in the measured lag are that values at the 15.9 and 84.1 percentile of the distribution.This is equivalent to 1σ error for a Gaussian distribution.
We also used JAVELIN to find the time delays between the optical and IR light curves and validated the ICCF measurements.In JAVELIN the driving optical continuum was first modeled using the DRW model with two free parameters; amplitude (σ d ) and time scale of variability (τ d ).To recover the IR continuum light curve a top-hat transfer function was convolved with the driving optical continuum light curve.We used the same search range as used for ICCF to find the lag, and adopted Markov Chain Monte Carlo approach to find the bestfitting model maximizing the likelihood.However, we observed that JAVELIN tends to produce strong aliases, particularly near the boundaries of the lag search window when a large lag search range was applied, leading to multiple peaks in the lag posterior distribution.This is a feature of aliasing in the light curves resulting from sparsely sampled multi-year data with seasonal gaps (Grier et al. 2019;Homayouni et al. 2020).Thus, we excluded the false peaks in the posterior distribution.First, we identified the most significant peak around the location of the ICCF peak, and subsequently excluded false peaks located outside of the ±1000 days range from the ICCF peak in the posterior distribution.The final lag was determined by calculating the median of the lag probability distribution around the most significant peak in the distribution.
Figure 3 shows examples of the results from ICCF and JAVELIN for the targets shown in Figure 2.For a few targets, along with a primary peak at smaller lag, there is a secondary peak at a larger lag (for examples, see the CCFs of OB44 and OB64 in Figure 3).While the primary peak represents the actual lag between the optical and W1-band light curves, the other peak arises mainly because of the shape and aliasing in the light curves (Guise et al. 2022).For example, in Figure 3, the CCF of the target OB44 shows a primary peak at around 110 days, and a secondary peak around 2200 days.The corresponding Auto Correlation Function (ACF) of the optical light curve has also two peaks, one around 0 lag and another secondary peak around 2090 days.Thus, the IR W1 light curve will also correlate with the optical light curve at the lag τ =109 days and at a latter time τ + 2090 days, resulting two peaks in the CCF.The same is applicable for the target OB64 shown in Figure 3. Thus, we calculated lags around the primary peaks in these cases.We obtained te lags between the optical and WISE W1 and W2-band light curves for all of those 729 (142+587) AGNs.

Quality assessment
In this section, we assess the dependability of the lag measurements obtained through cross-correlation analysis.The accuracy of the recovered time lags between the optical light curves, derived from various surveys, and sparsely sampled IR light curves is influenced by different sources of uncertainty.The presence of substantial data gaps in the IR light curves, coupled with a limited number of data points in the optical-IR light curve pairs, can lead to false positive peaks near the boundaries of the lag search window when using the ICCF method.The probability p(rmax) represents the likelihood of uncorrelated light curves producing a crosscorrelation coefficient equal to or exceeding rmax.This probability is presented as a function of the maximum crosscorrelation coefficient rmax in the left and right panels for W1 and W2-bands, respectively.The vertical dashed line represents rmax=0.6,whereas p(rmax)=0.2 is shown by the horizontal dashed line.The reliable time lags are defined with rmax ≥ 0.6, and p(rmax) < 0.2.
These factors collectively impact the reliability of the obtained time lags.Consequently, it is crucial to exercise caution when selecting targets with significant and reliable lags.
The maximum cross-correlation coefficient (r max ) is a key factor to assess the reliability in lag measurements.This coefficient signifies the strength of the correlation between two light curves.Previous DRM studies have often employed specific thresholds for r max , such as r max > 0.5 (Yang et al. 2020) or r max > 0.4 (Chen et al. 2023), as a means to determine the reliability of lag measurements.
However, it is crucial to acknowledge that relying solely on r max is insufficient for guaranteeing lag reliability.There are situations where r max might be relatively high even when the light curves possess only a limited number of data points but exhibit strong linear variation.Even small variability amplitude compared to the flux uncertainties in the light curves may result high r max as long as certain features in the light curves align well.Moreover, it is crucial to emphasize that the calculations for r max do not incorporate the uncertainties in flux.
To address these issues and evaluate the reliability of lag measurements, a new method called PyI 2 CCF6 has been introduced (Guo & Barth 2021;U et al. 2022, Guo et al., in prep).This method employs a null-hypothesis approach to assess whether the observed cross-correlation could be equally derived from two uncorrelated red-noise light curves.To put this method into practice, we utilized the publicly available PyI 2 CCF tool to investigate the reliability of our lag measurements.Our methodology involved generating 10 3 mock light curve pairs for each object, maintaining the same sampling and cadence as the observed light curves while utilizing DRW models.Subsequently, we calculated the lag significance indicator, referred to as p(r max ) τ >0 , by tallying the occurrences of positive lags (τ > 0) where the r max value exceeded the observed r max across all simulations.It is important to mention that we retained the same searching window that was initially used during the lag calculations.Finally, we apply a set of criteria to exclude measurements deemed unreliable, such as: 1. We define the lag measurements as reliable if r max ≥ 0.6 and p(r max ) τ >0 < 0.2.These thresholds are adopted following U et al. ( 2022) and Woo et al.
(2024) (See Figure 4).Using these criteria, we identify a total of 85 targets in the W1-band and 79 in the W2-band within Sample 1 as having reliable lags.Similarly, for Sample 2, out of 587 targets, we select 399 in the W1-band and 372 in the W2-band.
2. However, we discard 4 and 3 targets from Sample 1 in the W1 and W2-bands, respectively, due to negative lags within a 1σ uncertainty range, as negative lags are physically unacceptable (See Figure 5,left).Similarly, from Sample 2, 16 and 14 targets are rejected with negative lags within 1σ uncertainty in W1 and W2-bands, respectively.Additionally, we eliminate targets displaying a flat CCF spanning over thousands of days, as this condition results in an average lag computed over the entire search range (refer to Figure 5, right).Consequently, 18 targets from Sample 2 in both W1 and W2-bands are excluded from the analysis.
After applying the criteria, we narrowed down from 142 targets to 81 targets in the W1-band and 76 targets in the W2-band in Sample 1.In the case of Sample 2, we obtained 365 AGNs with W1-band light curves and 340 with W2-band light curves, out of the initial sample of 587 targets.In summary, our final sample consists of 446 (416) AGNs with a reliable and significant dust lag measurements with W1-(W2) band light curves.Note that if we ignore the systematic uncertainty of the light curve calibration between two different telescope data sets (see Section 3), the error of the photometry decreases by 30-40%.This reduction of the photometry errors in the light curves will slightly increase the quality of the cross-correlation.To evaluate the effect of the error treatment on the quality assessment, we remeasure the p(r max ) τ >0 for Sample 1 after removing the systematic error from the total error.We find that the change of the sample based on the quality assessment (see Figure 4) is insignificant.
In Figure 6, we present a comparison of the W1band lags in the observed-frame between the ICCF and JAVELIN methods.The final sample is represented by brown data points, while the rejected AGNs are shown in black.The majority of the rejected targets deviate from the expected 1:1 relation.In contrast, the lags obtained for the final sample demonstrate a good consistency, with a scatter of 0.17 dex.We present a few examples of the selected ICCF and JAVELIN analyses for Sample 1 and Sample 2 in Figures A3 and A4, respectively.
We checked that shifting the optical light curve by the median lag value (τ cent ) obtained from ICCF resulted in a good alignment with the IR light curves, as demonstrated in Figures 2, A1, and A2.Moreover, the error estimation in ICCF analysis is more conservative (Guo et al. 2022), and the ICCF measurements incorporate a cross-correlation reliability test.Consequently, in all subsequent discussions, our primary focus will be on the lags derived from the ICCF analysis.

Correcting for redshift effects on dust lags
The AGNs under investigation in this study span a wide range of redshifts as illustrated in Figure 1.Because of the wavelength-dependent nature of dust lag, objects at higher redshifts will exhibit shorter rest-frame time lags in specific wavelength bands, than the same luminosity AGNs at lower z (Sitko et al. 1993;Oknyanskij & Horne 2001;Yoshii et al. 2014;Minezaki et al. 2019;Sobrino Figaredo et al. 2020;Chen et al. 2023).This redshift effect must be corrected, especially when dealing with a sample over a large redshift range.
We used the correction factor in the form of (1+z) γ , as introduced by Yoshii et al. (2014) and Minezaki et al. (2019), to account for wavelength-dependent dust lags in order to remove the redshift effect.First, we determined the lag ratio, τ ratio = τ W 2 /τ W 1 , between The distribution of the lag-ratio (τratio = τW 2/τW 1) between the dust lags in W2 and W1-bands with and without AD contamination correction shown by red and black lines, respectively.The median values with their standard deviations are mentioned in the figure .dust lags in W2 and W1-bands for our selected final sample and computed γ using the formula, γ = log(τ W 2 /τ W 1 )/log(λ W 2 /λ W 1 ). Figure 7 represents the distribution of τ ratio , both with and without AD contamination correction (see Section 4.4 for details of AD contamination correction) obtained from ICCF analysis.Our analysis reveals median values of τ ratio to be 1.19 ± 0.43 (1.16 ± 0.55) and 1.22 ± 0.50 (1.26 ± 0.85) for the AD contamination corrected and uncorrected cases, respectively, obtained from ICCF (JAVELIN) analysis.These values closely align with the median values of τ ratio = 1.15 reported by Lyu et al. (2019) and 1.26 reported by Chen et al. (2023) within the 1σ range of the distribution.
By adopting a median τ ratio of 1.20 for both the AD contamination corrected and uncorrected cases, we find that γ ∼ 0.62, leading to a correction factor of (1 + z) 0.62 .
According to the dust radiation equilibrium model, the evaporation radius from the central ionizing source varies with the UV continuum luminosity and dust grain temperature, as R dust ∝ L 0.5 U V T −2.8 (Barvainis 1987).Since the dominant source of NIR emission from the torus is attributed to the thermal emission from dust with a temperature T, the wavelength dependency of the dust lag is given by τ ∝ λ 2.8 (Sitko et al. 1993;Oknyanskij & Horne 2001).The correction factor for the redshift effect is expected to be (1 + z) 2.8 .This the-Table 4. Details of the final sample and the results of their time lag measurements for W1 and W2-bands are provided.The complete table is available online in a machine-readable format, separately for Sample 1 and Sample 2. oretical prediction is significatnly different from our empirically derived factor.The discrepancy is partly due to the fact that the dust lag correlates with AGN luminosity with a much shallower slope than the expected value of 0.5 (Minezaki et al. 2019;Chen et al. 2023).Additionally, Oknyanskij & Horne (2001) reported that observed dust lags of certain AGNs are notably smaller than theoretical predictions.This discrepancy is attributed to the presence of thin ring or disk-like dust distributions, diverging from the spherically symmetric distribution of dust particles assumed in the theoretical framework.For empirical constraints, Minezaki et al. (2019) determined γ as 1.18 based on the K-and Hband lag measurements in six nearby Seyfert galaxies.
In the case of W1 and W2-band lag measurements, Chen et al. ( 2023) reported a γ ∼ 0.76 using a sample of 78 AGNs.These results indicate that the power-law index γ of the redshift correction factor may vary depending on the observed wavelength for dust lag measurements (Chen et al. 2023).
Combined with the cosmological time dilation correction factor of (1+z) −1 , the total correction factor for the observed-frame time lag becomes approximately (1+z) γ−1 ∼ (1+z) −0.38 .We uniformly applied this total correction factor to the measured time lag in the observed-frame, thereby obtaining the rest-frame time lags for the total sample.

Correction for the accretion disk (AD) contamination in the IR W1 and W2-band light curves
The observed IR fluxes in W1 and W2-bands contain contribution from the torus as well as from the AD (Tomita et al. 2006;Kishimoto et al. 2008;Lira et al. 2011), the latter of which makes the delay time between the optical and IR light curves shorter than the actual lag.Thus, the AD contamination needs to be removed from the observed IR fluxes to get the actual time lag.The AD emission in the IR can be estimated by considering a power-law spectrum of the accretion disk (Koshida et al. 2014) as follows where, F AD IR (t) and F OP (t) are the AD component of the IR flux and the optical flux in V-band (g-band for Sample 2) at time t, respectively.ν OP and ν IR are the effective frequencies of optical in V-band (g-band for Sample 2) and IR (W1, W2) bands, respectively, and a is the power-law index.We assumed a constant powerlaw index of a = 1/3 as suggested by the standard disk model (Shakura & Sunyaev 1973) to calculate F AD IR (t).Note that the polarized light spectra in the NIR strongly suggest intrinsic AD spectra with a = 1/3 (Kishimoto et al. 2008).In contrast, Koshida et al. (2014) utilized not only a = 1/3 but also a = 0.0, taking into account the potential maximal contribution of AD contamination in the IR bands, while Minezaki et al. (2019) used a = 0.1 following Yoshii et al. (2014).If a is smaller than a = 1/3, the time lag generally becomes larger, because the AD contamination in the IR-band flux is larger (i.e., a larger correction).To empirically test the effect, we performed the ADC with a = 0.0 using a representative subsample of 35 AGNs, which covers the entire range of luminosity and Eddington ratio of the total sample.We obtained insignificant effect, with a 0.03 dex increase in both W1 and W2 lag measurements (see also Koshida et al. 2014;Minezaki et al. 2019).
We estimated F AD IR (t) for each epoch of IR (W1 and W2) observations using the average flux values of the quasi-simultaneous epochs (within ± 5 days) in the optical light curve.If for any epoch in IR light curve, there is no quasi-simultaneous observations in the optical, we obtained the corresponding optical flux from the DRW modeled optical light curve using JAVELIN.It is worth noting that the influence of host galaxy flux contamination in the optical band is more evident in lower redshift AGNs.This may introduce a spurious redshift dependence in the dust lags.However, the effect of the host galaxy contribution to the ADC would be negligible since host galaxy flux remains constant and the hostgalaxy contribution to the total optical flux would not affect the lag measurements while host galaxy contribution reduces the variability amplitude.Finally, F AD IR (t) was subtracted from the observed IR (W1, W2) flux (F IR,obs ) to get the IR flux, which is coming from the torus (F IR,dust ) as Following this analysis, our next objective is to investigate how the AD contamination corrected and uncorrected dust lags are influenced by AGN luminosity and the Eddington ratio for the selected AGNs in our study.The AGN luminosity can be contaminated by stellar-light from the host galaxies.Liu et al. (2019) modeled and subtracted the stellar continua from the spectra and derived the continuum luminosity at restframe 5100 Å (L 5100 Å).We used their estimated L 5100 Å for Sample 1 (z < 0.33) and obtained the AGN bolometric luminosity (L bol ) by using a bolometric correction factor of 9.26 (Richards et al. 2006).The AGN bolometric luminosities for Sample 2 were collected from Yang et al. (2020).The same bolometric correction factor was used for Sample 2. In order to obtain the Eddington ratios for our own sample, we retrieved FWHM of Hβ values from Liu et al. (2019) and Shen et al. (2011).Subsequently, the Hβ BLR size (R Hβ ) were estimated by utilizing the R Hβ -L 5100 Å relation from Woo et al. (2024).Then, we employed the virial equation (Woo & Urry 2002), incorporating the FWHM of Hβ and a virial factor of f BLR = 1.12 (Woo et al. 2015), to derive the single epoch M BH .Finally, we obtained the Eddington ratio (L bol /L Edd ) using Eddington luminosity L Edd = 1.26 × 10 38 M BH /M ⊙ .These Eddington ratios were determined for approximately ∼72% of the AGNs in our final sample, specifically those with available FWHM values for Hβ.
In Figure 8, we present a comparison between the redshift-corrected rest-frame torus sizes, with and without AD contamination correction, denoted as  2023) did not take this corrective approach into account.Consequently, in our study, we implemented this correction for the W1 and W2band light curves.The key findings from this analysis are as follows: (1) correction of AD contamination leads to larger torus sizes in AGNs, (2) because of the positive correlation between torus size and AGN luminosity (see Section 5 for details), slopes in Equations 4 and 5 greater than unity suggest that AD contamination becomes more prominent as luminosity increases, (3) the size ratio D ratio (R IR,ADC /R IR,non−ADC ) between the size with/without the AD contamination correction shows a negligible correlation with the Eddington ratio: Spearman's rank correlation coefficient ρ = −0.003and a p-value of 0.955 (ρ = −0.041,p=0.492) based on W1 (W2)-band lags, implying no significant dependence on Eddington ratio, and (4) intercept values for the W1band are larger than those for the W2-band, exhibiting a similar slope of ∼ 1.12.This indicates that AD contamination is more significant in the W1-band, aligning with expectation as AD contaminates more at shorter IR wavelengths.All of these findings indicate that it is necessary to correct for AD contribution to the observed IR fluxes.
We compare our obtained redshift and AD contamination corrected torus sizes with that available in the literature in Section B, finding good agreement with Lyu et al. (2019), Chen et al. (2023), andYang et al. (2020) (See Figure B1).Table 4 describes the details of the finally selected targets and results of their lag-analysis.

The torus size and AGN luminosity correlation
In this section, we investigate the correlation between the measured torus size based on W1 and W2-band flux variations with AGN luminosity.First we will use bolometric luminosity of each AGN, and then we will also investigate the relation using the torus luminosity measured at W2-band.Note that we use the W1 and W2band lag measurements after correcting for AD contamination and redshift effect as described in Sections 4.4, and 4.3, respectively.
We use the following linear equation to parameterize where α and β are the slope and intercept, respectively.L 0 is the reference point, which is set close to the median of AGN luminosity L AGN used in the fitting.In fitting Equation 6, we employ the Bivariate Correlated Errors and intrinsic Scatter (BCES; Akritas & Bershady 1996).This method is particularly advantageous as it accommodates measurement errors with heterogeneity in variance, incorporates intrinsic scatter, and accounts for correlations within the measurement errors.

The torus size correlation with AGN bolometric luminosity
In Figure 9 we present the best-fit relation between torus size and AGN bolometric luminosity (L bol ) for W1(left) and W2-(right) bands, where L 0 is set to 10 45.7 erg s −1 .The torus size scales with the bolometric luminosity as, R dust,W 1 ∝ L 0.39 bol , and R dust,W 2 ∝ L 0.33 bol with an intrinsic scatter of 0.19 and 0.18 dex, in W1and W2-bands, respectively.Note that if we use the lag measurements without correcting for AD contamination, we obtain similar slopes of α =0.38 and 0.35 within the uncertainties for W1 and W2-bands, respectively, but with a smaller intercept (β) values and slightly larger intrinsic scatters.This result indicates that the presence of AD contamination in the observed IR light curves can cause an underestimation of the torus size.Therefore, correcting for AD contamination is essential to accurately recover the true torus sizes in AGNs.
To test the potential difference of the results due to the fitting method (Park et al. 2012), we employ an ad- ditional approach: Bayesian linear regression using the LINMIX ERR technique (Kelly 2007).This method excels in handling intrinsic scatter within the relationship, addressing measurement errors in both the independent and dependent variables, and effectively managing correlations among the measurement errors.We find that these different fitting methods yield consistent results with a similar slope and intercept within the uncertainties.We also note that two different lag analysis methods, i.e., ICCF and JAVELIN, provide consistent results.
Furthermore, depending on the scatter in the observed R dust,IR -L bol relation and the varied uncertainties in luminosities and torus sizes, we perform inverse linear regression fitting, i.e., AGN bolometric luminosity as a function of torus sizes.The inverse linear fitting results indicate a relationship: , and derived from the AD contamination corrected W1 and W2-band lags, respectively, with a considerably larger intrinsic scatter of σ ∼ 0.4 dex.From these findings, we derive dependencies of torus sizes on L bol as R W 1 ∝ L 0.52 bol and R W 2 ∝ L 0.49 bol , closely aligning with the expected value of 0.5 from the dustsublimation equilibrium model.It is important to note that the determination of the dust torus size occurs at a distance from the accretion disk, where dust particles evaporate because of the central ionizing luminosity.Hence, the AGN luminosity determines the torus size in AGN.Consequently, in alignment with our scien-tific objectives, we assert that our choice of the forward linear regression fitting between torus size as a function of AGN luminosity is more appropriate and is consequently employed for further analysis (refer to Isobe et al. (1990) and Feigelson & Babu (1992) for the choice of fitting method).
The slopes (α) we obtained in our analysis exhibit values ranging from 0.33 to 0.39.These slopes are notably shallower than the expected value of 0.5 predicted by the dust sublimation model.In contrast, Lyu et al. (2019) and Li & Shen (2023) bol ), and R dust,W 1 ∝ L 0.51 bol , respectively, in their findings, which differ considerably from our results.However, it is noteworthy that Chen et al. ( 2023) reported a value of α = 0.37 ± 0.03 for both W1-and W2-bands based on the study of DRM in AGNs with Hβ RM data.Our obtained slopes align with their findings within 1σ uncertainty.Furthermore, Minezaki et al. (2019) also reported a consistent slope of 0.42 from DRM based on K-band dust lags.On the other hand, recent work by Gravity Collaboration et al. ( 2023) mentioned a much shallower slope of approximately α = 0.4 from optical/NIR interferometry.Consequently, the observation of smaller slope values in the relation between torus size and luminosity in DRM is consistent with these findings and is not unexpected.
Our results indicate that the slope is likely shallower than the expected value of 0.5 from the dust radiation equilibrium model.This deviation could be attributed to the possibility that the optical luminosity at 5100 Å, which is used to obtain AGN bolometric luminosity, may not accurately represent the ionizing luminosity (Du & Wang 2019).

The torus size correlation with AGN IR luminosity
To further investigate the potential discrepancy in the observed R dust,IR -L bol relationship compared to the expected relation, we obtain a torus size-luminosity relationship (R dust,IR -L W 2 ) based on the luminosity of the infrared continuum emission in the W2-band, as it serves as luminosity indicators for AGNs, less obscured by dust.Thus, these relations can be applicable for obscured AGNs.We employ Equation 6 to characterize the relationship between R dust,IR and L W 2 , where L 0 is set to 10 45.0 erg s −1 .L W 2 is luminosity measured with W2band, which was derived from the average flux values of the W2-band light curves studied in this work, without correction for host contamination or galactic extinction.In Figure 10, we present the R dust,IR -L W 2 relationships for AD contamination corrected W1 and W2-band light curves in the left and right panels, respectively.All the fitting results are summarized in Table 5.
Interestingly, we observe similar slopes of ∼ 0.39 (W1) and 0.33 (W2) in the R dust,IR -L W 2 relations as those obtained from R dust,IR -L bol within the uncertainties, irrespective of whether we apply corrections for AD contamination or not.Our obtained slope differs considerably from the value of 0.5 reported in the torus size -IR luminosity relation at 12 µm, as per the Kband DRM measurements involving only 10 AGNs by Minezaki et al. (2019).This discrepancy is likely due to the limited quantity of data and luminosity range in their study.
The similarity in shallower slopes observed in the torus size -luminosity relations using both optical and IR luminosities can be attributed to the self-shadowing effect of the slim disk model (Kawaguchi & Mori 2011;Wang et al. 2014a;Du et al. 2018).Since dust particles within the torus reach thermal equilibrium with the IR radiation re-emitted by dust exposed to ionizing UV radiation, we would expect a slope of 0.5 in the R dust,IR -L W 2 relationship (Barvainis 1987).However, the self-shadowing effect of slim accretion disks leads to a contraction of the ionization front of the BLR and the torus, resulting in a reduction in torus size with increasing accretion rate.This effect yields the shallower slopes observed in both the R dust,IR -L bol and R dust,IR -L W 2 relationships.For further details, please refer to Section 6.4.3) name of the bands used to measure the lags, (4) slope of the best-fits with α as free parameter using lags from ICCF, (5) intercept of the best-fits with α as free parameter using lags from ICCF, (6) intrinsic scatter in dex, (7) intercept of the best-fits with α fixed to 0.5 using lags from ICCF, (8) slope of the best-fits with α as free parameter using lags from JAVELIN, (9) intercept of the best-fits with α as free parameter using lags from JAVELIN, and (10) intrinsic scatter in dex.

Comparison between torus sizes in W1 and W2bands
In Figure 11, we compare torus sizes in the W1 and W2-bands, revealing that R dust,W 2 exceeds that in the W1-band.This result aligns with our expectations, considering that lag generally increases with wavelength.Moreover, the mean torus size ratio between W2 and W1-bands is found to be ∼ 1.20.A linear fit to the torus sizes between W1 and W2 yields log(R dust,W 2 /400 light-day) = (0.90 ± 0.02) log(R dust,W 1 / 400 light-day) + (0.09 ± 0.01) (7) and after correcting for the AD contamination, we found log(R dust,W 2 /400 light-day) = (0.90 ± 0.02) log(R dust,W 1 / 400 light-day) + (0.08 ± 0.01) (8) Indeed, the fitting results indicate that R dust,W 2 is statistically larger than that in the W1-band.

Comparison between torus and BLR sizes
We have identified a total of 19 AGNs within our sample, all of which possess prior measurements of BLR size and M BH based on Hβ line obtained through BLR-RM analysis available in the existing literature (See Table 6).
Figure 12 illustrates a comparative analysis between the torus sizes observed in the W1 and W2-bands and their corresponding Hβ BLR sizes (R Hβ ) as determined by BLR-RM in the left and right panel, respectively for those 19 AGNs as shown by the red points.Notably, our examination reveals a discernible correlation between BLR sizes and our obtained torus sizes with a Spearman's rank correlation coefficient of ρ = 0.73 with a p-value of 2.56 × 10 −4 for the W1-band, and ρ = 0.71 with a p-value of 1.04 × 10 −3 for the W2-band.We employ the following linear equation to parameterize the relationship between R Hβ and R dust,IR log(R Hβ /1 light-day) = ϵ log(R dust,IR /1 light-day) + K (9) where ϵ denotes the slope and K is the intercept.We find that the BLR size scales with the torus size as, R Hβ ∝ R 1. 42±0.32  dust,W 1 with K = −2.18± 0.88, σ = 0.37 dex, and R Hβ ∝ R 1. 42±0.31  dust,W 2 with K = −2.29 ± 0.87, σ = 0.39 dex.The median size ratios between torus and BLR are approximately ∼ 9.5 and 15.1 for W1 and W2bands, respectively.In addition, Koshida et al. (2014) reported a significant four to five-fold disparity in the mean reverberation radii between the innermost dust torus based on K-band lag and the BLR.A median size ratio between the torus and BLR, R dust,K /R Hβ , was determined to be approximately ∼6.5 based on 10 objects (Gandhi et al. 2015).Du et al. (2015) reported R dust,K /R Hβ ∼4 for 10 low-accreting AGNs, and R dust,K /R Hβ ∼7 for high-accreting AGNs with a smaller sample size (4 objects).Kokubo & Minezaki (2020) also reported R dust,K /R Hβ ∼4 based on measurements from 15 Seyfert galaxies.Moreover, Mandal et al. (2018) found about 11.5 times larger torus size (based on K-band) than BLR size in H0507+164.Recently, Chen et al. (2023) found R dust,K /R Hβ ∼6.2, R dust,W 1 /R Hβ ∼9.2, and R dust,W 2 /R Hβ ∼11.2 based on a significantly larger sample of 78 AGNs with Hβ-RM measurements.Our findings indicate that the ratio of R dust,W 1 /R Hβ is ∼ 0.98 ± 0.36 dex, suggesting a torus size roughly 9.5 times larger than that of the BLR, which closely aligns with the results of Chen et al. (2023).In contrast, we observe a slightly higher ratio of R dust,W 2 /R Hβ at around 1.18 ± 0.36 dex, indicating a torus size approximately 15.1 times larger than the BLR size based on W2-band lags, derived from a sample of 19 AGNs.The discrepancy observed in W2-band data may be attributed to the smaller number of AGNs with Hβ-RM data in our sample compared to their sample.Nonetheless, taking into account the rms scatter of ap- Å marked with a is not corrected for host-starlight contamination.
proximately 0.36 dex, our measurements are consistent with previous findings reported in the literature.It is worth highlighting a consistent trend observed across these 19 AGNs: the torus sizes consistently surpass the corresponding BLR sizes.Since, the above R Hβ -R dust,IR relations are based on limited number of AGNs, to expand the sample size we incorporate measurements from Chen et al. (2023), who reported torus sizes for 78 AGNs with R Hβ measurements.Their sample does not cover the S82 region, and only 28 among the 76 targets from their sample overlap with our initial Sample 1 before applying the lag-quality assessment.After implementing our lag-quality assessment (see Section 4.2), only 12 targets remain, that overlap with their sample.Note that, Chen et al. (2023) did not correct for AD contamination in their measured lags.A direct comparison of the torus sizes based on W1-(W2) band lags for those 12 targets between our measurements and that from Chen et al. ( 2023) is discussed in Section B, which shows an offset of 0.16 (0.14) dex between the two sets of measurements.Additionally, we find that, among the remaining 66 (78-12) tar-gets, 60 have torus size measurements with r max ≥ 0.6.Thus, we include these additional 60 targets from Chen et al. (2023) in the R Hβ -R dust,IR plane represented by the grey points in Figure 12 and obtain the best fit relations using Equation 9. We use same redshift correction factor of (1 + z) −0.38 to get the rest-frame torus size, while disregarding the extra uncertainties caused by the AD contamination correction which was not performed for the additional 60 targets.For the total BLR-RM sample (19+60), R Hβ varies with R dust,IR as, R Hβ ∝ R 1.20±0.09dust,W 1 with K = −1.44 ± 0.23, σ = 0.34 dex, and R Hβ ∝ R 1.22±0.08dust,W 2 with K = −1.65 ± 0.21, σ = 0.33 dex.
While Chen et al. (2023) reported a linear relationship between R Hβ and R dust,IR , our findings reveal a superlinear relationship characterized by a slope exceeding unity.This aligns with the findings of Gravity Collaboration et al. (2023), who also identified a slope greater than unity of ∼ 1.12 and 1.10 in the R Hβ -R dust,IR relationship using K-band torus sizes derived from RM and optical/NIR interferometry, respectively.According to the photoionization model, the relationship be-tween R Hβ and R dust,IR can be expressed as (Chen et al. 2023): where, L ion and L OU V represent the ionizing and optical-UV luminosity, respectively.σ s , n e , k, T c , and T IR denote the Stefan-Boltzmann constant, electron number density, Boltzmann constant, temperature of the ionized clouds, and temperature emitting IR photons, respectively.The deviation from a linear relationship between R Hβ and R dust,IR may occur for two primary reasons.First, within AGNs, BLR gas clouds could exhibit a range of n e (Baldwin et al. 1995).Second, the optical luminosity measured at 5100 Å, utilized to derive L bol , might not precisely represent the ionizing luminosity, implying that L bol is not directly proportional to L ion (refer to the discussion by Czerny et al. 2019;Fonseca Alvarez et al. 2020).Consequently, the observed super-linear relationship between R Hβ and R dust,IR is not unexpected.Moreover, when considering the relationship between torus size in W1-(W2) band and AGN bolometric luminosity, we find that the BLR size scales with L bol as R Hβ ∝ L 0.47 bol (R Hβ ∝ L 0.40 bol ).Thus, relying on torus sizes in W1 and W2-bands, the slope of R Hβ -L bol relation ranges between 0.47 and 0.40 with a mean value of ∼ 0.43, which is smaller than the expected value of 0.5 from photoionization model (Davidson 1972).However, it closely aligns with the observed slope of 0.402 obtained from BLR-RM, as reported by Woo et al. (2024) within the uncertainty.
Consequently, these findings imply that the BLR may either constitute an inner extension of the torus or originate from infalling material associated with the torus.

Dependency of torus size on wavelength
The availability of nearly simultaneous monitoring observations across various IR-bands, coupled with optical data, enables us to determine torus sizes across distinct IR wavelengths and, in turn, examine whether there exists any wavelength-dependent trend in torus sizes.However, in our sample, only two targets have previous DRM measurements based on K-band data in the literature: ID2149795 (J213818.96+011222.3, with W1 and W2) and ID392158 (J232640.01−003041.4,with W1).Thus, we cannot establish a one-to-one relationship between the lags in the K-band and those in the other IR W1 and W2-bands.
Instead, we present the torus size -AGN bolometric luminosity relations based on our work using W1 and W2-band, along with the work by Minezaki et al. (2019) using K-band in Figure 13 (left).All these R dust -L bol relations have been adjusted for the AD contamination in their respective IR (W1, W2 and K-band) light curves.The R dust,K -L V relation is based on Kband lag and V-band luminosity of a sample of 36 quasars with relatively moderate to low luminosities (10 42.7 erg s −1 < L bol < 10 46.4 erg s −1 ), the majority of which have bolometric luminosities below 10 46 erg s −1 .We converted the V-band luminosity of the R dust,K -L V relation to bolometric luminosity using a correction factor of 10 (Yang et al. 2020).If we consider the same R dust,K -L bol relation to hold for high luminosity AGNs, Figure 13 (left) illustrates wavelength-dependent torus sizes: R dust,W 2 > R dust,W 1 > R dust,K for a given luminosity.We notice that as luminosity increases, the torus sizes in different IR bands gradually converge, suggesting that higher luminosity objects tend to exhibit a more compact torus structure, a notion in alignment with the findings of Kishimoto et al. (2011b) from interferometric observations.
The innermost radius of the torus composed of silicate grains is defined by (Barvainis 1987;Nenkova et al. 2008b) where T sub represents the dust sublimation temperature.Similarly, the sublimation radius for pure graphite grains is given by (Mor & Netzer 2012) The hot dust spectrum obtained by Mor & Netzer (2012) suggests that the existence of hot graphite dust makes a substantial contribution to the observed luminosity at 3.6 and 4.5 µm.On the other hand, using T sub = 1700 K and a dust grain size of ∼ 0.1 µm, Koshida et al. (2014) found that the sublimation radius is close to the K-band dust reverberation radius.For comparision we over-plot the theoretical prediction of the R dust -L bol relations as defined by the Equations 11 and 12 for two distinct sublimation temperatures T sub = 1500 K and 1800 K, respectively.
Using Equations 11 and 12, we estimate the sublimation radius of a low-luminosity AGN with L bol = 10 44 erg s −1 .For silicate dust with T sub = 1500 K, R dust,Si is 0.13 pc, while for pure graphite with T sub = 1800 K, R dust,C is 0.05 pc.Whereas, for a highluminosity AGN with L bol = 10 48 erg s −1 , the sublimation radius is significantly larger as ∼ 13.0 pc and 5.0 pc, respectively, for silicate with T sub = 1500 K and pure graphite grains with T sub = 1800 K.We also empirically determine the torus size from the observed R dust -L bol relationships (Equation 6).The derived radii are R dust,W 1 = 0.08 pc (R dust,W 2 = 0.13 pc), and R dust,W 1 = 3.0 pc (R dust,W 2 = 2.8 pc), respectively, for AGNs with L bol = 10 44 and 10 48 erg s −1 .By comparing the theoretical estimates and the empirically determined torus sizes, we find that for low-luminosity AGNs (∼ 10 43.4 − 10 45 erg s −1 ), the primary source of the 3.4 and 4.6 µm emission lies within the region bound by the graphite and silicate sublimation radii, respectively.These results align with the previous findings with the Spitzer Space Telescope data.For example, Vazquez et al. (2015) derived the torus size of NGC 6418 (L bol ≥ 10 43.3 erg s −1 ) as R dust,3.6µm= 0.031 pc with 3.6 µm data, and R dust,4.5µm= 0.040 pc with 4.5 µm data.
In contrast, for very high luminosity AGNs (L bol ≥ 10 47 erg s −1 ), W1 and W2-band emissions originate from a similar region, which is roughly consistent with the sublimation radius of graphite grains with T sub = 1800 K (See left panel of Figure 13).This alignment is well-supported by the computed emission spectrum of graphite dust in the model by Mor & Netzer (2012), which peaks in the 2-3 µm range and exhibits substantial emission in the 3.4−4.6µm range as well.Note that the observed R dust -L bol relations show shallower slopes than expected from theoretical relations, resulting in a deviation of the observed torus sizes from theoreti-cal predictions as luminosity increases.In addition, we present the Hβ BLR size-luminosity relation from Woo et al. (2024), demonstrating that the torus lies outside the BLR in AGNs.
By combining the best-fit parameters obtained by Minezaki et al. (2019) using K-band lags and our measurements in the W1 and W2-bands, we determine torus sizes from the best-fit torus size -AGN bolometric luminosity relations at K, W1 and W2-bands at luminosities L bol = 10 44 , 10 46 , and 10 48 erg s −1 spanning the entire dynamical range used in this study (see right panel of Figure 13).At L bol = 10 46 erg s −1 , nearly the median luminosity of our sample, we obtain the ratios of the dust torus sizes: R dust,K : R dust,W 1 : R dust,W 2 = 1.0:1.5:1.8.We find that the torus size scales with the wavelength as R dust ∝ λ 0.80 at L bol = 10 46 erg s −1 , indicating a stratified torus structure.Note that the dependence of R dust on wavelength becomes weaker as luminosity increases as indicated by the shallower slopes observed in the R dust -L bol relations.
Extensive efforts, comprising both observational studies (Tomita et al. 2006;Yoshii et al. 2014;Minezaki et al. 2019) and theoretical investigations (Barvainis 1987;Kawaguchi & Mori 2011), have been dedicated to gain insights into the structure of the dust torus in AGNs.The presence of a radial temperature gradient in the dust torus was reported by Tomita et al. (2006).Additionally, a clumpy dust torus model (Nenkova et al. 2008a;Kawaguchi & Mori 2011) proposes a vertical stratification of rotating dust cloud distributions within galactic nuclei (García-Burillo et al. 2019). Furthermore, Almeyda et al. (2017, 2020) through simulations of the reverberation response of a clumpy torus, found that the torus response exhibits a strong wavelength dependency.This phenomenon arises from the temperature gradient within the cloud surfaces of the torus and the highly anisotropic cloud emission at shorter wavelengths.Additionally, our findings support the idea that different IR emissions originate from distinct regions of the torus, thus lending support to the notion of a stratified torus structure.

Scatter in the R dust,IR -L bol
In this section, we explore the potential correlation between the scatter around the R dust,IR − L bol relations and various other AGN parameters, including the Eddington ratio, the flux ratios [O iii]/Hβ (R [OIII] ), Fe ii/Hβ (R F eII ), and redshift z.We define the deviations of the measured torus sizes from the corresponding best-fitting relationships as ∆R IR = logR IR -logR best−f it .
We compare ∆R IR with the Eddington ratio in Figure 14 (top left) considering the slope in the R dust,IR -L bol relation as free parameter.It is important to note that the Eddington ratio in this analysis is derived from a single-epoch M BH (see Section 4.4) using the BLR size -luminosity relation (Woo et al. 2024), which is based on AGNs in the luminosity range 10 42 < L 5100 Å < 10 46 erg s −1 with the majority of the AGNs having luminosities below 10 44.5 erg s −1 .This limitation imposes constraints on the precision of M BH measurements at higher redshift, where high-luminosity AGNs are more commonly observed.As a result, the canonical single-epoch M BH can lead to substantial over-estimations and, consequently, underestimations in the Eddington ratio, especially for AGNs with high accretion rates frequently found in the high-luminosity regime (also see section C).
Nonetheless, despite these limitations, we identify a mild negative correlation between ∆R IR and the Eddington ratio, with Spearman's rank correlation coefficients of ρ = −0.13 and a probability of p = 1.87 × 10 −2 in the W1-band, and ρ = −0.13 with a p-value of 2.61 × 10 −2 in the W2-band.
We also compare ∆R IR with the Eddington ratio in Figure 15 considering a fixed slope of α f ix = 0.5 in the R dust,IR − L bol relation as predicted by dust radiation equilibrium theory.This time, we find a moderate negative correlation between them, with ρ = −0.23,p = 2.36 × 10 −5 and ρ = −0.29,p = 3.66 × 10 −7 in W1 and W2-bands, respectively.
However, it is important to note that the Eddington ratio scales with L/R Hβ ∼ L/R dust,IR due to the linear correlation present between R Hβ and R dust,IR (Chen et al. 2023).Further, the change in ∆R IR varies as R dust,IR /L 0.5 .As a result, the negative correlation observed between ∆R IR and the Eddington ratio may be attributed to the self-correlation between the Eddington ratios and the deviation from the best-fit relation (Fonseca Alvarez et al. 2020).
Instead of Eddington ratio, we also utilized R [OIII] and R F eII , as indicators of Eddington ratio (Boroson & Green 1992;Sulentic et al. 2000;Ludwig et al. 2009;Marziani et al. 2003;Sun & Shen 2015;Du & Wang 2019).We collected Hβ, [O iii], and Fe ii luminosities from Liu et al. (2019) for Sample 1.For Sample 2, the Hβ, and [O iii] luminosities were retrieved from Shen et al. (2011), whereas Fe ii luminosities were obtained from Rakshit et al. (2020).Hence, we present ∆R IR as a function of R [OIII] and R F eII in Figure 14, and 15 for varying α and a fixed α f ix = 0.5, respectively.Notably, when α is treated as a free parameter, ∆R IR exhibits relatively weak correlations with R [OIII] and R F eII -positive and negative, respectively.However, these correlations become more pronounced when α is held constant at 0.5, as seen in Figures 14 and 15 for comparison.
In addition, we examine the dependence of ∆R IR with redshift z to investigate any residual redshift dependence.We find no significant dependence when α is treated as a free parameter.This suggests that the redshift dependence of dust lags is negligible after the z correction applied to our sample as detailed in Section 4.3.However, a noticeable correlation between ∆R IR and z emerges when we use α f ix = 0.5, indicating a reduction in torus sizes for AGNs at higher redshifts.However, this is mainly due to the luminosity effect since on average higher luminosity AGNs are selected at higher z.
Our findings are consistent with the previous studies by Chen et al. (2023).Notably, Yang et al. (2020) reported no dependency of ∆R W 1 on Eddington ratio.Note that they used Eddington ratio measured from single epoch M BH .In contrast, Chen et al. (2023) found that ∆R W 1 and ∆R W 2 systematically decrease with increasing accretion rate, suggesting the shortening of time lags for both BLRs and the dusty torus.They explained such behaviors by the self-shadowing effect of the slim disk model.Similarly, Fonseca Alvarez et al. (2020) found that AGNs have a shorter BLR size with lower R [OIII] .According to them, the deviation from the radius-luminosity relation depends on the distribution of the UV/optical SED and the relative contribution of ionizing radiation (for the case of Hβ BLR, see Woo et al. 2024).All of these findings align with our own results, implying that lag shortening of both BLRs and torus is caused by the same physical mechanism, i.e., the accretion rate.Additionally, the weakening of the overall correlation between ∆R IR and the Eddington ratio, R [OIII] , and R F eII as the slope becomes flatter suggests that the flattening of the R dust,IR -L bol relation from the dust radiation equilibrium model is influenced by the self-shadowing effect of the accretion disk as accretion rate increases.
However, the obtained R dust,IR -L AGN relations may suffer from several selection biases, such as, the IR lags may be underestimated for high luminosity sources because of the limited optical and IR baseline.Yang et al. (2020) mentioned that IR lags below 7 yr in the observed-frame can be detected within the base-line, whereas the lag detection probability will decrease toward longer lags in the observed-frame.
Utilizing high cadence light curve data in optical with strong features present in both optical and IR light curves, our analysis has allowed us to detect lags that are shorter than the cadence of the IR light curves (∼ 180 days) for a few targets.This lag detection was achieved through the application of both ICCF and JAVELIN techniques.It is worth noting that JAVELIN can detect lags shorter than the observing cadence (Li et al. 2019).However, it is important to consider that the IR lags might be prone to overestimation when the observedframe lags are significantly shorter than the cadence of the IR light curves, for instance, falling below 100 days.
To ensure the reliability of our lag measurements, we employed a comprehensive validation process in addition to assessing the quality of measured lags based on r max and p(r max ) τ >0 values, as described in Section 4.2.This validation procedure involved a systematic comparison of each lag measurement by shifting the optical light curves according to the calculated lag value and overlaying them onto the observed IR light curve.The results of this comparison demonstrated a consistent agreement between the shifted optical light curves and the observed IR light curves.As an illustrative example, please refer to Figure 2 for the target OB76 and Figure A1 for OB45 and OB60, all of which exhibit observed-frame W1 lags of less than 100 days.These figures provide clear visual evidence of the robustness of our lag measurements in such cases.The combined effect of these selection biases on the R dust,IR -L bol relation can be demonstrated properly by detail simulation, which is beyond the scope of this paper.

SUMMARY
We present DRM results for a total of ∼ 446 (416) AGNs with W1-(W2) band light curves, covering a large dynamic range in luminosities from 10 43.4 to 10 47.6 erg s −1 with redshift < 2 using 16-20 yr optical light curve data from different ground based optical surveys and IR data in W1 and W2-bands from WISE during 2010-2020.We explored the R dust,IR -L bol relations for different combinations of data.Our main findings are summarized below.
• We find that the torus size scales with the AGN bolometric luminosity as R dust,W 1 ∝ L 0.38 bol and R dust,W 2 ∝ L 0.35 bol in the W1 and W2-bands, respectively.Our best-fit slopes are much shallower than that expected from the thermal equilibrium model of R dust ∝ L 0.5 bol .However, our results are consistent with that obtained from the optical/IR interferometric observations, which also predict shallower slope than expected from the radiation equilibrium model.
• We remove the AD contamination from the observed IR (W1 and W2-bands) fluxes and obtained the R dust,IR -L bol relation.The correction of the AD contamination in IR fluxes results slightly different slope with R dust,W 1 ∝ L 0.39 bol and R dust,W 2 ∝ L 0.33 bol .As a result of this correction for AD contamination in the IR light curves, we observe larger torus sizes compared to cases where AD contamination is not accounted for in the analysis.
• Utilizing IR luminosity in the W2-band as the indicator for AGN luminosity, we identified analogous slopes in the R dust,IR -L W 2 relationships as those obtained for R dust,IR -L bol associations.
This observation implies that both unobscured and obscured AGNs exhibit a consistent trend within the torus size-luminosity plane.
• The torus sizes in W1 and W2-bands are found to be strongly correlated with R dust,W 2 larger than R dust,W 1 .We also compare our obtained IR R dust,IR -L bol relations in W1 and W2-bands with the R dust,K -L bol relation based on K-band lags available in the literature.We find a wavelength dependent torus size that increases with wavelengths as R dust,K : R dust,W 1 : R dust,W 2 = 1.0:1.5:1.8 (R dust ∝ λ 0.80 ) at L bol = 10 46 erg s −1 .This suggests a stratified dust torus.Additionally, a noteworthy observation is that objects with higher luminosities tend to display a more compact torus structure.
• The sizes of torus are consistently determined to be larger than that of the BLR, indicating that the dust torus exists outside of the BLR.Our analysis reveals a super-linear relationship between BLR size and torus size, as R Hβ ∝ R 1.20 dust,W 1 and R Hβ ∝ R 1.22 dust,W 2 , demonstrating a consistent slope observed in the R Hβ -L 5100 Å relationship available in the literature.
• We find that the scatter observed in the R dust,IR -L bol relation depends on Eddington ratio, while a mild correlation is evident with R [OIII] and R F eII .These findings collectively suggest that as the accretion rate increases, the torus size tends to shorten because of the self-shadowing effect of the slim disk, and the R dust,IR -L bol relation flattens compared to the predictions from the dust radiation equilibrium model.Nevertheless, we cannot definitively establish this correlation, except for the one arising due to the inherent relationship between AGN luminosity and Eddington ratio.We present examples of the light curves from a selection of 12 randomly chosen targets in Sample 1 and Sample 2, showcasing their wide range of bolometric luminosities (10 44.3 − 10 47.4 erg s −1 ).These examples, displayed in Figures A1 and A2, include the light curves in optical and IR W1 and W2-bands for Sample 1 and Sample 2, respectively.They serve as representative illustrations of our total sample.Furthermore, the lag measurements analyses between the optical and W1-band light curves for these targets, utilizing ICCF and JAVELIN, are depicted in Figure A3 and A4.

B. COMPARISON OF OUR LAG MEASUREMENTS WITH EARLIER WORK
In the context of our study, it is important to note that Lyu et al. (2019) previously reported lag measurements for 67 PG quasars.We identified 18 targets from our selected targets that overlap with the sample reported by Lyu et al. (2019).To assess the consistency between our new measurements and those obtained by Lyu et al. (2019) using their independent χ 2 −fit method, we compare the torus sizes in W1-band in Figure B1 (top left).With the exception of two targets, OB17 and OB32, our obtained torus sizes exhibit consistency (within 2σ) with those reported by Lyu et al. (2019), displaying an offset of 0.07 dex.
In Figure A1, we present the light curves of one of the most deviant targets, OB32, along with the corresponding results from the ICCF analysis shown in Figure A3.Notably, the optical light curve of OB32, shifted by our recovered time lag, aligns well with the W1-band light curve.For determining AGN bolometric luminosity, Lyu et al. (2019) employed optical-to-IR SED fitting and adopted an uncertainty of 0.3 dex in the luminosity.Figure B1 (top right) demonstrates the comparable bolometric luminosities we obtained for these 18 targets, in line with those derived from optical-to-IR SED fitting by Lyu et al. (2019).
In addition, from our analysis, we identified a subset of 12 targets that are common to both our sample and the sample reported by Chen et al. (2023).They used MICA (Li et al. 2016) in their time-lag measurements.Despite the robust agreement observed between the measurements obtained through our approach and those reported by Chen et al. (2023) (See the middle panels of Figure B1), our derived torus sizes consistently appear larger with an offset of 0.16 (0.14) dex for the W1-(W2) bands.This discrepancy arises from the fact that Chen et al. (2023) did not perform an AD contamination correction in their analysis, whereas our obtained torus sizes are corrected for AD contamination in the W1 and W2-band light curves.
Finally, we compare our measured torus sizes which are corrected for the AD contamination with those reported by Yang et al. (2020), who used JAVELIN without employing AD contamination correction for Sample 2. Despite most of our measurements lying above the 1:1 line due to the AD contamination correction, the two sets of measurements exhibit good consistency, showing an offset of 0.06 dex and a scatter of 0.15 dex.
Note that to ensure consistency in our comparison, we obtained the corresponding rest-frame time lags by multiplying the estimated observed-frame lags both from the literature and our own analysis, by the redshift correction factor (1+z) −0.38 .

C. COMPARISON BETWEEN EDDINGTON RATIOS DERIVED FROM SINGLE EPOCH AND REVERBERATION MAPPED M BH
We calculated the Eddington ratio for our sample by employing the single epoch M BH estimation method outlined in Section 4.4.In Figure C1, we present a comparative analysis of the Eddington ratios derived from the single epoch M BH estimates in our study (represented by red circles) and those obtained by Liu et al. (2019) (depicted by blue circles) with those obtained by using M BH measurements from BLR-RM.Since the black hole mass derived from RM is considered to be more accurate, Figure C1 clearly illustrates that the single epoch Eddington ratios tend to be underestimated as the accretion rate increases.This observation provides evidence for the need to account for the accurate measurement of M BH from BLR-RM to derive Eddington ratio.Baldwin, J., Ferland, G., Korista, K., & Verner, D. 1995, ApJL, 455, L119, doi: 10.1086/309827 Bao, D.-W., Brotherton, M. S., Du, P., et al. 2022, ApJS, 262, 14, doi: 10.3847/1538-4365/ac7beb Barvainis, R. 1987, ApJ, 320, 537, doi: 10.1086

Figure 1 .
Figure 1.Distribution of redshift (left) and bolometric luminosity (right) for different samples of AGNs with IR lag measurements.

Figure 2 .
Figure 2. From top left to bottom right: optical and IR W1 and W2-bands light curves of the targets (a) OB21, (b) OB44, (c) OB64, and (d) OB76 from Sample 1, as well as (e) ID24646, and (f) ID2497212 from Sample 2. The final calibrated optical light curves including data from CRTS, ASAS-SN, PTF and ZTF from Sample 1 are shown in the top panel, while for Sample 2 finally calibrated DES g-band optical light curves are presented.For Sample 1, the measurement errors (thick bars) and the total errors including the systematic uncertainty of the offset from PyCALI are presented.The bottom two panels represent the IR light curves in W1 and W2-bands.The W1 and W2-bands light curves shown in the figure are not corrected for the AD contamination.The shifted (by the lag from ICCF) and scaled version of the optical light curves (yellow circles) are also shown with the observed IR light curves, where the shade represents the 1σ uncertainty obtained in the scaling process.

Figure 3 .
Figure 3. From top left to bottom right: the ICCFs and JAVELIN lag analysis between optical and W1-band light curves (not corrected for AD contamination) for the targets (a) OB21, (b) OB44, (c) OB64, and (d) OB76, respectively, from Sample 1. Black solid line represents the ICCF between the optical and W1-band light curves and dashed green line is the ACF of the optical light curve.The CCCD obtained from ICCF is shown by the blue histogram, whereas the red histogram represents the lag probability distribution of the time delay from JAVELIN.The last two in the bottom panel are the ICCFs for the targets (e) ID24646, and (f) ID2497212 from Sample 2, where the vertical red dashed lines represent the JAVELIN lag measurements from Yang et al. (2020) with 1σ uncertainties shown by the red dotted lines.The observed-frame time lags obtained from ICCF with their uncertainties are also mentioned.
Figure 4.The probability p(rmax) represents the likelihood of uncorrelated light curves producing a crosscorrelation coefficient equal to or exceeding rmax.This probability is presented as a function of the maximum crosscorrelation coefficient rmax in the left and right panels for W1 and W2-bands, respectively.The vertical dashed line represents rmax=0.6,whereas p(rmax)=0.2 is shown by the horizontal dashed line.The reliable time lags are defined with rmax ≥ 0.6, and p(rmax) < 0.2.

Figure 5 .Figure 6 .
Figure 5. Examples of rejected targets based on ICCFs between the optical and W1-band light curves.The two targets ID7907462 (left) and ID954476 (right) from Sample 2 are the examples with negative lag and a flat CCF, respectively.The JAVELIN lag measurements from Yang et al. (2020) for these two targets are shown by the dashed red lines.The observed-frame time lags obtained from ICCF with their uncertainties are also mentioned at each panel.
Figure 7.The distribution of the lag-ratio (τratio = τW 2/τW 1) between the dust lags in W2 and W1-bands with and without AD contamination correction shown by red and black lines, respectively.The median values with their standard deviations are mentioned in the figure.

Figure 8 .
Figure 8. Comparisons of torus sizes with and without AD contamination correction are presented in W1 (left) and W2-(right) bands, color-coded by Eddington ratio.The AGNs which do not have Eddington ratio measurements are shown by the open grey circles.The best fits to all the data points are shown by the dashed brown lines, while the dashed black lines represent the 1:1 relation.The intrinsic scatter from the best fit is also mentioned at each panel.In the bottom panels, we represent the logarithm of the size ratio as a function of Eddington ratio.

Figure 9 .
Figure 9. Correlation between the torus size and AGN bolometric luminosity, measured in W1 (left) and W2-(right) bands after correcting for the AD contamination.The dashed blue line represents the best fit to the data points and the shaded blue region is the range driven by the uncertainties on the best fit.The intrinsic scatter is also shown at each panel in the figure.The dotted-dashed black colored line shows the best-fits for α f ix = 0.5.

Figure 10 .
Figure 10.Torus sizes in W1 and W2-bands after correcting for the AD contamination plotted against the AGN IR luminosity in W2-band in the left and right panels, respectively.The dashed blue line represents the best fit to the data points and the shaded blue region is the range driven by the uncertainties on the best fit.The intrinsic scatter is also shown at each panel in the figure.The dotted-dashed black colored line shows the best-fits for α f ix = 0.5.

Figure 11 .Figure 12 .
Figure 11.Comparison between the torus sizes in W1 and W2-bands obtained from ICCF in the rest-frame of the targets based on their AD contamination corrected IR light curves.The best-fit to the data points is shown by the dashed brown line, while the dashed black line represents the 1:1 relation.

Figure 13 .
Figure 13.Left: R dust -L bol relations based on W1, W2 and K-band lags, which are corrected for the AD contamination.Theoretical estimates are calculated with Equations 11 and 12, respectively, for silicate with a sublimation temperature T sub = 1500 K (black dashed line) and for graphite grains with T sub = 1800 K (black dotted-dashed line).The Hβ BLR size -luminosity relation adopted fromWoo et al. (2024) is also presented.Right: Torus size measured at K, W1, and W2-band for three different values of L bol = 10 44 , 10 46 and 10 48 erg s −1 , and the best-fit relation for each luminosity case.The error bars represent the intrinsic scatter in the corresponding torus size -luminosity relation.

Figure 14 .Figure 15 .
Figure 14.The deviations of the dust-torus sizes from the corresponding best-fitting relationships (∆RW 1, and ∆RW 2) as a function of Eddington ratio, R [OIII] , RF eII , and redshift z.The black points represent DRM sample corrected for AD contamination in the IR W1, W2 light curves.The Spearman's rank correlation coefficients and the corresponding probabilities are also mentioned in each panel.The blue dashed lines represent ∆RIR = 0.
been supported by the Basic Science Research Program through the National Research Foundation of the Korean Government (grant No. NRF-2021R1A2C3008486) and Samsung Science and Technology Foundation under Project Number SSTF -BA1501-05.S.W. acknowledges the support from the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2019R1A6A1A10073437).

Figure A2 .
Figure A2.Examples of light curves from Sample 2 in optical and IR W1-and W2-bands for the targets (a) ID47417, (b) ID1124174, (c) ID1189298, (d) ID1524735, (e) ID2023374, and (f) ID4065339.The IR W1-and W2-band light curves represented in the plot are not corrected for the AD contamination.Same notations are used as Figure 2.

Figure C1 .
Figure C1.Comparison between the Eddington ratios obtained from single epoch MBH and MBH from BLR-RM.The red and blue circles represent the objects with Eddington ratio single epoch calculated in this work (see Section 4.4) and retrieved from Liu et al. (2019), respectively, with the linear-fits shown by the dotted lines.The black dashed line is the 1:1 relation.

Table 1 .
Details of the surveys used in this work

Table 3 .
Optical and IR photometric light curves.The full table is available for Sample 1 and Sample 2 separately in a machine-readable form online.

Table 5 .
Linear regression fitting results of the torus size -luminosity relation.
Note: Columns are (1) luminosity used in fitting, (2) different cases, NON ADC: without AD contamination correction, ADC: correction for AD contamination applied, (

Table 6 .
Details of BLR-RM sample Note: Columns are (1) object name, (2) object ID number used in this work, (3) Hβ BLR size (rest-frame) obtained from BLR-RM, (4) host-corrected AGN optical luminosity at 5100 Å collected from the literature with BLR-RM measurements, (5) MBH obtained from RM, (6) Eddington ratio in log scale derived from reverberation-mapped MBH , and (7) references.L Kim et al. 2021017)can play an important role.The LSST+NGRST data set can be used to measure the size of the dust torus in local AGNs with redshift z < 0.5 since NGRST can observe up to ∼ 2 µm(Yang et al. 2020).In order to employ dust lag as a standard candle for cosmology and constrain cosmological parameters,Hönig et al. (2017)started a large DRM program, namely 'VEILS' (VISTA Extragalactic Infrared Legacy Survey) that will observe about 1350 Seyfert 1 galaxies in the redshift range of 0.1 < z < 1.2.The Spectro-Photometer for the History of the Universe, Epoch of Reionization and Ices Explorer (SPHEREx;Kim et al. 2021) mission aims to provide optical and NIR multiepoch spectroscopic data, which can be utilized in RM studies for bright AGNs.These RM studies can generate a unique data set that, when combined with complementary optical observations, can aid in comprehending the physical properties of the dust torus as well as the central structures of AGNs.ASAS-SN is supported by the Gordon and Betty Moore Foundation through grant GBMF5490 to the Ohio State University and NSF grant AST-1515927.The Development of ASAS-SN was supported by the NSF grant AST-0908816, the Mt.Cuba Astronomical Foundation, the Center for Cosmology and AstroParticle Physics at the Ohio State University, the Chinese Academy of Science South America Center for Astronomy (CASSACA), the Villum Foundation, and George Skestos.ZTF is funded by the National Science Foundation through grant number AST-1440341 and a partnership with Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches