STATISTICAL PROPERTIES OF MULTI-EPOCH SPECTRAL VARIABILITY OF SDSS STRIPE 82 QUASARS

We introduce the flux–flux correlation method (Choloniewski 1981; Winkler et al. 1992; Lyutyi & Doroshenko 1993) as an alternative to magnitude–magnitude or magnitude–color correlation analyses. Choloniewski (1981) first noticed that the quasars' UV–optical fluxes have a linear correlation in flux–flux space, and several authors have since confirmed this (e.g., Winkler 1997; Glass 2004; Sakata et al. 2011; Lira et al. 2011; Gu & Li 2013). In Figure 2, we show the Pearson correlation coefficient of the f_ν(r) − f_ν(i) plot for our sample (0 < z < 4) observed for more than 20 epochs (see Section 3.5 for more details of our sample). In terms of the Pearson correlation coefficient, almost all of the quasars have a strong positive flux–flux correlation. The strong correlation for a quasar light curve means that the flux difference spectrum, which is defined as the flux difference of the quasar spectra of any of the two observational epochs, keeps a nearly constant shape for at least several years.

Zoom In Zoom Out Reset image size

We express a linear regression line (y = b × x + a) in flux–flux space (f_ν(s) − f_ν(l) space) as

$$\begin{eqnarray} f_{\nu }(l)&=&b(s,l)\times f_{\nu }(s)+a(s,l)\\ (s&=&u,g,r,i ; \ l=g,r,i,z) \nonumber \end{eqnarray} \tag{ 1 }$$

where "s" ("short" wavelength, x-axis in flux–flux space), and "l" ("long" wavelength, y-axis in flux–flux space) indicate two of the five photometric bands, whose average wavelengths are λ(s) < λ(l) (where λ(u), λ(g), λ(r), λ(i), λ(z) = 3551, 4686, 6165, 7481, 8931 Å).⁸ Figure 3 shows an illustration of our method. When we fit a regression line in flux–flux space, the regression slope b(s, l) and the intercept a(s, l) are expressed as

$$\begin{eqnarray} b(s,l)&=&\frac{\Delta f_{\nu }^{{\rm var}}(l)}{\Delta f_{\nu }^{{\rm var}}(s)} \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} a(s,l)&=&\bar{f}_{\nu }(l)-b(s,l)\bar{f}_{\nu }(s) \end{eqnarray} \tag{ 3 }$$

where $\Delta f_{\nu }^{{\rm var}}$ is the varying broad-band flux range (flux difference) during the observations, and $\bar{f}_{\nu }$ is the baseline flux given as the mean flux of light curves for each bandpass. The effects of contaminations of baseline flux ($=\bar{f}_{\nu }$; e.g., host galaxy flux and time-averaged fluxes of emission lines) do not affect the regression slope (they only affect the regression intercepts), so the regression slope in the flux–flux space directly measures the ratio of the two broad-band fluxes of the flux difference spectrum of each quasar.

Zoom In Zoom Out Reset image size

Given the flux difference spectrum $f_{\nu }^{{\rm dif}}$ for a quasar, the regression slope in flux–flux space is expressed as (Richards et al. 2001),

$$\begin{equation} b(s,l)=\frac{\int _{}^{}f_{\nu }^{{\rm dif}}S_{\nu }(l)d\log \nu }{\int _{}^{}f_{\nu }^{{\rm dif}}S_{\nu }(s)d\log \nu }\times \left(\frac{\int _{}^{}S_{\nu }(l)d\log \nu }{\int _{}^{}S_{\nu }(s)d\log \nu }\right)^{-1} \end{equation} \tag{ 4 }$$

where S_ν is the throughput of the SDSS photometric system in each bandpass (Doi et al. 2010).

From Equation (4), the color of the flux difference spectrum m_s − m_l|_dif (e.g., r − i|_dif) in units of magnitude is expressed by b(s, l) as,

$$\begin{equation} m_s-m_l|_{{\rm dif}}=+2.5\log (b(s,l)). \end{equation} \tag{ 5 }$$

This means that the bluer color indicates a smaller value of Equation (5), and then a smaller value of b(s, l). In this way, the color of the flux difference spectrum is related directly to the observable b(s, l). In the present work, we refer to b(s, l) simply as "the color of the flux difference spectrum" and this can be converted to the conventional definition of color through Equation (5).

If we assume a power-law difference spectrum ($f_{\nu }^{{\rm dif}}\propto \nu ^{\alpha _{\nu }^{{\rm dif}}}$), we obtain (by Equation (4))

$$\begin{equation} b(s,l) \sim \left(\frac{\nu (l)}{\nu (s)} \right)^{\alpha _{\nu }^{{\rm dif}}}=\left(\frac{\lambda (s)}{\lambda (l)} \right)^{\alpha _{\nu }^{{\rm dif}}} \end{equation} \tag{ 6 }$$

where ν(s) and ν(l) indicate average frequencies for the SDSS photometric bands. Because we define b(s, l) as λ(l) > λ(s),

$$\begin{equation} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}b(s,l)>1 & ({\rm if} \ \alpha _{\nu }^{{\rm dif}}<0) \\ b(s,l)<1 & ({\rm if} \ \alpha _{\nu }^{{\rm dif}}>0). \end{array}\right. \end{equation} \tag{ 7 }$$

When the flux difference spectra of quasars have a power-law shape, the color of the flux difference spectra has no dependence on the quasar redshift because the wavelength-shift by cosmic expansion is canceled out in Equation (6). That is, b(s, l) with a power-law difference spectrum is constant as a function of redshift if $\alpha _{\nu }^{{\rm dif}}$ is fixed. Although the color of the difference spectrum is actually contaminated by BEL variability, we can infer $\alpha _{\nu }^{{\rm dif}}$ by investigating the mean redshift dependence of b(s, l) for sample quasars.

In flux–flux space, the "observed" color corresponds to the slope of a straight line that passes through the origin of the coordinates and an observed point at an epoch. The "observed" color is time-dependent if a(s, l) (the regression intercept in flux–flux space) is non-zero (Sakata et al. 2010, 2011; see also Appendix B).

We define the b^ratio(s, l) by taking the time-average of light curves in any combination of five photometric bands ("baseline flux" $\bar{f}_{\nu }(j), j = u,g,r,i,z$) and taking the ratio as (see Figure 3)

