The Time Domain Spectroscopic Survey: Changing-look Quasar Candidates from Multi-epoch Spectroscopy in SDSS-IV

Supermassive black holes (SMBHs) in the cores of galaxies have garnered a special fascination since the first quasars were discovered as strong, variable radio sources (e.g., Matthews & Sandage 1963; Sandage 1964) with bright optical quasi-stellar (point-source) counterparts (Oke & Schmidt 1963), with optical spectra showing a cosmological redshift and emission lines broadened by 10,000 km s⁻¹ or more (Schmidt 1963). Since then, massive efforts to study SMBHs have revealed that they span from 10⁴–10¹⁰ M_☉, found in the core of the Milky Way, in nearby dwarf galaxies (e.g., Baldassare et al. 2015; Afanasiev et al. 2018), in the dominant galaxies of galaxy clusters (e.g., McConnell et al. 2011), and out to redshifts above 7, when the universe was just 5% of its current age (e.g., Fan et al. 2006; Bañados et al. 2018; Wang et al. 2021). SMBHs may be luminous and easily detected in active galactic nuclei (AGN) if they are accreting significantly, or they may be dim and quiescent, as in our galaxy. SMBHs grow primarily by gas accretion or by BH-BH coalescence (e.g., Kauffmann & Haehnelt 2000), which should become detectable via gravitational waves in the near future with the pulsar timing arrays (e.g., Ransom et al. 2019; Arzoumanian et al. 2020 and NASA's Laser Interferometer Space Antenna mission; Danzmann 2018).

SMBH masses (M_BH) in AGN observationally correlate tightly with a variety of their host galaxy properties, such as the bulge mass and luminosity (e.g., Kormendy & Richstone 1995; McLure & Dunlop 2002; McConnell & Ma 2013), and the stellar velocity dispersion σ, (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Gültekin et al. 2009; McConnell & Ma 2013). Copious research has assumed or analyzed the possibility of feedback, where the observed M_BH–σ_* relation (i.e., the relationship between central BH mass and the stellar velocity dispersion σ_* in the bulge of a galaxy; Ferrarese & Merritt 2000; Gebhardt et al. 2000) is mediated by powerful outflows from the accreting SMBHs that regulate galaxy growth (e.g., Silk & Rees 1998; Fabian 2012). However, the reality may be more complex, e.g., the M_BH–σ_* relationship may simply be the result of stochastic averaging over galaxy merger histories (Jahnke & Macciò 2011), and AGN outflows may enhance rather than quench star formation (Zubovas & Bourne 2017). There are complex, high-scatter relationships between AGN activity, outflows, feedback, and galaxy properties (e.g., Fiore et al. 2017). Powerful and increasingly complex cosmological simulations such as IllustrisTNG (Pillepich et al. 2018) and FIRE (Hopkins et al. 2018) seek to tie together SMBH growth with cosmological structure formation.

All this makes clear that the accretion rate of SMBHs, i.e., quasar activity, is important to understand and characterize on all accessible timescales. Indeed, variability of 10%–20% observed on timescales of days to years is a fundamental characteristic of quasars that has often been used to identify samples with high efficiency (e.g., Palanque-Delabrouille et al. 2013). Larger amplitudes (a magnitude or more) of variability are also known, usually found on longer timescales and within large samples, but also sometimes serendipitously (e.g., Kynoch et al. 2019; Luo et al. 2020, and references therein). On cosmic timescales, variability may be crudely characterized by the total average time, or duty cycle spent in active accretion (Martini & Schneider 2003). The present-day SMBH mass function basically depends on quasar lifetime and duty cycle. We know that SMBH accretion activity varies significantly on timescales of millennia or longer. Kiloparsec-scale jets and cavities driven by SMBH activity are seen in massive radio galaxies, and piercing the hot plasma of galaxy clusters (e.g., Bîrzan et al. 2004). Double-double radio galaxies show jets that have stopped and restarted on such timescales (Mahatma et al. 2019). Extended emission line regions lacking a present-day ionizing source (Voorwerps; Lintott et al. 2009; Keel et al. 2012) are also being discovered that testify directly to long periods of AGN quiescence (up to ∼10⁵ yr so far; Schawinski et al. 2015).

We know that quasars are variable, and that emission lines respond to changes in the ionizing continuum, albeit with a less than linear relationship, known as the Baldwin effect (Beff), which is seen both in samples of quasars with single-epoch spectra (the ensemble Beff or eBeff, e.g., Baldwin et al. 1978; Green et al. 2001) and between spectral epochs of individual quasars (the intrinsic Beff or iBeff; Kinney et al. 1990; Goad et al. 2004). The temporal lag in the broad emission line (BEL) response to continuum changes continues to be used to measure SMBH masses via reverberation mapping (RM; Peterson 1993; Shen et al. 2016; Grier et al. 2017). On occasion, both the accretion luminosity (as represented by a power-law continuum component in the optical spectrum) and BELs are seen to vary quite strongly, and even more rarely, both may dim below detectability. Extreme variability observed in—particularly dis/appearance of—the BELs is called the changing-look phenomenon in AGN (CLAGN) or quasars (CLQs). ¹⁹ As the narrow emission lines (NELs), which originate from much larger regions (e.g., Bennert et al. 2002), tend to remain steady relative to the BELs, AGN may change their classification, based on which emission lines continue to show broad components (e.g., Osterbrock 1981) from Type 1 to Type 2, or by smaller increments, e.g., Types 1.8–1.5.

CLQs have been intensively studied as of late because of both their mystery and their utility. Depending on the precise definition of a CLQ (e.g., continuum and broad Hβ luminosity and their change), nearly 70 have been found just in the last few years, mostly identified from multi-epoch optical photometry and follow-up spectroscopy of previously known (spectroscopically identified) quasars (LaMassa et al. 2015; Ruan et al. 2016a; MacLeod et al. 2016; Runnoe et al. 2016; Gezari et al. 2017; Yang et al. 2018; Frederick et al. 2019; MacLeod et al. 2019; Ross et al. 2020; Sheng et al. 2020). The strong interest arises because for standard thin accretion disk models (Shakura & Sunyaev 1973), large luminosity changes in the optical-emitting regions of an SMBH accretion disk due to overall accretion rate changes are expected to occur on the viscous timescales, corresponding to thousands of years or more (Krolik 1999), whereas CLQs have been seen to change significantly on timescales of years or even as little as about a month (though the latter usually at lower Seyfert-like luminosities, e.g., Trakhtenbrot et al. 2019a). Therein lies the mystery. If instead, the optical continuum-emitting regions are reprocessing radiation from the inner UV and X-ray-emitting regions, observed timescales for CLQ variability may be compatible with predictions, as discussed in LaMassa et al. (2015, and references therein). Accretion disks supported by magnetic pressure may also allow powerful and rapid variability (Dexter et al. 2019; Scepi et al. 2021). However, there may be theoretical caveats to the timescale problem, as described elsewhere (e.g., in Lawrence 2012, 2016).

The utility of CLQs arises because we may witness quasars essentially turning on and/or off. This clearly demands revision of our understanding of the duty cycles and timescales of quasar activity, but also allows us to see clearly (without the contamination of a bright point-source nucleus) what a quasar host galaxy looks like. For instance, though the M_BH–σ_* relationship demands it, measuring the stellar velocity dispersion in a quasar host illuminated by a luminous active SMBH can be difficult, but becomes simple when the accretion activity ceases (Charlton et al. 2019; Dodd et al. 2021). A CLQ in a dim state also enables direct measurement of the host luminosity for accurate subtraction from the overall spectrum in brighter states (e.g., Bentz et al. 2006).

Tidal disruption events (TDEs) have been suggested as an explanation for intrinsic extreme variability. Debris from a star whose orbit brings it too close to the SMBH is pulled into a stream, of which half forms (or augments) an accretion disk. The resulting thermal radiation with ${T}_{\mathrm{eff}}\sim {({L}_{\mathrm{Edd}}/(4\pi \,{R}_{t}^{2}\sigma ))}^{1/4}$ (where R_t is the tidal radius of the BH) is expected to peak in the EUV or soft X-rays (0.25–2.5 × 10⁵ K) but also irradiate circumnuclear gas resulting in enhanced optical emission. The fallback of debris is expected to rise for about a month, and then decline as t^−5/3 (Rees 1988; Ulmer 1999). The CLQs found to date have light curves that show a wide variety of shapes. On rare occasions, what appears to be a CLQ may be a TDE (e.g., see the Merloni et al. 2015 analysis of the light curve of LaMassa et al. 2015). However, most CLQs last too long (∼years or decades) in the bright state for the emission to be primarily from a TDE (e.g., Runnoe et al. 2016; MacLeod et al. 2019) while others have a very different decay trend (e.g., Ruan et al. 2016a; Gezari et al. 2017). Still, rare, potentially longer-lasting light-curve patterns may emerge from unusual TDE systems involving, e.g., giant stars, or binaries.

Extreme variability may also occur due to strong changes in intrinsic absorption that could block the continuum source and shade the BEL region, but cloud crossing timescales (t_cross = $24\,{M}_{8}^{-1/2}\,{L}_{44}^{3/4}$ yr) are usually considered much too long. While modeling of dust extinction of the optical/UV quasar power-law continuum can sometimes reproduce the observed changes in the optical continuum emission, the observed BEL changes indicate a near disappearance of the ionizing flux (LaMassa et al. 2015; MacLeod et al. 2016). Obscuration of the continuum should also result in high linear polarization, which has not been found in dim state CLQs (Hutsemékers et al. 2019). Sheng et al. (2017) find that mid-infrared (MIR) variability of 10 CLQs appears to follow variability in the optical band with a time lag consistent with dust reprocessing. Their results argue that obscuration is not the cause of dimming in CLQs. Yang et al. (2018) found the MIR Wide-field Infrared Survey Explorer colors of CLQs are redder when brighter, indicating a strong hot dust contribution rather than reddening. In several cases, the IR tracks the optical variability (e.g., Nagoshi et al. 2021), again suggesting that obscuration is not the cause of dimming in CLQs.

TDEs and variable obscuration may contribute to samples of CLQ candidates, but larger samples are key to understanding the role and frequency of these phenomena relative to the predominant effects of accretion rate changes. The characterization of strong changes in accretion rate leading to changes in the broad line region (BLR) also demands larger samples, to study the dependence of the CLQ phenomenon on fundamental parameters such as luminosity, SMBH mass, Eddington ratio, and redshift.

For these reasons, the hunt for CLQs has intensified in the last few years, especially as large-area, multi-epoch surveys—both photometric and spectroscopic—are maturing. In this paper, we describe a search for CLQs among objects with multi-epoch spectroscopy from a dedicated program targeting point-source variables for spectroscopy, the Time Domain Spectroscopy Survey (TDSS) of the Sloan Digital Sky Survey (SDSS).

Throughout, we adopt a flat Λ cold dark matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.30.

Our parent sample consists of a subset of objects classified as quasars by the SDSS pipeline that appear in the SDSS Data Release 14 (DR14) quasar catalog (Pâris et al. 2018), which included 526,356 quasars. All the reduced and calibrated SDSS spectra analyzed herein were publicly released as part of SDSS DR14, which includes data through 2016 May 11. The SDSS-IV final quasar catalog DR16 (Lyke et al. 2020) describes an additional 225,082 new quasars, and additional epochs obtained for previously known SDSS quasars.

The Time Domain Spectroscopy Survey (TDSS) is a subprogram of the extended Baryon Oscillation Spectroscopic Survey (eBOSS; Dawson et al. 2016 within SDSS-IV (Blanton et al. 2017). A pilot survey for the TDSS dubbed SEQUELS actually started during SDSS-III (Eisenstein et al. 2011), and is described in Ruan et al. (2016b). The TDSS is the first major program intended to spectroscopically characterize generic variables, i.e., without any explicit selection criteria based on color or light-curve characteristics. The TDSS targeted photometric objects having morphology consistent with a point source, and that were imaged both by SDSS and Pan-STARRS1 (PS1; Kaiser et al. 2010), and showed significant variability within and/or between those two imaging surveys, as described by Morganson et al. (2015). The bulk of the TDSS sought single-epoch spectroscopy (SES) of such point-source, photometric variables. Results up through SDSS DR14 are described in S. F. Anderson et al. 2022, in preparation. However, the TDSS also included a smaller program for Few-Epoch Spectroscopy (FES), to examine the spectroscopic variability of certain select samples of objects with existing SDSS spectroscopy, totaling about 20,000 spectra, as described in MacLeod et al. (2018). These included hypervariable quasars, broad absorption line (BAL) quasars, quasars with possible double-peaked broad line profiles, white dwarf/main-sequence binaries, and dwarf carbon stars (Roulston et al. 2019). The TDSS is also described on the web. ²⁰

Another TDSS subprogram targeting objects with existing spectroscopy is the Repeat Quasar Spectroscopy (RQS) program, described in detail in MacLeod et al. (2018). The RQS program explicitly targeted known quasars for repeat spectroscopy so that quasar spectroscopic variability could be studied across a wide range of timescales. The majority (about 70%) of RQS targets are simply a bright magnitude-limited sample (RQS1; i < 19.1), but quasars with i < 20.5 and more than one spectrum in the archive were also prioritized. The RQS also includes fainter (i < 21) subsets with and without detected photometric variability above a tuned threshold ²¹ : RQS2, RQS2v, RQS3, and RQS3v. Unlike the TDSS SES or FES subprograms, the RQS sample also includes quasars with extended morphologies in SDSS imaging.

All TDSS DR14 spectra of quasars were reviewed in visual inspections for both quality assurance and potential science discoveries. We asked student coauthors at the University of Washington to study example spectra of normal quasars, and to expect such spectra for the great majority of new TDSS quasars. However, we asked the visual inspectors to be especially alert for new epochs of TDSS spectra that did not look like those of normal quasars, including, e.g., CLQs, quasars with odd emission line profiles, BALs, etc. Wherever an incoming new-epoch TDSS quasar spectrum was deemed potentially unusual by the inspector, that new-epoch spectrum was compared side by side with earlier epoch spectra contained in SDSS DR14. We note that this process strongly favors finding turn-off CLQs (in their fading or dim state). Furthermore, the method is both qualitative and incomplete, so that any statistical utility of the number of CLQs found here is at best limited to being a lower limit. However, human visual inspection can detect small changes that might be missed by a computer algorithm, and also discount a variety of apparent changes that are clearly caused by artifacts. There are also cases in between, where an apparently significant change may potentially be due to a problem, e.g., with spectral fiber plug placement in the plate, or glitches in the reduction. In the majority of such cases, we have verified whether the observed spectral change is likely real by checking against the behavior of the quasar's photometric light curve at the spectral epoch(s) in question (e.g., see Figure 1). However, several cases (see Table 2) called for follow-up spectroscopy with different telescopes and instruments to confirm the reality of the apparent changes first seen between the SDSS spectral epochs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (5 panels)

Tarball of figure set images

We restrict the TDSS CLQ candidate sample to redshift z < 0.9, for which the Hβ/O iii emission line complex remains below about 9500 Å in the observed frame, close to the red limit of early epoch SDSS spectra.

For clarity and statistical purposes, we define a bona fide CLQ to have a 3σ change in the flux of the broad component of Hβ, as described in detail in Section 3.5. As can be seen in Table 1, the largest share (23) of our CLQ candidates, of which four are bona fide CLQs, comes from the TDSS FES_HYPQSO sample. ²² This makes perfect sense, since they were prioritized for TDSS spectroscopy from the 2% most variable QSOs, outside the elliptical contours of a variability space defined both by variability within PS1 epochs, and between SDSS and PS1 epochs (Morganson et al. 2015; MacLeod et al. 2018). However, while large continuum changes enhanced the likelihood of their visual identification as CLQ candidates, we note (with a strong caveat for small number statistics) that the HYPQSOs do not represent the largest fraction of confirmed CLQs.

The next largest contributing TDSS selection is the TDSS RQS1 program, yielding 18 CLQ candidates and five confirmed CLQs. The RQS1 designation arises from a simple magnitude cut, with no variability criterion, so the similar final number of CLQs to the HYPQSO-selected sample is probably just a result of the larger RQS1 parent sample size. Indeed, given the small number of CLQs in Table 1, we can only say with confidence that a random spectroscopic resampling of quasars on these timescales (or longer) will yield a useful number of new CLQs.

As highlighted by Ruan et al. (2016a) and S. F. Anderson et al. 2022, in preparation, variability selection results in a broader quasar color distribution than traditional color selection, such as used to select most SDSS quasars, where the primary goal is to distinguish quasars from stars in photometric color space (e.g., Bovy et al. 2012), and the result is a bias against red quasars. The larger fraction of red quasars uncovered in variability-selected samples can be attributed to more high redshift objects, and/or redder continuum emission due, e.g., to their intrinsic continuum slope and/or internal dust reddening.

This work primarily focuses on (1) repeat spectroscopy, for which the initial target selection was most often via colors (Richards et al. 2002) and (2) redshifts where Hβ is visible, so our parent CLQ candidate sample has fairly typical initial (bright state) color distribution. Dimmer states have often lost much of their bright UV/blue continuum and are rather more dominated by the host galaxy, so the dim state population is indeed relatively red. In some quasars with redshift z < 0.35, strong broad Hα may persist even in the dim state, adding further flux to the red filter bands.

3.1. Spectroscopic Data

3.2. Visual Inspection and Rescaling

3.3. Light Curves and Spectra

3.4. Spectral Decomposition

We use PyQSOFit ²⁹ for spectral decomposition (Guo et al. 2018) of all of our SDSS spectra as well as for the follow-up spectra we obtained on other telescopes. We correct the spectra to the rest frame and correct for Galactic extinction using the extinction curve of Cardelli et al. (1989) and the dust map of Schlegel et al. (1998). We then perform a host galaxy decomposition using galaxy eigenspectra from (Yip et al. 2004a) as well as quasar eigenspectra from (Yip et al. 2004b). If more than half of the pixels from the resulting host galaxy fit are negative, then the host galaxy and quasar eigenspectral fits are not applied. We then fit the power-law, UV/optical Fe ii, and Balmer continuum models. The optical Fe ii emission template spans 3686–7484 Å, from Boroson & Green (1992), while the UV Fe ii template spans 1000–3500 Å, adopted from Vestergaard & Wilkes (2001), Tsuzuki et al. (2006), and Salviander et al. (2007). PyQSOFit fits these empirical Fe ii templates using a normalization, broadening, and wavelength shift. Next, we perform emission line fits, using Gaussian profiles as described in Shen et al. (2019). Depending on redshift and spectral coverage, we fit the following emission lines: Hαλ6564.6 broad and narrow, [N ii]λ6549,6585, [S ii]λ6718,6732, Hβλ4863 broad and narrow, O iii λ5007,4959, Mg ii λ2800 broad and narrow, C iii]λ1908, C iv λ1549 broad and narrow, Lyαλ1215 broad and narrow. We run all of these fits using Monte Carlo simulation based on the actual observed spectral error array, which in turn yields an error array for all our decomposition fits. An example of spectral decomposition is shown in Figure 2.

Zoom In Zoom Out Reset image size

The host galaxy fits used in PyQSOFit are limited to rest-frame wavelengths between 3450 and 8000 Å. Due to this limitation, to fit the Mg ii line complex, we also run PyQSOFit separately on all our spectra and epochs with host decomposition off.

We do not fit a polynomial continuum, as we found it often competed strongly with the power-law continuum, yielding unreasonable fits for both continuum components. At first, we experimented with fitting a host galaxy model only to the dim state spectrum, since it should be of the highest contrast there, and easiest to fit without contamination from the quasar continuum. However, we found that applying the best-fit dim state host model to all epochs often resulted in poor overall fits, likely due to differing combinations of seeing and spectral aperture for different epochs. For this reason, we fit a host galaxy component to every spectral epoch separately.

3.5. Defining and Identifying Bona Fide CLQs

For purposes of comparison between quasars, spectral epochs, or different studies of quasar samples, it is crucial to have a common definition of what we mean by the CLQ phenomenon. For instance, it is insufficient to merely say that "broad Hβ disappeared," since either a visual impression or a measurement can be strongly affected by the signal-to-noise ratio (S/N).

After running PyQSOFit on all of our spectra, we use the fits for the full continuum model (power law, the UV and optical Fe component, and Balmer continuum models) along with the galaxy host model if applicable. We first subtract the full continuum model and the masked ³⁰ galaxy host fits from each observed spectrum to arrive at a quasar line flux spectrum. For each such spectrum, we rebin both the flux spectrum and corresponding variance spectrum to 2 Å pixel⁻¹ so that different epoch spectra can be directly compared.

In the column marked "State" of Table 2, we denote our choice of dim and bright spectral epoch for each quasar by D and B, respectively. The dim state is chosen as that with the lowest 3240 Å continuum luminosity L₃₂₄₀, and marked with D in the State column. For all other spectra, we run our N_σ (Hβ) calculation on this quasar line flux array, at every pixel across the Hβ region from 4750–940 Å, as follows:

$$\begin{eqnarray}&&{N}_{\sigma }=({f}_{\mathrm{bright}}-{f}_{\dim })/\sqrt{{\sigma }_{\mathrm{bright}}^{2}+{\sigma }_{\dim }^{2}}.\end{eqnarray} \tag{ 1 }$$

Here, f is the flux in erg cm⁻² s⁻¹ Å⁻¹ in the pixel, whereas σ is the spectral variance, including as usual the propagated uncertainties from statistical and instrumental noise and sky background subtraction. The N_σ (Hβ) array is then smoothed with a median filter using a kernel size normally of 16 pixels (32 divided by the sampling rate of the spectrum in Å pixel⁻¹). We then subtract the N_σ (4750) value from this heavily smoothed array, and thereby find the maximum relative value of N_σ (Hβ), which is tabulated in Table 2. This is the same method used in MacLeod et al. (2019), as inspired by a similar usage in Filiz et al. (2012). An example of normalized dim and bright state spectra overplotted can be seen in Figure 3.

Zoom In Zoom Out Reset image size

We generally choose for a CLQ determination the bright state spectral epoch with the largest value of N_σ (Hβ). This is also the bright state we use wherever a comparison is made between just two epochs for each quasar, such as when plotting Hβ flux changes versus luminosity in Figure 8. Our selected bright state epochs are marked B in the State column of Table 2. In very few cases for quasars with more than two spectral epochs, rather than using N_σ (Hβ) to select the bright spectral epoch, we instead choose the spectral epoch with the highest L₃₂₄₀, if visual inspection of the N_σ (Hβ) spectral overplots or the light-curve trends present compelling evidence for that choice. For example, for J021359.79+004226.81, we chose MJD 51816 because while its N_σ (Hβ) value of 6.54 is not as high as that of 9.43 for MJD 57043, the enhanced blue spectral continuum is much stronger at MJD 51816.

Note that we do not use identical criteria to MacLeod et al. (2019) because they require a priori that ∣Δg∣ > 1 mag and ∣Δr∣ > 0.5 mag), and then verified CLQ status from pursuant follow-up spectroscopy. In comparison, our candidate sample begins with multi-epoch spectroscopy, and is, when necessary, verified by multi-epoch photometry. In Figure 4, we show an example where initially convincing spectroscopic changes are negated by comparison to photometry.

