Evolution of Quasar Stochastic Variability along Its Main Sequence

The structure function,⁹ SF(Δt), measures the statistical dispersion of two random variables (i.e., a magnitude pair) separated by time intervals, Δt. The structure function can be used to characterize the statistical dispersion of Δm for a sample of many similar quasars with the same (or close) Δt, where Δm is the magnitude difference between two observations. We adopted the interquartile range (i.e., IQR) to measure the statistical dispersion as it is robust against outliers or tails in the distribution. Therefore, we calculate the statistical dispersion as follows (MacLeod et al. 2010; Sun et al. 2015),

$$\begin{eqnarray}&&{\mathrm{SF}}_{\mathrm{IQR}}({\rm{\Delta }}t)=0.74\mathrm{IQR}({\rm{\Delta }}m),\end{eqnarray} \tag{ 1 }$$

where IQR(Δm) is the 25%–75% interquartile range of Δm. The constant 0.74 normalizes the IQR to be equivalent to the standard deviation of a Gaussian distribution. Therefore, 0.74 IQR is known as the normalized IQR (hereafter NIQR).

It should be noted that the measured statistical dispersion (i.e., Equation (1)) is a superposition of measurement errors and quasar variability. On very short timescales (e.g., days), the amplitude of quasar variability is small and the statistical dispersion is dominated by measurement errors. Therefore, we can estimate measurement errors from the statistical dispersion on timescales of a few days. On timescales of months to years, the contribution of measurement errors becomes negligible.

The CAR(1) process is often referred to as the damped random walk or the Ornstein–Uhlenbeck (OU) process; this process has proven to be effective in describing the light curves of quasar continuum emission (e.g., Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2010, 2012; Zu et al. 2013). The structure function of the CAR(1) process is given by

$$\begin{eqnarray}&&\mathrm{SF}({\rm{\Delta }}t| \tau ,\hat{\sigma })=\hat{\sigma }\sqrt{\tau (1-\exp (-{\rm{\Delta }}t/\tau ))},\end{eqnarray} \tag{ 2 }$$

where ${\rm{\Delta }}=| {t}_{i}-{t}_{j}|$ is the separation time between two observations. That is, the CAR(1) process is characterized by two parameters, $\hat{\sigma }$ and τ. The former, $\hat{\sigma }$, determines the short-term variability amplitude while the latter, τ, is the characteristic timescale.

It should be noted that quasar variability might be more complex than the CAR(1) process. Therefore, Kelly et al. (2014) proposed more flexible continuous-time autoregressive moving-average (i.e., CARMA(p, q)) models to describe quasar light curves; the CAR(1) process corresponds to the CARMA(1, 0) process. For each source in our parent sample, we use the Python CARMA package¹⁰ and adopt the Akaike information criterion (AIC; Akaike 1974) to choose the order of the CARMA(p, q) models (i.e., determining p and q that minimize AIC; see Section 3.5 of Kelly et al. 2014); we also calculated AIC for the CAR(1) process (hereafter AIC(1)). We found that, for most of our light curves (∼90%), the differences between the minimum AIC and AIC(1) is less than 10. Therefore, it seems that the data quality of our sample is insufficient to distinguish between the CAR(1) process and other more complex models. In Section 5.2, we model the structure functions as the CAR(1) process (i.e., Equation (2)); however, more complex models (i.e., Equation (3)) are also discussed. If quasar variability is indeed not driven by the stochastic models we assumed or the light curve is a nonstationary process, the uncertainties of our model parameters in Section 5.2 and Tables 3 and 4 might be inaccurate (or even underestimated; see e.g., White 1982).

We aim to explore the ensemble variability of quasar continuum as a function of Rfe. The ensemble structure functions for the three bins of the luminosity-matched sample are presented in Figure 2. Low-${R}_{\mathrm{Fe}{\rm{II}}}$ quasars tend to be more variable (for a statistical description of our conclusion, see Section 5.2). Our result is well expected if: (1) EV1 is indeed driven by Eddington ratio and (2) for fixed luminosity, high Eddington-ratio quasars are more stable. The former assumption is supported by independent tests (e.g., Shen & Ho 2014; Sun & Shen 2015). We discuss possible explanations of the second requirement in Section 5.3.

Zoom In Zoom Out Reset image size