$$\begin{equation} b^{{\rm ratio}}(s,l)\equiv \frac{\bar{f}_{\nu }(l)}{\bar{f}_{\nu }(s)}. \end{equation} \tag{ 8 }$$

Equation (5), which is the (time-averaged) color m₁ − m₂ of the quasar spectrum (e.g., r − i) and is expressed by b^ratio(s, l) as

$$\begin{equation} m_s-m_l=+2.5\log (b^{{\rm ratio}}(s,l)). \end{equation} \tag{ 9 }$$

We refer to b^ratio(s, l) simply as "the color of the quasar spectrum" in this paper. One can convert this to the conventional definition of color through Equation (9).

b^ratio(s, l) is expressed for a quasar with a spectrum f_ν as

$$\begin{equation} b^{{\rm ratio}}(s,l)=\frac{\int _{}^{}f_{\nu }S_{\nu }(l)d\log \nu }{\int _{}^{}f_{\nu }S_{\nu }(s)d\log \nu }\times \left(\frac{\int _{}^{}S_{\nu }(l)d\log \nu }{\int _{}^{}S_{\nu }(s)d\log \nu }\right)^{-1}. \end{equation} \tag{ 10 }$$

If we assume a power-law spectrum ($f_{\nu } \propto \nu ^{\alpha _{\nu }}$), we obtain

$$\begin{equation} b^{{\rm ratio}}(s,l) \sim \left(\frac{\nu (l)}{\nu (s)} \right)^{\alpha _{\nu }}=\left(\frac{\lambda (s)}{\lambda (l)} \right)^{\alpha _{\nu }}. \end{equation} \tag{ 11 }$$

Because we define b^ratio(s, l) as λ(l) > λ(s),

$$\begin{equation} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}b^{{\rm ratio}}(s,l)>1 & ({\rm if} \ \alpha _{\nu }<0) \\ b^{{\rm ratio}}(s,l)<1 & ({\rm if} \ \alpha _{\nu }>0). \end{array}\right. \end{equation} \tag{ 12 }$$

Richards et al. (2001) investigated the causes of the observed features in the color as a function of redshift for 2625 SDSS quasars by comparing the color–redshift relation model including the effects of emission lines. They confirmed that the average of the observed color of quasars as a function of redshift is well represented by a power-law continuum using α_ν ∼ −0.5 (similar to the result of Vanden Berk et al. 2001), and identified the effects of emission lines for the broad-band colors as deviations from power-law colors.

Because we are able to obtain the color of the flux difference spectrum for quasars using the flux–flux correlation method, the same kind of analyses used by Richards et al. (2001) for the flux difference spectrum can be applied. Currently, there are only a few studies related to the flux difference spectrum (Wilhite et al. 2005; Bian et al. 2012; Guo & Gu 2013). Among these, we refer to that of Wilhite et al. (2005) for the quasar composite flux difference spectrum in this paper. They constructed the composite residual (composite flux differences spectrum) from 315 quasars observed twice by SDSS spectroscopy. Figure 4 shows the histogram of the rest-frame time-lag between the two spectroscopic observations in Wilhite's sample, with a histogram of the time span of the photometric observations in our sample. The median rest-frame time lag of Wilhite et al.'s sample is 123 days, which is about one-tenth of ours (median ∼1200 days). Aside from the difference in the rest-frame time-lag, we can compare our results with the (redshifted) composite flux difference spectrum by Wilhite et al. (2005) using Equation (4) on the color–redshift space. Note that Wilhite et al. (2005) mainly focused on the ratio spectrum (the ratio of the composite difference spectrum to the composite spectrum), so they did not apply the correction of Galactic extinction to their composite difference spectrum. We discuss the effect introduced by this non-correction in Section 4.2.

Zoom In Zoom Out Reset image size

Here, we show the illustration of the effect of an emission line contamination for b(s, l) and b^ratio(s, l) on the redshift dependence (see also Richards et al. 2001). We can treat b(s, l) and b^ratio(s, l) in the same way, so here we show the case of b(s, l). We consider the case of s = r and l = i (i.e., b(r, i)) for example. When an emission line is redshifted and enters r band, the r-band flux becomes larger, and b(r, i) becomes smaller (bluer) than the case of a single power-law flux difference spectrum. At higher z, when the emission line enters the i-band, b(r, i) becomes larger (redder). This behavior is shown in Figure 5, assuming the continuum and line variability as

$$\begin{eqnarray} f^{{\rm dif}}_{\lambda } &=& f^{{\rm dif}}_{\lambda }({\rm continuum})+f^{{\rm dif}}_{\lambda }({\rm line}) \nonumber \\ &=&\left(\frac{\lambda }{3000{\rm \, \mathring{\rm{A}}}} \right)^{\alpha ^{{\rm dif}}_{\lambda }}+\frac{1}{2}\exp \left(-\frac{(\lambda -\lambda _{{\rm rest}})^2}{2\sigma _{\lambda }^2}\right) \end{eqnarray} \tag{ 13 }$$

where σ_λ is the emission line width (here σ_λ = 40 Å), and λ_rest is the central wavelength of the line (here λ_rest = 3000 Å). $\alpha _{\lambda }^{{\rm dif}}$ is as shown in Figure 5. The different power-law indices $\alpha ^{{\rm dif}}_{\lambda }$ lead to the different constant levels of b(s, l) following Equation (6).

Zoom In Zoom Out Reset image size

Schmidt et al. (2012) pointed out that previous analyses usually focused on the color variability of the quasar in color–magnitude space and suffered from covariance between the color and magnitude uncertainties. They avoided these error correlations by analyzing the color variability in magnitude–magnitude space, and concluded that the color variability has a bluer-when-brighter trend.

Although the magnitude–magnitude correlation method can avoid error correlations, it still has another problem. Schmidt et al. (2012) fitted a straight line for each quasar's magnitude–magnitude plot, and used the best fit slope as the indicator of color variability. Note that regression slopes on magnitude–magnitude space highly depend on the contaminations of baseline flux (e.g., host galaxy flux and time-averaged fluxes of emission lines) (e.g., Woo et al. 2007; Hawkins 2003). The regression slope in magnitude–magnitude space can be expressed as