Our use of a significance threshold as a criterion for determining CLQ status, while well defined and efficient, is not ideal in that it is dependent on the spectroscopic data quality rather than on intrinsic changes to the quasar spectra over time. We suggest a potential intrinsic definition later in Section 6.3.

In Table 2, we summarize data and measurements of our visually identified CLQ candidate sample. Each quasar has at least two spectral epochs. Under the SDSS ID of each quasar, we note the TDSS subprogram under which criteria it was targeted for repeat spectroscopy (see also Table 1). For each spectral epoch, we list the facility used to obtain the spectrum, the MJD, and luminosities and their uncertainties derived from our model fits for 3240 Å continuum, broad Hβ, 2800 Å continuum, and Mg ii emission, as available, in units of 10⁴² erg s⁻¹. The N_σ (Hβ) value we derive, as described above in Section 3.5 is listed for every epoch but that of the dim state.

In the final (Notes) column of Table 2, we confirm CLQ status for the bright state spectrum with the largest N_σ (Hβ) value above three, such that the CLQ designation appears at most once in the table for any given quasar. We find 19 CLQs in this TDSS sample, of which four were previously noted as CLQs in MacLeod et al. (2019) (J1055+2425, J1113+5313, J1434+5723) or MacLeod et al. (2016) (J0023+0035). In most cases, we have additional spectral epochs available.

