© The Authors 2025

1. Introduction

Active galactic nuclei (AGNs) are compact objects in which a substantial amount of gas accretes onto supermassive black holes (BHs) at the centers of galaxies; this forms an accretion disk that emits intense radiation, often comparable to or even exceeding the total emission of the host galaxy. AGNs produce powerful, multiwavelength light because their intricate central regions are governed by a broad range of physical processes that span multiple emission mechanisms (e.g., Elvis et al. 1994; Richards et al. 2006). In addition, the strength of the radiation from the AGNs fluctuates stochastically on various timescales from days to years (e.g., Ulrich et al. 1997). Notably, the UV and optical continuum, called the big blue bump, that arises from the accretion disk is inherently variable. This variability is echoed by other reprocessed light, such as broad emission lines and the near-infrared continuum (e.g., Peterson 1993; Kaspi et al. 2000; Koshida et al. 2014; Kim et al. 2024). While the physical cause of the variability is still under debate, various studies have claimed that it may be associated with instabilities within the accretion disk (e.g., Kawaguchi et al. 1998). According to this scenario, the variability timescale is expected to be correlated with the thermal timescale of the accretion disk (e.g., Czerny et al. 1999; Burke et al. 2021; Tang et al. 2023; Son et al. 2025). In the standard disk model, the thermal timescale is expected to be associated with the orbital timescale, which is proportional to the square of the wavelength and the square root of AGN luminosity.

The light curve of the UV and optical continuum has often been modeled as a damped random walk (DRW) process (Kelly et al. 2009; Kozlowski et al. 2010b, but see Kasliwal et al. 2015). In this framework, the power spectral density (PSD) is characterized by a power law whose index is inversely proportional to the square of the frequency at higher frequencies. Crucially, the PSD becomes flat at a frequency below a specific break frequency. However, the PSD can only be accurately estimated with multi-epoch photometric data with a regular cadence. For this reason, the structure function (SF) has been used as an alternative to study the AGN variability (e.g., Kozlowski et al. 2010b; Kozlowski 2017; Son et al. 2023). Computed as the root mean square of brightness at a given time difference (Δt), the SF can be easily estimated with sparsely sampled light curves obtained from ground-based telescopes. In the SF, the DRW model can be described by a single power law (i.e., SF ∝Δt) at a timescale lower than the characteristic timescale (τ) and equivalent to the break frequency in the PSD, above which the SF flattens above the characteristic timescale (e.g., Kozlowski et al. 2010a).

Previous studies based on multi-epoch broadband photometric data widely utilized the SF to characterize AGN variability (e.g., MacLeod et al. 2010; Kozlowski 2017; Li et al. 2018; Suberlak et al. 2021; Son et al. 2025). A variety of datasets have been used to infer the variability amplitude from the value of SF at a certain timescale and the characteristic timescale near the knee of the SF (e.g., Sánchez et al. 2017; Burke et al. 2023; Tang et al. 2023). Additionally, their dependence on AGN properties (BH mass, AGN luminosity (L), Eddington ratio, and rest-frame wavelength (λ)) was rigorously explored to reveal the mechanism of AGN variability and hence to investigate the structure of the accretion disk. For example, MacLeod et al. (2010) estimated the SFs of luminous AGNs in Stripe 82 based on the broadband photometric data of the Sloan Digital Sky Survey (SDSS) over a baseline of approximately ten years. They showed that the characteristic timescale (τ) directly estimated from the SFs is correlated with the wavelength and BH mass, but is nearly independent of luminosity. Using the PSD of a multi-epoch photometric dataset of ∼400 AGNs, spanning a wide range of BH mass (10⁴ ≤ M_BH/M_⊙ ≤ 10¹⁰), Burke et al. (2023) showed that τ is mostly governed by the BH mass. Intriguingly, these results are contradictory to the prediction from the standard accretion disk model (Shakura & Sunyaev 1973; see also Arévalo et al. 2024; Petrecca et al. 2024). More recently, Tang et al. (2023), using the photometric data from ATLAS (Tonry et al. 2018), demonstrated that the characteristic timescale is proportional to L^{0.539 ± 0.004}λ^{2.418 ± 0.023}, which is in broad agreement with the theoretical prediction from the standard accretion disk model (τ ∝ L^0.5λ²). These conflicting results may reflect the potential drawback of the broadband photometric data, which can be contaminated by strong emissions from the broad-line region (BLR) and continuum from neighboring wavelengths owing to the large width of the filter transmission (e.g., Patel et al. 2025).

