Understanding extreme quasar optical variability with CRTS: II. Changing-state quasars

We present the results of a systematic search for quasars in the Catalina Real-time Transient Survey exhibiting both strong photometric and spectroscopic variability over a decadal baseline. We identify 73 sources with specific patterns of optical and mid-IR photometric behavior and a defined spectroscopic change. These"Changing-State"quasars (CSQs) form a higher luminosity sample to complement existing sets of"Changing-Look"AGN and quasars in the literature. The CSQs (by selection) exhibit larger photometric variability than the CLQs. The spectroscopic variability is marginally stronger in the CSQs than CLQs as defined by the change in H$\beta$/[OIII] ratio. We find 36 sources with declining H$\beta$ flux, 37 sources with increasing H$\beta$ flux and discover seven sources with $z>0.8$, further extending the redshift arm. Our CSQ sample compares to the literature CLQ objects in similar distributions of H$\beta$ flux ratios and differential Eddington ratios between high (bright) and low (dim) states. Taken as a whole, we find that this population of extreme varying quasars is associated with changes in the Eddington ratio and the timescales imply cooling/heating fronts propagating through the disk.


INTRODUCTION
Quasar variability is generally regarded as a stochastic process. The summation of activity associated with accretion disk instabilities, ionizing continua, jets, stellar activity close to the core, and dust clouds are all potential contributors. Sparse studies, either in terms of sample size or temporal sampling, have produced a simple statistical model, the damped random walk (Kelly et al. 2009), which aims to describe this variability (but see Zu et al. 2013;Graham et al. 2014;Kasliwal et al. 2015;Kozlowski 2017;Smith et al. 2018;Moreno et al. 2018, for counter-discussions). The grow-E-mail:mjg@caltech.edu (MJG) ing availability of large collections of rich multiepoch data is, however, enabling a much more phenomenological approach. Systematic studies of the quasar population (or substantial fractions thereof) are now possible with different characterizations of variability in terms of discriminative features or statistical models. These aim to capture specific patterns of behavior associated with particular underlying physical processes. In this way, we have identified sources exhibiting periodic activity (Graham et al. 2015a,b), major flaring (Graham et al. 2017;Drake et al. 2019), and extreme broad line variability (Stern et al. 2017Ross et al. 2018). emission lines (BELs, prototypically Hβ) in their optical spectra and often large, order of magnitude changes in the optical photometry (LaMassa et al. 2015;Ruan et al. 2016;Runnoe et al. 2016;MacLeod et al. 2016;Gezari et al. 2017;Runco et al. 2016;Yang et al. 2018;Assef et al. 2018;Stern et al. 2018;Wang, Xu & Wei 2018;Ross et al. 2018;MacLeod et al. 2019). Such changing-look quasars (CLQs) are consistent with a change of spectral type (broad-lined to narrow-lined or vice versa) and may, in principle, be associated with a large change of obscuration, accretion rate, or accretion disk luminosity. Microlensing is also a potential cause of the CLQ phenomena. The term "changing-look quasar" is borrowed from X-ray astronomy where large changes in X-ray luminosity have been shown to be typically associated with varying absorption, e.g., Matt et al. (2003); Rivers et al. (2015a,b). Similar significant spectral variability has also been known for many years in a number of local low-luminosity AGN (e.g., Khachikian & Weedman 1971;Tohline & Osterbrock 1976;Penston & Perez 1984;Cohen et al. 1986;Bischoff & Kollatschny 1999;Aretxaga et al. 1999;Eracleous & Halpern 2001;Shappee et al. 2014;Denney et al. 2014;Li et al. 2015;Parker et al. 2016). We note that significant photometric and spectroscopic variability have also been detected in quasars with absorption lines systems, e.g., the extreme BAL QSO reported by Stern et al. (2017), but such objects are not normally considered as CLQs.
Although this type of behavior has so far seemed rare, it may well be that data sets are only now sufficient in size and temporal coverage to effectively detect such activity. Runco et al. (2016) studied 102 local (0.02 ≤ z ≤ 0.1) Seyfert galaxies with M BH > 10 7 M and found that ∼66% of objects showed variability (change in values) in width or flux in Hβ over a 3-9 year timeframe and Hβ completely disappeared in three sources. From a study of SDSS DR7 and DES Y3A1 data for 8,640 quasars, Rumbaugh et al. (2017, hereafter, R17) found that ∼ 10% exhibited extreme variability (|∆g| > 1 mag) sometime within a 15-year baseline (the actual distribution of the restframe baseline over which the maximum g-band variability was observed peaks around 1000 days but is largely insensitive to timescales beyond ∼1500 days). Correcting for selection incompleteness, they speculate that 30-50% of all quasars may actually show such behavior on these timescales.
The link between extreme spectroscopic and photometric variability is not clear, however. Despite R17's suggestion that extreme variable quasars (EVQs, |∆g| > 1) are good candidates for CLQs, a large amplitude photometric variation alone is not enough to identify them. On the one hand, we have previously reported extreme variability that is not associated with significant spectroscopic changes (Graham et al. 2017). On the other hand, while MacLeod et al. (2016) report on a sample of 10 CLQs with |∆g| > 1, Yang et al. (2018) present a sample of 21 CLQs, 15 of which have |∆g| < 0.5 (their full sample spans 0.03 ≤ |∆g| ≤ 1.89). Therefore, observations show that the spectroscopic CLQ and photometric EVQ phenomena are not directly correlated. What does seem to be indicated is that there are specific patterns of variable behavior that are likely associated with CLQs: for example, Lawrence et al. (2016, hereafter, L16) identify a sample of 15 quasars, including one CLQ, showing slow steady changes over several years, which they attribute to a mixture of changes in accretion state and microlensing. It is also unclear whether CLQs stand out from single epoch spectroscopy. L16 report that slow blue hypervariables have weaker Mg ii and [O iii] emission lines. However, R17 find that EVQs have lower Eddington ratios and larger Mg ii and [O iii] equivalent widths (EWs) than control quasars matched in redshift and luminosity. Although L16 and R17 form a superset containing some CLQs, they do suggest that quasars with lower accretion rates are more susceptible to changes in accretion rate and exhibiting more extreme behavior.
Strong correlations are reported as well between CLQ behavior and mid-infrared (MIR) variability (Sheng et al. 2017;Yang et al. 2018;Assef et al. 2018) with optical and MIR colors also changing with flux variation: a bluerwhen-brighter chromaticism in the optical and redder-whenbrighter in the mid-infrared. Given the pc-scale size of the MIR-emitting region, this clearly indicates that the strong variability is not due to an obscuring screen. Instead the chromatic trends are likely due to less host galaxy contributions and a stronger inner accretion disk contribution when quasars are more luminous.
Further evidence against obscuration being the primary cause of CLQs comes from optical polarimetric studies. If the disappearance of the broad emission lines originates from the obscuration of the quasar core by dusty clouds moving in the torus, high linear optical polarization would also be expected. Measurements of the polarization of CLQs (Hutsemekers et al. 2017(Hutsemekers et al. , 2019Marin 2017) are less than 1%, which suggests that the phenomenon is not due to obscuration but physical changes in the accretion disk and/or accretion rate. Such low polarization degrees indicate as well that these quasars are seen under inclinations close to the system axis. Finally, imaging of the host galaxies of four faded CLQs (Charlton et al. 2019) suggests that these are predominantly disrupted or merging galaxies that resemble AGN hosts, rather than inactive galaxies.
In this work, we present a search for quasars showing photometric and spectroscopic variability consistent with a change of state of activity. We will refer to these as changingstate quasars (CSQs) rather than changing-look quasars. Though the latter has been the conventional term in the literature to date, it is ill-defined with no clear phenomenology, be it photometric or spectroscopic, associated with it beyond "significant" variability. In the optical community, literature often ascribes variability to either changes in obscuration (as in the X-ray community) or changes in accretion rate, and it has been a challenge to identify physical mechanisms that would lead to accretion rate changes. The significant MIR variability associated with the optical variability in these sources largely rules out the obscuration scenario as already noted. More recently, several papers have noted that the variability events occur on thermal timescales Ross et al. 2018;Noda & Done 2018;Parker et al. 2019), implying that the luminosity changes are likely associated with rapid changes in the temperature of the accretion disk. In contrast, changes in the accretion rate would be expected on the viscous timescale, which is more than an order of magnitude longer (i.e., decades/centuries rather than years). CSQ is therefore the more appropriate term but we will continue to refer to those objects already identified in the literature as CLQs (see Table A1 and Fig. A1 for those covered by the data used in this work). This paper is structured as follows: in Section 2, we summarize the data sets used and in Section 3, we present the selection technique and criteria for identifying CSQs. We discuss our results in Section 4. Section 5 considers implications for the physical mechanisms behind the variability. We assume a standard WMAP 9-year cosmology (Ω Λ = 0.728, Ω M = 0.272, H 0 = 70.4 km s −1 Mpc −1 ; Jarosik et al. 2011) and magnitudes are approximately on the Vega system.

