The Structure Function of Mid-infrared Variability in Low-redshift Active Galactic Nuclei

Using the multi-epoch mid-infrared (MIR) photometry from the Wide-field Infrared Survey Explorer spanning a baseline of ∼10 yr, we extensively investigate the MIR variability of nearby active galactic nuclei (AGNs) at 0.15 < z < 0.4. We find that the ensemble structure function in the W1 band (3.4 μm) can be modeled with a broken power law. Type 1 AGNs tend to exhibit larger variability amplitudes than type 2 AGNs, possibly due to the extinction by the torus. The variability amplitude is inversely correlated with the AGN luminosity, consistent with a similar relation known in the optical. Meanwhile, the slope of the power law increases with AGN luminosity. This trend can be attributed to the fact that the inner radius of the torus is proportional to the AGN luminosity, as expected from the size−luminosity relation of the torus. Interestingly, low-luminosity type 2 AGNs, unlike low-luminosity type 1 AGNs, tend to exhibit smaller variability amplitude than do high-luminosity AGNs. We argue that either low-luminosity type 2 AGNs have distinctive central structures due to their low luminosity or their MIR brightness is contaminated by emission from the cold dust in the host galaxy. Our findings suggest that the AGN unification scheme may need to be revised. We find that the variability amplitude of dust-deficient AGNs is systematically larger than that of normal AGNs, supporting the notion that the hot and warm dust in dust-deficient AGNs may be destroyed and reformed according to the strength of the ultraviolet radiation from the accretion disk.


INTRODUCTION
Active galactic nuclei (AGNs) radiate strong multiwavelength continuum emission originating from complex structures around supermassive black holes (BHs), such as the accretion disk (AD), corona, jet, and dusty torus.One of the unique features of AGN emission is its variability over a wide range of time scales, indicating that the emission arises over a large range of physical scales (e.g., Matthews & Sandage 1963).Therefore, the variability can be used to probe the physical properties of central structures in AGNs, which, apart from interferometric observations (e.g., Swain et al. 2003;GRAVITY Collaboration et al. 2020a,b), are barely resolved with conventional imaging data.Variability has also been widely used for AGN selection (e.g., van den Bergh et al. 1973;Choi et al. 2014;Burke et al. 2023).While the physical origin of the variability remains under debate (e.g., Lyubarskii 1997;Ulrich et al. 1997;Kawaguchi et al. 1998;Livio et al. 2003;Dexter & Agol 2011;Kubota & Done 2018;Sun et al. 2020), the X-ray and UV/optical variability can be described by a stochastic process (e.g., Peterson 1997;Kelly et al. 2011).More specifically, the UV/optical continuum is well modeled by a damped random walk, for which the shape of the power spectral density is a broken power law (e.g., Kelly et al. 2009;MacLeod et al. 2010;Tang et al. 2023; but see Kasliwal et al. 2015).
The power spectral density of the light curves in the UV/optical band can be characterized by the variability amplitude, a power-law index of approximately −2.0 at high frequencies, and a characteristic frequency below which the power spectral density flattens.Various studies argued that these parameters may be related to the physical properties of AGNs (e.g., wavelength, BH mass, Eddington ratio, and AGN luminosity; Vanden Son et al. Berk et al. 2004;Li et al. 2018;Suberlak et al. 2021;Tang et al. 2023).However, different studies reached different conclusions regarding the dependence of the variability properties on the AGN properties, possibly attributed to the bias present in constructing the power spectral density introduced by sample selection, insufficient cadence, length of the time-series data, and method used to construct the power spectral density (Kozlowski 2017;Suberlak et al. 2021).For example, Li et al. (2018) demonstrated, using ground-based survey data and ensemble analysis, that the variability amplitude is inversely correlated with the AGN luminosity but is independent of the BH mass (see also Kelly et al. 2009).However, MacLeod et al. (2010) argued that the variability amplitude is correlated with the BH mass and inversely correlated with the Eddington ratio.In addition, the power-law slope at high frequencies is reported to be a constant (approximately −2) using datasets obtained from ground-based telescopes that may suffer from a sparse sampling rate.On the contrary, studies with high-cadence light curves obtained from the Kepler mission found that the power-law slope is not a constant and is correlated with the AGN properties (e.g., Mushotzky et al. 2011;Smith et al. 2018).
While the variability of the continuum from the innermost regions of AGNs has been extensively studied with multi-epoch X-ray, UV/optical, and radio observations, the variability in the mid-infrared (MIR) has yet to be investigated with a large dataset (e.g., Kozlowski et al. 2016;Son et al. 2022a;Li & Shen 2023).One of the obstacles toward realizing this goal is the lack of MIR time series owing to the relative difficulty of obtaining such data from ground-based telescopes.As the MIR continuum comprises AD emission reprocessed by the dust in the torus, it carries information not only on the intrinsic characteristics of the light from the AD but also helps us to constrain the structural properties of the torus (e.g., Kawaguchi & Mori 2011;Son et al. 2022a;Li & Shen 2023).In this regard, it is instructive to indirectly probe the continuum from the AD of type 2 sources to test the unification model of AGNs.Despite the importance of MIR variability in understanding the structure of the dusty torus, there have been few quantitative examinations of this variability (e.g., Kozlowski et al. 2010Kozlowski et al. , 2016;;Sánchez et al. 2017;Wang & Shi 2020).
Owing to the continuous survey conducted by the Wide-field Infrared Survey Explorer (WISE), ∼ 10 years of all-sky data are available, which permits systematic variability studies of AGNs to be conducted, albeit with a sparse cadence (∼ 6 months; Wright et al. 2010).The WISE database, in combination with complementary optical data, can be used to advance a variety of studies on AGN variability, including reverberation mapping, searching for changing-look AGNs or tidal disruption events, and investigating the physical properties of the hot dust components (e.g., Lyu et al. 2019;Sheng et al. 2020;Yang et al. 2020;Jiang et al. 2021;Son et al. 2022b;Li & Shen 2023).Indeed, these data are also ideal for studying the properties of the MIR variability for a large sample of AGNs.
To explore the characteristics of the MIR variability of nearby AGNs, we investigate the variability structure function (SF) of multi-epoch MIR data for a large sample of nearby AGNs obtained with WISE spanning a baseline of ∼ 10 yr, paying careful attention to treatment of noise estimation and host galaxy subtraction.In Section 2, we describe our sample selection and the method used to construct the MIR light curves.Section 3 introduces detailed methods for computing robust ensemble SFs.The SFs of various subsamples are presented in Section 4. The physical origins of the SF variations are discussed in Section 5. Finally, the results are summarized in Section 6.Throughout the paper, all magnitudes refer to the Vega system, and we adopt the cosmological parameters H 0 = 100h = 67.4km s −1 Mpc −1 , Ω m = 0.315, and Ω Λ = 0.685 (Planck Collaboration et al. 2020).