Multi-epoch spectroscopic data offer a better alternative to overcome the limitations of broadband photometric data in characterizing AGN variability. This method enables a more robust estimation of how AGN variability depends on the continuum wavelength. Acquiring time-series spectra of AGNs demands considerable resources and thereby limits the availability of sufficient datasets for comprehensive studies of AGN variability. Recently, however, several long-term (on approximate timescales of a few years) monitoring programs, primarily conducted for reverberation mapping (RM) experiments, offer opportunities to study the spectral variability for a large number of AGNs. For instance, Son et al. (2025) use the spectral data of a few hundred luminous AGNs from the SDSS RM project (Shen et al. 2019, 2024), which spans approximately seven years. Their study focuses on the SF across various rest-frame wavelengths. Carefully accounting for the uncertainty in the flux calibration of the observed spectra, Son et al. (2025) demonstrate that the variability timescale is proportional to L^{0.50 ± 0.03}λ^{1.42 ± 0.09}. This finding suggests that the dependence on AGN luminosity aligns with predictions from the standard accretion disk model. However, the observed dependence on wavelength significantly deviates from the model prediction, which implies a temperature profile that is significantly steeper. Although Son et al. (2025) clearly illustrate the utility of SF derived from the temporal spectral data in exploring the structure of accretion disks, their work is limited to a relatively small sample consisting predominantly of distant high-luminosity AGNs.

Conventional RM projects provide intensive monitoring data for a limited number of AGNs. While valuable, this approach restricts the sample size and the dynamic range of observable AGN properties. This study investigates spectral variability using an alternative method: analyzing multi-epoch spectroscopic observations of a large sample of AGNs from the main SDSS program. This method allows for a significant expansion of the sample size and extends the dynamic ranges of AGN properties, which facilitates more comprehensive studies of spectral variability. This work revisits AGN spectral variability, focusing on the SF methodology, by analyzing more than two million spectroscopic pairwise combinations generated from over 90 000 AGNs from the SDSS. Throughout the paper, we adopt the cosmological parameters from the 2018 Planck results (H₀ = 67.36 ± 0.54 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.6847 ± 0.0073, Ω_m = 0.3153 ± 0.0073; Planck Collaboration VI 2020).

2. Sample and data

2.1. Sample selection

The spectroscopic variability can be investigated in two main ways: (1) through continuous spectroscopic monitoring data of a limited sample of AGNs, primarily conducted for RM analyses (e.g., Son et al. 2025), and (2) through multi-epoch spectroscopic datasets, typically two observations for each object, across a large AGN sample. In this study, we adopted the latter method, as a statistically significant amount of data is available in the SDSS. We initially chose the sample from the SDSS quasi-stellar object (QSO) catalogs from data release 16 (DR16; Lyke et al. 2020) with a selection criterion of NSPEC ≥ 2, where NSPEC represents the number of duplicate observations. Starting from DR9, the SDSS spectra were acquired using the Baryon Oscillation Spectroscopic Survey (BOSS; Smee et al. 2013) spectrograph, which differs from the original SDSS spectrograph in both the fiber aperture and wavelength coverage. To minimize potential systematics in flux calibration due to the instrumental differences, we restrict our analysis to spectra obtained with the BOSS spectrograph from DR9, which results in a sample of 94 851 objects and 2 233 454 spectroscopic pairwise combinations.

2.2. AGN properties