DATA SETS
A systematic search for sources showing the particular behavioral patterns associated with CSQs requires a data set with a long temporal baseline and also a high sampling rate. Some candidates may be identified from relatively few epochs of data spread over a roughly decadal baseline, but such data sets will typically have insufficient resolution or sensitivity to detect specific forms. The Catalina Real-time Transient Survey (CRTS; Drake et al. 2009) represents the best data currently available with which to systematically define sets of quasars with particular temporal characteristics.

CRTS
The CRTS archive 1 contains the Catalina Sky Survey data streams from three telescopes -the 0.7 m Catalina Sky Survey (CSS) Schmidt and 1.5 m Mount Lemmon Survey (MLS) telescopes in Arizona and the 0.5 m Siding Springs Survey (SSS) Schmidt in Australia. These surveys, operated by the Lunar and Planetary Laboratory at the University of Arizona, were designed to search for near-Earth objects, but have proven extremely valuable for astrophysics topics ranging from Galactic transients (Drake et al. 2014) to distant quasars (Graham et al. 2014(Graham et al. , 2015b(Graham et al. , 2017. CRTS covers up to ∼2500 deg 2 per night, with 4 exposures per visit, separated by 10 min. The survey observes over 21 nights per lunation. The data are broadly calibrated to Johnson V (see Drake et al. 2013 for details) and the current CRTS data set contains time series for approximately 400 million sources to V ∼ 20 above Dec > −30 from 2003 to 2016 May (observed with CSS and MLS) and 100 million sources to V ∼ 19 in the southern sky from 2005 to 2013 (from SSS).
The error model used for CRTS DR2 is incorrect: errors at the brighter magnitudes are overestimated and those at fainter magnitudes (V > 18) are underestimated (Palaversa et al. 2013;Drake et al. 2014). In this analysis, we employ the improved error model derived in Graham et al. (2017); the actual CRTS error model will be fixed in a future release. We also note that none of the sources we consider have fewer than 10 observations in their light curve. We apply the same preprocessing steps described in Graham et al. (2015b) to all light curves, which remove outlier photometry points and combine all exposures for a given night to give a single weighted value for that night. We also remove sources associated with nearby bright stars or identifiable as blends 1 http://catalinadata.org from a combined multimodality in their magnitude and observation position, i.e., the spatial distribution of all points in a light curve is best described by n > 1 Gaussians.

WISE
In this paper we use MIR W1 (3.4µm) and W2 (4.6 µm) WISE data from the beginning of the mission in 2010 January through 2017 December, corresponding to the fourth year of NEOWISE operations 2 . Note that there is a gap between 2011 February and 2013 September when the satellite was in hibernation. For most sky positions, there are ∼12 observations of a source over a ∼1 day period with a sixmonth gap between repeat visits. We combine all exposures for a given 24 hour period with a signal-to-noise ratio greater than five to produce a single value using the same method as for CRTS data.

Spectroscopically confirmed quasars
The Million Quasars (MQ) catalogue 3 v5.2 contains all spectroscopically confirmed type 1 QSOs (577,146), AGN (30,062) and BL Lacs (1,615) in the literature up to 2017 August 5. Previous versions have formed the basis for the results of Graham et al. (2015b) and Graham et al. (2017). MQ (v5.2) also contains 1,297,111 photometric quasar candidates from SDSS or WISE. We crossmatched MQ against the CRTS data set with a 3 matching radius and find that 1,411,364 sources are covered by the full CRTS. Of these, 268,202 do not have enough observations (n < 10), leaving a data set of 1,143,162 quasars and quasar candidates. We also remove 3,724 known blazars based on the class designation in MQ and the BZCAT v5.0 catalog of blazars (Massaro et al. 2015). Table 1 gives a summary of this superset of MQ sources to which we now apply our selection criteria to identify CSQs.

Optical photometry selection
Photometrically, we expect to see a gradual change in magnitude associated with monotonically varying broad emission line strengths and/or continuum changes. This would show as a strong localized trend or enhanced variability strength over some timescale in the time series.
To identify local variability in a time series, we use a Bayesian blocks (BB) representation (Scargle et al. 2013) which provides an optimal segmentation of the data in terms of a set of discontinuous piecewise constant components (see Fig. 1). This approximation makes it easier to detect significant changes of behavior in the presence of irregular sampling, noise, and gaps. In particular, it is more sensitive to coherent magnitude variations over time characterized by the difference between the first and last piecewise segments of the BB than fitting a linear trend model to the data. We assume it unlikely that a quasar will undergo a transition and return to its initial state within the timeframe of our light curves. We also distinguish between this pattern of behavior and flaring as we identified in Graham et al. (2017).
To detect stronger variability on particular timescales, we use the Slepian wavelet variance (SWV; Graham et al. 2014) which provides a measure of the relative contributions of variability at specific timescales to the total variability in a time series. We have determined the median observed SWV for quasars in magnitude bins of width ∆m = 0.25. For a given source j, we then calculate the quantity: i.e., the sum of the differences in dyadic log-log space between the source SWV and the appropriate median SWV interpolated at the timescales of the source SWV. Note that dyadic logs are used since the wavelet bands are defined in terms of base-2 widths. Fig. 2 shows the distribution of the BB difference and SWV 1 for the quasars in our data set and for known CLQs. We note that using this characterization places the majority (82%) of the known CLQs within the 95th percentile Figure 2. The distribution of the magnitude difference between initial and final Bayesian block components and the similarity (aggregate distance) to the median Slepian wavelet variance for the quasars in our data set. Known CLQs are denoted by red circles, and blue circles indicate CSQs identified in this work. Open symbols indicate low luminosity sources with M V < −23. The contours indicate the 68th and 95th percentile levels respectively for a population of 128,000 spectroscopically confirmed quasars with V < 19. The green diamonds show quasars which have exhibited significant flaring activity (Graham et al. 2017). contour. These objects have predominantly been identified from their spectroscopic variability, e.g., from dual epoch SDSS spectra where the source is classified as a galaxy in one epoch and a quasar in the other. They will thus be in a quiescent galaxy state for at least a fraction of the time period covered by their CRTS light curve (see Fig. A1) and show less photometric variability over this period than the median quasar at the same magnitude. Since our focus is on extreme variability, we only consider sources outside the 95th percentile contour and with SWV 1 > 0 as candidate objects. The latter criterion selects objects with above median variability.

Mid-infrared Selection
Strong (> 0.4 mag) MIR variability has been shown to be a characterizing property of CLQs (Sheng et al. 2017;Assef et al. 2018;Stern et al. 2018) and so we also consider this as a discriminating feature. From the maximum ∆W1 (∆W2) distribution of MIR variability for quasars over the seven years of WISE observations, ∆m = 0.4 corresponds to the 82 nd (88 th ) percentile and 3% (3%) of sources vary by more than a factor of two in flux, or ∆m > 0.75, over the period. Yang et al. (2018) argue that a reasonable selection criterion is |∆(W1 − W2)| > 0.1 when |∆(W1)| > 0.2 but they define ∆(W1) as the magnitude difference between the brightest epoch in a WISE time series and either the first or last epoch, depending on whether the CLQ is turning on or off. This means that unless the MIR flux is monotonically varying, ∆(W1) will be some fraction of the total W1 variability. It also relies on prior knowledge of whether the CLQ is turning on or off to determine the appropriate magnitude difference to measure. We employ the absolute magnitude difference in a WISE time series, i.e., W max − W min -as a selection criterion and require either |∆W1| > 0.2 or |∆W2| > 0.2 (see Table 1). The optical variability constraints described above give 65,816 candidates from the initial 1.1 million source data set (rejecting blends) and the MIR variability constraint reduces this to 47,451 sources. Of these, 14,412 have z < 0.95 and are therefore suitable for spectroscopic confirmation (i.e., with optical spectroscopy Hβ falls within the wavelength coverage of the optical spectrum). A SDSS DR14 spectrum exists for 7,576 of these and multiepoch SDSS spectra with at least 100 (500) days between epochs are available for 466 (266) objects. For comparison, there are 213,358 SDSS DR14 (Abolfathi et al. 2017) sources classed as QSOs or AGN with z < 0.95, of which 8,213 (4,244) have additional spectra taken at least 100 (500) days after the initial epoch.

Spectroscopic Selection
Over the past three years we have obtained second epoch spectra (all at least > 500 days after the initial SDSS epoch) for an additional 172 candidates (and subsequent epoch spectra for another 35 sources) using either the Double Spectrograph (DBSP) on the Hale 200" telescope at Palomar Observatory, the Low Resolution Imaging Spectrometer (LRIS) spectrograph on the Keck I telescope at the W. M. Keck Observatory, or the Echellete Spectrograph and Imager (ESI), also on the Keck I telescope (see Table B1). All these spectra were processed using standard procedures and flux calibrated with observations of spectrophotometric standard stars from Massey & Gronwall (1990) observed on the same night. We note, though, that the photometric quality of the nights was not consistent across multiple observing runs. We fit all spectra with a single power law continuum and measure Hβ and [O iii] emission line profiles relative to this. We assume a two component Gaussian fit for Hβ to model a broad and narrow component and single Gaussians for the There is no objective definition of the spectral variability required to qualify as a CLQ/CSQ in the literature. Some authors (Ruan et al. 2016;Yang et al. 2018) rely on visual inspection to identify sources with obvious broad Hβ emission in one epoch and no detection in another. MacLeod et al. (2016) and MacLeod et al. (2019) compute the flux deviation between two spectra at any given wavelength to determine the significance of a broad emission line (BEL) change and assess the significance of a change relative to the underlying continuum at that wavelength. The earlier work considers a range of significance from < 2σ change for a faint spectrum to > 8σ, but generally an absence of Hβ at one epoch is still required. The more recent work assumes a significance in the flux deviation of Hβ greater than 3. Yang et al. (2018) describe two sources transitioning from type 1 to type 1.8 where Hβ does not entirely vanish. We can thus define a measure of this deviation to set a lower limit on the change in Hβ required.
The narrow [O iii] λ5007 emission line is not expected to vary on human timescales and so we can use the change in the flux ratio of Hβ to [O iii] between epochs as a robust indicator of Hβ change, i.e., the ratio should not be significantly affected by systematic errors due to observing conditions or spectral reduction. Since there is an expectation (almost by definition) that CSQs are associated with significant spectral variability, we will consider only those sources where the absolute value of the fractional change in Hβ/[O iii] > 0.3. We also reject all spectra with a signal-tonoise ratio less than 5 as defined in Stoehr et al. (2008) using: Fig. 3 shows the distribution of the ratio for all SDSS sources with multiepoch spectra and Hβ coverage, highlighting known CLQs (bottom panel). Table 2 gives details of the 73 objects that meet our full selection criteria.

RESULTS
The light curves and spectra of our final selected sample of 73 CSQs are shown in Fig. B1. We find 36 sources with declining Hβ and 37 sources with increasing Hβ. We reject 8 sources which meet the photometric selection criteria but have atmospheric O 2 A band absorption coincident with the Hβ -[O iii] complex in the observed frame. Fig. 4 shows an example of two sources which pass the photometric selection criteria but only one of which also shows spectroscopic variability. This demonstrates that there are a variety of phenomena which may produce quantitively similar photometric variability but are distinguishable with multiepoch spectra. Note that they might also be differentiated by other observables such as X-ray or radio behavior but this requires further investigation. Figure 5 presents the light curves and spectra of CSQs identified in this paper corresponding to the largest values in each of the selection parameters: SWV 1 , ∆BB, |∆W |, and ∆(Hβ / [O iii]. From Figs 5 and B1, we find that extreme Table 2. CSQs selected in CRTS and associated features including the median CRTS magnitude (V m ), the optical amplitude (Amp) , absolute change in W 1, Bayesian block change (∆BB), Slepian wavelet variance measure (SWV 1 ), and the change in flux ratio of Hβ to [O iii]. SMBH virial mass estimates are calculated as described in Sec. 5.1.    variable sources are a heterogeneous population reflecting the complexity of the physics of accretion with subsets dominated by particular processes. Different selection techniques drawing from this population may probe different physics and this needs to be borne in mind in any analysis.
We have identified a sample of CSQs on the basis of excess optical variability over specific timescales associated with a distinct change in levels of activity over a decade, strong MIR variability, and changes in the strength of Hβ emission. We are keen to compare this study to previous and contemporary analyses but note that this is not a straightforward exercise. Table A1 and Fig. A1 describe the collective set of CLQs from the literature with CRTS light curves and publicly available spectra. Since these are the result of a variety of selection techniques, we will consider in this section the two samples -CSQs and CLQs -in terms of the selection criteria used in this work. Fig. 6 shows that the CSQ sample spans a wider redshift range than the known CLQs. The distributions of these two samples as well as the general quasar population in terms of the photometric criteria compared to the spectroscopic are shown in Fig. 7. We also distinguish between low and high luminosity sources (taking M V < −23 as the fiducial boundary). The absolute magnitude of a source is given by: where A V is the Galactic extinction, DM is the distance modulus, K V is the K-correction, and m V is the median magnitude from the CRTS lightcurve. We obtain 4 extinction values at the source position from the Schlafly & Finkbeiner (2011) recalibration of the Schlegel et al. (1998) reddening maps. We assume a K-correction of: K = −2.5(α + 1) log(1 + z) for a power law SED of F ν ∝ ν α with α = −0.5.
It is clear that the CSQ sample represents a more extreme level of photometric variability than the known CLQ  set and that the lower luminosity CLQs have a less temporally characterizable variability (smaller SWV 1 values) than their brighter counterparts (Fig. 2). This is could be likely due to an overall stronger host contribution to the optical light curve from the generally lower redshift, lower luminosity CLQs as compared to CSQs. The spectroscopic variability also seems marginally stronger in CSQs with objects transitioning from a lower to a higher state of activity, although there is no distinction between high and low luminosity sources within this group.
The difference between the initial and final states in the Bayesian block representation of the light curve correlates with the fractional change in the Hβ/[O iii] ratio, although there is no difference between low and high luminosity sources. However, the flux change associated with just Hβ (dis)appearing is not sufficient to account for the full scale of the magnitude change seen in the light curves. Using a composite quasar spectrum, such as Selsing et al. (2016), the magnitude difference for a source with V ∼ 18 at z = 0.25 (which places Hβ roughly at the peak of the CRTS equivalent filter) with and without Hβ is only ∆m ∼ 0.07 mag. An accompanying reduction in continuum flux of a third is also required to give a magnitude change of ∆m ∼ 0.5, a typical value seen in CSQs. This suggests that there will be quasars experiencing the same physical changes as CLQ/CSQs but which may only exhibit a significant photometric change (corresponding to a change in continuum flux) without showing a corresponding spectral line change, i.e., a population of CLQs without the broad emission line variability.
There is no correlation, though, between the amplitudes of the MIR photometric variability and the optical spectral variability of these sources (or the general population of quasars). This argues that the physical mechanism underpinning the change of activity manifests differently at the AGN dust torus than at the broad emission line region, although it is also possible that a correlation does exist but that the different physical locations of the emitting regions within the quasar impose a several year temporal delay between MIR and spectral variability, (e.g., Jun et al. 2015) that simple amplitude measures or the available data do not Figure 7. The distributions of the photometric selection criteria -the distance to the median Slepian wavelet variance curve, SWV 1 , the difference between initial and final Bayesian block (BB) states, and the maximum W 1 difference -against the spectroscopic variability for the CSQs identified in this paper (blue), known CLQs in the literature (red) where the multiepoch spectra are available, and the general population of quasars with CRTS lightcurves (gray contours). Open circles indicate low luminosity sources with M V > −23. A positive BB difference indicates a transition from a higher (brighter) to a lower (fainter) state. capture. Fig. 8 shows a trend between the MIR and optical photometric variability for CLQ/CSQs but not for the general quasar population. Since the former are known to show spectral variability, this also supports the idea that continuum variability is the stronger component overall as compared to just spectral line changes.
Yang et al. (2018) reported a "redder when brighter" correlation between MIR color and magnitude amplitudes, indicative of a stronger contribution from the AGN dust torus when the AGN turns on. Fig. 9 shows that this correlation holds only for the known CLQs and possibly for the lower luminosity CSQs but not for the higher luminosity objects. CSQs are also a redder population at MIR wavelengths than CLQs (a large proportion of which would not pass the MIR selection criteria used in this paper). The higher luminosity, more variable CSQs may thus already have a strong contribution in their MIR flux from the AGN dust torus that the relative change associated with the changing-state mechanism is not so significant an effect as with the lower luminosity sources.
Our selection criteria are somewhat broader by design than those employed in other works so it is also worthwhile considering CSQs in the context of other criteria used. MacLeod et al. (2019) and other CLQ searches have employed a simple variability amplitude criterion (typically |∆g| > 1) to select sources exhibiting strong photometric variability over any of the available time baselines probed by surveys, such as SDSS and Pan-STARRS 1. We used the difference between the initial and final states of the Bayesian block representation of a light curve rather than its amplitude (i.e, the difference between its maximum and minimum states) to identify sources showing significant photometric change but avoiding objects showing flaring activity such as reported in Graham et al. (2017) and Lawrence et al. (2016). Fig. 10 shows that almost all CSQs lie outside the 95th percentile contour in the SWV 1 -BB amplitude plane but that flaring quasars from Graham et al. (2017) are the majority source at higher amplitudes (compare with Fig. 2). Known CLQs still form a predominantly less strongly variable population.
It is possible, however, that the long baseline of typically 4000 days between the initial and final states in the observed frame could bias our selection toward quasars with longer timescale variability and miss shorter timescale variability potentially associated with CLQs. The Slepian wavelet variance measure, SWV 1 , has a characteristic timescale associated with it (see Sec. 5.2) and Fig. 11 shows the distribution of BB amplitude and difference in relation to this timescale for the 1.1 million initial sources in MQ. It is clear that with increasing variability, measured either through the amplitude or the difference, there are fewer objects with longer characteristic timescales and no indication of a selection bias towards them.
The spectral variability constraint (fractional change in Hβ / [O iii] > 30%) that is employed here also differs from the ad hoc visual criteria that have defined CLQs in other searches. As M19 note, the application of such a quantitative definition could lead to a different set of ambiguities and make comparison with other CLQ samples difficult. Fig. 3 shows that our spectral criterion is not met by about 45% of the known CLQs for which we have multiepoch spectra. CSQs are associated with large continuum luminosity changes and MIR variability over long timescales but not necessarily with a complete (visual) absence of Hβ flux in their lower activity (fainter) state. M19 uses spectral flux ratios between high (bright) and low (dim) states to determine more clearly how the Hβ line varies relative to the continuum (measured at restframe 3460Å). Fig. 12 shows that the distribution of Hβ and continuum (measured at restframe 3240Å) flux ratios for CSQs is consistent with the known CLQs for which multiepoch spectra are available and the values reported by M19. Note that flux ratios measured for sources where Hβ has largely disappeared in the low state will be a lower limit.
A further complication can arise due to the timing of, e.g., spectroscopic follow-up observations for candidates. Given the data-taking, it is not always possible to obtain a spectrum and sample the light curve, say, at the largest amplitude of variability (in flux, BB, or SWV 1 parameter space). This could present a bias for/against a particular class of object depending on when the spectra are obtained relative to the light curve. We do not claim completeness in our sample: for example, we note that the changing source SDSS J2232-0806 (Kynoch et al. 2019) is a photometric candidate in this work but lacks a second epoch spectrum in our data set. Similarly 57 of the 262 CLQ candidates (22%) from MacLeod et al. (2019) also appear in our sample of CSQ candidates with single epoch spectra.
More generally, we can consider our sources in the context of the extreme variability quasars (EVQs) of R17. In a comparison with a control sample of SDSS quasars matched in redshift and g magnitude to their EVQ sample, R17 find that EVQs have a larger variability amplitude (from the structure function) than control quasars at all timescales from days to years. This is equivalent to the SWV 1 selection criteria we have employed with the source SWV 1 greater than the median SWV 1 at the source magnitude across timescales. However, we compare CSQ quasars to the general population of MQ sources with CRTS light curves rather than a magnitude and redshift-matched sample.

Luminosity
Quasar variability is known to be anticorrelated with luminosity in that low luminosity quasars have a larger probability of showing large amplitude variation over multi-year timescales (e.g., Hook et al. 1994, R17). Using the 5100Å continuum luminosity as a proxy for the intrinsic AGN luminosity, MacLeod et al. (2016) have shown that CLQs seem to preferentially be low-luminosity AGN. Fig. 13 shows the distribution of the peak [O iii] luminosity from multiepoch spectra as a function of redshift for the CSQs reported here, the known CLQs, and values measured for 30,000 SDSS quasar at z < 0.95. Our CSQs expands the range of changing AGN to more luminous sources and higher redshifts. Furthermore, this type of variability is not dependent on the luminosity of the object.
R17 find that there is a trend of decreasing Eddington ratio with variability. To determine the population distribution, we have calculated the bolometric luminosity for the 4,244 SDSS quasars with possible Hβ coverage (z < 0.95) and at least 500 days between multiepoch spectra using: where the solar constants for V-band are M ,V = 4.83 and L ,V = 4.64 × 10 32 erg s −1 and b V is the bolometric correction. A comparison with the bolometric luminosities calculated for DR12 quasars from SDSS by Kozlowski (2017) (hereafter K17) with L bol,V gives a mean bolometric correction of log 10 b V = 1.46, which we use as a fiducial value hereafter. We estimate the black hole virial mass using: log M M = a + b log L 5100 10 44 erg · s −1 + 2 log FWHM(Hβ) 1000 km · s −1 with a = 6.91 and b = 0.533 for Hβ (Ho & Kim 2015). We find that estimates for the same source from multiepoch spectra are typically consistent within 0.4 dex.
The Eddington ratio is defined as η Edd = L bol /L Edd where L Edd is the Eddington luminosity. Fig. 14 shows the distribution of the Eddington ratio for this sample as a function of variability amplitude as measured from the Bayesian block fits to the respective time series (note that here amplitude = 0.5 × (max value -min value)). It also shows this distribution for SDSS quasars with potential Hβ coverage (z < 0.95) with Eddington ratio estimates derived as above and also from K17 as a consistency check. The general population does not show any strong indication of decreasing Eddington ratio with increasing variability, as had been previously reported (e.g., R17). Instead, we find that low-amplitude (amp < 0.2) low-luminosity sources have a fractionally higher Eddington ratio (η Edd = 0.2) relative to an otherwise flat relationship between Eddington ratio and variability amplitude. We note that R17 used Eddington ratio values from Shen et al. (2011) whereas K17 derived their values from their MgII and CIV-based black hole virial mass  estimates. Our values are based on our own Hβ virial mass estimates and a bolometric luminosity correction derived from K17 data. Although virial mass estimators based on different lines can systematically disagree, this is not sufficient to explain this difference.
The discrepancy may lie in the R17 control sample having no upper redshift bound (the majority of their comparison sample is at z > 1) whereas both our CSQs and the control sample data is constrained to z < 0.95. There is a correlation between bolometric luminosity (Eddington ratio) and redshift in the R17 sample so that low variability amplitude bins in R17 will be biased toward higher Eddington ratio as a large fraction of objects will have z > 1. Higher amplitude bins have comparatively more sources with z < 1 and so are more consistent with our results. However, we agree that CLQ/CSQs do show the reported anticorrelation between Eddington ratio and amplitude of variability and that it is stronger for low luminosity (M V < −23) sources. This is consistent with attributing the changes seen to accretion physics occurring preferentially in lower activity systems but not necessarily just in low luminosity sources. We also note, though, that although the CSQ (this paper), CLQ (M19), and EVQ (R17) samples all have log(η Edd ) = −2 to −0.5, this does not imply the same physical mechanisms are necessarily involved across these samples.

Timescales
One of the selection criteria in Sec. 3 was an excess of variability relative to a median level for a magnitude range as measured by Slepian wavelet variance (SWV). We employed a cumulative measure looking for an overall significant signal rather than one at any specific timescale. However, CSQs  Fig. 7. The contour lines indicate the 50th, 70th, and 90th percentiles respectively. Figure 14. The distribution of Eddington ratio as a function of variability amplitude for SDSS quasars with z < 0.95 and multiple spectra separated by more than 500 days (contours and grey points; contours show 68th and 95th percentiles). The black points show the median value for this data set in each bin of variability and the green show the same quantity for all SDSS DR12 quasars from Kozlowski (2017) with z < 0.95. The values of known CLQs (red) and the CSQs (blue) reported here are also shown. The solid cyan line indicates the Eddington ratio trend from R17 Fig. 10 and the dashed cyan line shows a linear fit to the CLQ/CSQ sources. display a particular pattern of activity over the period covered by their light curve and therefore the timescale which contributes most to the variability of the source should be associated with (or even characteristic of) the physical mechanism driving the change. We have determined the restframe timescale for each CSQ in our data (and CLQ in the literature) at which the SWV of the source has its largest value relative to the median value (see Fig. 15) as well as the distribution of such timescales for 137,000 quasars with z < 1.1. The distribution of peak time values for the CSQ/CLQ sample is significantly different than the population distribution (>> 5σ according to the Anderson-Darling test) and so the timescales are indicative of process(es) associated with the observed variability.
As given in Stern et al. (2018), the relevant disk timescales for a black hole of mass M BH at R ∼ 150r g can be parameterized as: (1) where α is the disk viscosity parameter, h/R is the disk aspect ratio, R is the disk radius, and r g = GM BH /c 2 is the gravitational radius. The thermal timescale t th corresponds to the timescale on which the disk heats or cools with cooling and heating fronts crossing the disk of timescales of t front . The viscous disk timescale, t v , gives the characteristic timescale of mass flow.
(1)-(3), timescales associated with processes in AGN disks are expected to be functions of the disk aspect ratio (h/R), disk viscosity parameter (α), black hole mass M BH and disk radius (R(r g )). In Fig. 16 we show the loci of timescales from eqns. (2) and (3) as a function of both α and (h/R). Solid lines correspond to t front = 1yr and dash-dot lines correspond to t v = 1yr located at 50r g (red) and 150r g (black) respectively in a disk around a M BH = 10 8 M . We can read Fig. 16 as follows: if a CSQ is observed to change state on a one year timescale and we model the associated spectral change with the propagation of a front from the ISCO to 150r g , then disk properties (h/R, α) must live on the black solid line (e.g., both h/R ∼ 0.1, α ∼ 0.3 and h/R ∼ 0.3, α ∼ 0.1 are possible solutions for the model. Likewise if a CSQ varying on a timescale of one year is modeled in terms of a viscous change at 50r g , then the disk parameters (h/R, α) must live on the red dash-dot line (e.g., h/R ∼ 0.5, α ∼ 0.04 is a possible solution).
For a given value of R, each quasar will define a curve in the (h/R, α) plane and Fig. 17 shows the distribution of these for t front and t v for the CLQ/CSQ sample. There are clearly different preferred regions of the parameter space depending on whether the timescales are interpreted as front crossing or viscous. For example, King et al. (2007) argue that observations favor a typical range of α ∼ 0.1−0.4 and this would suggest disk scale ratios of ∼ 0.03−0.1 at R = 50r g for front crossing timescales but ratios of ∼ 0.2 − 0.3 for viscous timescales. Although there are degeneracies, broad constraints can be placed on viable disk geometries: a viscous timescale process favors a thicker disk and less change in disk thickness with increasing radius whereas a front crossing process can support a thinner disk but also one that expands more in height with increasing radius.

Physical mechanisms
Cooling fronts have been proposed as the mechanism for CLQ/CSQs, either as a result of a sudden change in torque applied by the magnetic field at the innermost stable circular orbit (ISCO) (Ross et al. 2018), or a drop in mass accretion rate causing an advection-dominated accretion flow (ADAF) in the inner disk and spectral state transition by disc evaporation (Noda & Done 2018). Sniegowska & Czerny (2019) also propose that for sources operating at a few percent of the Eddington limit, there is a radiation pressure instability in a narrow zone between the outer cold gas-dominated disk and an inner hot ADAF flow which can lead to outbursts producing changing look behavior. Noda & Done (2018) suggest that CLQ/CSQ sources can be placed in one of three groups, depending on which particular aspect of the process they exhibit: (1) a factor of two to four decrease in luminosity associated with disc evaporation/condensation; (2) large mass accretion rate change due to thermal front propagation; and (3) a variability amplitude of more than ten indicative of both phenomena. Most objects should be in either the first or the third group whereas the second group will contain sources showing large variability but not showing any significant spectral changes (although given the results of Sec. 4, MIR variability would be expected to be shown in addition to optical). Fig. 18 shows that most CSQs/CLQs are associated We assumed a fiducial timescale of one year in each case, at disk radii 50r g (red) and 150r g (black) around a black hole of mass M BH = 10 8 M . If timescales are shorter, the curves will shift in parallel towards the top right of the Figure. If the timescales are longer, the curves will shift in parallel towards the bottom left of the Figure. As the disk approaches a spherical configuration (H/r → 1), t front → t v as we see from eqns. (2) and (3)  with a change of Eddington ratio (accretion rate) of between 1% and 10% L edd , consistent with the observational predictions from Noda & Done (2018) and placing these sources in their third group. We note as well (see Fig. 19 and Table 3) that the magnitude of the change in η Edd shows a trend with (median) luminosity and also with the amplitude of variability: amp ∝ log(∆η Edd ) and ∆η Edd ∝ L/log(η Edd ), where the inverse relationship comes from Sec. 5.1. In other words, more extreme variability is associated with larger changes of η Edd in higher luminous systems but also with lower actual η Edd or, conversely, lower luminosity systems with higher η Edd but only able to support a smaller change in η Edd . If the magnitude of the change in η Edd correlates with either a change in torque at the ISCO or a change in mass accretion rate then larger systems show stronger fluctuations. This suggests that disk instabilities, e.g., magneto-hydrodynamical, may be a more likely cause than local perturbative events in the disk, e.g., an embedded supernova, since the latter should not scale with the size of the accreting system. Such instabilities may be driven by the larger environment: Charlton et al. (2019) reported that four CLQs are associated with galaxy mergers and Kim et al. (2018) have proposed that changing look activity in Mrk 1018 is due to a recoiling SMBH perturbing the accretion flow on a 29-year period. Alternatively, both disk instabilities and local perturbative events may be present but with a bias for the latter in lower luminosity systems. We expect that low luminosity AGN (LLAGN, should be more heterogeneous in origin than the high luminosity AGN population (filled circles in   . 7-20). This is because either the intrinsic SMBH mass is lower than for the higher luminosity population even if η Edd is comparable, or because η Edd is intrinsically lower than for the higher luminosity population, or some combination of these. So we should anticipate a lack of correlation between luminosity among the LLAGN and η Edd . This is confirmed in Fig. 20 where the LLAGN (open circles) form a scatterplot. Conversely, for the high luminosity AGN, there is an apparent correlation with η Edd above ∼ 10 44 erg s −1 in Fig. 20. While we expect LLAGN are more heterogeneous than the high luminosity AGN, a change in the accretion rate in LLAGN (effectively a change in η Edd ) should correlate with luminosity. This is indeed apparent in the left panel of Fig. 19. The LLAGN variability amplitude anti-correlates with the accretion rate (Fig. 14) which suggests that it is harder to significantly change the accretion rate in LLAGN.

CONCLUSIONS
We have identified 73 quasars exhibiting strong coherent optical and MIR photometric variability with significant contemporaneous spectroscopic variability that are comparable to the ∼60 existing changing-look quasars reported in the literature. Our sample, however, forms a higher luminosity (and higher redshift) counterpart to the known CLQs showing that this phenomenon is not restricted to low luminosity systems. The characterizing preference is rather for systems with low Eddington ratios and with the amplitude of the associated variability correlated with a change of Eddington    Figure 20. Eddington ratio as a function of luminosity. The same color scheme is used as in Fig. 7. ratio. Characteristic timescales of the photometric variability suggest that it most closely matches the timescale associated with a cooling/heating front propagating through the disk as has been proposed for individual sources Ross et al. 2018;Noda & Done 2018). The lack of large variability in smaller systems may also indicate disk instabilities associated with magnetic phenomena as the more likely physical cause for the fronts, particularly in larger systems.
The next generation of sky surveys, such as the Zwicky Transient Facility (ZTF; Bellm et al. 2019;Graham et al. 2019) and LSST, will effectively monitor all AGN in the sky ever few nights. Generative models of quasar variability can be learnt from archival data and predicted behavior compared to that observed with unexpected changes identified far more quickly -in weeks or months rather than years -than waiting for a significant magnitude change to be detected. In this way, these changes of the state of the accretion disk can be tracked with follow-up resources as they happen rather than serendipitously or after the fact. This will allow us to test more easily the theoretical explanations for this phenomena.
Future work will employ machine learning to identify potential further sources -the combination of known CLQs and the CSQs reported here ensures that there is now adequate coverage of the parameter space and a suitable training set can be defined. Further characterization of the sources will also aid this activity. We have also undertaken a program to find CSQs with z > 0.95, i.e, where Hβ does not fall into the optical spectral range. Candidates sharing the same photometric variability as their lower redshift counterparts have been identified and optical and near-IR spectra are being obtained, the latter to capture Hβ. Although multiepoch near-IR spectra are unlikely to exist, we will explore the possible correlations between Mg ii and Hβ variability for these sources relative to a more expected lack of correlation in the general population. We are interested as well in those objects which meet the photometric selection criteria but not the spectroscopic to understand if we are probing the same population but missing the spectral variability due to delayed discovery and followup or whether such objects are associated with a different phenomenon. Similar efforts are underway for sources with mid-IR variability but no associated optical change. Table A1. CLQs reported in the literature with CRTS coverage and associated features including the median CRTS magnitude (V m ), the optical amplitude (Amp), absolute change in W 1, Bayesian block change (∆BB), and Slepian wavelet variance measure (SWV 1 ). Sources marked with an asterisk have multiepoch spectra in the public domain. SMBH virial mass estimates are calculated as described in Sec. 5.1 except for sources where no spectra are available.

Name
Transition