The light curves of the 19 CLQs span as much as 20 yr, and show strong diversity in the character of variability. We sometimes find significant changes in brightness on surprisingly short timescales; J002311.06+003517.53 dims by Δg ∼ 1.5 mag within ∼1000 days. J105513.88+242553.69 shows instead a slow, steady dimming over 15 yr by Δg ∼ 2 mag. J143455.30+572345.10 dims by as much, on a similar timescale, but with sparser photometric coverage.

By contrast, some quasars show relatively minor shifts in photometry, yet still show changes in L₃₂₄₀ and N_σ (Hβ) that are significant, such as J024508.67+003710.68, which may have dimmed only about Δg ∼ 0.3 mag. However, its light curve, and those of many CLQ candidates, do not always sample epochs close to the spectroscopy. Sometimes, large, rapid changes are seen in photometry, but without nearby spectral epochs; J111329.68+531338.78 rebrightens in recent epochs by Δg ∼ 1.5 mag within just a year. The variability of individual CLQs is discussed in detail in the Appendix.

Most of the CLQs that we identify are turn-off CLQs, i.e., they satisfy our criteria by dimming with time. This is a consequence of our selection process, since our parent sample is a quasar catalog requiring identification of BELs. We nevertheless find that four of the 19 CLQs are turn-on CLQs, identified as such in the Notes column of Table 2: J100302.62+193251.28, J113651.66+445016.48 (also noted in MacLeod et al. 2019), J224113.54-012108.84, and J234623.42+010918.11. The latter is particularly striking, showing a brightening of Δg ∼ 1.5 mag early in the light curve. In such instances, the early epoch pipeline classification as quasar may be due to BELs other than Hβ, but in most cases, the continuum and broad Hβ emission were simply weaker, and have increased significantly during our monitoring.