SAMPLE AND DATA
For type 1 AGNs, the parent sample is drawn from the Data Release 14 quasar catalog of the Sloan Digital Sky Survey (SDSS; Pâris et al. 2018).Because the light contribution from the host galaxy is significant in the W1 band (3.4 µm), careful subtraction of the host component is necessary (Son et al. 2022a(Son et al. , 2023)).We constrain the host magnitude by fitting the spectral energy distribution (SED).We impose a low-redshift cut (z > 0.15) to mitigate significantly extended sources, whose photometric data among the different datasets can be inconsistent due to mismatches in the photometric method and spatial resolution.At z ≈ 0.15, the point-spread function (PSF) full-width at half maximum (FWHM) of 6 ′′ in the W1 band corresponds to ∼ 15 kpc, which ensures that most of the galaxy flux is included in the WISE photometry.To facilitate more effective estimation of the host galaxy contribution using restframe near-infrared data, we only choose type 1 AGNs with 2MASS (Skrutskie et al. 2006) counterparts with z < 0.4, which additionally ensures that the Hα region is covered by the SDSS spectra for AGN classification.These redshift cuts are also crucial for minimizing the effect of cosmic evolution on the dust properties of the torus and for neglecting possible dependences on rest- frame wavelength (e.g., Kawaguchi & Mori 2011).A total of 4295 type 1 AGNs are initially chosen.We use the same redshift cut (0.15 < z < 0.4) to select 6854 type 2 AGNs from the SDSS Data Release 8, based on spectral classification from the Baldwin, Philips & Terlevich diagram (i.e., "bptclass"=4;Baldwin et al. 1981;Veilleux & Osterbrock 1987) provided by the Max Planck Institute for Astrophysics and the Johns Hopkins University (MPA-JHU) value-added catalog 1 (Brinchmann et al. 2004).We exclude low-ionization nuclear emission-line regions (LINERs; Heckman 1980), whose extremely low luminosities and Eddington ratios (Ho 2008;Ho & Kim 2009) preclude accurate photometry of the nucleus based on the relatively low-resolution WISE images.For type 2 AGNs, matching with 2MASS is unnecessary because the total stellar masses from the MPA-JHU catalog are based on optical photometric data Kauffmann et al. (2003).
To investigate the SF of SDSS AGNs in the MIR, we use the multi-epoch data provided by the Near-Earth Object Wide-field Infrared Survey Explorer (NE-OWISE), which spans a baseline of ∼ 10 yr from December 23, 2013to December 13, 2022(Mainzer et al. 2011).In view of the typical positional uncertainty of the WISE dataset (∼ 0. ′′ 5; see also Assef et al. 2013; Son 1 https://www.sdss3.org/dr10/spectro/galaxy_mpajhu.phpet al. 2022a)2 , we adopt a radius of 2 ′′ to cross-match between SDSS and NEOWISE.We only use time-series data from the W1 band because its signal-to-noise ratio is substantially higher than that from the W2 band (4.6 µm).To secure a robust analysis of the SF, we only consider samples with more than 17 epochs in the W1 band.Finally, to exclude non-variable objects, which are not appropriate for further analysis, we employ variability probability criterion of P var > 0.95 (McLaughlin et al. 1996;Sánchez et al. 2017), where P var = 1 − Q and Q, the probability that the observed light curve originates from a non-variable object, is estimated using the χ 2 distribution at a given degree of freedom (N − 1, with N the number of observing epochs).A high P var (low Q) indicates that an object is likely variable.The final sample contains 3506 type 1 and 3074 type 2 AGNs.
As NEOWISE is designed to survey the entire sky over a six-month interval, any given field is observed ∼ 13−14 times on average over a few days during each visit.We retrieve the photometric measurements from the single exposures of each visit and compute a mean value with a 3σ clipping.To discard suspicious measurements, we only use photometric data flagged with cc_f lags = 0, qual_f rame > 0, qi_f act > 0, saa_sep > 0, and moon_masked = 0.The method from Lyu et al. (2019) ∆t is adopted initially to calculate the uncertainty of the W1 magnitude: (1) where N s is the number of single exposures in each epoch, m i and σ i,pho represent the magnitude and uncertainty of a single exposure, m epoch is the mean magnitude in each epoch, and σ sys is the systematic uncertainty (∼ 0.016 mag) due to the instability of the system.However, from extensive testing of the SF from non-variable sources, we find that the initial uncertainty can be substantially overestimated (see Section 3.3.2),and the final uncertainty is estimated better by dividing σ in by √ N s (σ e = σ in / √ N s ).We extract various AGN properties (Figure 1) derived from the spectral measurements of Rakshit et al. (2020) to explore the physical connection between the SF and properties of the AGN.The BH mass (M BH ) for type 1 AGNs is derived using the virial method: M BH = f Rv 2 /G, where R is the size of the broadline region inferred from the continuum luminosity at 5100 Å (L 5100 ) using the size−luminosity relation (Kaspi et al. 2000;Bentz et al. 2013, v the velocity width of broad Hβ emission, and f is a scaling factor determined by the geometry and kinematics of the line-emitting region.We adopt the BH mass estimator from Ho & Kim (2015) appropriate for both bulge types (see Ho & Kim 2014 on the systematic difference of f for classical and pseudo bulges): M BH = 10 6.91 (FWHM/1000 km s −1 ) 2 (L 5100 /10 44 erg s −1 ) 0.533 M ⊙ .For type 2 AGNs, we utilize the empirical relation between BH mass and the total stellar mass of the host (M * ), as calibrated by Greene et al. (2020) for all (early and late) galaxy types: M BH = 10 7.43 (M * /3 × 10 10 M ⊙ ) 1.61 M ⊙ .The stellar masses, extracted from the MPA-JHU catalog, are derived by fitting the spectral energy distribution from the SDSS photometry with the stellar population model.The BH mass measurements have uncertainties of ∼ 0.5 dex and ∼ 0.8 dex for type 1 and type 2 AGNs, respectively.We use the [O III] luminosity to trace the strength of both AGNs types, and, as an additional check, L 5100 for type 1 sources alone.Converting the [O III] luminosity to the bolometric luminosity requires consideration of possible correction for dust extinction (see Kong & Ho 2018, for an extensive discussion of the bolometric conversion).For our sample of type 1 AGNs, L 5100 correlates more strongly with the observed [O III] emission than with the line luminosity after extinction correction based on the Balmer decrement.We therefore choose the bolometric correction of Heckman et al. (2004) Ensemble SF of non-variable sources.The dashed line denotes the mean error of the sources.The blue circles and red stars denote the mean of the SF and the square root of the mean SF, respectively.

Definition
AGN variability can be modeled with the damped random walk model, wherein the variability amplitude of the light curve as a function of a time lag obeys a power law on short time scales and flattens on longer time scales (e.g., Kelly et al. 2009).Alternatively, the variability can be modeled non-parametrically with the SF, which is defined as the mean magnitude difference (SF) as a function of time lag (∆t; e.g., Simonetti et al. 1985;Kawaguchi et al. 1998;Kozlowski 2016): where N ∆t,pair denotes the number of pairs associated with ∆t, m is the observed magnitude, and σ e represents the uncertainty in each epoch.As expected from the damped random walk model, the SF in the optical band is divided into two power laws, SF ∝ ∆t γ at a break time scale t break , such that: γ ≈ 1 for ∆t ≤ t break and γ ≈ 0 for ∆t < t break .
Ideally, σ e can be estimated from the standard deviation of the multi-epoch photometric data [i.e., SF(0)] for the non-variable objects and is known to be independent of ∆t.Alternatively, the dependence on ∆t can be computed from the SF derived from the light curves of non-variable sources as SF 2 non−variable (∆t) = σ 2 e (t) + σ 2 e (t + ∆t) (Kozlowski et al. 2016).As de-scribed in Section 3.3.2,we evaluate the ensemble SF using the WISE light curves from non-variable sources and find that the SF depends on ∆t as SF(∆t) = SF(0)(0.999+ 0.055 × ∆t/1000).We use this relation to estimate the SF.Note that the observed ∆t obs is converted to the rest frame, namely ∆t = ∆t obs /(1 + z).
Owing to the sparse cadence of the NEOWISE data (∼ 6 months), it may not be ideal to use the SF from individual objects to explore the statistical properties of the MIR variability.Instead, we utilize an ensemble SF, which is computed by averaging the SFs for a given ∆t obtained from subsamples with similar AGN properties (e.g., AGN type, AGN luminosity, BH mass, and dust properties; Almaini et al. 2000;Vanden Berk et al. 2004;Li et al. 2018).

Host Subtraction
The contribution of the host galaxy can be significant in the W1 band and should be carefully removed to yield a robust estimate of the AGN luminosity.The SF, therefore, is highly sensitive to host subtraction.For example, over-subtraction of the host light can lead to an overestimation of the SF.For type 1 AGNs, we perform a fit to the SED spanning from the optical to the MIR, using integrated photometry derived from SDSS, 2MASS, and WISE, as detailed in Son et al. (2023).We consider two components, one for the host and the other for the AGN.The AGN component is represented by three SEDs-hot dust-deficient (HDD), warm dustdeficient (WDD), and normal AGNs-empirically determined based on the presence of hot dust near 1 − 3 µm and warm dust around 3−10 µm.While the HDD AGNs are defined as lacking both hot and warm components in their SED, the SED of normal AGNs is well fit with the template SED of bright QSOs from Elvis et al. (1994).The host component is modeled with seven templates: an old stellar population with 7 Gyr old stars and empirical templates of inactive galaxies (Hubble type E, S0, Sa, Sb, Sc, and Sd) from the Spitzer Wide-Area Infrared Extragalactic Survey library (Polletta et al. 2007).
For type 2 AGNs, we find that SED fitting does not provide satisfactory results, possibly because the extended morphology of these sources may introduce systematic offsets in the photometry among the different datasets.Li et al. (2023) demonstrate, for instance, that the photometry provided by the 2MASS and WISE catalogs can be systematically underestimated if the sources are extended.Therefore, we instead estimate the luminosity expected in the W1 band based on the stellar mass derived from the SDSS photometry, in combination with the mass-to-light ratio in the W1 band adopted from Kettlety et al. (2018).The average flux fraction of the host in the W1 band is 0.20 ± 0.12 for type 1 AGNs and 0.47 ± 0.24 for type 2 AGNs.Because of their lower bolometric luminosity, type 2 AGNs have significantly more dominant hosts than type 1 AGNs.

Systematic Uncertainties
As the sparse cadence, instrumental characteristics, and photometric errors can introduce systematic uncertainties into the SF, we perform extensive simulations to evaluate them.

Robust Structure Function Measurements
There are two ways to estimate an ensemble SF: (1) the mean of √ SF 2 (SF) or (2) the square root of the mean of SF 2 ( SF 2 ).Ideally, the two methods would yield similar results.However, the sparse cadence of the WISE dataset may introduce uncertainties in the ensemble SFs.We perform simulations with artificially generated light curves to test which method better reproduces the true SF.For simplicity, we assume that the light curves obey the damped random walk model with a break time scale of 650 days.The variability amplitude is assumed to be 10.6% of the mean flux (SF ∞ = 0.15 mag), comparable to our sample.In total, 10 4 light curves are generated randomly, considering the actual cadences and redshift distributions of our sample.We subsequently estimate the ensemble SF using the mock light curves.The photometric error is not included in this initial simulation.This experiment reveals that SF 2 successfully reproduces the input SF, whereas the SF tends to be underestimated by SF (Figure 2).More interestingly, for large ∆t (≳ 1000 days), the discrepancy between the input SF and SF becomes as severe as 20%.This clearly indicates that SF 2 needs to be used to trace the true SF of our sample.To account for the effect of photometric noise, we perform the same simulation by adding to the light curves a photometric noise of 0.017 mag, which is similar to the typical error in our sample.We again find that SF 2 is more suitable than SF, but it deviates from the input value for ∆t ≳ 2500 days.Unless otherwise noted, SF refers to SF 2 throughout this paper.

Photometric Errors
As a particular field is visited every ∼ 6 months by NEOWISE, the same field is observed multiple times with 14 ± 7 single exposures in each visit.The photometry measurements are performed in a single exposure and combined into the representative magnitude in each epoch.Therefore, estimating the photometric noise (σ e ) in this calculation is not straightforward.One way to compute the photometric noise is to utilize the SF of non-variable sources (e.g., Kozlowski et al. 2010).For this purpose, we generate 10 4 light curves with a constant flux associated with constant photometric noise (σ p ), following the distributions of the cadences and redshifts of our sample.From the ensemble SF taken from these light curves, we find SF ≈ √ 2 × σ p , with SF independent of ∆t.The √ 2 comes from the fact that σ 2 p is subtracted twice, from t and t + ∆t (see Equation 2).
To mimic realistic observations, we perform the simulation using the NEOWISE multi-epoch data from nonvariable sources.For this purpose, we randomly choose 10 4 galaxies drawn from the set of inactive galaxies with neither AGN nor SF activity in the optical spectra contained in MPA-JHU catalogs by imposing the same redshift cut (0.15 < z < 0.4) as the AGN sample.As described in Section 2, we initially calculate the photometric noise (σ in ) in each epoch using Equation (1).From this simulation, we find that the photometric uncertainty based on σ in is significantly larger than the standard deviation of the magnitudes (σ rms ) in the light curves of non-variable sources.Our experiment yields σ rms ≈ σ in / √ N s , where N s denotes the number of single exposures in each epoch.From the SF measurements, we independently find that SF can be approximated by √ 2σ in / √ N s (Figure 3), although it depends weakly on ∆t, such that SF(∆t) ≈ √ 2σ in / √ N s × (0.999 + 0.055∆t/1000.We adopt this relation when estimating the SF of our sample, noting that the dependence on ∆t barely affects the SF measurements. We note that the MIR continuum of inactive galaxies may be variable because of several possible origins, including supernova explosions, weak nuclear activity, or tidal disruption events, which can naturally introduce bias in the error estimation, even if their occurrence rates are relatively low (e.g., Jiang et al. 2021;Son et al. 2022b;Sun et al. 2022).If the inactive galaxies vary significantly in the MIR, our error estimation should be taken as an upper limit.To quantify this, we assume that the photometric noise can be described as σ in / √ N s and rigorously estimate the ensemble SF of the inactive targets by subtracting this photometric noise using Equation (2).From this experiment, we found that the SF is ∼ 0.02 mag regardless of ∆t, which shows that the intrinsic variability of the inactive galaxies is almost negligible.Alternatively, if the inactive galaxies significantly vary in the MIR, the photometric errors are systematically smaller than our estimation (i.e., σ < σ in / √ N s ).Even in this case, the photometric noise is substantially smaller than the SFs in our sample AGNs; therefore, its impact on the SF measurements is minimal.

Uncertainties of the Ensemble Structure Function
Finally, the uncertainty of the SF for an individual object is estimated from the standard deviation of each ∆t bin, while that of the ensemble SF is estimated by bootstrap resampling of the SFs for the individual galaxies.We create 10 3 realizations of the SF for each object, and the final uncertainty is taken to be half of the difference between the 14th and 86th percentiles of the ensemble SF in each bin.

Accretion Disk Contribution
Although the W1 flux is generally dominated by the reprocessed emission from the hot and warm dust in the torus, the flux contribution from the AD may not be negligible for type 1 AGNs (e.g., Asmus et al. 2011).In addition, variability in the continuum emission from the AD can bias the SF.To properly remove the effect of the AD in the SF measurements, it would be ideal to subtract the contribution from the AD inferred from the corresponding optical photometric data in each epoch.Unfortunately, this is unavailable for our sample.Instead, we compute the expected W1 flux inferred from the SDSS g-band flux by assuming F ν ∝ ν 1/3 , as predicted from theory and verified by polarization measurements of quasars (Vanden Berk et al. 2001;Kishimoto et al. 2008;Li & Shen 2023).Even after the host sub- traction, the average flux contribution from the AD in the W1 band is 0.09 ± 0.05 (Figure 4).To estimate the expected SF at the W1 band solely due to the AD, we adopt the g-band SF of AGNs from (Li & Shen 2023) that have a median L bol ≈ 10 45.3 erg s −1 at z < 0.7, similar to that of our sample.This SF is modeled with a variability amplitude SF ∞ = 0.544 mag, break time τ = 853.1 days, and a power-law index β ≈ 2γ = 1 for small values of ∆t.We also account for the dependence of SF on wavelength, SF ∝ λ −0.30±0.05(Ivezic et al. 2004), when converting the SF from g to W1.We conservatively adopt an AD flux fraction of 0.17, which is the 95th percentile value of the distribution of the AD fraction in the W1 band for our sample.The expected SF for ∆t = 1 yr is ∼ 0.032, which indicates that the SF due to the AD is almost negligible.Therefore, we neglect this effect throughout the study.

RESULTS
We estimate the ensemble SFs for various subsamples according to AGN type, AGN property, and dust property.Previous studies have often modeled ensemble SFs of sparsely sampled light curves with a single power-law function (e.g., Kozlowski 2016).However, such a function is insufficient for describing the SFs in our sample.Instead, we fit the ensemble SF with a broken power-law function  Ensemble SF for subsamples binned by the continuum luminosity at 5100 Å, in units of erg s −1 .The dashed lines denote the fits with a broken power law.
(3) where ∆t is the time lag in days, A is the variability amplitude at 1 yr, t break is the break time scale, γ 1 is the power-law slope for ∆t ≤ t break , and γ 2 is the powerlaw slope at ∆t > t break .For direct comparison between subsamples, we primarily consider the variability amplitude at 1 yr (A) and the slope for the lag at low ∆t (γ 1 ).We estimate the 1σ uncertainties of the fitted parameters through bootstrapping the ensemble SFs with their 1σ errors.On account of the sparse cadence of the WISE light curves, t break is found to be sensitive to the initial conditions in the fitting of the ensemble SF, although its uncertainty appears to be small.Therefore, we decide not to overinterpret the trend of t break in this study.The fitted parameters of the various subsamples are summarized in Table 1.

Type 1 vs. Type 2 AGNs
According to the AGN unification model, type 1 and type 2 AGNs are determined by the viewing angle with respect to the opening angle of the dusty torus (Antonucci 1993; Urry & Padovani 1995).We use the MIR SFs of type 1 and type 2 AGNs to test this model, which cannot be performed in the optical because the continuum emission from the AD is heavily obscured in type 2 AGNs.We find that the ensemble SFs of type 1 and type 2 AGNs are well fitted with a broken power-law model and are consistent in terms of their overall shape (γ), although type 2 AGNs exhibit a slightly shallower power law than their type 1 counterparts (Figure 5).On the other hand, the variability amplitude of type 1 AGNs is systematically larger than that of type 2 sources.These trends remain the same when the sample is divided into subsamples based on AGN properties.Interestingly, the similar power-law indices of type 1 and type 2 AGNs conflict with the results derived from optical SFs, for which the ensemble SFs of type 2 AGNs are significantly flatter than those of type 1 AGNs, possibly owing to obscuration (De Cicco et al. 2022).Finally, our values of γ 1 (0.47 − 0.51) broadly agree with γ ≈ 0.45 derived from the MIR (3.6 and 4.5 µm) light curves of high-z quasars studied by Kozlowski et al. (2016).

Dependence on AGN Properties
The physical connection between the variability properties and AGN properties has been extensively studied with optical light curves.We revisit this topic using MIR light curves.One of the most striking features in our sample is the dependence of the SF on the AGN luminosity.Using the monochromatic luminosity at 5100 Å as a tracer of the bolometric luminosity for type 1 AGNs, Figure 6 clearly shows that low-luminosity AGNs tend to exhibit a larger variability amplitude (A) and a shallower slope (γ 1 ) than high-luminosity AGNs.The inverse correlation between A and the AGN luminosity is consistent with that found for the optical SF (e.g., Kelly et al. 2009;Caplar et al. 2017;Sun et al. 2018;Tang et al. 2023), which suggests that MIR variability may be directly related to the intrinsic variability of the optical light.The dependence of the slope on the luminosity is further discussed in Section 5.1.
To directly compare type 1 and type 2 AGNs, we adopt the [O III] luminosity as a proxy for the AGN luminosity.Although the overall trends remain the same, the dependence on the AGN luminosity as inferred from L [O III] is slightly weaker compared to that derived using L 5100 (Figure 7).This can be attributed to the less direct link between [O III] luminosity and the bolometric luminosity.Intriguingly, the ensemble SFs of type 2 AGNs of low-luminosity, defined here as L [O III] ≤ 10 41.3 erg s −1 , significantly deviates from those of their type 1 counterparts, in the sense that the SF of low-luminosity type 2 AGNs tend to have lower amplitude and shallower slope.
We also examine how the ensemble SFs change with the BH mass and Eddington ratio (Figures 8 and 9).Although our results do not show any clear trend, γ 1 increases with increasing BH mass in type 1 AGNs.How-  ever, this trend is not clearly detected in type 2 AGNs, possibly owing to the relatively large uncertainties in the BH mass estimation.In addition, A is weakly correlated with the Eddington ratio only in type 1 AGNs, presumably reflecting the dependence of the SF on AGN luminosity.

Dependence on Dust Properties
Type 1 AGNs exhibit diverse MIR SEDs depending on the strength of the emission from their hot and warm dust components, which radiate mostly in the 1 − 3 µm and 3 − 10 µm regions, respectively.AGNs can be classified into three categories according to their SED  Ensemble SF for HDD (red circles), WDD (green squares), and normal (blue diamonds) AGNs.The fits with a broken power law are denoted by corresponding dashed lines.
shapes: normal, HDD, and WDD.We classify our sample of type 1 AGNs through SED fitting with empirically determined template spectra of the nearby quasars investigated by Lyu et al. (2017), who identified HDD and WDD AGNs as those whose MIR SEDs deviate from that of normal AGNs by more than ≥ 0.3 dex in the wavelength range ∼ 1 − 3 µm and ∼ 3 − 10 µm, re-spectively.The fiducial reference SED for normal AGNs were taken from Elvis et al. (1994).As described in detail in Son et al. (2023), our fitting uses photometric data covering the ugriz bands from SDSS, the JHK S bands from 2MASS, and the W1, W2, W3, and W4 bands from WISE.Using the Bayesian Information Criterion to determine the goodness of fit, we obtain best-fit models for 449 normal, 819 WDD, and 451 HDD AGNs.
To test whether the dust properties are physically related to the MIR variability, we estimate the ensemble SFs for normal, WDD, and HDD AGNs. Figure 10 clearly demonstrates that the WDD and HDD AGNs exhibit larger variability amplitudes and steeper powerlaw slopes than normal AGNs.One caveat is that these findings could be a by-product of the dependence on AGN properties, as the dust properties are known to be closely connected with AGN parameters such as bolometric luminosity, Eddington ratio, and BH mass (e.g., Hao et al. 2010;Jiang et al. 2010;Jun & Im 2013;Lyu et al. 2017;Son et al. 2023).To mitigate against this effect, we reestimate the ensemble SFs at a fixed, narrow range of [O III] luminosity and BH mass, but the dependence on dust properties still persists (Figure 11).We conclude that dust-deficient AGNs vary more drastically in the MIR than normal AGNs.which is heated by the UV/optical photons from the AD, two separate components determine the MIR variability.On the one hand, the intrinsic variability of the flux from the AD should be reflected in the MIR variability.On the other hand, the MIR light curves are heavily smoothed compared to the optical light curves because the variability time scale associated with the torus is larger than that of the UV/optical continuum (e.g., Suganuma et al. 2006;Hönig & Kishimoto 2011;Koshida et al. 2014).Consequently, the MIR SF can also be sensitive to the structure and physical properties of the torus.Here, we discuss the main drivers of the MIR variability based on the ensemble SFs of various subsamples.
The most striking finding is the dependence of A and γ 1 on the AGN luminosity.In particular, more luminous AGNs have systematically lower variability amplitudes (Figure 12), which agrees well with trends known in the optical (e.g., Vanden Berk et al. 2004;Kelly et al. 2009;Tang et al. 2023; but see Sánchez-Sáez et al. 2018).This may imply that there is a direct connection between the optical and MIR variability.By contrast, the dependence of the slope γ 1 for low ∆t on the AGN luminosity has not been identified in most previous studies involving optical light curves, suggesting that this trend may be attributed to the geometry of the torus.
Interestingly, Li & Shen (2023) demonstrated that the ensemble SFs of the MIR flux, in combination with those of the optical flux, can be used to estimate the inner radius of the dusty torus (R in ).The MIR light curves can be represented by the convolution of the optical light curves with a transfer function of the torus, which depends on the torus geometry.In a toy model considering the torus geometry, Li & Shen (2023) showed that the slope of the SF for low ∆t increases with increasing inner radius because the short-term variability in the optical light curves can be easily smoothed out with a large R in (see Figure 4 in Li & Shen 2023).In addition, the MIR variability can be substantially suppressed for large R in , as the torus transfer function is more widely spread over time at larger R in .Therefore, the inverse correlation between [O III] luminosity and A (γ 1 ) may suggest that R in is proportional to the AGN luminosity.Interestingly, this finding agrees well with the prediction from the size−luminosity relation, in which R in is tightly correlated with the AGN luminosity (e.g., Koshida et al. 2014;Yang et al. 2020;Li & Shen 2023).This finding is also in good agreement with the receding torus model (e.g., Barvainis 1987;Lawrence 1991;Simpson 2005), in which Rin is expected to be proportional to the AGN luminosity.
Meanwhile, some studies have argued that the slope γ of the SF (SF ∝ ∆t γ ) is weakly correlated with the AGN luminosity even for the optical light curves (e.g., Kozlowski 2016;Smith et al. 2018).However, the variation of γ in the optical studies is significantly smaller than that in the MIR, hinting that γ derived from MIR variability is physically connected with the torus geometry.The distribution of γ 1 as a function of [O III] luminosity in type 1 and type 2 AGNs is almost identical (Figure 12).This suggests that both AGN types share a similar torus geometry, in line with expectation in the AGN unification paradigm.This result, however, seems to be at odds with the finding that type 2 AGNs have systematically lower variability amplitudes than type 1s.Extinction by the torus might be to blame, an explanation that resonates with previous observational and theoretical studies (Nenkova et al. 2008;Ramos Almeida et al. 2011;Hickox et al. 2017).Intriguingly, low-luminosity type 2 AGNs exhibit relatively low variability amplitudes, in apparent contradiction with type 1 AGNs and previous studies.The discrepancy is more evident over long-time scales.
While the physical origin of the variability behavior is unclear, we suspect that the W1 flux is contaminated by non-variable sources (e.g., warm and cold dust, or hot nanometer-sized grains in the host galaxies ;Tran 2001;Xie et al. 2018).Alternatively, the physical properties of the torus are distinctive from those of normal AGNs.In particular, the correlation between the variability properties and the Eddington ratio for type 1 and type 2 AGNs are clearly different (Figure 12).Some have argued that AGNs with very low Eddington ratios could not effectively form dense clouds around the AD (e.g., Elitzur & Ho 2009;Elitzur et al. 2014).Such type 2 AGNs, so-called unobscured or intrinsic type 2 AGNs (Panessa & Bassani 2002;Tran et al. 2011;Ho et al. 2012), naturally lack a broad-line region and dusty torus.This may account for the lower detection rate of the broad-line region in LINERs than in normal AGNs (e.g.Ho et al. 1997;Ho & Kim 2009).According to this model, the low variability amplitude of low-luminosity type 2 AGNs can be attributed to the small torus covering factor, which poses a challenge to the traditional AGN unification model (Ichikawa et al. 2015).This may also reveal that the low A in low-luminosity type 2 AGNs can be a by-product of the dependence on the Eddington ratio, as the bolometric luminosity is correlated with the Eddington ratio for type 2 AGNs (Figure 1).To test this scenario, we compare the [O III] luminosity and W1 luminosity (L W1 ) solely from the AGN, which is obtained from the SED fit described in Section 4.3 (Figure 13).For type 2 AGNs, we adopt the AGN templates with a large extinction (A v = 20) to account for the intrinsic  extinction by the torus.We find that not only type 2 AGNs exhibit a lower fraction of L W1 to [O III] luminosity compared to type 1 AGNs, possibly because of the extinction, but the fraction also dramatically decreases toward low [O III] luminosity.This again implies that the covering factor of the torus in low-luminosity type 2 AGNs may be significantly lower than in type 1 AGNs of similar AGN luminosity.Dust-deficient AGNs systematically exhibit larger A, γ 1 , and γ 2 compared to normal AGNs.As expected from the correlation between R in and γ 1 inferred from the torus model of Li & Shen (2023), the steeper slopes for dust-deficient AGNs may indicate that the dust, which radiates the MIR continuum in dust-deficient AGNs, is likely distributed over larger radii than in normal AGNs (Li & Shen 2023).This interpretation is consistent with the fact that dust-deficient AGNs lack hot and warm dust, which is more centrally located within the torus.In addition, a larger variability amplitude in dust-deficient AGNs than normal AGNs regardless of AGN properties may imply that the hot and warm dust is physically unstable in dust-deficient AGNs.Perhaps it can be easily destroyed and reformed, depending on the strength of the high-energy photons from the AD, which would naturally lead to a high variability amplitude in the MIR.Alternatively, the systematic differences between normal and dust-deficient AGNs in terms of their MIR SF may arise from the fact that the two types intrinsically exhibit different variability character-istics in the UV/optical continuum from the AD.This can be further addressed in future studies with optical light curves.

SUMMARY
We assemble a 10-year long multi-epoch MIR dataset for type 1 and type 2 AGNs selected from SDSS that are also observed by NEOWISE in the W1 band.With these data, we examine the MIR variability by rigorously computing the ensemble SF of subsamples according to various AGN properties.We also perform an extensive simulation to consider the sparse cadence of the WISE dataset, investigate the systematic uncertainty involved, and evaluate the uncertainties in the SF measurements.This analysis yields the following results.
• The overall shapes of the ensemble SFs for both type 1 and type 2 AGNs are nearly identical, except that the MIR variability amplitude of type 1 AGNs is systematically larger than that of type 2 sources.This may be attributed to the large extinction in type 2 AGNs.
• The variability amplitude is inversely correlated with the AGN luminosity, likely a reflection of the intrinsic variability of the light from the AD.The same trend holds for the slope of the SF for low ∆t (γ 1 ), which we interpret as the positive correlation between the inner radius of the torus and the AGN luminosity.
• Dust-deficient AGNs tend to have larger variability amplitudes and slopes in the ensemble SF than do normal AGNs.While the physical origin of the larger variability amplitudes and slopes remains unclear, one factor may be that the covering factor of the centrally located hot and warm dust sensitively responds to the variation of UV photons from the AD.
• Low-luminosity type 2 AGNs tend to have lower variability amplitudes than their type 1 counterpart or high-luminosity type 2 AGNs.We suggest that low-luminosity type 2 AGNs have a characteristically low torus covering factor.Alternatively, the W1 flux is significantly contaminated by the cold dust in the host galaxy.
The above results lead us to conclude that MIR variability can be a powerful tool for probing the torus properties of AGNs.For example, time-series nearinfrared spectroscopic data obtained through the upcoming space mission Spectro-Photometer for the History of the Universe, Epoch of Reionization and Ices Ex-Son et al. plorer (SPHEREx) will play a crucial role in this effort (e.g., Doré et al. 2018;Kim et al. 2021).
We thank the anonymous referee for constructive comments that greatly helped to improve the manuscript.This work was supported by the National Science Foundation of China (11721303, 11991052, 12011540375, 12233001), the China Manned Space Project (CMS-CSST-2021-A04, CMS-CSST-2021-A06), and the National Research Foundation of Korea (NRF), through a grant funded by the Korean government (MSIT) (No. 2022R1A4A3031306 and 2023R1A2C1006261).

Figure 1 .
Figure 1.Distribution of Eddington ratio versus (a) bolometric luminosity and (b) BH mass for our sample of AGNs.The blue circles and open histogram denote type 1 AGNs, while filled red circles and filled histogram represent type 2 AGNs.

Figure 2 .
Figure 2. Ensemble SF estimated from light curves simulated (a) with no noise added and (b) with noise added.The dashed line is the input SF.The blue circles and red stars denote the mean of the SF and the square root of the mean SF, respectively.
Figure 3.Ensemble SF of non-variable sources.The dashed line denotes the mean error of the sources.The blue circles and red stars denote the mean of the SF and the square root of the mean SF, respectively.

Figure 4 .
Figure 4. Distribution of the flux fraction of the continuum from the accretion disk in the W1 band, after subtracting host contribution, for type 1 AGNs.The expected flux from the accretion disk (FW1,AD) is calculated from the g-band flux.The dashed and dotted lines denote the mean value and standard deviation, respectively.

Figure 5 .
Figure 5. Ensemble SF of type 1 (open blue circles) and type 2 (filled red circles) AGNs.The fitting results with a broken power law are denoted by the blue and red dashed lines.
Figure 6.Ensemble SF for subsamples binned by the continuum luminosity at 5100 Å, in units of erg s −1 .The dashed lines denote the fits with a broken power law.

Figure 8 .
Figure 8. Ensemble SF of subsamples binned by BH mass, in units of M⊙.The open and filled symbols denote (a) type 1 and (b) type 2 AGNs, respectively.The dashed lines represent the fits with a broken power law.
Figure 9. Ensemble SF of subsamples binned by Eddington ratio (λ Edd ), using bolometric luminosities calculated from the [O III] luminosities.The open and filled symbols denote (a) type 1 and (b) type 2 AGNs, respectively.The dashed lines represent the fits with a broken power law.

Figure 11 .
Figure 11.Ensemble SF for HDD (red circles), WDD (green squares), and normal (blue diamonds) AGNs at fixed (a) [O III] luminosity and (b) BH mass.The fits with a broken power law are denoted by corresponding dashed lines.

Figure 12 .
Figure 12.Dependence of A (upper panels) and γ1 (lower panels) on the [O III] luminosity (left), BH mass (middle), and Eddington ratio (right).The open blue and filled red circles denote type 1 and type 2 AGNs, respectively.

Figure 13 .
Figure 13.Correlation between [O III] luminosity and W1 luminosity, in units of erg s −1 .The symbols are the same as in Figure 1.The solid and dashed lines represent the linear fits for type 1 and type 2 AGNs, respectively.
Figure 7. Ensemble SF of subsamples binned by [O III] luminosity, in units of erg s −1 .The open and filled symbols denote (a) type 1 and (b) type 2 AGNs, respectively.The dashed lines represent the fits with a broken power law.
Son et al.