Extreme Variability Quasars in Their Various States. I: The sample selection and composite SDSS spectra

We start from the Sloan Digital Sky Surveys (SDSS) date release 14 (DR14) quasar catalog (DR14Q, Pâris et al., 2018), which contains 526,356 spectroscopically confirmed quasars with luminosity $M_{i}[z=2]<-20.5$ over 9376 $\rm deg^{2}$. Note the DR14Q catalog only provides single epoch photometry (i.e., the primary SDSS magnitude) for each source. In order to select EVQs, we construct the $g$- and $r$-band light curves for each quasar through further gathering all archive photometric observations from both SDSS¹¹1https://skyserver.sdss.org/casjobs/ and Pan-STARRS1²²2https://mastweb.stsci.edu/ps1casjobs/home.aspx (PS1) databases.

The SDSS photometric observations were taken by the drift-scan camera (30 2k $\times$ 2k CCDs) installed on the 2.5m Sloan telescope (Gunn et al., 2006). In the SDSS DR14Q catalog, the primary photometric observations were obtained with SDSS-I/II, lasting from 2000 to 2007 (covering 11,663 $\rm deg^{2}$), and SDSS-III (before 2009) on $\rm\sim 3000\ deg^{2}$ new sky area. Using a matching radius of 1″, we gather all available SDSS $g$-and $r$-band photometry for DR14Q quasars. Referring to the recommendation on the use of photometric processing flags from SDSS³³3https://www.sdss.org/dr16/algorithms/photo_flags_recommend/, we reserve detections with $\rm{mode=1}$ or 2, “clean” flags and PSF magnitude error $<$ 0.2 mag in both $g$- and $r$-band. About 45% quasars (237,099 quasars, specifically) have multi-epoch SDSS photometry (and 28,672 of them had been observed over more than five epochs).

We also collect the $g$- and $r$-band photometry from the PS1 $3{\rm\pi}$ survey, with up to four exposures in each band per year, conducted from 2009 to 2013. In total, within a matching radius of 1″, 510,838 SDSS DR14 quasars have counterparts in the PS1 archive (with mean epochs of around 10 and 12 in $g$- and $r$-band, respectively). Similar to the processing of SDSS data, we also filter the matched PS1 detections according to their photometric info flags. We rule out PS1 detections with PSF magnitude error $>$ 0.2 mag or flagged as:

1. 1.

   Peak lands on diffraction spike, ghost or glint;

2. 2.

   Poor moments for small radius, try large radius or could not measure the moments;

3. 3.

   Source fit failed or succeeds but with low S/N, high Chi square or too large for PSF;

4. 4.

   Source model peak is above saturation;

5. 5.

   Size could not be determined;

6. 6.

   Source has crNsigma above limit;

7. 7.

   Source is thought to be a defect;

8. 8.

   Failed to get good estimate of object’s PSF;

9. 9.

   Detection is astrometry outlier.

After conversion between PS1 and SDSS photometry (see §2.2), we could merge photometric data points from both SDSS and PS1, and obtain for every quasar relatively long-term $g-$ & $r-$ band light curve with a length of $4\sim 15$ years. We further remove from the light curves a small fraction of photometric measurements which could be unreliable according to the consistency check between $g$- and $r$-band light curves (see §2.3).

A mean magnitude $g_{\rm mean}$ is then derived from the clean light curve for each quasar. We need such a mean magnitude to represent the average brightness of a quasar over a long duration, and to be compared with the synthetic magnitude $g_{\rm spec}$ derived from the SDSS spectrum to determine if the spectrum was captured during a bright, normal or dim state. In order to avoid the mean magnitude being over-dominated by a few measurements with very small uncertainties, we simply calculate the mean without error weighting. The light curves of most quasars contain far more data points from PS1 than from SDSS (and contrarily for quasars in SDSS Stripe 82). In order to avoid the mean magnitudes being over-dominated by data from one instrument or by data from a single observing season with a large number of epochs, we first calculate the yearly mean magnitudes and then the final mean from the yearly mean values.

In Fig. 1 we plot the distributions of $\Delta g_{\rm fix}$ =$g_{\rm PS1}-g_{\rm SDSS}$ (and $\Delta r_{\rm fix}$ =$r_{\rm PS1}-r_{\rm SDSS}$) versus redshift for SDSS-PS1 matched quasars, where $g_{\rm PS1}$ ($r_{\rm PS1}$) and $g_{\rm SDSS}$ ($r_{\rm SDSS}$) are the mean PS1 magnitude and the primary SDSS magnitude, respectively. Subtle but clear offsets between PS1 and SDSS magnitudes (see the dark-blue lines which plot the mean offsets between PS1 and SDSS photometries averaged with a redshift bin of $\Delta z=0.05$) are seen within certain redshift ranges. This is owing to the slight difference of the filter transmission curves between PS1 and SDSS. For example, given that the red side transmission of the PS1 $g$-band is considerably higher than that of SDSS $g$ at $5200\sim 5500$Å, where locates the Mg II emission line for quasars at $z\sim 0.9$, PS1 $g$-band would receive more Mg II emission line photons and result in a brighter magnitude compared with SDSS $g$. Moreover, the general decreasing of the $\Delta g_{\rm fix}$ beyond $z\sim 2.6$ could be due to the fact that the spectral slope of quasars is redder at the wavelength shorter than 1216Å (intrinsically and/or because of intergalactic medium absorption), which makes the energy distribution of the photons that $g$-band received concentrates to the red side and magnifies the effect of filter transmission difference near $5200\sim 5500$ Å we mentioned before. The dependence on the redshift of the mean offset could indeed be well recovered through comparing the synthetic PS1 and SDSS magnitudes of the mean quasar spectrum of Yip et al. (2004) (the dot-dashed lines in Fig. 1). Note the mean spectrum extends down to 900 Å  in the rest frame thus results in a cut at $z\sim 3.4$ (for $g$ band).

To eliminate the systematic offsets between PS1 and SDSS magnitudes for our quasars, a correction is applied to the PS1 magnitudes of each quasar using the mean $\Delta g_{\rm fix}$ ($\Delta r_{\rm fix}$) of its closest 1000 neighbors in redshift space.