Zoom In Zoom Out Reset image size

There are three quasars that distinctly appear by visual inspection to be CLQs based on both spectra and light curves, but we measure only 2 ≤ N_σ (Hβ) < 3. We mark these in the Notes column of Table 2 as VIc. Two of them were previously published as bona fide CLQs—J1002+4509 (MacLeod et al. 2016) and J1150+3632 (Yang et al. 2018). These would likely have passed our CLQ criteria, had the spectral S/N been higher. We show the light curves and spectra for these objects in Figure 5.

Zoom In Zoom Out Reset image size

While our candidate CLQs extend to redshift 0.9, the CLQs we identify are all at z < 0.5. Figure 6 shows histograms of the full DR16 quasar sample between 0 < z < 0.9 (black spikes), our TDSS FES+RQS targets (blue), CLQ candidates (green), and CLQs (red). For clarity, each histogram is normalized so that it peaks at one. We also plot the absolute and apparent i-band magnitudes against redshift for the same samples.

Zoom In Zoom Out Reset image size

Despite the small number of statistics, these plots indicate that CLQs are likely from a typical quasar population, and even when compared to our TDSS FES+RQS target sample with its brighter magnitude limits and/or variability criteria, tend to be found at lower redshifts and brighter magnitudes, likely due to our S/N criterion. The i-band magnitudes shown in Figure 6 are taken from SDSS, so are typically from an earlier, brighter phase. The i band, with effective wavelength 7628 Å (Fukugita et al. 1996) samples rest-frame wavelengths near 5900 Å (5330, 4737, 4487 Å) for the median redshifts of the CLQs (0.3; candidates 0.43; FES+RQS 0.61, full DR16 sample 0.7). This i band is a reasonable choice because these rest-frame wavelengths are normally redward of ∼4000 Å, where the quasar continuum changes strongly between dim and bright states for quasars.

4.1. Luminosity Changes, Timescales, and Flux Ratios

Figure 7 shows the rest-frame time interval between spectral epochs ΔMJD versus luminosity change for our full CLQ candidate sample, both for the 5100 Å power-law continuum and broad Hβ. The only clear trend is that a larger fraction of candidates are bona fide CLQs at longer epoch separations, for turn-off CLQs (ΔMJD > 0). Statistics are inadequate in our sample to judge such trends in turn-on CLQs. Interestingly, it is not readily apparent that larger multi-year time spans lead to larger luminosity changes in our sample. Luo et al. (2020) found that even among EVQs (i.e., photometrically selected to have Δg > 1 mag), only a small fraction (<4%) showed monotonic variability over a ∼16 yr time span ³¹