As the ultimate goal of this study is to investigate how the AGN variability depends on the physical properties of AGNs, it is essential to derive the AGN properties accurately. We directly measured the continuum fluxes (f₁₃₆₇ or f₃₀₅₀) at 1367 Å or 3050 Å from the mean spectrum of each object in the rest frame, which is free from strong emission lines. They are corrected for Galactic extinction using E(B − V) values from Schlafly & Finkbeiner (2011) and the extinction curve of Cardelli et al. (1989) and converted to the conventional luminosities (L₁₃₅₀ and L₃₀₀₀) at 1350 Å and 3000 Å assuming f_λ ∝ λ^−1.56 (Vanden Berk et al. 2001). Figure 1 compares L₁₃₅₀ and L₃₀₀₀ for objects with both measurements available. As the slope of the correlation clearly deviates from unity, we used an ordinary least-squares linear regression to derive the conversion between the two luminosities: log (L₃₀₀₀/erg s⁻¹) = 0.94 log L₁₃₅₀/erg s⁻¹) + 2.59. Finally, the bolometric luminosity (L_bol) was computed with a bolometric correction of (L_bol = 3.81L₁₃₅₀) adopted from Richards et al. (2006). We used L₁₃₅₀ whenever possible. If L₁₃₅₀ was not covered by the SDSS spectrum, L₃₀₀₀ served as an alternative.

Fig. 1. Density distributions of luminosities at 1350 Å and 3000 Å for our sample. The solid and dashed red lines denote the best regression and its 1σ scatter, respectively. The best-fit result displayed in the panel was used to convert between the two luminosities.

For the BH mass (M_BH), we used the measurements from Wu & Shen (2022), based on broad C IVλ1549 and Mg IIλ2800 emission from the single-epoch spectra¹. Priority was given to the Mg II because C IV-based BH masses can suffer from systematic biases (e.g., Ho et al. 2012; Mejía-Restrepo et al. 2016). Subsequently, the Eddington ratio (L_bol/L_Edd) was computed based on the Eddington luminosity L_Edd = 1.26 × 10³⁸M_BH/M_⊙ erg s⁻¹. Radio-loud AGNs are often associated with extreme variability, which makes them unsuitable for our analysis and thus need to be excluded (e.g., Vanden Berk et al. 2004; MacLeod et al. 2010). To identify radio-loud AGNs, we calculated the radio-to-optical flux ratio determined at 2500 Å (f₂₅₀₀) and 6 cm (f_6 cm; Kellermann et al. 1989). We initially measured f₂₅₀₀ at the nearest continuum windows (2230, 3050, 3500, or 3800 Å) in the rest frame, depending on the redshift of the target, and subsequently extrapolated to 2500 Å by assuming a power-law continuum derived from the SDSS composite spectrum (f_λ ∝ λ^−1.56; Vanden Berk et al. 2001). We estimated f_6 cm from the 20 cm flux density measurements of radio counterparts matched within 2″ from the FIRST survey (Lyke et al. 2020). For the flux conversion from 20 cm to 6 cm and k-correction, we assumed f_ν ∝ ν^−0.5 (Kellermann et al. 1989; but see Ho & Ulvestad 2001). Objects with a radio loudness (R = f_6 cm/f₂₅₀₀) greater than 10 (2954 objects) were classified as radio-loud objects and discarded for further analysis. The distributions of BH mass, bolometric luminosity, and Eddington ratio of the final sample of 91 897 sources are shown in Figure 2.

Fig. 2. Distributions of BH mass and bolometric luminosity for our sample. The solid, dashed, and dotted lines correspond to Eddington ratios of 1, 0.1, and 0.01, respectively.

2.3. Continuum fluxes and uncertainties

To minimize the effect of the emission lines, including Fe II multiples, and starlight from the host galaxy, we estimated the continuum fluxes at 1367, 1479, 1746, 2230, 3050, 3500, and 3800 Å in the rest frame (Son et al. 2025). At each continuum wavelength, the flux was estimated using the mean flux value, and its uncertainty (σ_in) was calculated as $\sqrt{\frac{\sum _{n=1}^N{\sigma }_i^2}{N}}$, where N is the number of spectral elements within a continuum window² of 40 Å and σ_i is the uncertainty associated with each spectral element. However, the flux calibration of the SDSS spectra can introduce systematic uncertainty, which should be considered for accurate SF estimates (e.g., Margala et al. 2016; Son et al. 2025). While in principle this uncertainty can be estimated using multi-epoch spectra of standard stars, the standard stars are systematically brighter than our sample, which could possibly lead to an underestimate of the uncertainty.

Instead, we utilized the SF of the spectroscopic pairs in the AGN sample with short time intervals (Δt < 2 days in the observed frame), assuming that the AGN variability on such timescales is negligible. The SF of an observed pair was calculated as