A commonly used criterion for selecting EVQs is to ask for a change in magnitude of $|\Delta g|>1$ mag between any two epochs in the light curve (e.g., MacLeod et al., 2016; Rumbaugh et al., 2018; Guo et al., 2020a). By applying $|\Delta g|>1$ on our sample, we can get $\sim$ 56,000 candidates, however over 1/6 of them have $|\Delta r|<0.5$ which are highly suspected and are likely due to defects in photometry measurements, such as by ghost, glint, cosmic ray, CCD problem or unknown instrumental problems. Such defects should be rejected before we could build a reliable sample of EVQs.

In Fig. 2 and Fig. 3, we present for example the $g$- and $r$-light curves for two sources, each showing a clear defect in $g$-band photometry during an individual epoch. In Fig. 2 the SDSS $g$-band cutout of the problematic epoch (marked in red) shows much blurry and diffuse signal comparing with the corresponding $r$-band data and another $g$-band image obtained 116 days later. In Fig. 3 the red dot mark an epoch during which the PS1 $g$-band signal is clearly contaminated by an artificial CCD feature.

Such epochs could be identified through checking the consistency between $g$- and $r$-band variability revealed in the light curves. For SDSS, simultaneous $g$- and $r$-band photometry are always available. However, this is not true for PS1. For each PS1 $g$-band photometry we identify an $r$-band counterpart which was obtained closest in time. In Fig. 4 we plot the distribution of the time intervals between $g$-band and the corresponding (closest in time) $r$-band exposures. For most PS1 $g$-band exposure we could pair it with an $r$-band exposure obtained within 30 days, and use this quasi-simultaneous $r$-band data to examine the reliability of the $g$-band photometry (assuming a quasar does not significantly vary within a month). We simply drop those PS1 $g$-band exposures without corresponding $r$-band exposure obtained within 30 days.

Then for an individual quasar, we could calculate $\Delta g$ and $\Delta r$ (a later measurement minus an earlier measurement) between any pair of epochs, and we plot $\Delta g$ vs $\Delta r$ for all quasars in Fig. 5. We expect that the variability revealed with good photometry should follow the general trend seen in the $\Delta g$ and $\Delta r$ plot (orthogonal distance regression yields $\Delta r=0.878\Delta g$). In Fig. 5 we define a green zone within which the variability well follows the general trend:

and consider those pairs out of the green zone potentially affected by defects in photometry. For a single quasar with N epochs in the light curve, an individual epoch contributes to N$-$1 pairs of epochs. If more than a half of the N$-$1 pairs of epochs locate out of the green zone, we mark this epoch problematic and then drop it from the light curves. A total of 14,899 epochs (0.27% out of 5,559,539) are marked “problematic” with this approach and excluded. In Fig. 5, pairs with the “problematic” epoch involved are plotted in gray, yielding a butterfly-shaped distribution and demonstrating the inconsistency between $\Delta g$ and $\Delta r$ caused by photometry defects.

We consider all available SDSS spectra for the DR14Q quasars and compare their synthetic magnitudes (spectrophotometry) of the spectra with the aforementioned photometric mean magnitudes (see §2.1) to determine their spectral states. Same as the photometric magnitude, we demand the g-band synthetic magnitude error $<$ 0.2 mag. In total, about 1/7 of the spectra were taken during SDSS-III or earlier, about a half taken after PS1.

We fit the quasar spectra mainly following Shen et al. (2011, hereafter S11) and using the PyQSOFit code (Guo et al., 2018). For each spectrum, after correcting for the Galactic extinction adopting the dust map of Schlegel et al. (1998) and the Fitzpatrick (1999) extinction law assuming ${\rm RV}=3.1$, we shift the spectrum to the rest frame using the redshift given in the DR14Q catalog. A global continuum including a power-law and an iron emission template pseudo continuum (Boroson & Green, 1992; Vestergaard & Wilkes, 2001; Salviander et al., 2007) is fitted with tens of separated line-free spectral windows, and the monochromatic luminosity (at $\lambda=$ 5100Å ($L_{5100}$), 3000Å ($L_{3000}$), and 1350Å ($L_{1350}$) for quasars at various redshifts) is then derived. The host galaxy contamination is not considered as most of our sources are high redshift luminous quasars (see S11).

Following S11, we fit various emission lines separately. For H$\beta$ (of objects with $z\leq 0.89$), we use a power-law continuum with iron template to fit within the wavelength windows of [4435, 4700]Å and [5100, 5535]Å. The emission lines are fitted within [4700, 5100]Å with 8 Gaussians, three for the broad components of H$\beta$ with Full Width at Half Maximum (FWHM) $>$ 1200 \textkm s^-1, one for the narrow component with FWHM $<$ 1200 \textkm s^-1, and the rest four for the narrow [O III] $\lambda\lambda$4959,5007 doublet (one core and one wing component for each line). The FWHMs and the velocity offsets of H$\beta$ narrow line and the [O III] core component are tied up. The same restriction is also applied on the wing component of [O III]. For Mg II ($0.35\leq z$ $\leq$ 2.25), we fit the continuum spectra utilizing the same continuum model as above, but within the spectral windows of [2200, 2700]Å and [2900, 3090]Å, and the line over [2700, 2900]Å with three Gaussians for the broad component and one for narrow. For C IV ($z\geq 1.5$), only a power-law continuum is used over windows of [1445, 1465]Å and [1700, 1705]Å, and only 3 Gaussians for broad component over the spectral range of [1465, 1700]Å. We note that, for simplicity and to be easy to reproduce, the numbers of broad Gaussians we used for the lines are fixed rather than variable in S11. To ensure the reliability of the fitting, we only adopt the results of lines with a median spectral S/N (signal-to-noise ratio) around the line-fitting region $>$ 3, which is roughly corresponded to a $\pm 20\%$ fitting bias of FWHM and EW for high-EW objects (Shen et al., 2011). We finally stress that hereafter, unless otherwise stated, the H$\beta$ or Mg II line refers to the broad component we derive from the spectral fitting.