$$\begin{eqnarray*} &&b^{{\rm mag}}(s,l) \nonumber\\ &&\quad \approx \frac{\log [(\bar{f}_{\nu }(l)+\Delta f_{\nu }^{{\rm var}}(l)/2)/(\bar{f}_{\nu }(l)-\Delta f_{\nu }^{{\rm var}}(l)/2)]}{\log [(\bar{f}_{\nu }(s)+\Delta f_{\nu }^{{\rm var}}(s)/2)/(\bar{f}_{\nu }(s)-\Delta f_{\nu }^{{\rm var}}(s)/2)]}\nonumber \end{eqnarray*}$$

where the notation is the same as Equations (2) and (3). Unlike in flux–flux space (Equation (2)), the regression slope in magnitude–magnitude space is dependent not only on the flux difference $\Delta f^{{\rm var}}_{\nu }$, but also on the time-averaged flux levels $\bar{f}_{\nu }$ ("baseline flux," Equation (8)). Thus, b^mag(s, l) only tells whether the "observed" color variability is bluer-when-brighter (b^mag < 1), redder-when-brighter (b^mag > 1), or achromatic (b^mag = 1). To derive the flux difference spectrum from b^mag, one needs to employ some models or estimate the contamination $\bar{f}_{\nu }$. On the other hand, the flux–flux correlation method enables us not only to avoid error correlations, but also to derive the flux difference spectrum model independently (Equation (2)).

Another statistical method to derive the flux difference spectrum exists. Meusinger et al. (2011) investigated the rest-frame wavelength dependence of the quasar variability using the structure functions V^mag of five SDSS photometric bands as a function of rest-frame wavelength. However, the structure function analyses, which are done in units of magnitude, are not appropriate to study the spectral variability or the flux difference spectrum because the structure function can be expressed as

$$\begin{eqnarray*} V^{{\rm mag}} &\approx & [\Delta m]^2 \nonumber \\ &\approx & \left[ -2.5\log \left(\frac{\bar{f}_{\nu }(s)+\Delta f^{{\rm var}}_{\nu }(s)/2}{\bar{f}_{\nu }(s)-\Delta f^{{\rm var}}_{\nu }(s)/2} \right) \right]^2.\nonumber \end{eqnarray*}$$

Then again, this is dependent not only on the flux difference $\Delta f^{{\rm var}}_{\nu }$, but also on the time-averaged flux levels $\bar{f}_{\nu }$.

In summary:

• 1.  
  The flux–flux correlation method can avoid the error correlations that exist in magnitude–color or flux–color correlation methods as was pointed out by Schmidt et al. (2012).
• 2.  
  The regression slope b(s, l) derived by the flux–flux correlation of a quasar multi-band light curve is not affected by baseline fluxes (e.g., host galaxy flux or time-averaged fluxes of emission lines).
• 3.  
  Regression slopes represent the color of the flux difference spectrum for each quasar, so we are able to infer the (mean) flux difference spectrum using redshift dependence of the derived regression slopes.

The linear regression for the data that have intrinsic scatter (σ_int) and error bars on both axes (e.g., flux–flux space linear regression) is complex and there are several different methods (see Park et al. 2012; Feigelson & Jogesh Babu 2012; Cappellari et al. 2013, and references therein). We adopt the MPFITEXY IDL routine (Williams et al. 2010) to fit a linear relation to the flux–flux correlation for each quasar. According to Park et al. (2012), the FITEXY method generally provides the least-biased result compared to the other methods. The MPFITEXY routine depends on the MPFIT package (Markwardt 2009). This routine is able to cope with intrinsic scatter, which is automatically adjusted to ensure χ²/(degrees or freedom) ∼1 (see Tremaine et al. 2002; Novak et al. 2006; Park et al. 2012, for more detail and justifications).

We fit a straight line in flux–flux space for each band pair light curve, excluding the data points with a "Bad observations" flag (MacLeod et al. 2012) or fainter than limiting magnitudes (22.0, 22.2, 22.2, 21.3, and 20.5 mag for u, g, r, i, and z bands, respectively).⁹ We focused on the properties of long-term and multi-epoch-averaged variability, so we exclude quasars that have less than 20 photometric epochs. A linear relation y = b(x − x₀) + a' is assumed in the fitting procedure, where x₀ is taken to be the average of x values (i.e., $\bar{f}_{\nu }(s)$). The reference value x₀ is necessary to minimize the uncertainty in the estimate of a' and the correlation between a' and b (Tremaine et al. 2002). Regression intercept a is calculated as a = a' − bx₀. Linear regressions for several quasars did not converge because of the low signal-to-noise ratio (S/N). Thus, we eliminated them from our sample. Furthermore, we confirmed that all of the regression slopes with a negative value were the result of a bad fitting (because of the small variability or low S/N) and we also eliminated them from our sample. In this paper we examine the modal redshift trends, so these eliminations had little effect on our analyses.

As an example, we show the histogram of the number of observational epochs of the (r, i) band pair light curves without a "Bad observations" flag (MacLeod et al. 2012) for our quasar sample in Figure 6. Resulting sample sizes for each band pair are shown in Table 1. Redshift distribution is not significantly modified from Figure 1 by our sample selection.

Zoom In Zoom Out Reset image size

Figure 7 shows the regression slopes b(s, l) as a function of redshift for each quasar (dots). The curves in Figure 7 are:

• •  
  Solid and dashed curves. The modal color (b(s, l)) of the flux difference spectrum (solid) with 1σ outer quartiles (dashed), derived by linear regression analyses in flux–flux space (i.e., mode and 1σ outer quartiles of the dots in the figure).
• •  
  Dash–dotted and double-dot–dashed curves. The modal color (b^ratio(s, l)) of the time-averaged spectrum (dash-dotted) with 1σ outer quartiles (double-dot–dashed), given by Equation (8).

Zoom In Zoom Out Reset image size

The modal value is defined following Hopkins et al. (2004), and the bin size is taken to be 0.2 in redshift.