Zoom In Zoom Out Reset image size

For quasars in general, stochastic variability, with a structure function described by an asymptotic long-timescale rms variability SF_∞ ∼ 0.2 mag and a rest-frame damping timescale τ ∼ 200 days provides a reasonable description of quasar variability (MacLeod et al. 2010). EVQs are better represented by SF_∞ ∼ 0.4 mag and timescale τ ∼ 600 days, but the relative fractions of longer-term variability trends (monotonic versus single peak/dip versus complex) are not well reproduced by simple stochastic variability models (Luo et al. 2020).

Our results indicate that large multi-epoch spectroscopic samples (or large time domain photometric samples with prompt spectroscopic follow-up of strong variables) are likely to find a small fraction of quasars to be CLQs, but that CLQs may be found spanning a wide range of timescales sampled. Whether epoch separations less than a few years will reveal substantial CLQ behavior is not clear, since our sample is sparsely populated in that range (see Figure 7). However, the larger RM program planned for SDSS-V could detect and track perhaps 2000 quasars over more than about 150 epochs to detect a larger number of short-timescale CLQs, such as SDSS J141324+530527, reported by Dexter et al. (2019). The medium tier of the All Quasar Multi-Epoch Spectroscopy (AQMES) survey of SDSS-V will obtain hundreds more spectra with rest-frame intervals of months (see Section 6.3).

4.2. Broad Line Variability

In Figure 8, we plot broad Hβ line strength against 3240 Å continuum strength, showing a strong correlation in both luminosity and flux space. The correlation is not surprising, since Hβ is the optical emission line most used for RM studies. The success of RM confirms the strong influence of photoionization in the BLR. The ionizing source is assumed to be small (relative to the BLR) and quasi-isotropic, Hβ is indeed seen to react in many quasars to L₅₁₀₀, with a delay from about 1–100 days. Thus in the great majority of cases here, the broad Hβ line flux has had time to react to continuum changes between the bright and dim states measured.

Zoom In Zoom Out Reset image size

The Mg ii emission line plays an important role in quasar studies, as it is used increasingly for RM estimates of M_BH (e.g., Shen et al. 2016), but also for virial (single-epoch) estimates when Hβ becomes inaccessible from the ground (z > 1; e.g., Vestergaard & Osmer 2009; Wang et al. 2009; Shen et al. 2011).

During our visual inspection of these multi-epoch spectra of strongly variable quasars, we noted that the variation in the strength of the Mg ii emission line is generally much lower than that of Hβ. For instance, the Mg ii line may remain relatively strong even in the dim state, after the continuum and broad Hβ emission have significantly faded. The relatively weak response of Mg ii to continuum changes has been noted repeatedly in the literature (e.g., Clavel et al. 1991; Kokubo et al. 2014; Cackett et al. 2015; Zhu et al. 2017).

Homan et al. (2020) presented a large study of Mg ii variability in extremely variable (∣Δg∣ > 1) quasars, and found considerable complexity in different quasars' Mg ii response to continuum changes, both for line strength and profile. The response of Mg ii to continuum variations was different not only between quasars, but even between epochs in individual quasars.

In an RM study, Sun et al. (2015) found that the Mg ii line tends to vary less than Hβ, with a broader response function, and that no clear Mg ii radius-luminosity relation may exist for Mg ii. This could be because while the Balmer lines are recombination dominated, at BLR densities, Mg ii is predominantly collisionally excited. Furthermore, the Mg ii BLR may not necessarily end where the photoionization conditions are no longer optimal, but rather at the dust sublimation radius Guo et al. (2020).

Interestingly, Roig et al. (2014) searched among some 250,000 luminous galaxies with SDSS/BOSS spectra in the redshift range 0.35 < z < 1.1 and found 293 examples (∼0.1%) with strong, broad Mg ii lines, but lacking a blue continuum or any broad Balmer line emission. If not for their broad Mg ii emission, these spectra would normally be classified as Seyfert 2 or LINER galaxies. If the Mg ii BLR is larger than the Hβ BLR, these could be dim state CLAGN where the Mg ii region remains illuminated. Indeed, Mg ii CLAGN have subsequently been discovered by Guo et al. (2019).

Significant broad Hα emission sometimes remains after Hβ disappears in many CLQs, as evident in some of the z < 0.4 CLQs shown in Figure 1 and also in, e.g., LaMassa et al. (2015), Yang et al. (2018), and Sheng et al. (2020). By contrast, in the luminous SDSS/BOSS sample of galaxies with broad Mg ii emission (Roig et al. 2014), there is little evidence for broad Balmer line emission. They stacked 162 broad Mg ii galaxies with 0.35 < z < 0.57 (covering both Mg ii and Hα) and detected only narrow Balmer line components (250 km s⁻¹ for Hα and 150 km s⁻¹ for Hβ). This implies that the primary broad line emitting region size increases from Hβ to Hα, to Mg ii. Grier et al. (2017) studied Hα and Hβ lags in the SDSS RM project and found the Hα lags to be ∼40% longer than for Hβ (see their Figure 10). Homayouni et al. (2020) measured Mg ii lags and also found they exceed Hβ lags by a similar amount.

In Figure 9 we plot the change in Mg ii against continuum, in both luminosity and flux; evidently, the trend is indeed weak.

Zoom In Zoom Out Reset image size

In Figure 10, we show the best-fit broad line equivalent width W_λ versus the underlying continuum luminosity for both Hβ and Mg ii. ³² Treating every spectrum as if it were an individual QSO, there is a weak, possibly positive trend in the overall Hβ eBeff, which has been noted before for Balmer lines in QSO samples (Greene & Ho 2005; Rakić et al. 2017). The linear equation in this sample of spectra for the Hβ eBeff is $\mathrm{log}\,{W}_{\lambda }({\rm{H}}\beta )=0.29\pm 0.06\,\mathrm{log}({L}_{5100})-10.9\pm 2.5$.

Zoom In Zoom Out Reset image size

This is expected, since broad Hβ emission reacts (albeit with a lag) directly and almost linearly to continuum changes. The Hβ plot shows extreme scatter, largely because the primary selection for our CLQ candidate sample is to have extreme variability in continuum and broad Hβ emission. Expected sources of eBeff scatter in quasar samples can likely be attributed to differences between quasars in their inclination angles, BLR, absorber geometries, SMBH masses, and accretion rates. In contrast, for the iBeff, only the last factor changes substantially in any given QSO. Yet here, the Hβ intrinsic slopes, shown in Figure 10 by dashed lines between our designated dim and bright states (see Table 2), seem to show at best a weak trend, and many of the lines are nearly vertical. This is again likely due to our a priori selection of extremely variable Hβ strength. Indeed, the average slope of the Hβ iBeff is much steeper, at 2.05 ± 0.07.