$$\begin{array}{c}\hfill {\mathrm{SF}}_{2\mathrm{days}}^2={(m(t)-m(t+\mathrm{\Delta }t))}^2-{\sigma }_{\mathrm{in}}^2(t)-{\sigma }_{\mathrm{in}}^2(t+\mathrm{\Delta }t)(\mathrm{\Delta }t<2\mathrm{days}),\end{array}$$(1)

where Δt is the time lag between the observed pair, m is the observed magnitude at a given continuum wavelength, and σ_in is the initial uncertainty of magnitude in the single-epoch spectrum. We found that the average SF² at Δt < 2 days (SF_2days²) is mostly greater than 0 across all the continuum wavelengths, which indicates that the original uncertainties (σ_in) were significantly underestimated. Notably, the average SF_2days² depends on the continuum brightness and is well described by a double power-law model (Fig. 3). Therefore, we incorporate this excess, which depends on the flux density and continuum wavelengths, when estimating the structure function over longer timescales (Sect. 3.1). Specifically, the SF_2 day² was modeled using the form

$$\begin{array}{c}\hfill {\mathrm{SF}}_{2\mathrm{days}}^2=A{\left(\frac{f}{f_b}\right)}^{-{\alpha }_1}{\left\{\frac{1}{2}[1+{\left(\frac{f}{f_b}\right)}^{1/\mathrm{\Delta }}]\right\}}^{({\alpha }_1-{\alpha }_2)\mathrm{\Delta }},\end{array}$$(2)

Fig. 3. Ensemble SF at Δt < 2 days (SF2days) plotted against continuum brightness across different wavelengths. The solid line represents the fitting results with a smooth double power-law model.

where f is the continuum flux density, f_b is the break flux density, A is the SF_2days² at the low-luminosity end, α₁ is the power-law slope at f_b ≫ f, α₂ is the power-law slope when f ≫ f_b, and Δ is the smoothness parameter. As α₁ was found to be consistent with 0 across wavelengths, we fixed its value to 0 in the fit. As previous studies showed that the flux calibration is highly uncertain at the edges of the BOSS spectrum (Son et al. 2025), we restrict our analysis to spectra within the 3820–9150 Å range to minimize this effect.

3. Method and results

3.1. Analysis of spectral variability with ensemble SF

To analyze the AGN spectral variability, we utilized the ensemble SF of the spectral dataset. While SFs are typically calculated for individual objects, we grouped observed pairs of AGNs with similar physical properties instead (e.g., BH mass, bolometric luminosity, and Eddington ratio). We then computed the ensemble SF by averaging the individual SFs at each time lag. Therefore, the square of the ensemble SF was calculated as

$$\begin{array}{cc}\hfill {\mathrm{SF}}^2(\mathrm{\Delta }t)& =\frac{1}{N_{\mathrm{\Delta }t,\mathrm{pair}}}{\displaystyle \sum _{i=1}^{N_{\mathrm{\Delta }t,\mathrm{pair}}}}{(m(t)-m(t+\mathrm{\Delta }t))}^2\hfill \\ \hfill & -{\sigma }_{\mathrm{in}}^2(t)-{\sigma }_{\mathrm{in}}^2(t+\mathrm{\Delta }t)-{\mathrm{SF}}_{2\mathrm{day}}^2,\hfill \end{array}$$(3)

where N_Δt, pair is the number of observed pairwise combinations at the given Δt across the AGN subsample, which shares similar AGN properties, and SF_2day² is estimated from the mean magnitude of the observed pairs. It is important to note that Δt was expressed in the rest frame for this calculation. We performed this analysis in each continuum window. The distributions of the number of observed pairs as a function of Δt across different continuum windows are shown in Figure 4. This indicates that a sufficient number of pairs exist up to Δt ≤ 1000 days. The extreme AGN variability can occur in a minority of AGNs, manifesting as changing-look AGNs, episodic flares, and tidal disruption events. These phenomena can substantially increase the variability amplitude and thus strongly influence the estimation of the ensemble SFs. Therefore, we applied the 3σ clipping method when averaging the SFs to minimize this effect.

Fig. 4. Distribution of observed spectroscopic pairwise combinations at each Δt for our sample. Colors denote continuum wavelengths.

3.2. Dependence of SF on AGN properties and wavelength