The tendency between ${R}_{\mathrm{Fe}{\rm{II}}}$ and the quasar variability amplitude might be induced by FWHM of Hβ¹¹ since there might be an anticorrelation between FWHM and ${R}_{\mathrm{Fe}{\rm{II}}}$ (see Table 1) and M_BH ∝ FWHM². In order to verify this speculation, we explore quasar variability as a function of FWHM after controlling ${R}_{\mathrm{Fe}{\rm{II}}}$, L_bol , and z. Therefore, we construct samples as follows. First, we select quasars within ${10}^{45.3}\ \mathrm{erg}\,{{\rm{s}}}^{-1}\leqslant {L}_{\mathrm{bol}}\leqslant {10}^{45.6}\ \mathrm{erg}\,{{\rm{s}}}^{-1}$ and 0.4 < ${R}_{\mathrm{Fe}{\rm{II}}}$ < 1. We now choose a slightly wider luminosity bin to increase the statistic. Second, these sources are divided into two bins according to FWHM, i.e., the low- (high-) FWHM bin with FWHM being smaller (larger) than $\mathrm{Median}(\mathrm{FWHM})$. Third, we ensure L_bol and Fe ii strength of the two samples are matched via the methodology in Section 2. The number of quasars in the low- (high-) FWHM bin is 75 (74). The median values of FWHM for the two subsamples are $3127\ \mathrm{km}\,{{\rm{s}}}^{-1}$ and $6278\ \mathrm{km}\,{{\rm{s}}}^{-1}$, respectively. As shown in Figure 3, the ensemble structure functions for the two subsamples are quite similar. Therefore, it seems that the relation between quasar variability and the virial M_BH is rather weak or absent since, for fixed quasar luminosity, M_BH ∝ FWHM².

Zoom In Zoom Out Reset image size

We also control FWHM, L_bol, and z, and divide sources into two ${R}_{\mathrm{Fe}{\rm{II}}}$ bins following the method mentioned above. That is, we select quasars within ${10}^{45.3}\ \mathrm{erg}\,{{\rm{s}}}^{-1}\leqslant {L}_{\mathrm{bol}}\leqslant {10}^{45.6}\ \mathrm{erg}\,{{\rm{s}}}^{-1}$ and $3000\ \mathrm{km}\,{{\rm{s}}}^{-1}\lt \mathrm{FWHM}\lt 5000\ \mathrm{km}\,{{\rm{s}}}^{-1}$ and divide them into two bins according to ${R}_{\mathrm{Fe}{\rm{II}}}$. We calculate the structure functions for the two bins. We again find that sources with larger ${R}_{\mathrm{Fe}{\rm{II}}}$ tend to be less variable (see Figure 4). These conclusions provide additional evidence supporting the claim that orientation determines the dispersion of FWHM (e.g., Shen & Ho 2014). We discuss this idea in Section 5.1.

Zoom In Zoom Out Reset image size

In the previous section, we demonstrate the relation between quasar variability and ${R}_{\mathrm{Fe}{\rm{II}}}$. To examine whether there is an additional dependence on L_bol, we compare the ensemble structure functions of the ${R}_{\mathrm{Fe}{\rm{II}}}$-matched sample (see Figure 5). On short timescales (i.e., 1 ≲ Δt ≲ 100 days), there is a clear anticorrelation between quasar variability and L_bol (for a statistical description of our conclusion, see Section 5.2). This tendency diminishes on long timescales (i.e., Δt ≳ 100 days). Therefore, it seems that: (1) L_bol controls the short-term (1 ≲ Δt ≲ 100 days) quasar variability and (2) ${R}_{\mathrm{Fe}{\rm{II}}}$ drives quasar variability on timescales of Δt ≳ 10 days.

Zoom In Zoom Out Reset image size

According to our inspection of the structure function described in Section 4, quasar variability at a given wavelength in the UV/optical bands and on timescales from weeks to years can be characterized by L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$. There is no additional correlation between quasar variability and FWHM. Our results can be well explained in the framework that the Eddington ratio and orientation govern most of the quasar diversity (Shen & Ho 2014). According to this scenario, the EV1 is driven by the Eddington ratio; high-Fe ii-strength sources have high Eddington ratios and are less variable; FWHM is a tracer of orientation and does not correlate with quasar variability.

To test whether FWHM traces orientation, we compare the $r-W1$ color of the low-FWHM sample with that of the high-FWHM sample, where W1 refers to the WISE 3.4 μm band. To obtain W1, we cross-match our quasars with the ALLWISE catalog¹² (Wright et al. 2010; Mainzer et al. 2011) with the maximum matching radius of 2''. The left panel of Figure 6 presents our results. Indeed, sources in the broad-FWHM bin tend to have redder SED than those in the narrow-FWHM sample (the p value of the Anderson–Darling test is <0.01). Similar results have been obtained by Shen & Ho (2014). Therefore, broad- (narrow-) FWHM sources are consistent with being viewed more edge-(face-) on. If so, the geometry of BLR is disk-like rather than spherical, which is consistent with other observations (e.g., Jarvis & McLure 2006; Pancoast et al. 2014; Grier et al. 2017a; Storchi-Bergmann et al. 2017; Xiao et al. 2018). The orientation scenario also naturally explains the lack of correlation between quasar variability and FWHM (Figure 3).