At high redshift, u-, g-, and r-band fluxes are affected by Lyα forests or damped Lyα systems, which absorb UV flux and decrease u-, g-, and r-band variability (i.e., regression slopes become larger at higher redshifts). This makes it difficult to estimate intrinsic quasar spectrum and spectral variability. The shortest transmission wavelengths for u, g, and r band are 3000 Å, 3700 Å, and 5400 Å, respectively, so the Lyα emission line (1216 Å) enters in these bands at z = 1.47, 2.04, 3.44. The arrows with the label "Lyα" in Figure 7 indicate the redshift ranges in which Lyα forests change the observed color. In the later sections, we focus on redshift ranges below z = 1.47, 2.04, and 3.44 for band pairs containing u, g, and r bands, respectively, for clarifying the discussion about intrinsic spectral variability.

The dash-dotted curves (b^ratio(s, l)) in Figure 7 are generally larger than 1. This is consistent with the result of the geometric mean composite spectrum for SDSS quasars presented by Vanden Berk et al. (2001); the spectral index of the composite spectrum $\alpha ^{{\rm com}}_{\nu }\sim -0.44 < 0$ ($\alpha ^{{\rm com}}_{\lambda }\sim -1.56$) indicates b(s, l) > 1 (Equation (12)). The upper thin solid horizontal line in each panel of Figure 7 indicates a power-law spectrum model with α_ν = −0.44 calculated by Equation (4).

On the other hand, the solid curves (b(s, l)) in Figure 7 are generally less than 1 which implies $\alpha _{\nu }^{{\rm dif}}>0$ (Equation (7)). The composite flux difference spectrum presented by Wilhite et al. (2005) was reported to have a spectral index $\alpha _{\nu }^{{\rm dif}} = 0$, which implies b(s, l) = 1 (Equation (6)). However, the solid curves in Figure 7 are less than 1 within the 1σ ranges (dashed curves). We attempt to employ the spectral index $\alpha _{\nu }^{{\rm dif}}=+1/3$ suggested by Sakata et al. (2010) (the lower thin, solid horizontal line in each panel of Figure 7). This spectral index seems to reproduce the mean value of the b(s, l) for all the band pairs well. This inconsistency between Wilhite et al.'s (2005) result and ours is discussed in Section 4.2.

b^ratio(s, l) > b(s, l) is the general trend for quasar variability as shown in Figure 7 (see also Appendix B). This indicates that the quasar composite flux difference spectrum is flatter (bluer) than the composite spectrum. In other words, quasar UV–optical spectra tend to become flatter (bluer) when in brighter phases. To clarify the "observed" bluer-when-brighter trend for individual quasars, we show the ratio of b(s, l) to b^ratio(s, l) given as

$$\begin{equation} c(s,l) \equiv \frac{b(s,l)}{b^{{\rm ratio}}(s,l)} \sim \left(\frac{\lambda (s)}{\lambda (l)} \right)^{\Delta \alpha _{\nu }} \end{equation} \tag{ 14 }$$

as a function of redshift in Figure 8, where $\Delta \alpha _{\nu }=\alpha ^{{\rm dif}}_{\nu }-\alpha _{\nu }$ (approximate value is from Equations (6) and (11)). By definition, as in Equations (7) and (12),

$$\begin{equation} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}c(s,l)>1 & ({\rm if} \ \alpha _{\nu }^{{\rm dif}}<\alpha _{\nu }) \\ c(s,l)<1 & ({\rm if} \ \alpha _{\nu }^{{\rm dif}}>\alpha _{\nu }). \end{array}\right. \end{equation} \tag{ 15 }$$

Dashed horizontal lines in Figure 8 indicate $\Delta \alpha _{\nu }=\alpha _{\nu }^{{\rm dif}}-\alpha _{\nu } =-\alpha _{\lambda }^{{\rm dif}}+\alpha _{\lambda }=7/3 -1.56 \sim 0.77$ (difference of the spectral indices shown in Figure 7), and solid horizontal lines indicate the case of $\alpha _{\nu }^{{\rm dif}}=\alpha _{\nu }$. Figure 8 shows that almost all of quasars have $\alpha _{\nu }^{{\rm dif}}>\alpha _{\nu }$. This means the "observed" bluer-when-brighter spectral variability is a common trend in quasars.

Zoom In Zoom Out Reset image size

Besides the continuum (power-law) variability, we are able to see several spectral features in Figure 7 and also in Figure 8. These features are mainly made up of the emission line contributions (e.g., Wilhite et al. 2005). It is known that the emission line variability is smaller compared with the continuum variability (intrinsic Baldwin effect; Wilhite et al. 2005). As expected from the intrinsic Baldwin effect, the solid curves in Figure 7 show less features compared to the dash-dotted curves. Nonetheless, we are still able to identify the coincidence and the difference of features (the redshift dependence) between the solid and dash-dotted curves. The most apparent difference between them, which corresponds to the difference between the composite and composite difference spectrum, is the bump in intermediate redshifts (e.g., the redshift range z ∼ 1.0–1.5 of g − r band pair in Figure 7). This difference makes several sharp transitions on c(s, l) as a function of redshift shown in Figure 8. The peaks of the bump of dash-dotted curves correspond to the spectral region of Mg ii (2800 Å) emission line and Fe ii pseudo-continuum, and this peak is not seen in the solid curves. This implies that the variability of these emissions are weak. This point is discussed in Sections 5 and 6.

There have been several statistical studies about spectral shape and spectral variability for quasars selected and spectroscopically confirmed within SDSS. Here we show the consistency of these previous spectroscopic results with our photometric results.

Figure 9 shows the comparison of our results (modal color as a function of redshift in Figure 7) with the composite (Vanden Berk et al. 2001) and composite difference (Wilhite et al. 2005) spectrum using Equations (10) and (4). In Figure 9 we plot not only the modal value, but also the (weighted) mean value as a function of redshift to check the effects of the asymmetric color distribution, which has a red wing at a given redshift due to intrinsic dust reddening (Richards et al. 2003). We are able to see in Figure 9 that the mean values are slightly larger (redder) than the modal value, but the difference is small. This indicates that our results are not affected heavily by this method of taking the average.

Zoom In Zoom Out Reset image size

The composite spectra composed by taking the arithmetic or geometric means (Vanden Berk et al. 2001; Wilhite et al. 2005) are compared with our result through Equations (4) and (10). The geometric mean preserves the average power-law slope, insofar as quasar spectra can be accurately described by a power-law. On the other hand, the arithmetic mean retains the relative strength of the non-power-law features, such as emission lines (Vanden Berk et al. 2001).