To explore how the SF varies with AGN properties and continuum wavelengths, we estimated the SFs for subsamples categorized by their AGN properties, such as BH mass, bolometric luminosity, and Eddington ratio. By accounting for the distributions of the physical parameters of the entire sample, we divided the sample into subgroups of five bins covering 10^44.75 − 10^46.30 erg s⁻¹, 10^7.3 − 10^9.6 M_⊙, and 10^−2.2 − 10^0.3 for the AGN luminosity, BH mass, and Eddington ratio, respectively. On average, there are approximately 2500 pairwise combinations per bin, given wavelength, and Δt, ranging from 65 to about 19 000.

As illustrated in Figure 5, SFs are highly correlated with the bolometric luminosity and wavelengths. Notably, AGNs with higher luminosity or those at longer wavelengths tend to have lower SFs or longer variability timescales than their lower-luminosity or shorter-wavelength counterparts.

Fig. 5. Ensemble SFs for subgroups categorized according to bolometric luminosity (a), BH mass (b), and Eddington ratio (c). The size of the symbols denotes the continuum wavelength.

The standard accretion disk model theoretically predicts that the thermal timescale is intrinsically tied to the orbital timescale. Previous observational studies also proposed that AGN variability timescales are governed by the thermal timescale (e.g., Collier & Peterson 2001; Kelly et al. 2009). This implies that the observed dependence of SFs on the AGN luminosity and wavelength can be interpreted through variations in the thermal timescale and the radial structure of the accretion disk. Accordingly, this model predicts that the variability timescale is strongly correlated with both AGN luminosity and wavelengths. Motivated by this, we quantified the dependence of the variability timescale on AGN luminosity and wavelength (τ ∝ L^αλ^β), assuming that all AGNs share a universal ensemble SF when the time lag is normalized by the characteristic timescale (τ) (e.g., Tang et al. 2023; Son et al. 2025).

We assumed that SF ∝ a(Δt/τ)^b, with τ scaling as τ ∝ L^αλ^β, where a is the SF at Δt = τ and b is the power-law index for the entire sample. We then performed the fit using the ensemble SFs computed from subsamples divided by AGN luminosity and wavelength. In this calculation, we only included the ensemble SFs derived from more than 500 spectroscopic pairs to minimize a possible systematic error due to small sample sizes (Fig. 6). This analysis, which is based on the χ² minimization, yields that τ ∝ L^{0.62 ± 0.07}λ^{1.74 ± 0.23}. At longer wavelengths (λ > 3050 Å), the continuum may be contaminated by the Balmer continuum from the BLR and by stellar light from the host galaxy. In particular, the effect of the Balmer continuum can be complex, as it may be less variable than the underlying continuum and its variability may be delayed due to the extended geometry of the BLR. To mitigate this, we repeated the experiment excluding the wavelengths greater than 3050 Å, which resulted in τ ∝ L^{0.68 ± 0.07}λ^{1.12 ± 0.24}. Specifically, the power-law index (β) for the wavelength changes dramatically from ∼1.74 to ∼1.12, which suggests a substantial contribution from either the Balmer continuum or the host light. Additionally, to mitigate the potential impact of outliers in the SFs, we applied a 3σ clipping during the fitting process. The fitting result remained unchanged within the uncertainties.

Fig. 6. Ensemble SFs with Δt normalized by the characteristic timescale (τ). We estimated τ assuming that all AGNs share a universal SF, with τ scaling as τ ∝ Lαλβ.

For the Eddington ratio, similar trends were observed in the sense that AGNs with higher Eddington ratios tend to have smaller SFs or longer timescales, possibly due to their correlation with the AGN luminosity (Fig. 5a). With an assumption of the universal ensemble SF, we found τ ∝ (L_bol/L_Edd)^{0.70 ± 0.08}λ^{1.99 ± 0.24} when the fitting was done with the entire wavelength range up to 3800 Å. The fitting results also change significantly to τ ∝ (L_bol/L_Edd)^{0.67 ± 0.08}λ^{1.43 ± 0.31} when the wavelength over 3050 Å is excluded. However, intriguingly, the dependence of SFs on BH mass is not significant.

4. Discussion

4.1. Primary parameter that determines AGN variability