Zoom In Zoom Out Reset image size

[O iii] EW has also been proposed as a tracer of orientation (e.g., Risaliti et al. 2011). We also show the distributions of [O iii] EW for the broad and narrow FWHM samples in the right panel of Figure 6. Contrary to our expectation, we cannot reject the null hypothesis that the two distributions of [O iii] EW are drawn from the same population (the p value of the Anderson–Darling test is 0.4). Therefore, we conclude that [O iii] EW is driven by ${R}_{\mathrm{Fe}{\rm{II}}}$(i.e., the EV1; Boroson & Green 1992) or the maximum disk temperature (Panda et al. 2017a; Panda et al. 2017b; Panda et al. 2018) rather than orientation.

Previous works (e.g., MacLeod et al. 2010, 2012; Kozłowski 2016) aimed to find correlation between quasar variability as a function of L_bol, and M_BH. Often in these works, M_BH is estimated via the single-epoch virial black hole mass estimators, i.e., $\mathrm{log}{M}_{\mathrm{BH}}={p}_{0}+{p}_{1}\mathrm{log}L+{p}_{2}\mathrm{logFWHM}$, where p₀, p₁, and p₃ are constants (e.g., Vestergaard 2002; Vestergaard & Peterson 2006; Shen et al. 2011). However, as we demonstrate in Section 4.1 and Figures 3 and 4, there is no clear relation between quasar variability and FWHM. Therefore, we relate quasar variability to L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$.

The main purpose of this section is to provide new empirical relations for future variability modeling. Therefore, for simplicity, we assume quasar variability is a CAR(1) process (which can, in practice, describe the light curves well; see, e.g., Kelly et al. 2009; MacLeod et al. 2010, 2012; Zu et al. 2013).

We aim to explore the correlations between the CAR(1) parameters (i.e., τ and $\hat{\sigma }$) and quasar properties (i.e., L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$). Following Kozłowski (2016), we constrain $\hat{\sigma }$ and τ by modeling the ensemble structure function with

$$\begin{eqnarray}&&{\mathrm{SF}}^{2}({\rm{\Delta }}t| \tau ,\hat{\sigma })={\hat{\sigma }}^{2}\tau (1-\exp {(-{\rm{\Delta }}t/\tau )}^{\beta })+{\sigma }_{{\rm{p}}}^{2},\end{eqnarray} \tag{ 3 }$$

where σ_p is the uncertainty of the magnitude difference between two observations separated by Δt. We fix β = 1 (i.e., the CAR(1) process, see Equation (2)) in our subsequent analysis (we try to set β as a free parameter in Section 5.3).

To explore the dependence of τ on L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$, we perform the following analysis. For each bin of the ${R}_{\mathrm{Fe}{\rm{II}}}$-matched sample, we assume that the ensemble CAR(1) parameters are determined by

$$\begin{eqnarray}&&\mathrm{log}\tau ={c}_{11}+{\kappa }_{11}\mathrm{log}({\bar{L}}_{\mathrm{bol}}/{10}^{45.5}\ \mathrm{erg}\,{{\rm{s}}}^{-1}),\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&\mathrm{log}\hat{\sigma }={c}_{21}+{\kappa }_{21}\mathrm{log}({\bar{L}}_{\mathrm{bol}}/{10}^{45.5}\ \mathrm{erg}\,{{\rm{s}}}^{-1}),\end{eqnarray} \tag{ 5 }$$

where ${\bar{L}}_{\mathrm{bol}}$ is the average of L_bol in each bin. We also try to remove galaxy contamination to ${\bar{L}}_{\mathrm{bol}}$ by applying the empirical relation of Equation (1) in Shen et al. (2011). We then calculate the theoretical structure function from these two equations and Equation (3).

We fit the theoretical structure functions to the three ensemble structure functions of the ${R}_{\mathrm{Fe}{\rm{II}}}$-matched sample via a Bayesian approach. The likelihood function is

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}p(f| x,\mathrm{pms},{\sigma }_{\mathrm{int}})\\ \quad =\,-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{i=3}\displaystyle \sum _{n}\left[\displaystyle \frac{{({f}_{n,i}-{f}_{\mathrm{model},n,i})}^{2}}{{s}_{n,i}^{2}}+\mathrm{ln}(2\pi {s}_{n,i}^{2})\right],\end{array}\end{eqnarray} \tag{ 6 }$$