As shown in the top panel of Figure 9, b^ratio(r, i) of our result and of the composite spectrum (Vanden Berk et al. 2001) are generally consistent. Vanden Berk et al. (2001) presented the geometric and arithmetic composite spectra, but their difference is very small (their spectral index is α_ν = −0.44 and −0.46, respectively, with the uncertainty ∼0.1 due to the spectrophotometric calibration) and they have almost the same b^ratio(r, i)–redshift relation. A power-law index α_ν = −0.44 (shown as dotted line in upper panel of Figure 9) is consistent with b^ratio(r, i), and as we have shown in Figure 7, all b^ratio(s, l) are also consistent with α_ν = −0.44.

On the other hand, b(r, i) of the regression results and the geometric composite difference spectrum (Wilhite et al. 2005) (middle panel of Figure 9) seem to be inconsistent in that the b(r, i) of the regression results is generally lower than that of the composite difference spectrum for any redshift. This indicates that their composite geometric difference spectrum is steeper than our sample mean color of the difference spectra. We note that their difference spectra, composed by taking the geometric mean and the arithmetic mean, are largely different. The authors did not specify the reasons for this large discrepancy. They rely on the geometric composite spectrum when inferring the continuum spectral index $\alpha _{\nu }^{{\rm dif}}=0.0$, and we should compare our result with the geometric one when inferring the continuum variability. Also, we note that their spectra are not corrected for Galactic extinction (see Abazajian et al. 2004; Wilhite et al. 2005). The composite spectra with no correction for Galactic extinction obtained by Wilhite et al. (2005) has the spectral index $\alpha _{\lambda }^{{\rm com}} = -1.35$, which is steeper (redder) than the dereddened composite spectrum by Vanden Berk et al. (2001) ($\alpha _{\lambda }^{{\rm com}} = -1.56$). So, the composite flux difference spectrum with $\alpha _{\lambda }^{{\rm dif}} = -2.0$, with no correction for Galactic extinction obtained by Wilhite et al. (2005), means that the intrinsic composite flux difference spectrum should have a flatter (bluer) spectral index (i.e., $\alpha _{\lambda }^{{\rm dif}} < -2.0$), which is consistent with our result $\alpha _{\lambda }^{{\rm dif}} \sim -7/3$ or $\alpha _{\nu }^{{\rm dif}} \sim +1/3$. In the bottom panel of Figure 9 we show b(r, i) with no correction for Galactic extinction as a function of redshift. In this plot, although bumpy features in the b(r, i)–redshift relation (which correspond to the emission lines variability and are discussed in Sections 5 and 6) of the regression results are still not in good agreement with those of the geometric composite difference spectrum (Wilhite et al. 2005), it is clear that our result and that of the geometric composite difference spectrum (Wilhite et al. 2005) are approximately consistent in that they have similar average power-law values. Thus, a large fraction of the discrepancy between our result and that of Wilhite et al. (2005) seen in the middle panel of Figure 9 can be attributed to a correction or no correction for the Galactic extinction. It also indicates that the composite flux difference spectrum is not so different in spite of the large difference in time scale between ours and that of Wilhite et al. (2005; see Figure 4).

Hereafter, we refer to the approximate form of the power-law continuum flux difference spectrum as

$$\begin{equation} f_{\nu }^{{\rm dif}} \propto \nu ^{+1/3} \end{equation} \tag{ 16 }$$

(i.e., $\alpha _{\nu }^{{\rm dif}}\sim +1/3$). This power-law index is well known as the standard disk model prediction (Lynden-Bell 1969; Shakura & Sunyaev 1973). Tomita et al. (2006) and Lira et al. (2011) concluded that the flux difference spectra in near-infrared (NIR) wavelengths for their AGN samples are well represented by a combination of the dusty torus component (blackbody spectrum) and the accretion disk component ($\alpha _{\nu }^{{\rm dif}}\sim +1/3$) on the assumption that the flux difference spectrum conserves the spectral shape (spectral index) of the underlying accretion disk spectrum (see also Paltani & Walter 1996; Collier et al. 1999). Note that the observations of the polarized flux spectra of AGNs—which is assumed to be produced by electron scattering interior to the dust sublimation radius and is regarded as a copy of the spectrum originating in the accretion disk—have also revealed the underlying continuum spectral shape consistent with α_ν = +1/3 (α_ν = +0.44 ± 0.11, Kishimoto et al. 2008). Our result can be considered to be a UV–optical extension of the previous results on the flux difference spectrum obtained in the optical–NIR wavelengths (Tomita et al. 2006; Lira et al. 2011) because the dust in the nuclear region cannot become higher temperature than the dust sublimation temperature (T_sub ∼ 1500 K) and the accretion disk component dominate in UV–optical wavelengths.

In Section 4.3, we discuss the validity of the accretion disk model in UV–optical wavelengths.

Pereyra et al. (2006) fitted the composite flux difference spectrum presented by Wilhite et al. (2005) with a standard thermal accretion disk model (Shakura & Sunyaev 1973) in which the luminosity varies following the changes in accretion rate from one epoch to the next, and concluded that most of the UV–optical variability observed in quasars may be due to processes involving changes in mass accretion rates. Sakata et al. (2011) investigated this model for 10 SDSS quasars individually, and also concluded that the model can explain the UV–optical variability quantitatively. Definitions of the model adopted by Pereyra et al. (2006) and Sakata et al. (2011) are slightly different in that the former depends on the first order Taylor expansion about the average characteristic temperature T* (see Pereyra et al. 2006), but the predicted flux difference spectra are identical for the modest variation amplitude, which is valid for most quasars.

In the formalization of Pereyra et al. (2006), the model flux difference spectrum has two free parameters: a normalization constant and the average characteristic temperature T*. One of the advantages of focusing on the color–redshift relation is that it does not depend on the inclination angle (including limb-darkening effect, e.g., Phillips & Meszaros 1986). When inferring the color of the flux difference spectrum (observable b(s, l) in the present work), this model requires only one parameter T* because the normalization constant does not affect the color of the difference spectrum. T* is defined as

$$\begin{eqnarray} T^*&=&\left(\frac{3GM_{{\rm BH}}\dot{M}}{8\pi \sigma _S (3R_S)^{3} }\right)^{1/4}\qquad \qquad \qquad \, \, \qquad \qquad \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} &=&90{,}600\, {\rm K} \left(\frac{\epsilon }{1/12}\right)^{-1/4}\left(\frac{M_{{\rm BH}}}{10^{9}\,M_{\odot }} \right)^{-1/4}\left(\frac{{\rm ER}}{0.1} \right)^{1/4} \end{eqnarray} \tag{ 18 }$$