Using the best-fit broad Gaussian models, we measure for each line the EW and FWHM with PyQSOFit, and the line asymmetry with Pearson’s skewness coefficient: ${\rm skewness}=3(\lambda_{\rm mean}-\lambda_{\rm median})/\sigma_{\lambda}$ (Vanden Berk et al., 2001). The $\lambda_{\rm median}$ is where the wavelength bisects the area of the emission line model while the $\lambda_{\rm mean}$ is defined as

We stress that the skewness under such definition can reveal the shape of the emission line model only, regardless of the systematical shift.

We estimate the black hole (BH) mass of quasars based on single-epoch spectrum assuming virialized BLR (S11). With the continuum luminosity as a proxy for the BLR radius and the FWHM of broad line as a proxy for virial velocity, the virial BH mass can be given following the expression as

where the $L_{\lambda}=L_{5100}$ for H$\beta$, $L_{\lambda}=L_{3000}$ for Mg II  and $L_{\lambda}=L_{1350}$ for C IV. We adopt the calibration parameters in S11 (cf. their Equations 5, 8, and 6, respectively):

We also adopt the same bolometric corrections (BCs) as S11 to estimate the bolometric luminosity ($L_{\rm bol}$) where ${\rm BC}_{5100}=9.26$, ${\rm BC}_{3000}=5.15$, and ${\rm BC}_{1350}=3.81$, respectively.

We note that a recent work has compiled spectral properties as well as SMBH mass and bolometric luminosity measurements for the DR14Q quasars (Rakshit et al., 2020, hereafter R20) using the same PyQSOFit code and similar parameters. In R20, the continuum components are fitted to the whole spectrum with various templates, while in this work, we prefer to fit the continuum with nearby line-free spectral windows. Comparing to the global fit in R20, the local fit perform better in fitting continuum flux with simple models within a limited region, so that it could provide more precise emission line spectra. Besides, as we will state below, in this work we are more interested in spectra in the most extreme state for EVQs, 1,363 of which are not the primary thus not included in R20. In addition, in this work, detailed spectral measurements (such as line skewness, bisectional line center and the properties of composed spectra) are required. Therefore we choose our independent spectral fitting results in this work for the following analyses.

Nevertheless, we present in Fig. 6 the comparison of our measurements of luminosity, line EW and FWHM with those of R20 and as well as S11, showing rather negligible deviations and small scatter between the measurements.

We will not go deeper discussing the details of the weak deviation and smaller scatter between the measurements, however, adopting the measurements from R20 does not alter the results presented in this work.

Using the clean and paired $g$- and $r$-band light curve we derived in §2.3, we consider a source as an EVQ if any of its two photometry pairs satisfied $|\Delta g_{max}|>1$ mag and $|\Delta r_{max}|>0.8$ mag quasi-simultaneously. By this criterion, 14,012 EVQs are selected, a catalog of them will be released in a future paper in this series. We plot $|\Delta g_{max}|$ versus $g_{mean}$ of the sample in the upper panel of Fig 7. Though most EVQs are faint sources, most of their $|\Delta g_{max}|$ are statistically significant (with S/N $>$ 5, lower panel of Fig 7), as we have dropped photometric data points with magnitude error $>$ 0.2 mag (see §2).

We note that a simultaneous $|\Delta r_{max}|>0.8$ mag adopted in this work is a strong and conservative request, rejecting which would yield 22,740 candidates instead⁴⁴4This number is still considerably smaller than $\sim$ 56000 aforementioned if applying $|\Delta g_{max}|>1$ mag on the $g$-band light curves before the “cleaning” process described in §2.3. This is because the “cleaning” process has excluded potentially problematic epochs which could yield spuriously large $\Delta g$ (see Fig. 5), and also a portion of $g$-band data points without paired (obtained within 30 days) $r$-band PS1 exposures (see Fig.4).. Practically, EVQs with $|\Delta g_{max}|>1$ mag and $|\Delta r_{max}|>0.8$ mag simultaneously have more extreme variability than those selected with a single $|\Delta g_{max}|>1$ mag criterion. This could be clearly seen in Fig. 8.

Before we determine the state of the spectra of EVQs, we first exclude 3,991 low quality spectra ($g$-band synthetic magnitude error $<0.2\ {\rm mag}$), leaving 2,447 EVQs with no available spectra. To parameterize the states of the spectra, we calculate the magnitude difference between their synthetic photometry $g_{\rm spec}$ derived from SDSS spectra (the AB magnitude evaluated from the spectroFlux given by SDSS which was derived through convolving the spectrum with the corresponding filter) and the mean photometric magnitude ($g_{\rm mean}$, see §2.1). To derive the stacked spectra and explore possible variation of the spectral feature in different states, re-binning the sample is necessary. According to the commonly used EVQs criterion ($|\Delta g_{max}|>1$), we class the spectra with $|g_{\rm spec}-g_{\rm mean}|>0.5$ into extreme states (symmetrically into extremely bright state or extremely dim state). We further divide the rest 8,376 spectra into three classes (see below for details), so that, the whole sample is divided into five classes, which enable us to explore the gradual variation of spectral features from extremely dim to extremely bright states. On the whole, the criterion can be expressed as follow:

1. 1.

   $-\Delta g_{max}<g_{\rm spec}-g_{\rm mean}<-0.5$ : Extremely Bright State (EBS);

2. 2.

   $-0.5\leq g_{\rm spec}-g_{\rm mean}<-0.2$ : Bright State (BS);

3. 3.

   $-0.2\leq g_{\rm spec}-g_{\rm mean}<0.2$ : Median State (MS);

4. 4.

   $0.2\leq g_{\rm spec}-g_{\rm mean}<0.5$ : Dim State (DS);

5. 5.

   $0.5\leq g_{\rm spec}-g_{\rm mean}<\Delta g_{max}$ : Extremely Dim State (EDS).