where x represents a set of quasar parameters (e.g., L_bol, ${R}_{\mathrm{Fe}{\rm{II}}}$); pms is a collection of parameters c₁₁, c₂₁, κ₁₁, and κ₂₁; ${f}_{\mathrm{model},n,i}$ and ${f}_{n,i}$ are the theoretical and observational structure functions, respectively; i = 1, 2, 3 represents the three bins; and n indicates each Δt(n). The vector s_n,i denotes the summation of the measurement uncertainty of f and the (possible) intrinsic scatter (which is assumed to be σ_int f_model,n,i).¹³ That is, ${s}_{n,i}^{2}={({\sigma }_{\mathrm{int}}{f}_{\mathrm{model},n,i})}^{2}+{f}_{\mathrm{err},n,i}^{2}$, where ${f}_{\mathrm{err},n,i}$ is the bootstrap uncertainty of ${f}_{n,i}$. The priors are summarized in Table 2. We use the MCMC code emcee to sample the posterior distributions of the parameters.

Table 2.  Priors of the Parameters

	Parameter	Min	Max	Distribution
	c1i	0.0	4.0	Uniform
Equations (4) and (5) (i = 1)	c2i	−10.0	10.0	Uniform
or	${\kappa }_{1i}$	−2.0	2.0	Uniform
Equations (7) and (8) (i = 2)	${\kappa }_{2i}$	−2.0	2.0	Uniform
	$\mathrm{ln}{\sigma }_{\mathrm{int}}$	−10.0	10.0	Uniform
	a	−20.0	6.0	Uniform
Equation (11)	b	−5.0	5.0	Uniform
	c1	−5.0	5.0	Uniform
	$\mathrm{ln}{\sigma }_{\mathrm{int}}$	−10.0	10.0	Uniform
	a	−20.0	6.0	Uniform
	b	−5.0	5.0	Uniform
Equations (14)	Ωm	0.0	1.0	Uniform
	ΩΛ	0.0	1.0	Uniform
	$\mathrm{ln}{\sigma }_{\mathrm{int}}$	−10.0	10.0	Uniform

Download table as:  ASCIITypeset image

The best-fitting structure functions are the solid lines in Figure 5. The posterior distributions of c₁₁, c₂₁, κ₁₁, and κ₂₁ are shown in Figure 7 and are summarized in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Statistical Properties of the Parameters for the Ensemble Structure Function as a Function of L_bol or ${R}_{\mathrm{Fe}{\rm{II}}}$

	Parameter	Median ± NIQR
	c11	2.55 ± 0.03
	c21	−1.85 ± 0.01
Equations (4) and (5)	κ11	0.50 ± 0.08
	κ21	−0.26 ± 0.02
	$\mathrm{ln}{\sigma }_{\mathrm{int}}$	−2.95 ± 0.10
	c12	2.38 ± 0.07
	c22	−1.72 ± 0.02
Equations (7) and (8)	κ12	0.001 ± 0.070
	κ22	−0.08 ± 0.02
	$\mathrm{ln}{\sigma }_{\mathrm{int}}$	−2.41 ± 0.10
	c12	2.38 ± 0.04
	c22	−1.72 ± 0.02
Equations (7) and (8) (with κ12 ≡ 0)	κ12	0 (fixed)
	κ22	−0.08 ± 0.01
	$\mathrm{ln}{\sigma }_{\mathrm{int}}$	−2.42 ± 0.10

Download table as:  ASCIITypeset image

The correlation (i.e., the slope κ₁₁) between τ and L_bol for fixed ${R}_{\mathrm{Fe}{\rm{II}}}$ (or fixed Eddington ratio) might simply reflect the dependence of τ on M_BH (e.g., Kelly et al. 2009; MacLeod et al. 2010, 2012; Kozłowski 2016). If so, we expect a strong correlation between τ and Eddington ratio (or ${R}_{\mathrm{Fe}{\rm{II}}}$) for fixed L_bol. To test this argument and explore quasar variability as a function of ${R}_{\mathrm{Fe}{\rm{II}}}$, we fit the ensemble structure functions of the luminosity-matched sample with

$$\begin{eqnarray}&&\mathrm{log}\tau ={c}_{12}+{\kappa }_{12}{\bar{R}}_{\mathrm{Fe}\mathrm{II}},\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&\mathrm{log}\hat{\sigma }={c}_{22}+{\kappa }_{22}{\bar{R}}_{\mathrm{Fe}\mathrm{II}}.\end{eqnarray} \tag{ 8 }$$