Based on the SF analysis for subgroups of the AGN sample, we found that the variability timescales strongly depend on both AGN luminosity and Eddington ratio. As noted in the previous section, the two parameters are closely linked to each other, which makes it difficult to identify the primary driver of AGN variability. To disentangle this degeneracy, we examined the effect of each parameter by fixing one parameter as a constant while varying the other. First, we selected a subset of AGNs with a narrow range of Eddington ratio (−1 ≤ log (L_bol/L_Edd)≤ − 0.7) and recalculated the ensemble SFs for subsamples divided by AGN luminosity. To investigate the dependence of τ on AGN luminosity, we performed the same analysis under the assumption of a universal SF. This yielded a scaling relation of τ ∝ L^{0.95 ± 0.10}λ^{4.15 ± 0.36}. This analysis revealed that the variability timescale retains a significant dependence on AGN luminosity, which suggests that AGN luminosity plays an important role (Fig. 7). However, given the limited sample size³, the extracted scaling parameters should be interpreted with caution, as they may not reflect true physical relationships.

Fig. 7. Ensemble SFs of AGNs with a fixed Eddington ratio, grouped by bolometric luminosity.

As a counterexample, we fixed the AGN luminosity within the range of 10^45.8 − 10^46.2 erg s⁻¹, and assessed the ensemble SFs for subsamples categorized by Eddington ratio. Intriguingly, under these conditions, the dependence on the Eddington ratio vanished, as indicated by the relation τ ∝ (L_bol/L_Edd)^{−0.01 ± 0.06} (Fig. 8). This result implies that the Eddington ratio has a minimal impact on the variability timescale when AGN luminosity is held nearly constant. Taken together, AGN luminosity is the main parameter that governs the variability timescales, rather than the Eddington ratio. We note, however, that the marginal correlation between the Eddington ratio and BH mass suggests that secondary effects resulting from this relationship warrant consideration. While this hypothesis could, in principle, be tested by additionally fixing the BH mass, the limited sample size precludes such an analysis within the scope of this study. Nevertheless, given that the BH mass appears to exert the least significant influence on the SF across the entire sample (Fig. 5), its impact is likely not dominant.

Fig. 8. Ensemble SFs of AGNs with a fixed bolometric luminosity, grouped by Eddington ratio.

4.2. Physical interpretations

In theoretical studies based on the standard accretion disk framework, the thermal timescale scales as τ_th ∝ M_BHR^3/2, where R is the radial distance in the accretion disk. Additionally, the radial temperature profile of the disk follows $T(R)\propto {(M_{\mathrm{BH}}\dot{M}{R}^{-3})}^{1/4}$, where Ṁ represents the mass accretion rate. If the AGN luminosity is linearly correlated with Ṁ, the thermal timescale should be proportional to L^0.5λ² MacLeod et al. (2010), Tang et al. (2023), Son et al. (2025). Interestingly, this prediction is broadly consistent with our observational findings, which show that the variability timescale, inferred from the SFs, is primarily dependent on the AGN luminosity and wavelength.

Quantitatively, however, the best-fit power indices derived from our analysis (α = 0.62 ± 0.07 and β = 1.74 ± 0.23) exhibit marginal deviations from the theoretical prediction within a 2σ uncertainty. We note that the deviation in β dramatically increases when the wavelengths longer than 3050 Å, where substantial contamination from the Balmer continuum and stellar light is expected, are excluded from the fit (β = 1.12 ± 0.24). The fitted value of β falls below the theoretical prediction with a significance of ∼4σ. This result underscores the critical role of accurately accounting for host galaxy starlight and Balmer continuum contamination in studies of AGN variability.

Interestingly, the disagreement in β with the standard disk model has also been reported in previous observational studies (e.g., Son et al. 2025, but see also Tang et al. 2023). It is worthwhile noting that β is directly linked to the radial temperature profile under the assumption that the flux at a given wavelength primarily originates from the radius where the local temperature corresponds to that wavelength via Wien’s displacement law (λ ∝ 1/T). The significantly low value of β suggests that the radial temperature profile in the accretion disk may be steeper than predicted by the standard disk model (e.g., Son et al. 2025). This effect can result in a harder spectral energy distribution (SED) of the UV continuum, which is typically contrary to previous findings that report a softer-than-expected UV SED in AGNs (e.g., Zheng et al. 1997; Shang et al. 2005). The larger disk size, relative to predictions from the standard accretion disk model, can also be attributed to a shallower temperature profile (e.g., Sun et al. 2019). Conversely, an intrinsically harder UV continuum is required to explain the strength of broad emission lines in luminous AGNs (e.g., Lawrence 2012), which is consistent with a steeper temperature profile.