where M_BH is the black hole mass, $\dot{M}$ is the mass accretion rate, σ_S is the Stefan–Boltzmann constant, and K denotes Kelvin. The disk inner edge is assumed to be at the inner stable circular orbit (R_in = 3R_S, where R_S is Schwarzschild radius). The Eddington luminosity is defined as L_Edd = 4πGM_BHm_pc/σ_T (where m_p is the proton mass and σ_T is the Thomson cross section) and the Eddington ratio, denoted "ER," is given as ER ≡ L_bol/L_Edd. Radiative efficiency , defined as $L_{{\rm bol}}=\epsilon \dot{M}c^2$, is fixed¹⁰ to = 1/12. Changes in $\dot{M}$ lead to changes in T* as ${T}^*\propto \dot{M}^{1/4}$ and then to changes in flux. Pereyra et al. (2006) fitted their continuum variability model to the geometric mean flux difference spectrum of Wilhite et al. (2005), and obtained T* = 92, 700 K. The upper panel of Figure 10 shows the characteristic temperature of our quasar samples calculated by Equation (18) using the fiducial virial black hole mass and the Eddington ratio listed in the quasar property catalog (Shen et al. 2011). T* = 92, 700 K, obtained by Pereyra et al. (2006), is very reasonable for SDSS quasars.

Zoom In Zoom Out Reset image size

This model predicts $\alpha _{\nu }^{{\rm dif}}=1/3$ around optical to near-infrared wavelengths (e.g., Sakata et al. 2010). This spectral index is well known as the long wavelength limit for a standard thermal accretion disk model (Shakura & Sunyaev 1973). At the long wavelength limit, the accretion disk model spectrum with varying accretion rate turns out to be simple scaling, thus the flux difference spectrum predicted in this model also has $\alpha _{\nu }^{{\rm dif}}\sim +1/3$. Because we obtained $\alpha _{\nu }^{{\rm dif}}\sim +1/3$ as shown in Figure 7, this model is approximately valid. However, we should be careful because:

• 1.  
  The long wavelength limit approximation is not valid for the wavelength range of SDSS photometry for high-redshift quasars. In UV wavelengths, $f_{\nu }^{{\rm dif}}$ predicted by this model becomes gradually steeper (redder) because the UV spectral range corresponds to the turnover spectral region of the blackbody radiation from the accretion disk inner boundary at the inner stable circular orbit (i.e., the hottest part of the disk), as seen in the lower panel of Figure 10.
• 2.  
  Pereyra et al. (2006) used the composite difference spectrum not corrected for Galactic extinction (Wilhite et al. 2005), and the effect of the non-correction is not investigated in their paper. The same is true for Sakata et al. (2011), in which the correction for the Galactic extinction was not properly applied.¹¹ It is possible that the non-correction for Galactic extinction leads to the wrong conclusion.

So, it is worth comparing our (Galactic-extinction corrected) result with the model prediction to clarify the validity of the standard disk model. We compare our results with the model proposed by Pereyra et al. (2006) in b(r, i)–redshift space as shown in Figure 11. In Figure 11, Pereyra's model (with their best-fit value of T* = 92,700 K) is actually in agreement with Wilhite et al.'s composite flux difference spectrum. However, the model prediction is apparently steeper (b(r, i) is larger) than our result. This discrepancy is mainly due to whether or not Galactic extinction is corrected (see the middle and bottom panels in Figure 9). In Figure 11, we also plot the model prediction range in b(s, l)–redshift space, which is from T* = 50, 000 K (upper boundary) to 200, 000 K (lower boundary). The model prediction is apparently not consistent with the b(s, l)–redshift region obtained in our analyses in that the model predicts larger b(s, l) values and larger spread in b(s, l)–redshift space than observed values, particularly in the high redshift range.

Zoom In Zoom Out Reset image size

These indicate that the intrinsic composite flux difference spectrum of quasars is actually flatter (bluer) than the prediction of the accretion disk model, particularly in the UV wavelength range. In other words, we cannot see the model-predicted "UV turnover" feature (e.g., Shang et al. 2005) in the flux difference spectrum. This result is qualitatively consistent with Schmidt et al. (2012), but our method enables us to examine the difference between observed variability and the model prediction quantitatively.

We are able to make the model flux difference spectrum flatter (up to $f_{\nu }^{{\rm dif}}\propto \nu ^{1/3}$ at all UV–optical wavelengths) if we assume the disk characteristic temperature to be, for example, T* > 200, 000 K, but this is not valid for most of our quasar samples as shown in the upper panel of Figure 10. Beyond the standard accretion disk model (Shakura & Sunyaev 1973), several more sophisticated accretion disk models are proposed (see, e.g., Hubeny et al. 2000; Davis et al. 2007). However, accretion disk model spectra based on non-LTE atmosphere calculations (Davis et al. 2007) usually have redder color than the standard thermal accretion disk model because of opacity effects (e.g., Hubeny et al. 2000) and this leads to the redder color of the flux difference spectrum as shown in Schmidt et al. (2012). Although radiative transfer calculations for a disk+corona configuration (rather than those for a bare disk above) result in a weaker "UV turnover" (Kawaguchi et al. 2001), the expected spectra are redder than the long wavelength limit of the standard disk as well. In addition, the accretion disk of quasars or luminous AGNs is thought to be seen almost face-on, so the relativistic effects for the color of the observed spectrum is probably negligible (Hubeny et al. 2000). In short, it is difficult to explain the flux difference spectrum with $\alpha _{\nu }^{{\rm dif}} \sim 1/3$ by any of the existing accretion disk models.

However, it is true that the reasonable range of physical parameters related to the accretion disk model (M_BH and Eddington ratio) can qualitatively explain the absolute continuum flux and flux variation amplitude as indicated in Sakata et al. (2011) and Gu & Li (2013), and $\alpha _{\nu }^{{\rm dif}}$ obtained here can be naturally attributed to the well-known long wavelength limit value predicted in standard accretion disk models. Thus, we can only conclude that the continuum spectral variability in quasars cannot be explained by the accretion disk model with varying mass accretion rate, but it is probably related to the thermal accretion disk itself. Finally, we note that it seems to be difficult to explain the large coherent flux variation within the UV–optical wavelength range (e.g., a strong linear flux–flux correlation in each individual quasar) by the accretion disk local fluctuations or hot spot scenarios (e.g., Dexter & Agol 2011; Schmidt et al. 2012; Meusinger & Weiss 2013).