The priors are summarized in Table 2. The statistical properties of the distributions are summarized in Table 3. To our surprise, the correlation between τ and ${R}_{\mathrm{Fe}{\rm{II}}}$ is statistically insignificant as κ₁₂ = 0.001 ± 0.070. Therefore, we conclude that τ depends mostly on L_bol.

We then refit the ensemble structure functions of the luminosity-matched sample with Equations (7) and (8) but fix κ₁₂ ≡ 0 (i.e., we assume τ does not depend on ${R}_{\mathrm{Fe}{\rm{II}}}$). The statistical properties of the distributions are summarized in Table 3. The best-fitting structure functions are the solid lines in Figure 2. By fixing κ₁₂ ≡ 0, the intrinsic scatter of the fit ($\mathrm{ln}{\sigma }_{\mathrm{int}}=-2.41$) is similar to that of the previous fit ($\mathrm{ln}{\sigma }_{\mathrm{int}}=-2.42$). That is, τ and ${R}_{\mathrm{Fe}{\rm{II}}}$ are not tightly correlated.

Combining the best-fitting relations for the ${R}_{\mathrm{Fe}{\rm{II}}}$- and luminosity-matched samples, we can derive quasar variability as a function of L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$, i.e.,

$$\begin{eqnarray}&&\mathrm{log}\tau =2.49+0.50(\mathrm{log}{L}_{\mathrm{bol}}-45.50),\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\mathrm{log}\hat{\sigma }=-1.788\mbox{--}0.26(\mathrm{log}{L}_{\mathrm{bol}}-45.50)-0.08{R}_{\mathrm{Fe}{\rm{II}}}.\end{eqnarray} \tag{ 10 }$$

For each S82 quasar with "good" data (e.g., at least 10 epochs and small measurement errors), MacLeod et al. (2010) fit the CAR(1) process to the light curve and constrained $\hat{\sigma }$ and τ. In principle, we can adopt their data and fit the best-fitting parameters as a function of quasar properties. However, Kozłowski (2017) recently demonstrated that if the baseline is not ∼5–10 times larger than τ, the best-fitting CAR(1) parameters are biased. The biases are negligible for $\hat{\sigma }$ but are rather strong for τ. Therefore, we should only focus on $\hat{\sigma }$.

The function we use to relate $\hat{\sigma }$, L_bol, and ${R}_{\mathrm{Fe}{\rm{II}}}$ is

$$\begin{eqnarray}&&\mathrm{log}\hat{\sigma }=a+b(\mathrm{log}{L}_{\mathrm{bol}}-45)+{c}_{1}{R}_{\mathrm{Fe}\mathrm{II}}.\end{eqnarray} \tag{ 11 }$$

For comparison, we also try to fit the following function:

$$\begin{eqnarray}&&\mathrm{log}\hat{\sigma }=a+b(\mathrm{log}{L}_{\mathrm{bol}}-45)+{c}_{2}\mathrm{log}\,\mathrm{FWHM}.\end{eqnarray} \tag{ 12 }$$

We fit the functions of Equations (11) and (12) to quasars in the range 0.5 < z < 0.89 (i.e., a narrow range of redshift) via a Bayesian approach. The likelihood function is

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{ln}p(\mathrm{log}\hat{\sigma }| x,{\sigma }_{x},\mathrm{pms},{\sigma }_{\mathrm{int}})\\ \quad =\,-\displaystyle \frac{1}{2}\displaystyle \sum _{n}\left[\displaystyle \frac{{(\mathrm{log}{\hat{\sigma }}_{n}-\mathrm{log}{\hat{\sigma }}_{\mathrm{model},n})}^{2}}{{s}_{n}^{2}}+\mathrm{ln}(2\pi {s}_{n}^{2})\right],\end{array}\end{eqnarray} \tag{ 13 }$$

where x is $[\mathrm{log}({L}_{\mathrm{bol}}),{R}_{\mathrm{Fe}\mathrm{II}}]$ (or $[\mathrm{log}({L}_{\mathrm{bol}}),\mathrm{FWHM}$]); σ_x is the uncertainty of x; pms represents parameters a, b, and c; ${\hat{\sigma }}_{\mathrm{model}}$ is given by Equation (11) or (12); and σ_int is a summation of the measurement uncertainty of $\hat{\sigma }$ and the intrinsic scatter. Furthermore, ${s}_{n}^{2}={\sigma }_{\mathrm{int}}^{2}+{(b{\sigma }_{L})}^{2}+{(c{\sigma }_{c})}^{2}$, where σ_L and σ_c are the measurement errors of L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$ (or FWHM), respectively, and σ_int represents the statistical dispersion due to either measurement errors of $\mathrm{log}\hat{\sigma }$ or the intrinsic scatter. The priors are summarized in Table 2.