The line best fitting the ensemble Mg ii Beff in our full candidate CLQ sample is $\mathrm{log}\,{W}_{\lambda }(\mathrm{Mg}\,\mathrm{II})=-0.32\,\pm 0.05\mathrm{log}({L}_{2798})+16.0\pm 2.4$. In contrast, the average slope best fitting the iBeff is much steeper, at −0.828 ± 0.001. Steeper slopes for the iBeff compared to the eBeff are well known, and indicate a subdued response to continuum changes. This is especially true, and clearly visible in Figure 10 for Mg ii. Yang et al. (2020) investigated the variability of the Mg ii line in a sample of EVQs, finding that, in contrast to the Balmer lines, the FWHM of broad Mg ii does not react strongly to continuum changes. In Figure 11, we demonstrate a slight positive trend of broad Mg ii width with luminosity, indicating that the Mg ii broad lines react differently than does Hβ.

MacLeod et al. (2019) and Rumbaugh et al. (2018) have argued that the EVQ and (even more so) CLQ populations have lower Eddington ratios L_bol/L_Edd than do normal quasars. However, since their CLQ samples were selected starting from criteria of photometric variability, it is interesting to compare here with our primarily ³³ spectroscopically selected sample.

Zoom In Zoom Out Reset image size

The Eddington luminosity is L_Edd = 1.26 × 10³⁸ M_BH in erg s⁻¹, with M_BH in solar units. For our TDSS sample, we derive M_BH separately from our PyQSOFit model fits to broad Hβ in all spectra. Assuming that the BLR is virialized, we use

$$\begin{eqnarray*}&&{M}_{\mathrm{BH}}={f}_{\mathrm{BLR}}\displaystyle \frac{{R}_{\mathrm{BLR}}\,{\rm{\Delta }}{V}^{2}}{G},\end{eqnarray*}$$

where f_BLR is a virial factor to characterize the kinematics, geometry, and inclination of the BLR clouds, R_BLR is the characteristic radius of the BLR, ΔV is the velocity of the BLR clouds, usually traced by the FWHM or the line dispersion σ_Hβ of the broad Hβ line, and G is the gravitational constant. Specifically, we use FWHM from the PyQSOFit output Hb_whole_br_fwhm. Errors in single-epoch virial estimates are larger than for the more robust RM measurements (Peterson 1993; Shen et al. 2016; Grier et al. 2017), but can be reduced by taking into account the mass accretion rate, as in Yu et al. (2020). Their best fit for the virial factor, which incorporates the effects of BLR geometry, orientation, and kinematics, is

$$\begin{eqnarray*}&&\mathrm{log}\,{f}_{\mathrm{BLR}}=-1.1\,\mathrm{log}\,[\mathrm{FWHM}({\rm{H}}\beta )/2000]+0.48.\end{eqnarray*}$$

For the characteristic radius of the BLR, we use their best fit:

$$\begin{eqnarray*}&&\mathrm{log}\,{R}_{\mathrm{BLR}}=0.42\,\mathrm{log}\,{L}_{44}-0.28{R}_{\mathrm{Fe}}+1.53,\end{eqnarray*}$$

where L₄₄ is our L₅₁₀₀ continuum fit values from Table 2 in units of 10⁴⁴ erg s⁻¹, R_Fe is the ratio of iron to Hβ flux (taken from PyQSOFit outputs Hb_whole_br_area and fe_op), and R_BLR is in lt-days. The resulting M_BH estimates are shown plotted against redshift in Figure 12 for the designated bright and dim state for each QSO in our sample (red and blue points, respectively), with those values connected by a dashed line. The apparent overall correlation between logM_BH and redshift is due to a combination of the SDSS flux limit (which due to volume sampled and the quasar luminosity function induces a correlation between luminosity and redshift) and the known correlation between quasar luminosity and M_BH. The sometimes significant differences in mass estimates for a single quasar glaringly highlight the uncertainties in single-epoch virial mass estimates, expected to be particularly severe in a CLQ sample defined by large changes in both luminosity and broad Hβ line emission. However, the epochal differences in virial mass estimates do not generally exceed the individual uncertainties.

Zoom In Zoom Out Reset image size

To calculate bolometric luminosity, we use L_bol =8.1λ L _5100Å from Runnoe et al. (2012). Figure 13 shows the distribution of Eddington ratios as logL/L_Edd for several different samples. Relative to most previous CLQ samples initially selected explicitly on large amplitude variability (e.g., MacLeod et al. 2016, 2019), our TDSS sample can be thought of as primarily spectroscopically selected; candidates are identified by visual inspection of multi-epoch spectra. While the initial TDSS selection for spectroscopic targeting was indeed by variability, the selection thresholds were at levels associated with typical quasar (∼0.2 mag; Morganson et al. 2015) rather than more extreme (∼1 mag) variability. However, as demonstrated in Table 1, the TDSS repeat spectroscopy programs provide a good portion of our CLQ candidates, especially the Disk Emitter (DE) and hypervariable QSO (HypQSO) FES subprograms. The latter objects generally required a variability measure ≳0.5 mag (MacLeod et al. 2018, Figure 5.).

Zoom In Zoom Out Reset image size

Eddington ratio distributions for several quasar samples in Figure 13 are shown as histograms, normalized to have the same total area. To similar redshifts as our CLQ sample (z ≤ 0.9), the general quasar population shows a broad distribution peaking at log L/L_Edd ∼ −0.8 or about 10%–15% of L_Edd. Extremely variable (∣Δg∣ > 1) quasars from Rumbaugh et al. (2018) show an even broader distribution, but centered at lower L/L_Edd values near about 4%. The bright state L/L_Edd values for the sample of (29) CLQs from MacLeod et al. (2019) share a similar range and mean as the EVQs, and to our TDSS CLQ sample, all centered near Eddington ratios of a few percent, about twice as low as L/L_Edd for the broader quasar distribution. This may be evidence that the accretion rate is more variable at lower L/L_Edd, yielding stronger variability and occasional darkening of the BLR. However, it may also be that the FWHM of broad Hβ emission, used to derive M_BH, is typically larger in EVQs and CLQs for other reasons, thereby yielding lower estimates of L/L_Edd. Indeed, Ren et al. (2021) find that EVQs show a subtle excess in the very broad line component compared with control samples. Strong turbulence in the inner accretion disk may be the reason for the continuum variability (e.g., Kelly et al. 2009; Cai et al. 2020) and may also launch more gas into the inner BLR, where the broadest line component is formed. We also plot in Figure 13 histograms of these same samples, but mapped to the disk wind parameter of Elitzur & Shlosman (2006) and Elitzur et al. (2014) that determines whether or not a BLR can form. The red vertical dashed line indicates the predicted critical value above which BELs should be observable in this model. All the dim state CLQs lie below it, as the model predicts.

6.1. Diversity

6.2. Classification

6.3. Future Work