In summary, our conclusions in this section are:

• 1.  
  We show that the continuum component of the flux difference spectra is (in average) well approximated by a power-law shape with $\alpha _{\nu }^{{\rm dif}}\sim +1/3$, and it is bluer than the spectroscopic composite flux difference spectrum ($\alpha _{\nu }^{{\rm dif}}=0.0$) obtained by Wilhite et al. (2005). We conclude that the discrepancy is due to the Galactic extinction; the light curves used in our analyses are corrected for Galactic extinction but the spectra used in Wilhite et al. (2005) are not. If we derive the color of the flux difference spectrum without the correction, our result shows consistent color and rest-frame wavelength dependence of the flux difference spectrum with Wilhite et al. (2005), in spite of the large difference in time scale (Figure 4).
• 2.  
  We compare the b(s, l)–redshift relation with "the standard disk with varying mass accretion rate model" (Pereyra et al. 2006). One of the advantages of our method is that it is based on the direct comparison of observable color of the difference spectrum (b(s, l)) with the model predicted difference spectrum, and does not require additional models or assumptions for non-variable spectral components. We concluded that the flux difference spectrum is flatter (bluer) than the model prediction. In particular, we confirm that the model-predicted "UV turnover" is not seen in the flux difference spectrum. This is surprising because not only the standard accretion disk model (Shakura & Sunyaev 1973) but also any accretion disk models cannot produce such a flat spectrum as $\alpha _{\nu }^{{\rm dif}}\sim +1/3$. This result is qualitatively consistent with Schmidt et al. (2012), but our method enables us to examine the difference between observed variability and the model prediction quantitatively.

As seen in Figure 7 and discussed in the previous section, the most apparent difference between b^ratio(s, l) (color of the time-averaged spectrum) and b(s, l) (color of the flux difference spectrum) exists in features of the Fe ii pseudo-continuum and the Mg ii emission line. The spectral region containing the Fe ii pseudo-continuum, Balmer continuum emission, and the Mg ii emission line is called SBB (Wills et al. 1985). Note that in the actual quasar spectra, higher order Balmer lines are merging to a pseudo-continuum and yielding a smooth rise to the Balmer edge (Dietrich et al. 2003), which is also contained in SBB. Here we discuss the SBB variability.

To clarify the difference between b^ratio(s, l) and b(s, l) (as a function of redshift) arising from the complex emission line composition, we decompose the SBB in the composite spectrum (Vanden Berk et al. 2001) into power-law continuum, Mg ii emission line, Fe ii pseudo-continuum, and Balmer continuum and then investigate the color (corresponding to b^ratio(s, l) and b(s, l)) as a function of redshift for the composite spectra with and without some SBB components.

Kurk et al. (2007) decomposed the Vanden Berk et al. composite spectrum into a power-law continuum, Balmer continuum, and residual. They confirmed that the residual was well described by the Fe ii emission template (Vestergaard & Wilkes 2001) and a Mg ii emission line. We follow their decomposition procedure, which is done in the 2200–3646 Å region. First, the power-law continuum with α_λ = −1.54 is normalized to the values in the 4200–4230 Å region. Second, the Balmer continuum is created by Equation (A2) and normalized to the value at 3675 Å. The residual spectrum (the composite spectrum minus a power-law minus Balmer continuum) is composed of Fe ii pseudo-continuum and a Mg ii emission line. The Mg ii emission line is decomposed by simply assuming the Fe ii pseudo-continuum in 2700 Å to 2900 Å region is described by linear interpolation of the residual spectrum from 2700 Å to 2900 Å.

Using the decomposed components, we define four spectra of the SBB region f_λ(SBB) as:

• •  
  VdB: the composite spectrum of Vanden Berk et al. (2001).

  $$\begin{equation} f_{\lambda }({\rm SBB})=f_{\lambda }({\rm PL})+f_{\lambda }({\rm BaC})+f_{\lambda }({\rm Fe\,{\scriptsize {II}}})+f_{\lambda }({\rm Mg\,{\scriptsize {II}}}) \end{equation} \tag{ 21 }$$

  where "PL" denotes a power-law component.
• •  
  VdB w/o Fe ii: VdB minus Fe ii pseudo-continuum

  $$\begin{equation} f_{\lambda }({\rm SBB})=f_{\lambda }({\rm PL})+f_{\lambda }({\rm BaC})+f_{\lambda }({\rm Mg\,{\scriptsize {II}}}). \end{equation} \tag{ 22 }$$
• •  
  VdB w/o Mg ii: VdB minus Mg ii emission line

  $$\begin{equation} f_{\lambda }({\rm SBB})=f_{\lambda }({\rm PL})+f_{\lambda }({\rm BaC})+f_{\lambda }({\rm Fe\,{\scriptsize {II}}}). \end{equation} \tag{ 23 }$$
• •  
  VdB w/o Fe ii, Mg ii: VdB minus Fe ii pseudo-continuum minus Mg ii emission line

  $$\begin{equation} f_{\lambda }({\rm SBB})=f_{\lambda }({\rm PL})+f_{\lambda }({\rm BaC}). \end{equation} \tag{ 24 }$$

These spectra are shown in the upper panel of Figure 15. Then, we calculate the color of these spectra by Equation (4) or (10) as a function of redshift. Results are shown in the bottom panel of Figure 15. We can see that the Balmer continuum, Fe ii pseudo-continuum, and the Mg ii emission line make distinguishable bump features at different redshift ranges. The peak in b(r, i)[VdB] at z ∼ 1.5–1.8 are made of a Mg ii emission line (z ∼ 1.5) and Fe ii pseudo-continuum (z ∼ 1.8) (in agreement with Richards et al. 2001), and this feature is not seen in b(r, i)[mode]. The b(r, i)[mode] curve resembles the case of "VdB w/o Fe ii, Mg ii," (i.e., f_λ(SBB) = f_λ(PL) + f_λ(BaC)). This indicates that the variability of Mg ii emission line and Fe ii pseudo-continuum are relatively weaker compared to the Balmer continuum emission.