The statistical properties of the parameters a, b, c₁, and σ_int for $\hat{\sigma }$ as a function of L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$ (i.e., Equation (11)) are presented in Table 4. Our results indicate that while the short-term variability is mainly driven by L_bol, an additional dependence on ${R}_{\mathrm{Fe}{\rm{II}}}$ (or Eddington ratio) is also statistically significant.

Table 4.  Statistical Properties of the Parameters for the CAR(1) Parameter $\hat{\sigma }$ as a Function of Quasar Properties

	Parameter	Median ± NIQR
	a	−1.70 ± 0.02
Equation (11)	b	−0.29 ± 0.02
	c1	−0.05 ± 0.01
	$\mathrm{ln}{\sigma }_{\mathrm{int}}$	−1.81 ± 0.03
	a	−1.86 ± 0.11
Equation (12)	b	−0.30 ± 0.02
	c2	−0.03 ± 0.02
	$\mathrm{ln}{\sigma }_{\mathrm{int}}$	−1.78 ± 0.02

Download table as:  ASCIITypeset image

In the works of MacLeod et al. (2010) and Kozłowski (2016), they explored the dependencies of ${\mathrm{SF}}_{\infty }$ and τ on L_bol and M_BH. Using their best-fitting relations, we can also obtain the relation between $\hat{\sigma }$, L_bol, and Eddington ratio. In both works, the dependence of $\hat{\sigma }$ on L_bol is close to our result. Kelly et al. (2009, 2013) also obtained a similar relation using light curves of the international AGN Watch projects.¹⁴ However, the correlation between $\hat{\sigma }$ and Eddington ratio is statistically insignificant in these works.

It is quite possible that in previous works the correlation between $\hat{\sigma }$ and Eddington ratio was diluted by the large uncertainty in M_BH due to orientation. Indeed, after controlling L_bol and Rfe, the ensemble structure function does not depend on FWHM (see Figure 3). To confirm our guess, we explore the dependence of $\hat{\sigma }$ on L_bol and FWHM by fitting Equation (12). The priors are summarized in Table 2. The statistical properties of the distributions are summarized in Table 4. As we expected, there is indeed no correlation between $\hat{\sigma }$ and FWHM (the slope, c₂, is statistically consistent with 0). Therefore, the additional dependence of $\hat{\sigma }$ on Eddington ratio is missed in previous works.

In this work, we find a dependency of the variability parameters on L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$. Therefore, it is likely that the optical/UV variability is produced in the quasar central engine. Several models have been proposed to explain the connection between the optical/UV variability and quasar properties. For instance, Li & Cao (2008) proposed that variations in the global accretion rate drive quasar optical/UV variability (see also Zuo et al. 2012). However, such models failed to explain timescale-dependent color variability (e.g., Sun et al. 2014; Cai et al. 2016; Zhu et al. 2016). Instead, models with local fluctuations (possibly regulated by some common variations; see Cai et al. 2018) in the accretion disk are more compatible with observations. The local fluctuation model can also produce the CAR(1) process (Lin et al. 2012). Meanwhile, X-ray reprocess might also play a role (e.g., Czerny et al. 1999) although no significant correlation between X-ray and UV/optical variations has been found (Kelly et al. 2011, 2013) and the color variability might not be explained by X-ray reprocess (Zhu et al. 2017).

Kelly et al. (2013) proposed that the variance of the short-term variability per τ_TH is a constant. If so, for a fixed observational timescale, ${\hat{\sigma }}^{2}\propto 1/{\tau }_{\mathrm{TH}};$ from the accretion disk theory, we expect τ_TH scales with ${L}_{\mathrm{bol}}^{1/2}$. Therefore, this scenario predicts $\hat{\sigma }\propto {L}_{\mathrm{bol}}^{-1/4}$. This scenario can explain our best-fitting relation between $\hat{\sigma }$ and L_bol (see Tables 3 or 4).