As could be said for many astrophysical experiments in the time domain, the ideal campaign to understand the CLQ phenomenon would consist of rapid cadence, multi-epoch, multiband observations of a large sample of quasars, observed over a long time span. Even well short of the ideal, spectroscopy in the time domain is already proving to have profound effects on our interpretation of quasar classes. Concerted time domain spectroscopic surveys for CLQs have certain advantages over surveys where follow-up spectroscopy is triggered by photometric or color variability. Spectroscopy allows a focus directly on the BLR emission, and so can even detect CLQs where the observed continuum does not change substantially (e.g., in cases where the BLR may be subject to a different ionizing continuum than evident along the observer's sightline) or is mired in host galaxy emission. Repeat spectroscopy surveys, also can find CLQs where the AGN luminosity never dominates the host galaxy emission.

Many more spectroscopically identified CLQs will be discovered during SDSS-V, which is just beginning, and is expected to last about 5 yr. The BH Mapper program of SDSS-V (Kollmeier et al. 2017) includes a significant program of repeat quasar spectroscopy dubbed "AQMES" for All Quasar Multi-epoch Spectroscopy. The AQMES Medium program targets about 2000 known QSOs to a given magnitude limit (of i = 19.1) for about 12 epochs of spectroscopy over about 300 deg² of sky. AQMES Wide will similarly target about 20,000 known quasars to a similar magnitude limit, but for only a few (1–3) epochs over about 3000 deg² of sky. Combining SDSS-V quasar spectra with prior SDSS archival spectra of the same objects will afford a broad survey of continuum and emission line changes spanning the full parameter space of redshift, luminosity, SMBH mass, Eddington ratios, and rest-frame timescales for a variety of emission lines. Many more CLQs will be discovered but, perhaps more importantly, they will be placed in the broader context of quasar spectroscopic variability. Furthermore, better statistics along with the availability of detailed light curves across the sky should allow for enhanced constraints on the distribution of these parameters, to determine whether CLQs represent a uniquely variable species, the extreme tail of a smooth distribution of variability, or simply normal quasars in a highly variable phase. The first interpretation might require, e.g., evidence for distinct properties apart from variability itself, such as distinct host galaxy features, or star formation histories. The power to discriminate between any of these interpretations benefits from increased length and cadence of photometric and spectroscopic monitoring.

Our adoption of the N_σ (Hβ) > 3 criteria has been to be conservative in the presence of large samples of variable S/N spectroscopic data, and to create a numerical criterion that might be adopted by the community to facilitate easy aggregation and/or comparison of samples. However, if for example, we had spectra of S/N = 100, a 3σ change could admittedly represent quite a minor change in broad Hβ luminosity. Ideally, criteria to declare identification of a bona fide CLQ should be based on intrinsic properties, measured with adequate significance. Based on our Figure 8 (and similarly, MacLeod, Figure 3), among CLQ candidate samples, there does not appear to be a clear lower limit to Δ L(3240) and Δ L(Hβ) that could provide intrinsic criteria. Rather, criteria involving fractional changes might make the most sense. Figure 14 shows the fractional change in both Δ L(3240) from the power-law fit and Δ L(Hβ) for our candidate CLQ sample, including all spectral epochs contrasted against the designated dim epoch for each CLQ candidate. Contrasted spectral pairs with N_σ (Hβ) > 3 are shown in red. The blue dashed lines indicate a 30% change in each axis. A sample of typical quasars with multi-epoch spectral analysis would be expected to cluster near the origin of this plot. Since the full sample shown here was selected to be candidate CLQs by visual inspection, the majority of fractional changes between spectral pairs are found in the upper right corner, where they may meet our suggested intrinsic criteria for CLQ if the S/N of the spectra in the region of Hβ is also sufficient.

Zoom In Zoom Out Reset image size

Nearly all the CLQs designated as such only by the N_σ (Hβ) criterion in this paper also meet the other two criteria (fractional change in both L(3240) and broad Hβ > 30%). One exception (left of the blue dashed box) is J224113.54-012108.84, comparing MJDs 58367 and 55824. As we remark in the Appendix, this QSO has unusual continuum properties and perhaps a spectroscopic calibration issue. Another exception is J100302.62+193251.28, comparing MJDs 57817 and 53762, which appears below the blue dashed box, so a fractional change in broad Hβ than prescribed, although we find N_σ (Hβ) = 4.31. We note that both of these exceptions are turn-on CLQs, which remain underrepresented in the literature and deserve further study.

The spectral energy distributions (SEDs) of quasars are known to change significantly with Eddington ratio. More luminous objects have a softer ionizing continuum. The impact of SED changes on line-forming efficiency and even the structure of the BLR is well documented (Zheng & Malkan 1993; Green 1998; Korista et al. 1998). Clearly, the extremes of accretion rate occurring in CLQs may also correspond to significant changes in SED. Optical, UV, and X-ray monitoring of CLQs during transitions is key to understanding these changes, and their propagation from the central engine to the BLR. A few nearby Seyferts have been monitored closely throughout state changes, e.g., Krumpe et al. (2017) for Mkn1018 (z = 0.043) and Trakhtenbrot et al. (2019b); Ricci et al. (2020) for 1ES 1927+65 (z = 0.019), but these have 10–100 times lower SMBH masses, and high Eddington ratios compared to our CLQ sample. Some efforts to chart SED changes in CLQs, including X-rays have been made, but most have relied on asynchronous observations (e.g., Ruan et al. 2019). A Chandra program (PI: P. Green) is in progress, catching CLQs in optical and X-ray in both bright and dim states. We may expect many more CLQs to afford real-time multiwavelength monitoring in coming years via eROSITA (Predehl et al. 2021) in X-rays and the Rubin Observatory/LSST in the optical (Ivezić et al. 2019).

Many thanks to Qian Yang for comments and suggestions.

S.F.A., P.J.G., C.L.M, B.R., and J.J.R. are supported by the National Science Foundation under grant Nos. AST- 1715121, and AST-1715763. C.L.M. and Q.Y. were partially supported for spectroscopic observations and analysis by the National Aeronautics and Space Administration through Chandra Award Nos. GO8-19089X and GO9-20086X, issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract No. NAS8-03060.

Photometric light curves used herein are based on observations obtained with the Samuel Oschin 48 inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under grant No. AST-1440341 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and the Lawrence Berkeley National Laboratory. Operations are conducted by COO, IPAC, and UW.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics ∣ Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

The authors acknowledge the use of the Catalog Archive Server Jobs System (CasJobs) service at http://casjobs.sdss.org/CasJobs, developed by the JHU/SDSS team.

This work uses data obtained with the 1.2 m Samuel Oschin Telescope at Palomar Observatory as part of the Palomar Transient Factory project, a scientific collaboration among the California Institute of Technology, Columbia University, Las Cumbres Observatory, the Lawrence Berkeley National Laboratory, the National Energy Research Scientific Computing Center, the University of Oxford, and the Weizmann Institute of Science. The Zwicky Transient Facility is supported by the NSF under grant No. AST-1440341 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW.

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE).

This work made use of data products from the Catalina Sky Survey. The CSS survey is funded by the National Aeronautics and Space Administration under grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the U.S. National Science Foundation under grant Nos. AST-0909182 and AST-1313422.

Facilities: Sloan - Sloan Digital Sky Survey Telescope, ING:Herschel - , Magellan:Baade (LDSS2 imaging spectrograph - , Boller & Chivens spectrograph) - , MMT (Blue Channel spectrograph - , BinoSpec) - , PS1 - , PO:1.2 m - , SO:Schmidt - , SO:1.5 m - , SO:1m - .

Software: Astropy (Astropy Collaboration et al. 2013, 2018), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), TOPCAT (Taylor 2005).