Alternatively, our results are consistent with predictions from the Corona-Heated Accretion-disk Reprocessing (CHAR) model, in which a steep temperature profile can be reproduced if the flux at a given wavelength is significantly contaminated by emission from neighboring radii (Sun et al. 2020). Within this framework, Zhou et al. (2024) demonstrate that the characteristic timescale depends on AGN luminosity and wavelength, with power-law indices of α = 0.65 and β = 1.19, which are in excellent agreement with our results. The results of Zhou et al. (2024) were derived from simulated light curves constructed using the AGN parameters of QSOs in Stripe 82, as adopted from Stone et al. (2022).

4.3. Comparison with previous studies

Mostly, AGN variability has been studied from a photometric dataset obtained with broadband filters (e.g., MacLeod et al. 2010; Sánchez-Sáez et al. 2018; Arévalo et al. 2023). Photometric data, however, can be severely affected by the emission from BLR, which can lead to systematic biases in characterizing the AGN variability. For example, based on the multi-epoch photometry of QSOs in Stripe 82, MacLeod et al. (2010) showed that τ exhibits moderate correlation with BH mass and wavelengths but little correlation with AGN luminosity. More recently, Arévalo et al. (2023) analyzed multi-epoch g-band photometry from the Zwicky Transient Facility (ZTF; Masci et al. 2019). Burke et al. (2021) analyzed broadband light curves from AGNs across a wide range of BH masses (10⁴ M_⊙ < M_BH ≤ 10¹⁰ M_⊙), collected from surveys like SDSS and various RM projects. Both studies discovered a correlation between τ and BH mass, contrary to the results in this study. Although the direct cause of this discrepancy with our findings remains unclear, it is worth noting that, unlike this study, Tang et al. (2023) estimated τ directly from the SFs or PSDs of multi-year broadband photometric light curves. Accurate measurements of τ are challenging, even when using light curves with baselines exceeding ten times the value of τ (Kozlowski 2017).

On the other hand, Tang et al. (2023) analyzed the approximately five-year multi-epoch broadband photometry of ∼5000 QSOs from the Asteroid Terrestrial-impact Last Alert System (ATLAS) survey using the same method in this study and found that τ is primarily governed by AGN luminosity and wavelength. While this result aligns with our finding, their reported value of β ∼ 2.4 is systematically larger than that obtained in this study. To examine the impact of broadband filters, we artificially generated the broadband photometric light curves using the spectroscopic datasets employed in this study. To mimic broadband conditions, we set the width of the continuum windows to 1000 Å. When we applied the same analysis method, we found α = 0.67 ± 0.06 and β = 2.13 ± 0.22, which agrees well with the results of Tang et al. (2023). More importantly, excluding continua with central wavelengths longer than 3050 Å yielded a similar trend (α = 0.68 ± 0.07 and β = 1.77 ± 0.28), with no significant change in the parameters derived from including the entire wavelength range. These findings demonstrate that, in broadband photometric data, contamination from strong emissions and the Balmer continuum originating from the BLR is non-negligible. Additionally, the range of AGN luminosity also differs among the observational studies (∼10^{46 − 47} erg s⁻¹ for Tang et al. 2023 and ∼10^{44.8 − 46.2} erg s⁻¹ for our study).

On the other hand, Son et al. (2025) utilized the multi-epoch spectroscopic dataset obtained through the SDSS RM project and analyzed the SFs of ∼800 QSOs, which, to our knowledge, was the first effort to specifically construct the SF for a large sample. Applying the same method in this study, they show that τ is mostly determined by AGN luminosity and wavelengths with α = 0.50 ± 0.03 and β = 1.42 ± 0.09. It is noteworthy that their results are consistent with ours at the ∼2σ level, despite the differences in the spectroscopic datasets: Son et al. (2025) employed an approximately seven-year time series of spectra for hundreds of QSOs, whereas this study used 2 million spectroscopic pairwise combinations for tens of thousands of QSOs. This confirms the reliability of our approach to explore AGN variability. Son et al. (2025) show that the fitting result for α and β is insensitive to the wavelengths of the continuum windows, revealing that the contribution from the Balmer continuum may not be crucial. However, on the contrary, the best-fitting value for β changes significantly with and without longer wavelengths above 3050 Å. This may be attributed to the fact that our sample comprises less luminous AGNs compared to Son et al. (2025), and hence the contribution from the host galaxy is more substantial in our sample.