In contrast to previous works, we find a correlation between $\hat{\sigma }$ and ${R}_{\mathrm{Fe}{\rm{II}}}$. The additional dependence of $\hat{\sigma }$ on ${R}_{\mathrm{Fe}{\rm{II}}}$ might be induced by X-ray reprocessing. High- and low-${R}_{\mathrm{Fe}{\rm{II}}}$ (Eddington ratio) quasars tend to have weaker and stronger X-ray emission, respectively (Lusso et al. 2012). As a result, X-ray reprocessing is more efficient and can induce more variations in UV/optical bands for low-${R}_{\mathrm{Fe}{\rm{II}}}$ sources. A promising alternative explanation is that Eddington ratio might correlate with gas metallicity (e.g., Matsuoka et al. 2011). If so, high-${R}_{\mathrm{Fe}{\rm{II}}}$ quasars are iron-overabundant, and their accretion disks are more stable (Jiang et al. 2016).

The scatter of $\hat{\sigma }$ as a function of L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$ is slightly smaller than that of the relation between $\hat{\sigma }$, L_bol , and FWHM. These scatters are caused by measurement errors (which is 0.088 dex) and intrinsic scatter. Guo et al. (2017) argued that the intrinsic scatter is caused by the deviation from the CAR(1) process on long timescales (see their Figure 9). Based on this hypothesis, they constrained the PSD of quasar variability on long timescales to be steeper than f^−1.3. According to our best-fitting results, the intrinsic scatter in their work is slightly overestimated since they related $\hat{\sigma }$ to FWHM. Therefore, the PSD of quasar variability on long timescales approaches the 1/f relation. Such a PSD is expected from the local variations of accretion rate (Lyubarskii 1997; Noble & Krolik 2009).

We find that the characteristic timescale, τ, is mostly driven by L_bol (see Section 5.2; Table 3). This solo dependence and the normalization encourage us to link τ with the thermal timescale (τ_TH). The best-fitting slope (0.50 ± 0.08) is remarkably consistent with the theoretical expectation (i.e., the thermal timescale ${\tau }_{\mathrm{TH}}\sim {\text{}}{L}_{\mathrm{bol}}^{0.5}$). It should be noted that even if τ is the thermal timescale, there might still be an anticorrelation between τ and M_BH for fixed L_bol. This is because the thermal timescale of an accretion disk depends positively on iron abundance (Jiang et al. 2016); quasars with high Eddington ratios might be more metal-rich than those with low Eddington ratios (Matsuoka et al. 2011). However, such a correlation is not found in our results. It is possible that this correlation is weak and is unable to be revealed in our data.

Recent works suggested that significant deviations occur on very short timescales (i.e., ∼days; see, e.g., Mushotzky et al. 2011; Kasliwal et al. 2015; Smith et al. 2018). However, on timescales we consider here (months to years), this deviation should not be very important. Kozłowski (2016) revealed a positive correlation between β and L_bol by studying the S82 quasars. We then refit Equation (3) to the ${R}_{\mathrm{Fe}{\rm{II}}}$-matched samples via the same Bayesian approach but set β as a free parameter. We do not find a significant correlation between β and L_bol. The discrepancy might be caused by the following factors. First, our selected S82 quasars have much lower luminosity than those of Kozłowski (2016). Second, we use ${R}_{\mathrm{Fe}{\rm{II}}}$ rather than the ratio of L_bol to the virial M_BH (which is likely biased by orientation) to trace the unknown Eddington ratio.

The strong correlation between τ and L_bol is also found by Caplar et al. (2017). We note, however, that they adopted a different method to constrain τ. In some other previous works (MacLeod et al. 2010; Kozłowski 2016), τ is found to be insensitive to L_bol but depends on the virial M_BH. The differences between our results and those of Kozłowski (2016) might also be caused by the factors mentioned above.¹⁵

However, it should be noted that quasar variability on long timescales is likely not consistent with the CAR(1) process (MacLeod et al. 2012; Guo et al. 2017). If this is found to be true, it is unclear whether or not we can directly compare τ with physical timescales. The forthcoming era of time domain astronomy will be key to understanding the physical nature of τ.

Our work and many previous works (e.g., MacLeod et al. 2010; Kozłowski 2016; Caplar et al. 2017) indicate that the short-term UV/optical variability amplitude (or $\hat{\sigma }$) critically depends on L_bol. In this work, we also find an additional dependence of $\hat{\sigma }$ on ${R}_{\mathrm{Fe}{\rm{II}}}$. This additional dependence is statistically significant but rather weak since the slope is −0.05 ± 0.01 (see Table 4). In practice, we can ignore this additional dependence and fit $\hat{\sigma }$ as a function of L_bol only (i.e., the parameter c₁ in Equation (11) is fixed at zero). The best-fitting parameters are $\hat{\sigma }=(-1.74\pm 0.012)-(0.30\,\pm 0.018){L}_{\mathrm{bol}};$ the scatter (i.e., σ_int which is a combination of measurement errors and the intrinsic scatter) is exp(−1.81 ± 0.026) which is the same as that of Equation (11). We can, in principle, estimate L_bol from the $\hat{\sigma }$–L_bol relation without assuming any cosmological models. Therefore, it is possible to use the short-term UV/optical variability of quasars as a probe of cosmology parameters.