Zoom In Zoom Out Reset image size

Mg ii line variability has been examined in several papers (see Woo 2008; Hryniewicz et al. 2013 and references therein). Goad et al. (1999a, 1999b) showed that continuum flux and the broad high-ionization emission lines of Lyα, C iv, N v(1240 Å), He ii(1640 Å), and the semi-forbidden lines of the Si iii]+C iii] blend varied significantly (a factor of five at 1365 Å and three at 2230 Å for continuum, and a factor of two for Lyα and C iv), while the low-ionization lines of Mg ii and Fe ii do not vary significantly (less than 7%) during 11 months of observations for the Seyfert 1 galaxy NGC 3516. Goad et al. (1999a) listed three possible physical reasons explaining the absence of a response in the Mg ii emission line: (1) the continuum band responsible for the Mg ii emission line (i.e., 600 to 800 eV, Krolik & Kallman 1988) is invariant, (2) the emission line is insensitive to continuum variations, and (3) the line-emitting region is physically extended and thus has yet to respond to the observed continuum variations. Situation (1) was disproved in Goad et al. (1999a) by observational evidence that showed that the amplitude of the continuum variations generally increases toward shorter wavelengths, at least up to EUV wavelength range for several AGNs (e.g., Marshall et al. 1997; Cagnoni et al. 1999; Halpern et al. 2003). Strong soft X-ray variability has been observed in many AGNs, so (1) is definitely not true. Situation (3) is not true, especially in our sample, with three to four years (in quasar rest-frame) observational duration. Thus, we can expect situation (2) to be responsible for an absence of variability in the Mg ii emission line. Goad et al. (1993) defined the line responsivity η as a linear coefficient relating line emissivity with ionization parameter, assuming that locally the emission-line gas responds linearly to continuum variations (Goad et al. 1999a; Korista & Goad 2004; Goad et al. 2012). The larger η is, the more responsive the line is to continuum variations. Although η is model dependent, high-ionization emission lines generally have η ∼ 1, Balmer lines have η ∼ 0.6, and Mg ii emission line has η ∼ 0.2 (Goad et al. 1999a; Korista & Goad 2004). We are able to interpret our result with respect to the weak variability of Mg ii as the consequence of the small responsivity η for the low-ionization lines.

The Fe ii emission line is also a member of the low-ionization lines. In the same way as the Mg ii emission line, the weak variability of Fe ii can be the consequence of the small responsivity η. However, the situation may be more complicated with regard to Fe ii.

The origin of Fe ii emissions is still controversial (e.g., Dong et al. 2011). Baldwin et al. (2004) showed that Fe ii emissions from a photoionized broad emission line region (BELR) cannot explain all the Fe ii emission features, and suggested that a large amount of Fe ii emission originates in either a separate collisionally ionized component or BELR with considerable velocity structure (microturbulence). The collisionally ionized component would not vary in response to continuum variations, thus Fe ii pseudo-continuum variability would be weak as a whole. On the other hand, several authors suggest that Fe ii is emitted in infalling photoionized gas (Hu et al. 2008; Ferland et al. 2009), and much recent observational evidence supports this picture (Barth et al. 2013, and references therein). When this is the case, the small amplitude of Fe ii variability can be explained as the consequence of the different origin of Fe ii, (i.e., Fe ii emission does not come from usual BELR, so it cannot vary in the same way as other BELs). However, more detailed investigation of the Fe ii emission is beyond the scope of this paper.

Figure 15 indicates that the Balmer continuum is a dominant variable component in SBB region. We can expect that the Balmer continuum and Balmer emission lines are formed coherently by some photoionization processes in the broad line region, as we can see in the fourth global eigenspectrum presented by Yip et al. (2004). There is no wonder that the Balmer continuum is variable like the Balmer emission lines (e.g., Hα and Hβ). The variability of the Balmer continuum has been discussed observationally in Wamsteker et al. (1990) and Maoz et al. (1993), and theoretically in Korista & Goad (2001). Our result is the first study clearly showing the variability of the Balmer continuum emission in large quasar samples.

Wilhite et al. (2005) showed that the ratio spectrum of the composite difference spectrum to the composite spectrum was depressed around 3000 Å (see also Bian et al. 2012). Photometric variability using structure function formalization (Vanden Berk et al. 2004; MacLeod et al. 2012) also shows the depression around 3000 Å. This behavior can be interpreted by our result concerning the variability and non-variability of SBB components. Because the Fe ii pseudo-continuum (∼2200–3000 Å) and Mg ii emission line (∼2800 Å) exist in the composite spectrum, although not in the composite difference spectrum, the ratio is significantly depressed. Moreover, as discussed in Bian et al. (2012), the Balmer continuum should obey the intrinsic Baldwin effect in the same way as Balmer emission lines, so the ratio spectrum around the Balmer continuum region should be slightly depressed.

In summary, our conclusions in this section are:

• 1.  
  We first decompose the SBB spectral region of the high resolution SDSS quasar composite spectrum (Vanden Berk et al. 2001) into a power-law continuum, Balmer continuum, Fe ii, and Mg ii following Kurk et al. (2007). Then, we investigate the color–redshift relation of the composite spectrum with and without each component. We confirm that the composite spectrum without Fe ii and Mg ii emission line components has similar color–redshift relations with the color of the flux difference spectrum. This indicates that Fe ii and Mg ii have (on average) small variability, and the Balmer continuum is variable in broad spectral region. This is the first time that the strong variability of Balmer continuum and small variability of Fe ii and Mg ii are confirmed as the statistical property of large quasar samples.
• 2.  
  Small variability of the Mg ii emission line is naturally interpreted by the small responsivity η of the low-ionization line to the ionizing continuum variation (Goad et al. 1999a).
• 3.  
  Small variability of the Fe ii emission line is also interpreted by the small responsivity η of the low-ionization line, but the origin of the Fe ii emission is controversial, so it is possible that the small variability of Fe ii is due to the different origin of Fe ii from normal photoionization BELR.
• 4.  
  It is no wonder that the Balmer continuum is variable like the other Balmer emission lines because we can expect that they are formed coherently by some photoionization processes in the broad line region.