To mitigate potential biases introduced by broadband photometric data and/or the limited range of AGN physical parameters resulting from a relatively small sample size, simultaneous monitoring with both photometric and spectroscopic datasets over a large sample is ideal. Such an approach will be enabled by the Time Domain Extragalactic Survey (Frohmaier et al. 2025) conducted with the 4 m Multi-Object Spectroscopic Telescope, in coordination with the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST; Ivezić et al. 2019).

5. Summary

We utilized multi-epoch spectroscopic data of ∼90 000 luminous AGNs from SDSS DR16 and examined the continuum ensemble SFs for subsamples categorized by various AGN properties, such as BH mass, AGN luminosity, and Eddington ratio. Based on this analysis, we reached the following conclusions:

• Overall, the UV SFs are strongly dependent on wavelength, AGN luminosity, and Eddington ratio. However, they are weakly or barely correlated with the BH mass. To identify the primary parameter determining the SFs of AGNs, we re-examined the SFs by fixing either the AGN luminosity or the Eddington ratio. This revealed that the SFs are solely dependent on wavelength and AGN luminosity.

• Under the assumption of a universal SF shared by all the AGNs in which the variability timescale (τ) is determined by wavelength (λ) and AGN luminosity (L), τ is found to scale with α = 0.62 ± 0.07 and β = 1.74 ± 0.23 following the equation τ ∝ L^αλ^β. This finding aligns with the prediction (α = 0.5 and β = 2) from the standard disk model to within a 2σ level.

• Excluding continua above 3050 Å allowed us to minimize the flux contribution from the host galaxy and Balmer continuum. We found β = 1.12 ± 0.24, which indicates a significant ∼4σ deviation from the standard disk model. This outcome may suggest that the temperature profile of the accretion disk is steeper than that predicted by the standard model.

• The findings of this study significantly deviate from those of previous studies based on broadband photometric data. To robustly examine the dependence of SFs on AGN properties, particularly wavelength, it is essential to use spectroscopic datasets. This underscores the importance of spectroscopic monitoring of AGNs.

Acknowledgments

We are deeply grateful to the anonymous referee for insightful and constructive comments that have substantially enhanced the manuscript. LCH was supported by the National Science Foundation of China (12233001) and the China Manned Space Program (CMS-CSST-2025-A09). This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (Nos. RS-2024-00347548 and RS-2025-16066624).

References

All Figures

	Fig. 1. Density distributions of luminosities at 1350 Å and 3000 Å for our sample. The solid and dashed red lines denote the best regression and its 1σ scatter, respectively. The best-fit result displayed in the panel was used to convert between the two luminosities.
In the text

	Fig. 2. Distributions of BH mass and bolometric luminosity for our sample. The solid, dashed, and dotted lines correspond to Eddington ratios of 1, 0.1, and 0.01, respectively.
In the text

	Fig. 3. Ensemble SF at Δt < 2 days (SF2days) plotted against continuum brightness across different wavelengths. The solid line represents the fitting results with a smooth double power-law model.
In the text

	Fig. 4. Distribution of observed spectroscopic pairwise combinations at each Δt for our sample. Colors denote continuum wavelengths.
In the text

	Fig. 5. Ensemble SFs for subgroups categorized according to bolometric luminosity (a), BH mass (b), and Eddington ratio (c). The size of the symbols denotes the continuum wavelength.
In the text

	Fig. 6. Ensemble SFs with Δt normalized by the characteristic timescale (τ). We estimated τ assuming that all AGNs share a universal SF, with τ scaling as τ ∝ Lαλβ.
In the text

	Fig. 7. Ensemble SFs of AGNs with a fixed Eddington ratio, grouped by bolometric luminosity.
In the text

	Fig. 8. Ensemble SFs of AGNs with a fixed bolometric luminosity, grouped by Eddington ratio.
In the text