The $\hat{\sigma }$–L_bol relation (i.e., Equation (11) with c₁ fixed at zero) can be revised as

$$\begin{eqnarray}&&\mathrm{log}\hat{\sigma }=a+b\mathrm{log}(4\pi f/{10}^{45}\ \mathrm{erg}\,{{\rm{s}}}^{-1})+2b\mathrm{log}({D}_{L}),\end{eqnarray} \tag{ 14 }$$

where f is the observed flux and D_L, the luminosity distance, is a function of the cosmological model and can be independently measured if we know $\hat{\sigma }$, f, a, and b. Given the intrinsic scatter of the $\hat{\sigma }$–L_bol relation, such constraints can only be made with a large sample of quasars that span over a wide range of cosmic history.

To illustrate this idea, we perform the following simulation of 10⁵ quasars. For each quasar, the intrinsic L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$ and their measurement errors are assigned according to the randomly selected quasar from our parent sample. We then calculate ${\hat{\sigma }}_{\mathrm{th}}$ from our best-fitting Equation (11); a Gaussian noise with standard deviation of exp(−1.81) (see Section 5.2) is added to ${\hat{\sigma }}_{\mathrm{th}}$ to generate the observed $\hat{\sigma }$. We also assign galaxy contamination according to Equation (1) of Shen et al. (2011). The observed $\hat{\sigma }$ is diluted by the nonvariable galaxy emission. The observed L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$ are generated by perturbing intrinsic L_bol and ${R}_{\mathrm{Fe}{\rm{II}}}$ with their measurement errors. In addition, the galaxy emission is added to the observed L_bol. To calculate the observed flux, we assume a flat ΛCDM cosmology with h₀ = 0.7 and Ω_m = 0.3; redshift is randomly assigned from a uniform distribution within [0.1, 0.89]. We then fit Equation (14) to the simulated mock sample by considering the ΛCDM cosmology with h₀ = 0.7 via a Bayesian approach. Both Ω_m (i.e., the matter density fraction) and Ω_Λ (i.e., the dark energy fraction) are free parameters. The likelihood function is the same as Equation (13) and the priors are summarized in Table 2.

The posterior distributions of the model parameters are presented in Figure 8. Even if Ω_Λ is not constrained, the recovered Ω_m = 0.28 ± 0.03 is accurate. We note that if the sample size is limited to 10⁴, the recovered Ω_m = 0.36 ± 0.12. Therefore, the large sample size is one of the key factors.

Zoom In Zoom Out Reset image size

Our simulated sample might be available in the era of time-domain astronomy (e.g., with the Large Synoptic Survey Telescope). However, it remains unclear whether the scatter of the $\hat{\sigma }$–L_bol relation depends on the sample size/redshift or not. In order to test this hypothesis, we select sources with 0.96 < z < 1.48 or 1.48 < z < 2.03. Their r or i bands correspond to the rest-frame of 0.5 < z < 0.9 g band. We calculate the differences between their $\hat{\sigma }$ and the expectated values from our best-fitting $\hat{\sigma }$–L_bol relation. Some of the differences are due to a combination of the scatter of the relation and the measurement error of L_bol, σ_L (i.e., the total scatter is $\sqrt{{\sigma }_{\mathrm{int}}^{2}+{(b{\sigma }_{L})}^{2}}$). We then calculate the ratio of the differences to this total scatter. We find that for 95% of sources the ratio is less than 3. Many of the remaining 5% of sources are highly variable (20% have $\hat{\sigma }\gt 1$). Such sources might be "changing-look" AGN candidates (MacLeod et al. 2016); the origin of such variability could be different. Therefore, it is unlikely that the scatter of the $\hat{\sigma }$–L_bol relation significantly depends on the sample size/redshift. We could also use high-redshift (z ∼ 2) quasars to constrain cosmological parameters; the accuracy would be further improved.

In addition to our method, it is also proposed that the BLR and dust reverberation (Watson et al. 2011; Yoshii et al. 2014), the nonlinear relation between the UV and X-ray luminosities (Risaliti & Lusso 2015), the X-ray variability and broad line width (La Franca et al. 2014), and the saturated luminosity of super-Eddington AGNs (Wang et al. 2013) could also be adopted as distance measurements. In conclusion, AGNs will play a more important role in measuring the universe (for a recent review, see Czerny et al. 2018).