Here, we discuss the individual CLQ light curves and spectra shown in Figure 1. All the CLQs but four are turn-offs, dimming with time, unless noted otherwise below. We include a quantitative estimate of AGN spectral type, an adaptation of the Lick classification scheme (Osterbrock 1981), offered by Winkler (1992) using the flux ratio R = flux(total Hβ)/flux(total [O iii]λ5007) to define Types 1, 1.2, 1.5, and 1.8 with R values R > 5, 2.0 < R < 5.0, 0.333 < R < 2.0, and R < 0.333, respectively. For Type 1, permitted lines are all broad. Type 1.5 has broad components about equal in flux to narrow. Type 1.9 has no detectable Hβ BEL, while Type 2 has neither Hβ nor Hα BELs detectable, i.e., permitted and forbidden lines have similar widths near ∼1000 km s⁻¹ . Types 1.2 and 1.8 correspond simply to finer divisions in the relative strength of broad and narrow lines. Note that Hα is generally not within the SDSS spectral range for z > 0.4.

002311.06+003517.53: This CLQ presents several strong luminosity changes over a nearly 20 yr time span in the observed frame. Brightening of 1.5 mag in the g band occurs within a few months. Later, steady dimming of 1.5 mag occurs over about a 3 yr timeframe. Six spectral epochs confirm strong accompanying changes in the continuum and broad line emission. From bright to dim state, the spectrum effectively changes from a Type 1.2 to a Type 1.5.

020514.77-045639.74: The observed photometric epochs do not coincide well with the dimmest spectral epoch, but dimming by nearly a magnitude in the g band is evident in the light curves. The continuum shape does not change as strongly as most CLQs, and Hα weakens much less than Hβ. From bright to dim state, the spectrum effectively changes from a Type 1.0 to a Type 1.5.

021259.59-003029.43: This quasar was targeted by TDSS both in the HYPQSO (highly variable) and DE (disk emitter) samples. It dimmed by about 1.2 mag in the r band, and likely more than that in the g band. The two latter, dim state spectra are essentially identical. There appears to have been a re-brightening by nearly a magnitude after the second dim state spectrum, on a timescale of about 6 months or less. In all the spectra, several of the higher order Balmer lines are double peaked. In the bright state spectrum, there is a pronounced red wing to the broad Hβ emission. This quasar's classification remains at Type 1.5.

021359.79+004226.81: This quasar dimmed in the g band by about 1.5 mag, and about half that in the r band. Again, all the higher order Balmer lines are double peaked. In the intermediate state spectrum, the Hβ broad emission has broad shoulders that disappear entirely in the dim state spectrum. From bright to dim state, the spectrum effectively changes from a Type 1.5 to a Type 1.8.

024508.67+003710.68: Here the g band dims by about 0.5 mag, while the r-band continuum remains steady. There are likely flux calibration issues with the later epoch spectrum. Nevertheless, it appears that essentially all broad line emission disappears. From bright to dim state, the spectrum effectively changes from a Type 1.2 to a Type 1.9.

024932.01+002248.35: This quasar suffers dimming of about Δg = 1.5 mag, whereas the r-band continuum dims by about half that amount. The blue continuum can be seen to almost entirely disappear in the spectra. Some broad Hβ remains. Details of these spectra and model fitting are presented separately in Figures 2 and 3. From bright to dim state, the spectrum effectively changes from a Type 1.0 to a Type 1.5.

100302.62+193251.28: Judging from the spectra, this CLQ appears to have turned on, though it is difficult to gauge the change in magnitude, since the later brighter state spectrum has no nearby photometric epoch. The strong O iii lines in this object indicate a very bright narrow line region, likely illuminated from past bright epochs. Broad emission from Hβ declines markedly, yet Mg ii remains strong. From bright to dim state, the spectrum stays at Type 1.5.

105058.42+241351.18: Here, the g band likely dims by about half a magnitude. While Hα decreases slightly, all the higher order Balmer lines seem to disappear. From bright to dim state, the spectrum effectively changes from a Type 1.2 to a Type 1.5.

105325.40+302419.34: The g band dims monotonically by more than a magnitude over about a dozen years. Broad Hβ nearly disappears; broad Hα weakens significantly as well, but the final state may be Type 1.9. From bright to dim state, the spectrum also effectively changes from a Type 1.2 to a Type 1.5.

105513.88+242553.69: The g band dims monotonically by about 2 mag, with an accompanying striking change in the blue spectroscopic continuum. Hβ, Hγ, and Hδ all disappear. From bright to dim state, the spectrum effectively changes from a Type 1.0 to a Type 1.9. It may end up as a Type 2 quasar, but at z = 0.496, Hα is not visible in the spectrum to make that determination.

111329.68+531338.78: The g band likely dimmed by about a magnitude between bright and dim states, but spectral epochs fall within gaps in photometric coverage. All the visible Balmer lines show a marked decrease in flux. From bright to dim state, the spectrum effectively changes from a Type 1.0 to a Type 1.2.

113651.66+445016.48: In this turn-on CLQ, the g-band magnitude brightens by about 0.5 mag with a strong upturn in the blue spectral continuum. Broad Hβ and Hγ strengthen considerably, as does Fe ii emission, but with little if any change in Hα. Changing from a dim to bright state, the spectrum effectively changes from a Type 1.5 to a Type 1.2.

113706.93+481943.68: This CLQ seems to have dimmed steadily by about a magnitude in the g band over about a decade, during which broad Hβ emission disappeared. From bright to dim state, the spectrum effectively changes from a Type 1.2 to a Type 1.5.

135415.54+515925.77: This CLQ seems to show large rapid variability even in the r band superimposed on a steady overall dimming. The g band dimmed very significantly, by about 2 mag. The higher order Balmer lines appear to show unusual double narrow peaks in the dim state. The spectral classification remains at Type 1.5.

143455.30+572345.10: This CLQ dimmed by a full 2 mag in the g band, showing a spectacular disappearance of the BEL components, which began in the bright state with strong red wings. From bright to dim state, the spectrum effectively changes from a Type 1.2 to a Type 1.8.

163620.38+475838.36: There is a complete disappearance of broad Hβ emission in this CLQ, while the continuum dimming is rather quite modest. Judging from the photometric points, the early epoch spectrum would likely have been much more spectacular had it been obtained about a year earlier. From bright to dim state, the spectrum effectively changes from a Type 1.5 to a Type 1.8.

224113.54-012108.84: This is a turn-on CLQ that brightened by a full magnitude in the g band, as is clear in the <4000 Å continuum. Four of the five spectra show N_σ (Hβ) > 3 relative to the dim state. Hα appears to maintain its broad emission throughout, and retains a Type 1 classification for all its spectral epochs. The continuum changes are quite unusual, since the reddest spectrum MJD 55824, which we classify as the dim state, is actually brighter in the red than in the bright bluer states. There may be some calibration issues with the SDSS spectra.

231625.39-002225.50: For this CLQ, the continuum dimming between spectral epochs is quite modest—likely no more the Δg ∼ 0.2 mag. With N_σ (Hβ) = 3.16, these epochs barely pass our CLQ criterion, but broad Hβ clearly diminishes significantly by visual inspection. The spectrum remains a Type 1.5 throughout.

234623.42+010918.11: This turn-on CLQ brightened in the g band by at least 1.5 mag, within about a year in the observed frame, and remained bright. However, if the flux calibrations are correct, there were relatively rapid and significant changes in the spectral continuum blueward of about 4000 Å in the six later epochs of spectroscopy. The spectrum remains a Type 1.5 throughout.