where the $\Delta g_{max}$ is the greatest magnitude change in each $g$-band photometric light curve. A sketch of the above criteria is shown in Fig. 9. For those 2,341 EVQs with multiple SDSS spectra, we only keep the most extreme spectrum the $g$-band synthetic magnitude of which departs most from the mean photometric magnitude from SDSS and PS1. We defer the study of spectral variability in these individual EVQs with multiple spectra to a future work in this series. We note that dividing the parent sample into more classes would reduce the number of quasars in each class, and adopting different number of classes or using boundaries different from what we adopt would not alter the results in this work. Furthermore, there could be other strategies to re-bin the EVQs into various states. For instance, one may choose to re-bin according to ($g_{\rm spec}-g_{\rm mean}$)/$\Delta g_{max}$ (instead of $g_{\rm spec}-g_{\rm mean}$, see §A), i.e., to normalize the magnitude deviation $g_{\rm spec}-g_{\rm mean}$ by the maximum variability amplitude even seen. However, $\Delta g_{max}$ only represents the maximum variability amplitude of a quasar seen by SDSS and PS1 photometric survey (i.e., would be significantly affected by the sampling), but not necessary an intrinsic physical property of the quasar. Meanwhile, this alternative strategy would not alter the results in this work either.

There are a few EVQs with too bright (13) or too dim (32) synthetic magnitude ($|g_{\rm spec}-g_{\rm mean}|>\Delta g_{max}$). We note that the SDSS spectra taken by recent BOSS campaigns have a smaller fiber diameter (2″ vs 3″) than the former SDSS campaigns. The smaller fiber might increase the possibility of the fiber-drop that results in biased spectra with significantly low flux density which has been reported in literature (e.g., Shen et al., 2015; Sun et al., 2015; Guo et al., 2020a). More information about fiber-drop can be found in Dawson et al. (2012). We suspect the majority of the too dim spectra of 32 EVQs (4 from SDSS-I and II and 28 from SDSS-III and IV) were due to fiber-drop, and exclude them from this study. The too bright spectra of 13 EVQs may be physically real signals or spurious (such as due to natural or artificial solar system objects accidentally passing through the line of sights during spectroscopic observations). Presently, we also exclude them from this paper, and defer studies on these individual sources to a future work.

In total, 11,520 EVQs remain, including 2,152 in EBS, 3,142 in BS, 3,615 in MS, 1,619 in DS, and 992 in EDS. In the following analyses, in order to avoid confusion caused by plotting too many subsamples in a single plot, when necessary we also merge the EBS/BS subsamples into ABS (all bright states) and the EDS/DS into ADS (all dim states). We note that there are considerably fewer sources in DS (EDS) compared with BS (EBS). That is because a clear portion of sources that should belong to DS (EDS) are excluded for their poor synthetic magnitude. Nevertheless, even if we keep those sources with too faint spectra, the dim samples are still considerably smaller. This could be a selection bias⁵⁵5See Shen & Burke (2021) for a different manifestation of this kind of bias caused by AGN variability. that if a quasar was in a dim state, it may not have been spectroscopically identified, or even may not have been selected for spectroscopic observation, and would be missed by the quasar catalogs. However such biases would not affect the following analyses in this work, as the control samples (see §3.2 below) are built to have matched redshift and spectra-derived monochromatic luminosity (thus matched brightness in the spectra) with our EVQ samples.

We then build control samples of subtly-variable quasars with essentially the same physical parameters including redshift, bolometric luminosity and SMBH mass with our EVQs, for the comparison of the spectra properties between EVQs and the control samples. The matching in redshift, bolometric luminosity and SMBH mass is critical as the spectral properties of quasars could significantly evolve (intrinsically or due to intricate observational effects) with such parameters.

The subtly-variable quasars samples are selected out of quasars with $|\Delta g_{max}|<0.4$ and also $|g_{\rm spec}-g_{\rm mean}|<0.4$ (around 140,000 quasars satisfy such criteria). We build the first set of control samples using the $L_{\rm bol}$ and $M_{\rm BH}$ of EVQs measured directly from their single-epoch spectra (§2.4), by selecting the most alike subtly-variable quasar for each EVQ in the space of $L_{\rm bol}$, $M_{\rm BH}$ and redshift (hereafter direct control sample, or DCS). Naturally, each DCS control sample has the same size and matched redshift, $L_{\rm bol}$, $M_{\rm BH}$, broad line FWHM, and also the Eddington ratio ($\eta_{\rm Edd}$) as those of its corresponding EVQ sample.

However, EVQs are experiencing extreme variations, thus the bolometric luminosity derived from the single-epoch spectrum must have been biased by such variations, particularly for those quasars with spectra captured in their extreme states. To better represent the long-term average brightness of a quasar we apply a correction to the aforementioned single epoch bolometric luminosity ($L_{\rm bol}$, and also the monochromatic luminosity $L_{\lambda}$) for each quasar:

We then build the 2nd set of control samples (luminosity-“corrected” control samples, hereafter LCS) with matched redshift, $L_{\rm bol,``cor"}$ and single-epoch SMBH mass. Note during the correction we simply assume the variability amplitude of $L_{\rm bol}$ is the same as that of $g$ band luminosity. Alternatively, we may simply request matching in average $g$ band luminosity (effectively in $g_{\rm mean}$ since redshift is also matched), which however would not alter the results in this work.

Besides, the single-epoch virial black hole mass estimates may also biased by variability (the breathing of broad line region, i.e., the change of line width with luminosity in individual AGNs). It was found while H$\beta$ line display normal breathing expected from the virial relation (Gibson et al., 2008; Denney et al., 2008; Park et al., 2011; Barth et al., 2015; Runco et al., 2016), Mg II shows much weaker breathing (e.g., Shen, 2013; Yang et al., 2019; Guo et al., 2020b; Homan et al., 2020), and C IV exhibits even anti-breathing (e.g., Richards et al., 2002; Wilhite et al., 2005a; Shen et al., 2008; Sun et al., 2018; Wang et al., 2020). Simply assuming no breathing of Mg II and C IV (i.e., the line width does not vary with luminosity)⁶⁶6The broad line breathing in EVQs could be explored using individual EVQs with multiple spectra. We would defer this to a future work in this series., we could further apply a correction to the single-epoch SMBH mass through replacing $L_{\lambda}$ with $L_{\lambda,``cor"}$ during the calculation of mass in Equation 5:

The 3rd set of control samples (for Mg II and C IV samples only) are selected to have matched redshift, $L_{\rm bol,``cor"}$ and $M_{\rm BH,``cor"}$ (Mass-“corrected” control samples, MCS hereafter). Further note that it has been shown $M_{\rm BH}$ estimated using high-ionization emission lines like C IV (e.g., Sulentic et al., 2007; Shen & Liu, 2012; Runnoe et al., 2013; Coatman et al., 2016, 2017) should be used with cautious as C IV may not be virialized. In this case, the control samples could be interpreted more precisely as matched in C IV line FWHM (but not necessarily in SMBH mass).

Note neither the LCS nor MCS samples are perfect control samples of our EVQs, also because the broad line breathing, the bolometric correction factors and the calibration parameters for the virial mass might be different for EVQs and normal quasars, or vary from source to source. However, further meticulous fix is out of our scope and seems unnecessary as we will show later our results in this work are insensitive to the choice of the three control samples we built.

To display the overall spectra of EVQs and the control samples, we construct composite spectra for each sample aforementioned. Before stacking, each spectrum is corrected for Galactic extinction and shifted to the rest frame. For each sample, we normalize the spectra from various quasars at the wavelength of which we derive the monochromatic luminosity (see §2.4, also marked in Fig. 12) and derive a geometric mean spectrum. To avoid confusion in the plots, we merge EVQs in Bright State and Extremely Bright State into EVQs(ABS), and those in Dim State and Extremely Dim State into EVQs(ADS). Their corresponding control samples are equally treated. The composite spectra of EVQs and control samples and plotted in Fig. 12.

Compared with the control samples, EVQs in brighter states exhibit clearly bluer spectra, and reversely redder spectra in dim states. The trend is consistent with the so-called “bluer-when-brighter” pattern widely seen in AGNs and quasars (see §5.1 for discussion).

From Fig. 12 we also see clear line residuals in the spectra ratios of EVQs and the control samples, particularly for EVQ(ADS) and EVQ(Median), showing EVQs tend to have stronger emission lines compared with their control samples. We present detailed comparison of the line EW in §4.2 and line profile in §4.3.

Note in Fig. 12 we only plot the DCS control samples (to avoid confusion). Replacing DCS samples with LCS/MCS samples will not alter the results. A minor note is that the DCS control samples for ABS, MS, and ADS exhibit somehow slightly different spectral slopes between themselves, but only significant in the low redshift bin (H$\beta$ sample). This is likely because these direct control samples for EVQs have lower bolometric luminosity in dim states compared with those in bright states, that the relatively stronger host contamination yields a redder spectrum⁷⁷7The difference is much weaker or disappears if we instead plot LCS/MCS samples which show indistinguishable luminosity distributions.. The effect of host contamination is much weaker at shorter rest-frame wavelength thus the difference almost disappears in higher redshift bins, which is also consistent with previous studies which found that near UV spectral slopes of quasars show little dependence on luminosity (e.g. Telfer et al., 2002; Bonning et al., 2007).

Following the procedures we adopt to fit the individual spectra, we also fit the composite spectra to derive the line parameters. To illustrate the difference of the line EW between EVQs and their control samples, we plot in Fig. 13 the best-fit line EWs (of [O III], C IV, broad Mg II and H$\beta$) derived from the stacked spectra of EVQs in various states (and of their control samples).

Clearly, using the control samples as references (the main results we present below are indeed insensitive to the choice of control samples), we find that EVQs in their extremely dim and dim states tend to have larger line EWs and contrarily smaller EWs in their extremely bright and bright states. This could be primarily be attributed to the so-called intrinsic Baldwin effect (iBeff), i.e., emission line EW in individual AGNs often decreases when AGNs brighten (e.g. Pogge & Peterson, 1992; Kinney et al., 1990; Homan et al., 2020).

Puzzlingly, the broad H$\beta$ line behaves differently. Though the iBeff has been clearly detected in individual AGNs (Goad et al., 2004; Rakić et al., 2017), the stacked spectra of our EVQs exhibit no iBeff of broad H$\beta$.

We further find that, from the overall trend, EVQs tend to have larger emission line EWs compared with the control samples though the extent varies from dim to bright states due to the iBeff. By comparing the EW in MS which is free from the iBeff, we find that the EW of broad Mg II in EVQs is higher by $\sim$ 47% compared with the control sample. For broad H$\beta$, [O III] and C IV, the numbers are $\sim$ 27%, $\sim$ 25% and $\sim$ 25% respectively. Spectra ratio plots in Fig. 12 also reveal clear line residuals around [O III], Mg II and C IV line, further demonstrating that EVQs have systematically larger line EWs compared with control samples. A minor note is that, in Fig. 12 we could barely see line residuals around H$\beta$. The larger best-fit EW of H$\beta$ in EVQs shown in Fig. 13 might be due to the excess of the very broad component of H$\beta$ in EVQs (see §4.3).

We plot in Fig. 14 the stacked emission line spectra of our EVQs in various states, in comparison with the control samples. Following Vanden Berk et al. (2001), the stacked emission line spectra were obtained through normalizing each spectrum by the continuum flux density at the corresponding line center, median stacking the spectra, fitting the median spectra and subtracting the continuum models. From Fig. 14 we could barely see differences in the line profiles between EVQs and the control samples. This is actually expected as the control samples were built to have matched luminosity and SMBH mass with EVQs (thus matched broad line FWHM as the SMBH mass is a virial product of continuum luminosity and line FWHM). Note the situation is slightly different for LCS. LCS samples were built to have matched long-term averaged luminosity and single-epoch spectroscopy-based SMBH mass (compared with EVQs), thus the matching in line FWHM is not guaranteed. In fact, the LCS samples for EVQs in bright (dim) states tend to have slightly larger (smaller) broad line FWHMs. However as we will show below, such an effect is weak and negligible. For MCS, since we use the long-term averaged luminosity to derived SMBH mass, matching in broad line FWHM is also guaranteed. We again stress that the results we provide below are insensitive to the choice of the control samples.

To further explore potential subtle differences in the line profiles between EVQs and the control samples, we plot in Fig. 15 the broad line spectra normalized by accumulated line fluxes. The stacked line profile of broad H$\beta$ and Mg II seem symmetric, and the C IV line exhibit clear redward skewness. Through plotting the spectra in logarithm space in Fig. 15, we find EVQs tend to have stronger broad line wings compared with the control samples. The excesses in such broad-base components are generally seen in all three broad lines we concern (H$\beta$, Mg II, and C IV) and all states of EVQs. They seem symmetric in Mg II (with excesses seen in both the blue and red wings), but redward asymmetric in C IV (only seen in the red wing) and probably also in the dim states of H$\beta$.

We also plot in Fig. 16 the distribution of Pearson skewness measured from individual sources for our EVQs and control samples. The patterns illustrated in Fig. 16 are consistent with what we have revealed from the stacked line profiles (Fig. 15). The median skewness values of broad H$\beta$ and Mg II are much closer to zero than that of C IV  showing H$\beta$ and Mg II are mainly symmetric while C IV exhibits redward asymmetry. Because of the clear excess of the redshifted broad-base component, the C IV line of EVQs is dramatically more redward skewed compared with the control samples, and the K-S test gives $p=1.1\times 10^{-32}$ between EVQs and DCS and $p=2.1\times 10^{-37}$ between EVQs and LCS; The K-S test (see Fig. 16) also reveals statistical difference between the skewness parameter distributions between EVQs and their control samples for H$\beta$ and Mg II  suggesting H$\beta$ and Mg II are also slightly more redward in EVQs, though their median skewness values (of EVQs and control samples) are very close.

To explore the contribution of the excess quantitatively, the broad lines are further fitted with core and wing components. We use two Gaussian with FWHM $<$ 6000 \textkm s^-1 to represent the core component and one Gaussian with FWHM $>$ 6000 \textkm s^-1 for the very broad component (wing). The stacked spectra of EVQs in each state and of their corresponding control samples are fitted together with the line widths and center linked and normalization free to vary. As a result, we find the very broad component in EVQs accounts for higher fraction of the total line flux compared with their control samples. The fraction contributed by the very broad component in EVQs (and in their control samples) to total line flux is $\sim$ 55.0% (39.5%) for H$\beta$, $\sim$ 47.6% (36.1%) for Mg II and $\sim$ 55.2% (47.6%) for C IV, respectively.

Through comparing the composite spectra of EVQs at various states, we find that EVQs follow the general “bluer-when-brighter” pattern widely seen in quasars and AGNs (e.g. Cutri et al., 1985; Wamsteker et al., 1990; Clavel et al., 1991; Giveon et al., 1999; Webb & Malkan, 2000; Trevese & Vagnetti, 2001; Vanden Berk et al., 2004; Wilhite et al., 2005b; Meusinger et al., 2010; Sakata et al., 2011; Schmidt et al., 2012; Bian et al., 2012; Zuo et al., 2012; Sun et al., 2014; Ruan et al., 2014; Cai et al., 2016; Guo & Gu, 2016; Cai et al., 2019; Guo et al., 2020a). The color variation is actually also visible in Fig. 5, where we could see that EVQs (which contribute blue dots in the green zones) follow the same $\Delta r$ vs $\Delta g$ trend of normal quasars. The best orthogonal regression slope ($\Delta r$ = 0.878 $\Delta g$) is less than unity, also demonstrating a “bluer-when-brighter” pattern (similarly, see Fig. 1 in Rumbaugh et al. 2018, but from a much smaller sample).

The similar “bluer-when-brighter” trend seen in EVQs and normal quasars suggests a common physical origin of the variation. Historically, the “bluer-when-brighter” behaviour had been attributed to host galaxy contamination (e.g. Choloniewski, 1981; Hawkins, 2003), or changes in global accretion rate (e.g. Pereyra et al., 2006). Such models have been refuted by recent observations (e.g. Schmidt et al., 2012; Sun et al., 2014; Zhu et al., 2016). It is now more widely accepted that UV/optical variations are due to magnetic turbulence in the accretion disc (e.g. Kelly et al., 2009), and the flux and color variations could be modelled with thermal fluctuations (Dexter & Agol, 2011; Ruan et al., 2014; Cai et al., 2016, 2018, 2020). The same mechanism may also be responsible for the extreme variability seen in EVQs.

The iBeff can be expressed by a simple formula $\rm{EW}\propto L_{cont}^{\beta}$, where $L_{cont}$ is the continuum luminosity (e.g., Kinney et al., 1990; Pogge & Peterson, 1992). In this work, using the control samples of EVQs as reference, we estimate the mean iBeff slope ($\beta$) of EVQs by linear fitting the data with the following relation:

where $EW_{EVQs}(i)/EW_{ctrl}(i)$ is the ratio of the line EW of an EVQ to that of the corresponding normal quasar in the control sample, and ${g_{spec}(i)-g_{mean}(i)}$ is the deviation of its synthetic spectroscopy magnitude from its mean photometric magnitude. The derived iBeff slopes are presented in Table 1 and Fig. 17. In Table 1 we also present the Spearman’s rank correlation coefficients between the two quantities on each side of Equation 14, which yield correlation patterns consistent with the results from the linear regression.

We note that the iBeff slopes of EVQs derived using the DCS as reference are considerably different from those using LCS or MCS as reference. This is particularly prominent for Mg II and C IV that using DCS as reference yields much flatter iBeff slopes. This is because EVQs in dim and extreme dim states have systematically lower luminosity compared with EVQs in bright and extreme bright states (see Fig. 11), thus their DCS control samples would suffer the ensemble Baldwin effect (see Fig. 13). Clearly here using LCS or MCS as reference is less biased to probe the iBeff of EVQs here.

We see significant iBeff of Mg II and C IV in our EVQs, and the iBeff slopes are comparable to those reported in literature for individual AGNs with multiple spectra. Statistically, the slope of Mg II iBeff is slightly steeper than that of C IV, which is contrary to the trend seen in eBeff that the eBeff is stronger for lines with higher ionization energy (Dietrich et al., 2002), which might indicate that the physical origins between the iBeff and eBeff are not connected (Rakić et al., 2017). One possibility is that Mg II line is produced at larger distance and does not respond to continuum variation as fast as C IV thus stronger iBeff is expected.

However, it is rather puzzling that the broad H$\beta$ line of EVQs does not show clear iBeff. This is not only contrary to Mg II and C IV of EVQs, but also to the clear iBeff of H$\beta$ commonly detected in individual AGNs. For instance, a study on a few (six) long-term monitored AGNs revealed the iBeff of H$\beta$ in all subsets of type 1 AGNs (i.e. Seyfert 1, narrow line Seyfert 1 or high-luminosity quasars) with $\beta\gtrsim-0.4$ (Rakić et al., 2017). But note the reported iBeff slopes in literature vary from source to source, or even from year to year for the same source (Rakić et al., 2017; Goad et al., 2004). We notice that the host galaxy contamination could reduce the continuum variability amplitude thus alter the iBeff slope (see §3.6 in S11 for relevant discussion of the role of host contamination on the ensemble Baldwin effect). However, even if after we restrict to EVQs at $\log{L_{bol}}>45.8$ (the most luminous $\sim$10% EVQs in the H$\beta$ sample), we still get rather weak iBeff (with slopes well above -0.12).

H$\beta$ line is previously known to be exceptional in the ensemble Baldwin effect. The eBeff always occurs in high-ionization emission lines whose correlation will be steeper when ionization energy goes higher (Dietrich et al., 2002). However, the Balmer lines, like H$\alpha$ and H$\beta$, exhibit no correlation (e.g., Kovačević et al., 2010; Popović & Kovačević, 2011) or even a weak positive correlation (Croom et al., 2002; Greene & Ho, 2005) between the EW and luminosity of continuum, although Ly$\alpha$ line does exhibit eBeff (e.g. Dietrich et al., 2002). Meanwhile, a recent work of Kang et al. (2021) found that while more variable quasars have clearly stronger (with larger EW) Mg II, [O III] and C IV lines (after correcting the effects of bolometric luminosity, black hole mass and redshift), broad H$\beta$ shows rather weaker correlation between EW and UV/optical variability amplitude. This trend is also visible in Fig. 13 and 14 in this work, that while EVQs have systematically stronger Mg II and C IV line compared with the control samples, the difference is less prominent for broad H$\beta$. It is yet unclear why H$\beta$ behaves exceptionally. Dietrich et al. (2002) proposed that the lack of eBeff of Balmer lines might be due to the intricate and unclear radiation processes of Balmer lines (Netzer, 2020).

We finally note that it would be intriguing to examine the iBeff of individual EVQs with multiple spectra. We would defer this to a future paper of this series.

In Fig. 13 we see EVQs have systematically larger line EWs compared with their control samples. The difference is highly prominent for EVQs in extremely dim state because of the iBeff, remains statistically significant in medium state, and is even visible in the extremely bright state for Mg II. A similar conclusion, free from our sample division, can also be found in Table 1 where we fit the correlation between the $\log{EW_{EVQs}/EW_{ctrl}}$ and the ${g_{spec}-g_{mean}}$ of our EVQs. Comparing with all three control samples, the intercepts are well above zero in the three emission line samples, further support our findings in stacked spectra.

Rumbaugh et al. (2018) has reported that EVQs seem to have systematically larger EW in UV emission lines (compared with normal quasars with matched redshift and luminosity), and attributed such phenomenon to the overall lower Eddington ratio of EVQs. However, in this work, we find EVQs have stronger emission lines, even compared with a group of control sample with matched luminosity and SMBH mass thus matched $\eta_{\rm Edd}$. The results are unaltered if we rectify the measurements of bolometric luminosity and the virial black hole mass which might be biased by the extreme variability (see the comparison with LCS and MCS in Fig. 13).

The discovery that EVQs have systematically stronger emission lines is in good agreement with a recent study of Kang et al. (2021), who found the UV/optical variation amplitude of quasars in SDSS Stripe positively correlate with emission line EWs⁸⁸8Kang et al. (2021) has shown that the correlation between broad H$\beta$ line EW and variability amplitude is however much weaker. This is also consistent with the pattern shown in our Fig. 12 and 13. , after controlling the effects of redshift, luminosity, and SMBH mass. One possible explanation for such correlation is that stronger disc fluctuations could lead to harder quasars SED (Cai et al., 2018), thus provide relatively more ionizing photons. Alternatively, stronger disc turbulence may be able to launch BLR clouds with larger sky coverage.

Notably, Ross et al. (2020) reported several C IV changing look quasars, which show intrinsic Baldwin effect of C IV line, and have larger C IV line EWs comparing with sources with matched SMBH mass (see Fig. 3 of Ross et al. 2020). Therefore, we can see that normal quasars, EVQs, and CLQs exhibit similar intrinsic Baldwin effect, and they follow the same trend that more variable quasars tend to have systematically larger line EWs. These facts suggest common physical processes behind these various populations.

Comparing the broad line profiles of EVQs with their control samples, we find EVQs exhibit subtle excess of the very broad line component (VBC, see Fig. 15). The excess is similarly visible in all states of EVQs.

Statistical studies on line profiles had suggested the broad line region of quasars consists of two components: a very broad line region (VBLR) closer to the central SMBH and an intermediate line region (ILR) at larger distance (Wills et al., 1993; Brotherton et al., 1994; Sulentic et al., 2000). The existence of multiple BLR components is further supported by variation studies which revealed distinct variation patterns of the two line components (e.g. Sulentic et al., 2000; Hu et al., 2020; Guo et al., 2020a). Since the VBLR lies closer to the SMBH, the VBC could easily reverberate (respond to the variation of the central ionizing continuum) in short time. Contrarily, the ILR may appear non-reverberating or reverberating at much longer time. The observed “anti-breathing” of C IV line (i.e., the line broadens when luminosity increases) could also be attributed to the combination the two components (e.g., Denney, 2012; Wang et al., 2020).

The likely physical origin of the VBLR is optically thin gas located near the black hole (Popovic et al., 1995; Corbin, 1997a, b). The Keplerian velocity of that gas could lead to a very broad line width and the SMBH gravity could yield systematical line redshift. The excess of the VBC we discover in EVQs suggests that the strong disc turbulence associated with the extreme variability could launch more gas into VBLR from the accretion disk.

The excess of the VBC is clearly redward skewed in C IV, but not in Mg II and H$\beta$. This is likely because the dominant ILR C IV flux comes from the accretion disk wind which is significantly outflowing thus yielding systematically blue-shifted emission (see also §5.6), while the VBLR is not outflowing. Furthermore, the gravitational redshift of the VBLR region could be more prominent for C IV which could locate at smaller radii compared to Mg II and H$\beta$.

It’s interesting to note that radio loud AGNs and blazars tend to show excess of redshifted very broad components of broad emission lines, particularly in high ionization line C IV (Punsly et al., 2020), and the C IV red wing luminosity excess was found to correlate with radio loudness (the spectral index from 10 GHz to 1350Å, Punsly 2010). Such redshifted very broad component could be produced by gas lying deep within the gravitational potential of the central SMBH, and for the nearly face-on orientation in blazars, the gravitational redshift could be comparatively large (Punsly et al., 2020). The redward excess is somehow similar to what we find in the C IV profile in EVQs. However, limiting our study to radio quiet quasars (quasars with $f_{6cm}/f_{2500}>10$ based on FIRST detections, assuming a radio spectral index of $\alpha=-0.5$, are defined as radio loud, Jiang et al. 2007) does not alter the results in this work. Furthermore, the radio loud fraction in our EVQ sample ($\sim 5.06\pm 0.18\%$) is also comparable to that in the control samples ($\sim 3.95\pm 0.16\%$)⁹⁹9Two fractions are slightly different likely because the photometric selection of radio detected and non-detected quasar candidates for SDSS spectroscopic observations were not uniform (Richards et al., 2002). Thus the discoveries in this work are not due to a small fraction of radio loud quasars in our samples, but represent the properties of the general population of EVQs.

In Fig. 18 we further plot the Mg II and C IV line skewness of radio loud EVQs, compared with the full sample. The median skewness parameters and the K-S test indicate that radio loud EVQs do show (but slightly) more redward skewness (statistically marginal for C IV).

Punsly et al. (2020) also find that blazars with lower Eddington ratio tend to show strong redward asymmetry which can be explained as the Eddington ratio influenced the distance of the most efficient BLR to the center black hole. However, in this work since the control samples were selected to have matched Eddington ratios, thus the stronger very broad component and more redward asymmetric C IV line profile we found in EVQs can’t not be attributed to lower Eddington ratios.

It has been widely reported that the high-ionization BELs in luminous quasars often significantly shifted blueward with respect to low-ionization lines (e.g., Gaskell, 1982; Wilkes, 1986; Corbin, 1990; Sulentic et al., 2000, 2007; Richards et al., 2002, 2011; Baskin & Laor, 2005; Shen et al., 2008; Wang et al., 2011; Denney, 2012; Shen & Liu, 2012; Coatman et al., 2016, 2017; Sun et al., 2018). Such blueshift is commonly contributed to accretion-disk winds (e.g., Richards et al., 2011; Denney, 2012).

We test whether the C IV blueshift in EVQs differs from that in the control samples, using quasars with redshift 1.5 $<$ z $<$ 2.25 for which the low ionization Mg II line is available to derive the systematical redshift. We measure the velocity offset of C IV line in each quasar, for both the line peak and the bisectional line center ($\lambda_{b}$ that bisect the total line flux, Sun et al. 2018), with respect to the systematical redshift determined from the Mg II line peak. We find that while both EVQs and their control samples show clear C IV blueshift on average (with negative median values of the offset), the median peak blueshift is marginally ($<$ 2$\sigma$) weaker in EVQs (see Fig. 19). Note the K-S test shows that the distribution of C IV line peak velocity offset in EVQs is different from that of the control sample DCS with a p-value of $6\times 10^{-10}$. This indicates that, compared with the median line peak blueshift, the line peak velocity offset distribution better reveals the difference between EVQs and the control sample.

We further find the C IV bisectional line center of EVQs shows a redward-skewed distribution (with a positive median velocity offset), clearly different from their control samples (with negative offsets), and the K-S test yields a p-value of $5\times 10^{-65}$ between EVQs and DCS. That is, because of the more redward-skewed line profile of C IV in EVQs, though the line peak is blueshifted, the bisectional line center is contrarily redshifted. This pattern is consistent with what we have revealed from the stacked line profile (see §4.3). The difference in the line profile may also partially account for the weaker C IV peak blueshift in EVQs. We note that the series weak emission lines (He ii $\lambda$1640, O iii] $\lambda$1663 and He ii $\lambda$1671) reside in the red wing of C IV could also affect the bisectional line center. However, considering that those lines are relatively weak and submerged in the flux of C IV, it is very hard to measured those lines precisely with current low SNR spectra. Nevertheless, given their line centers which are $\sim 100\AA$ away from that of C IV, they are not likely to affect the peak of C IV nor to contribute the red excess of C IV.

It has been theoretically proposed that the disc wind in quasars could be driven by either radiation pressure (continuum and/or UV line) or magneto-centrifugal forces, or some combination thereof (Murray et al., 1995; de Kool & Begelman, 1995; Proga et al., 2000). As EVQs have stronger emission line compared with their control samples with matched redshift, UV monochromatic luminosity and SMBH mass, they may have relatively stronger ionization continuum (see §5.3), thus stronger continuum radiation pressure is expected. The fact that the C IV line peak blueshift in EVQs is weaker than that in the control samples disfavors the scenario that the wind acceleration is dominated by continuum radiation pressure.