# Extreme Variability Quasars in Their Various States. II: Spectral Variation
Revealed with Multi-epoch Spectra

Wenke Ren CAS Key Laboratory for Research in Galaxies and Cosmology,
Department of Astronomy, University of Science and Technology of China, Hefei,
Anhui 230026, China  School of Astronomy and Space Science, University of
Science and Technology of China, Hefei 230026, China  Junxian Wang CAS Key
Laboratory for Research in Galaxies and Cosmology, Department of Astronomy,
University of Science and Technology of China, Hefei, Anhui 230026, China
School of Astronomy and Space Science, University of Science and Technology of
China, Hefei 230026, China  Zhenyi Cai CAS Key Laboratory for Research in
Galaxies and Cosmology, Department of Astronomy, University of Science and
Technology of China, Hefei, Anhui 230026, China  School of Astronomy and Space
Science, University of Science and Technology of China, Hefei 230026, China
Xufan Hu  CAS Key Laboratory for Research in Galaxies and Cosmology,
Department of Astronomy, University of Science and Technology of China, Hefei,
Anhui 230026, China  School of Astronomy and Space Science, University of
Science and Technology of China, Hefei 230026, China

(Received September 10, 2023; Revised November 26, 2023; Accepted December 14,
2023)

###### Abstract

We previously built a sample of 14,012 extremely variable quasars (EVQs) based
on SDSS and Pan-STARRS1 photometric observations. In this work we present the
spectral fitting to their SDSS spectra, and study the spectral variation in
1,259 EVQs with multi-epoch SDSS spectra (after prudently excluding spectra
with potentially unreliable spectroscopic photometry). We find clear “bluer-
when-brighter” trend in EVQs, consistent with previous findings of normal
quasars and AGNs. We detect significant intrinsic Baldwin effect (iBeff, i.e.,
smaller line EW at higher continuum flux in individual AGNs) in the broad Mg
II and C IV lines of EVQs. Meanwhile, no systematical iBeff is found for the
broad H$\beta$ line, which could be attributed to strong host contamination at
longer wavelengths. Remarkably, by comparing the iBeff slope of EVQs with
archived changing-look quasars (CLQs), we show that the CLQs identified in
literature are mostly likely a biased (due to its definition) sub-population
of EVQs, rather than a distinct population of quasars. We also found no
significant broad line breathing of either H$\beta$, Mg II or C IV, suggesting
the broad line breathing in quasars may disappear at longer timescales ($\sim$
3000 days).

Quasars (1319), Time domain astronomy (2109), Active galactic nuclei (16),
Black hole physics (159)

††journal: ApJ††facilities: Sloan, PS1, PO:1.2m††software: PyQSOFit (Guo et
al., 2018), astropy (Astropy Collaboration et al., 2013; Price-Whelan et al.,
2018; Astropy Collaboration et al., 2022), pandas (pandas development team,
2020; Wes McKinney, 2010), scipy (Virtanen et al., 2020), numpy (Harris et
al., 2020),

##  1 Introduction

One of the hallmarks of quasars is their universal variability across all
wavelengths. Quasar light curves typically exhibit variations with an
amplitude of $\sim$0.2 mag in the rest-frame UV/optical band, occurring on
timescales of months to years (e.g., Vanden Berk et al., 2004; Sesar et al.,
2007; MacLeod et al., 2010). These variations are often described by a damped
random walk (DRW) model first proposed by Kelly et al. (2009). Reverberation
mapping (RM) studies have revealed strong correlations between the variation
of the continuum and the broad emission lines (BELs) (e.g., Grier et al.,
2017, 2019; Homayouni et al., 2020), suggesting that the optical-UV
variability of quasars is driven by the innermost region of the accretion disk
(e.g., Ross et al., 2018). However, due to the stochastic and complex nature
of quasar variability (e.g. Cai et al., 2018; Sun et al., 2020), the physical
processes underlying these variations are still a topic of debate.

In the cutting-edge field of quasar variability research, the phenomenon of
changing-look (CL) is particularly intriguing. It has long been noted that a
few AGNs exhibit drastic changes in their broad emission lines, with their
spectral types transitioning between Type 1 and Type 2 (e.g., Cohen et al.,
1986; Goodrich, 1995). Initially, due to limited data sizes and time spans,
the discovered CLQs are very rare and usually transient. The observed changes
in broad lines were often attributed to dust extinction (Risaliti et al.,
2007), tidal disruption events, or supernovae (Aretxaga et al., 1999).
However, with the accumulation of data, an increasing number of CLQs
($\sim$100) have been identified (e.g., Ruan et al., 2016; MacLeod et al.,
2019; Green et al., 2022). Most CL phenomenon are accompanied by extreme
continuum variations, and their new spectral states can keep for months to
years. These observational results indicate that the CL phenomenon should be
driven by intrinsic but violent, accretion dominated activity. What are the
underlying mechanisms for such violent process? Are they due to possible
accretion transition in analogy to X-ray binaries (e.g. Ruan et al., 2019;
Yang et al., 2023), or just more violent version of normal variability? Much
larger samples of sources with violent variability are desired to distinguish
such possibilities.

Since the identification of CLQs relies on repeated spectral observations
showing emerging/disappearing of broad emission lines, which are relatively
rare, some related studies initially focus on quasars with extremely strong
continuum variability (MacLeod et al., 2016; Rumbaugh et al., 2018; Sheng et
al., 2020). Compared with CLQs, much larger samples of such quasars could be
easily built solely based on broadband photometric variability, and such
samples could provide essential clues to the physical nature of violent
variability and its relation with changing-look and normal (less extreme)
variability.

Through systematic research on extremely variable quasars (EVQs), Rumbaugh et
al. (2018) found that the equivalent width (EW) of the broad emission lines of
EVQs are larger than that of normal quasars. Note Kang et al. (2021) further
discovered a general correlation between the variability and emission line EW
in large sample of quasars. It is clear that quasar variability has an
underlying influence on the broad-line regions (BLRs) and thus is worth
investing more efforts. Furthermore, studies with small samples of AGNs have
illustrated that emission lines change along with continuum variations, i.e.,
the EW of broad lines decreases as the continuum increases (intrinsic Baldwin
effect, iBeff, e.g. Goad et al., 2004; Rakić et al., 2017); the full width at
half maximum (FWHM) of different lines also exhibits varying correlations with
continuum luminosity (e.g. Wang et al., 2020). However, restricted by the
signal-to-noise ratio (SNR) of available survey spectra, such changes are
often not eminent enough to distinguish from noise when extending to quasars,
as most of which exhibit low intrinsic variability. In this sense, studying
the emission line variability in large samples of EVQs could also shed new
light on the BLR physics.

In our previous work Ren et al. (2022) (hereafter, Paper I, ), we constructed
a sample of EVQs consisting of 14,012 sources using photometric data from the
Sloan Digital Sky Survey (SDSS) and Pan-STARRS (PS1). We obtained the
composite spectra of the EVQs in various states (according to the relative
brightness of each EVQ during SDSS spectroscopic observation compared with the
mean brightness from multi epoch photometric observations) and of control
samples with matched redshift, luminosity and black hole mass. We verified
that EVQs tend to have stronger emission lines in H$\beta$, [O III], Mg II,
and C IV. The EW changes between different states shows a clear iBeff in Mg II
and C IV, but not in H$\beta$. Furthermore, we found that EVQs have stronger
broad line wings compared to normal quasars, suggesting that EVQs could launch
more gas from the inner disk into the very broad line region (VBLR).

In this work, we focus on the EVQs with multi-epoch SDSS spectra to explore
the spectral variations in individual EVQs. Based on the EVQ sample we
constructed in Paper I, we select a sub-sample with repeated spectral
observations in SDSS. The structure of this paper is as follows. In §2, we
describe the sample selection and methods of spectral measurement. The results
are presented in §3. We then discuss our results in §4 and provide a summary
in §5. Throughout this paper, we adopt a flat $\rm{\Lambda}$CDM cosmology with
$\Omega_{\Lambda}=0.7$, $\Omega_{m}=0.3$, and
$H_{0}=70~{}\rm{km~{}s^{-1}~{}Mpc^{-1}}$.

##  2 The Data and Sample

###  2.1 The EVQ Sample

We start from the EVQ catalog constructed by Paper I, which is selected from
the SDSS data release 14 (DR14) quasar catalog (DR14Q, Pâris et al., 2018).
Here we provide a brief overview of the selection criteria and spectra fitting
procedure. Readers are advised to refer to Paper I for detailed description.
Based on all photometric observations ($g$ and $r$) from SDSS and PS1 database
for each quasar in DR14Q, we amend the minor differences in filter
transmission between two surveys, reject epochs with potentially unreliable
photometric measurements, and build clean light curves for every source.
Sources with robust detection of extreme variability in both $g$ and $r$ light
curves (with $\Delta g_{max}>1$ and $\Delta r_{max}>0.8$) are selected as
EVQs. In total, we select 14,012 EVQs with 20,069 archived SDSS spectra. Each
spectrum is processed using PyQSOFit (Guo et al., 2018) following the
procedure described in Paper I. The catalog of the EVQs and the spectral
fitting results are now presented in Table 1 of Appendix A.

Among the 14,012 EVQs, 1354 sources have repeated spectral observations in
SDSS archive with a continuum SNR $>3$. In Appendix §B, we develop a method to
identify spectra with potentially unreliable spectroscopic photometry (such as
due to fiber-drop) by comparing the observed spectroscopic photometry with the
expected value derived through modeling the photometric light curves with a
damped random walk (DRW) process. After excluding such possibly dubious
spectra, 1259 multi-spectra EVQs are remained for further analyses.

###  2.2 Spectral measurements

The spectra are fitted following Shen et al. (2011) using PyQSOFit (Guo et
al., 2018). For each spectrum, we fit with a power-law
($f_{\lambda}\propto\lambda^{\alpha}$) and a broadened Fe II template for
continuum within tens of separated emission line-free windows and give the
measurement of continuum and Fe II properties. Then, we cut the spectrum
around three luminous broad emission lines: H$\beta$, Mg II, and C IV if they
are within the spectral coverage and then fit the lines locally. The continuum
is re-fitted using the windows on both sides of the emission lines. After
detracting the continuum, we use a few groups of Gaussian to fit the
residuals: 4 Gaussians (2 cores + 2 wings) for [O III] doublet, 1 for narrow
H$\beta$ and three for broad H$\beta$ in H$\beta$ region; 1 Gaussian for
narrow and 3 for broad Mg II region; and 3 Gaussians for only broad C IV
component in C IV region. The detailed procedure and parameter settings can be
found in Paper I.

For each emission line component, we calculate the fitted peak wavelength
(PEAK), integrated flux (FLUX), equivalent width (EW), line dispersion (second
moment, SIGMA), and full width at half maximum (FWHM). Since the broad
component of H$\beta$, Mg II, and C IV (consisting of triple Gaussians) may
have asymmetric shapes indicating complex structures of BLRs, we further
measure the following properties. The bisectional center (BISECT, Sun et al.,
2018), where the wavelength splits the line into two equal areas in flux, has
been used in Paper I to qualitatively reflect the wing structures in a broad
line. To further parameterize the shape and asymmetry of the broad component,
we adopt the: i) full width at half maximum (FWHM), ii) full width at quarter
maximum (FWQM), and iii) full width at 10% maximum (FW10M) to denote the BLRs
from different distance to the SMBH. We also measure the corresponding shifts
of the centers of the above three widths, Z50, Z25, and Z10 to show the
intrinsic gravitational redshift or inflow/outflow of BLRs (cf. Fig. 6 in
Rakić 2022).

We employ a Monte Carlo approach to assess the statistical uncertainties of
the parameters we measured. This is done through generating 50 mock spectra by
adding Gaussian noise to each original spectrum using the reported flux
density errors, fitting the mock spectra with the same routines, and deriving
the scatter of each parameter.

We adopt the bolometric corrections (BCs) and the calibration parameters in
Shen et al. (2011) to estimate the bolometric luminosity $L_{bol}$ and the BH
mass $M_{BH}$ of a single-epoch spectrum. We take fiducial $L_{bol}$ and
$M_{BH}$ calculated from $L_{5100}$ and H$\beta$ for sources at $0.08\leq
z<0.7$, from $L_{3000}$ and Mg II at $0.7\leq z<1.9$, and from $L_{3000}$ and
Mg II at $1.9\leq z<4.0$ respectively. Our following results will be based on
such three sub-samples. The Eddington ratio is also given in our catalog. We
adopt the mean value from multiple-epoch measurements for each target and plot
the general physical property distribution of multi-spectra EVQs in Fig. 1.

Figure 1: The redshift, $L_{bol}$, $M_{BH}$, and $\eta_{Edd}$ distribution of
multi-spectra EVQs. The mean values from multiple spectra are adopted for each
EVQ. We split the sample into three groups according to redshift. The source
number of each group is noted in the parenthesis in the legend. As a
comparison, we also plot the re-normalized distribution of SDSS DR14Q (Rakić,
2022), divided by the same redshift criterion and drawn with corresponding
colored shades.

##  3 Results

Figure 2: The correlation between the variability of spectral properties and
that of continuum luminosity. The variation of continuum power-law slope
($\alpha$), broad line EW and FWHM are shown in rows sequentially. The
variability for each parameter is derived through subtracting the mean value
of multi-epoch spectral measurements for the same EVQ from a single-epoch
measurement. The Spearman correlation coefficients and the p-values (derived
through bootstrapping the sample, hereafter the same, see the 2nd parapha of
§3 for details) are shown in the upper right of each plot. For correlations
with p-value$<$0.01, we use a black solid line to show the orthogonal distance
regression result and present the derived regression slope $k$ with 1$\sigma$
uncertainties derived through bootstrapping (hereafter the same).

In this section, we explore how the spectral properties of EVQs vary along
with the continuum flux, in three distinct samples split by redshift to
guarantee the emission line measurements. In total, we have 329 sources with
920 spectra for the H$\beta$ sample, 676 sources with 2206 spectra for the Mg
II sample, and 254 sources with 799 spectra for the C IV sample. About $2/3$
of the objects have only two spectra, and there are 15 sources with more than
20 observations which are all from the SDSS-RM program (Shen et al., 2014).

To present sources with different luminosities together and explore their
average variability, our results are primarily expressed in terms of
variation. To fully leverage each observation, we derive the deviation of each
single-epoch measurement from the mean value for each EVQ. In this case, one
single-epoch spectrum will contribute one data point to our following plots.
We note that the various data points from the same EVQ are thus not completely
independent to each other (since calculating the mean value consumes one
degree of freedom for each EVQ) and the traditional p-value of the Spearman
correlation evaluated from Student’s t-test is no longer valid. Therefore, we
estimate the p-value by bootstrapping the quasar sample instead. Note since we
run bootstrap 100,000 times, thus we are unable to give exact p-value smaller
than 1.e-5.

###  3.1 Color variation

The so-called “bluer-when-brighter” (BWB) pattern has been widely seen in
quasars and AGNs (e.g. Trevese & Vagnetti, 2002; Sun et al., 2014; Guo & Gu,
2016). In Paper I, the composite spectra of EVQs indirectly support such a
pattern. We quantitatively test this conclusion with the multi-spectra EVQs by
exploring the spectral slope variability in individual sources.

The first row in Fig. 2 shows the changes of spectral slope versus the changes
of log-luminosity. The anti-correlation between $\Delta\alpha$ and
$\Delta\log{L}$ is found in all three samples, denoting a consistent BWB
pattern up to $z\sim 4$. We use the orthogonal distance regression (Boggs &
Rogers, 1989) method provided by scipy.odr package to fit the data with the
data errors considered (hereafter the same, unless otherwise stated). The
regression slope $k$ for the three samples are $-1.6\pm 0.3$, $-0.90\pm 0.12$,
and $-1.4\pm 0.3$ ((with the 1$\sigma$ statistical uncertainties derived from
bootstrapping the samples, hereafter the same), respectively, roughly
consistent with the results from PG QSOs (e.g. Pu et al., 2006).

The slightly steeper slope from the H$\beta$ sample could be due to the
contamination of the host galaxy as we do not decompose the host galaxy
component while fitting the continuum. As a result, for those low-z sources
for which the host contamination could be more prominent, the
$\Delta\log{L_{5100}}$ will be underestimated, and the $\alpha$ at dim states
will be overestimated. Besides, the relative flux of the host galaxy component
may also vary depending on the observation conditions (Zhang, 2013; Hu et al.,
2015). Nevertheless, given their high luminosity as quasars (see Fig. 1) and
the strong correlation at higher redshift, the host alone should not be fully
responsible for the observed spectral changes. Indeed, it has been well
demonstrated that the host contamination can not dominate the observed color
variation in SDSS quasars, as host contamination can not reproduce the
observed timescale-dependence of the color variation (Sun et al., 2014; Cai et
al., 2016; Zhu et al., 2016).

###  3.2 Variation of line EW and FWHM

In the second row of Fig. 2, we explore the iBeff of the three most luminous
broad lines. Strong anti-correlations between the variation of line EW and
that of the continuum luminosity are found in Mg II and C IV while only
marginal anti-correlation in H$\beta$. The averaged iBeff slope for Mg II and
C IV are comparable but slightly steeper than the results revealed with the
stacked spectra in Paper I and those from literature (Pogge & Peterson, 1992;
Kong et al., 2006; Homan et al., 2020; Zajaček et al., 2020). The steeper
iBeff of Mg II and C IV could likely be attributed to the large variability of
our sample, as Homan et al. (2020) also found a steeper slope in their more
variable sample.

We note that the missing anti-correlation between $EW_{H\beta}$ and $L_{5100}$
does not means the absence of iBeff in H$\beta$, as the host component has not
been properly considered in our spectral fitting. The blend of the host galaxy
in the spectra causes significant underestimation of $EW_{H\beta}$,
particularly in dim states, thereby weakening the underlying anti-
correlations. Subtracting the host component from such low SNR individual
spectra, which is beyond the scope of this work, would be highly unreliable
and would introduce great uncertainties to spectral measurements (e.g. Wu &
Shen, 2022).

Figure 3: The broad H$\beta$ iBeff is visible in a sub-sample of more luminous
EVQs, suggesting the absence of broad H$\beta$ iBeff shown in the lower left
panel of Fig. 2 is due to host contamination. The labels are similar to those
in Fig. 2.

To alleviate the host contamination, we build a luminous sub-sample from our
EVQs which contains about 1/3 of the total sample with
$\log\overline{L}_{5100}>44.5$. The broad H$\beta$ iBeff of this sub-sample is
shown in Fig. 3. A robust negative correlation is observed in this sub-sample,
indicating that the disappearance of the iBeff of H$\beta$ line in the full
sample is likely due to the host contamination. We would further discuss the
effects of host contamination in §4.1, where the connection between the EVQs
and CLQs are investigated.

We also explore the so-called ”breathing” of the broad lines, e.g., how the
line width changes with the luminosity (see the third row of Fig. 2). Only
weak and marginal anti-breathing (positive correlation) is found in the C IV
sample, and no clear sign of breathing is observed in the Mg II sample. These
results are generally consistent to those reported in Wang et al. (2020).
However, we find that the H$\beta$ sample also shows no breathing according to
Spearman test, contradicting previous findings (e.g., Denney et al., 2008;
Park et al., 2011; Wang et al., 2020). This paradox will be discussed in §4.2.

###  3.3 Line profile and asymmetry

Figure 4: The variation of broad line profile and asymmetry (see text for
description) versus that of continuum luminosity. The labels are similar to
those in Fig 2

In Paper I, we found that the EVQs have slightly stronger broad wings compared
with control samples. Here we investigate whether the broad line profiles vary
with continuum flux in EVQs with multi-epoch spectra. We use the line width
ratio (FWQM/FWHM and FW10M/FWHM) to denote the line profile and the shift of
line center (Z50, Z25, and Z10) to quantify the line asymmetry. The changes of
above parameters with the continuum luminosity are shown in Fig. 4.

According to the Spearman correlation coefficient and its p-value, we do not
find any strong dependency between line profiles and continuum flux. This
suggests the BLR structure does not systematically change with continuum
variation. Such a result also indicates that the EVQs are not an unusual
population with exceptionally large variability.

##  4 Discussion

###  4.1 EVQs and CLQs

In §3.2 we found the H$\beta$ EW of EVQs does not vary with continuum flux
(but with considerable scatter), i.e., showing no clear iBeff on average. In
such cases, there will naturally be some EVQs that have lower H$\beta$ EW in
their noisy dim-state spectra which may yield the so-called “vanishment” of
broad Balmer lines (i.e., CLQs). In this section, through this entry point, we
will delve deeper into whether CLQ is a distinct population or simply a biased
sub-population of EVQ.

####  4.1.1 Archived CLQs

The study of CLQs has been heating up in recent years, decades after the
discovery of Changing-look AGNs (e.g., Tholine & Osterbrock, 1976; Cohen et
al., 1986). So far, nearly a hundred of CLQs have been claimed, but a unique
classification is still lacking. The term “Changing-look” is an iconic
description (highly qualitative but not quantitative) firstly used to denote
the appearance or disappearance of broad emission lines. Since the CL
phenomenon is relatively rare at the early stage, it was feasible to use AGN
types or merely by visual checks to identify CL-AGNs. Such subjective
identification method is inherited to the research of CLQs.

Since the CL phenomenon is usually accompanied by dramatic continuum changes,
one developed approach to efficiently search for CLQs is to find candidates
from quasars with strong continuum variation, and then pick CLQs through
visual check (e.g. MacLeod et al., 2016). Obviously, the process of visual
check can be easily influenced by individual subjective judgements.

With the accumulation of multi-epoch spectra, researchers have explored
another method for hunting CLQs. For the spectra of CLQs in dim states usually
have no broad lines, they were often classified as galaxies by the initial
pipelines. Therefore, one can search for CLQs by directly examining sources
with repeated spectra but with different classifications (e.g. Ruan et al.,
2016; Yang et al., 2018). Such semi-quantitatively search for CLQs is highly
dependent to the classification pipeline, and in many cases still dependent on
visual inspection (e.g., 19 confirmed CLQs out of 10,204 candidates by visual
check in Yang et al., 2018).

In order to further standardize the selection of CLQs, a few recent works
start to focus on the change of emission line flux itself. MacLeod et al.
(2019) proposed that a quasar could be deemed as a CLQ if the SNR (signal-to-
noise ratio) of the broad line flux change between two epochs is greater than
3. Meanwhile Yang et al. (2018) and Sheng et al. (2020) directly use the SNR
of broad line to decide whether a broad component is apparent or not. The most
recent work by Green et al. (2022) combined the quantitative criteria and
visual inspection. These efforts indeed make the search for CLQs more
objective and consistent. However, the criterion above highly rely on the SNR
of the spectra, and they are unstable when applying to spectra observed in
different conditions and could be sensitive to the quasar luminosity. For
example, as MacLeod et al. (2019) have noted that their methods may select
quasars with minor variation in H$\beta$ line flux when their spectra have
very high SNR. Furthermore, applying quantitative criteria on the broad line
flux is not straightforward, as the broad line flux is expected to vary with
continuum in all quasars and AGNs (the cornerstone of reverberation mapping
studies).

To explore whether there are intrinsic differences between currently archived
CLQs and EVQs, we collected all archived H$\beta$ CLQs and refitted their
spectra using the same method we applied to our EVQs. Note in literature there
are only 3 MgII CLQs and 3 C IV CLQs claimed (Guo et al., 2019, 2020; Ross et
al., 2020, see also §4.1.3), thus here we focus on Hb select CLQs only. We
gathered a total of 79 securely-identified CLQs from LaMassa et al. (2015);
Ruan et al. (2016); MacLeod et al. (2016); Yang et al. (2018); MacLeod et al.
(2019); Sheng et al. (2020); Green et al. (2022); López-Navas et al. (2022).
To ensure the consistency of analyses, we only retained CLQs with repeated
SDSS spectra in both dim and bright states. As a result, 36 sources are
preserved, all with $z<0.65$. Among those CLQs, 14 are included in our EVQ
sample. The other 22 CLQs do not satisfy our strict selection criteria of EVQs
(i.e., $\Delta g_{max}$ $>$ 1 mag and $\Delta r_{max}$ $>$ 0.8 mag, in SDSS &
PS1 photometric light curves well cleaned against potential photometric
defects, see Ren et al. 2022).

Indeed, our independent spectral fitting reveal broad H$\beta$ lines in the
dim states of 33 CLQs, indicating in most archived CLQs the broad H$\beta$
indeed does not completely disappear in the dim state. In the remaining three
cases, two of them have their H$\beta$ lines in the dim states completely
submerged by the host galaxy, leaving an obvious absorption feature. The dim-
state spectrum of the third CLQ is flawed without enough wavelength coverage
around H$\beta$. These 3 dim-state spectra are therefore excluded from further
analysis below.

####  4.1.2 Comparing the iBeff of EVQs with CLQs

Figure 5: The orthogonal distance regression correlation between the changes
of continuum luminosity and the changes of line EW. The variation of continuum
luminosity and line EW are calculated by subtracting the mean value of each
object from the measurement of each epoch, which is identical to that in Fig.
2. Upper: the iBeff of H$\beta$ of our EVQs and CLQs collected from
literature. Lower: the Mg II iBeff for the sources plotted in the upper panel
but with Mg II line measurements. The fitting results and the Spearman
correlation coefficient $r$ (with p-values) of each sample is listed. For
correlations with p-value$<$0.01, we use a black solid line to show the
orthogonal distance regression fitting result.  Figure 6: The distribution of
H$\beta$ EW and iBeff slope versus the mean continuum luminosity. The iBeff
slope is given by the $\Delta EW_{Hb}$/$\Delta L_{5100}$ between the brightest
and dimmest states. The histogram distributions of the parameters are shown
aside the scatter plot. The median of each distribution is shown with
vertical/horizon lines in accordance style. From the EVQ sample, the luminous
sub-sample with $\log\overline{L}_{5100}>44.5$ are plotted with gray lines in
the histograms.  Figure 7: Upper: the correlation between the total luminosity
at 5100$\AA$ and the host galaxy fraction from the RM spectra of Shen et al.
(2014). The dashed line represents a simple linear regression (adopting
$L_{5100}$ as the independent variable), which is then utilized to estimate
the host fraction at given input $L_{5100}$ for our EVQs/CLQs plotted in the
lower panel. Lower: the iBeff of H$\beta$ of the EVQs and CLQs with the host
contamination corrected. The orthogonal distance regression fitting results
and the Spearman correlation coefficient $r$ (with p-values) of each sample
are listed.

In Fig. 5, we plot the variation of line EW along with the continuum of our
H$\beta$ EVQ sample, in comparison with the 33 archived H$\beta$-selected
CLQs. As already shown in Fig. 2, there is no clear correlation between
$\Delta EW_{Hb}$ and $\Delta L_{5100}$ in the EVQ sample, i.e., no iBeff.
Meanwhile a prominent positive correlation is observed in the CLQ sample
(i.e., a reverse iBeff). This is actually expected as the definition of
changing-look automatically requires much weaker H$\beta$ broad line in dim-
state CLQs. However, when it comes to the Mg II line for the same
H$\beta$-selected CLQs (i.e., lower panel of Fig. 5), CLQs also exhibit iBeff
(suggesting CLQs are physically similar to normal quasars), but with a slope
marginally flatter than EVQs.

We further examine the mean broad H$\beta$ EW versus the mean 5100Å
luminosity, and H$\beta$ iBeff slope versus the mean 5100Å luminosity, for the
EVQs and H$\beta$-selected CLQs (Fig. 6). The iBeff slope is calculated using
the brightest and dimmest epochs of each source. Similar to the results in
Fig. 5, the majority of CLQs have a reverse iBeff slope, with a median of
$0.60\pm 0.18$ 111We note that the only study in literature discussing the
change in EW of CLQs is presented in Green et al. (2022), which reports a
notably higher average iBeff slope of 2.05 compared to ours. However, their
result could be severely influenced by an individual target
J091234.00+262828.32, which has an iBeff slope of $\sim$107 due to its minor
change in $L_{5100}$, primarily affected by the host galaxy. The average slope
of the rest CLQs of Green et al. (2022) is comparable with our result. , while
the median for EVQs is $-0.11\pm 0.05$. The overall luminosity and mean
H$\beta$ EW of CLQs are lower than those of EVQs. However we do not see clear
correlation between mean broad H$\beta$ EW and 5100Å luminosity, or between
iBeff slope and 5100Å luminosity, for our uniformly selected EVQ sample. The
smaller mean H$\beta$ EW and 5100Å luminosity of CLQs (compared with EVQs) are
likely caused by observational biases, that EVQs with weaker broad H$\beta$
line and lower 5100Å luminosity are more likely be classified as CLQs (due to
limited spectral SNR).

Could the different iBeff slope of CLQs (compared with EVQs, as shown in Fig.
5 and 6) be due to potentially stronger host contamination in CLQs? Given the
overall low SNR of the SDSS spectra, directly decomposing the host component
in individual sources through spectral fitting is yet unreliable.
Nevertheless, by referencing to the results from Shen et al. (2014), we can
estimate the host fraction as a simple function of the luminosity at 5100Å and
give a rough correction for host contamination. In the upper panel of Fig. 7,
we present the correlation between the total luminosity (AGN plus host) at
5100Å, $L_{5100,tot}$, and the host fraction at 5100Å, $f_{host,5100}$,
decomposed by Shen et al. (2014) utilizing their reverberation mapping data.
Through simple linear regression we yield the expected average $f_{host,5100}$
as a function of $L_{5100,tot}$. For each EVQ/CLQ, we conservatively estimate
its host contribution using the $L_{5100,tot}$ from its dimmest epoch and
subtract the host contribution from the spectra of all epochs. The corrected
result is shown in the lower panel of Fig. 7, in which we could see clear
H$\beta$ iBeff in our EVQs, with a slope similar to that from the luminous
sub-sample (Fig. 3) and closer to the results in literature (Goad et al.,
2004; Rakić et al., 2017). This further suggests the absence of H$\beta$ iBeff
in EVQs shown in Fig. 2 is due to host contamination. However, after the
correction, the average iBeff slope of CLQs remains statistically different
from that of EVQs, which could be an expected consequence of the definition of
“changing-look”.

####  4.1.3 Are CLQs a distinct population or just a biased sub-population of
EVQs?

As we have shown above, the initial definition of “changing-look” requires the
disappearance of broad H$\beta$ line in the spectrum of dim state, but note
the definition is qualitative but not quantitative. In reality, the broad
H$\beta$ lines of CLQs selected in literature do weaken in the dim state, but
not “completely disappear” in most of them.

Meanwhile, the broad H$\beta$ line EW of our EVQs does not systematically vary
with continuum, i.e., on average with no iBeff, but with large scatter in the
iBeff slope in individual sources. Thus individual EVQs with reverse iBeff
slope (smaller broad H$\beta$ line EW in dimmer state) could more likely be
identified as “changing-look” due to the limited SNR of the spectra of dim
state. Meanwhile the CLQs selected in literature do exhibit reverse H$\beta$
iBeff as expected.

In this sense, CLQs are not necessarily a distinct population, but a sub-
population (biased by its definition) of EVQs. In this scheme, the fundamental
driven mechanism for both CLQs and EVQs is the continuum variation.

Previous researches implied that the lower Eddington ratio should be
responsible for the dramatic change of the BELs in CLQs (e.g. MacLeod et al.,
2016). Similarly, as shown in Fig. 6, we show CLQs have systematically lower
luminosity and lower H$\beta$ EW, compared with our EVQs. However, all these
facts could be attributed to selection biases: 1) the stronger host
contamination in sources with lower luminosity/Eddington ratio could hinder
the detection of broad H$\beta$ in the noisy dim state spectra; 2)
intrinsically weaker H$\beta$ lines are more difficult to be detected in the
noisy dim state spectra.

We note in literature while most CLQs were selected based on H$\beta$ line, Mg
II or C IV CLQs were rarely found (Guo et al., 2020; Ross et al., 2020). This
could also be naturally be attributed to the different iBeff slopes of
H$\beta$, Mg II and C IV we shown in Fig. 2. Since both Mg II and C IV of EVQs
exhibit strong iBeff, sources in the dim state tend to have larger EW of Mg II
and C IV, thus they are less likely be classified as changing-look. Meanwhile,
as EVQs exhibit no apparent iBeff on average in H$\beta$ line (like due to
host contamination), because of the large scatter in the iBeff slope from
source to source, some EVQs could exhibit reverse iBeff in H$\beta$ and are
more likely to identified as changing-look in their dim states.

The large scatter in the observed iBeff slope may partially be attributed to
the lag between broad line and continuum variation. Meanwhile, the fact that
the observed optical continuum variation does not necessarily coordinate with
that of the UV ionizing continuum (Xin et al., 2020; Gaskell et al., 2021)
could also play a significant role. Note such non-coordinated variations in
different bands (see also Sou et al. 2022) are natural subsequence of the
inhomogeneous disc fluctuation model (Cai et al., 2018, 2020).

To conclude, though we can not completely rule out dramatic vanishment of BLR
in some CLQs, most CLQs selected in literature appear consistent with a biased
sub-population of normal quasars with extreme continuum variability, i.e., the
essence of the “changing-look” phenomena is the extreme continuum variability.

###  4.2 Line Breathing

As shown in Fig. 2, the breathing of H$\beta$ is missing in our EVQ sample,
contradicting the results of previous works on AGNs and quasars (e.g., Denney
et al., 2008; Park et al., 2011; Wang et al., 2020). In addition, the
breathing relation in the Hb line within CLQ sample is also insignificant,
with a Spearman test $r=-0.29$ and $p=1.0\times 10^{-2}$.

It is not uncommon for AGNs to be caught undergoing changes in their breathing
relations. For example, Lu et al. (2022) reveal that the H$\beta$ of NGC 5548
well breathes in the most recent 5 years, however, previous decades of data
deviate significantly. Another interesting example is Mrk 50, in which a
breathing mode transition (from anti-breathing to normal breathing) was
captured during an $\sim$ 80-day RM observation (cf. Fig. 16 and Fig. 21 in
Barth et al. 2015). A more interesting fact is that the breathing period of
NGC 5548 coincides with its luminous stage, while Mrk 50 restored its
breathing when it was at its dimmest state, thus it seems the breathing mode
transition may not be simply one-way related to the flux level. Furthermore,
although the majority of sources in Wang et al. (2020) exhibit H$\beta$
breathing, a considerable proportion of their sources exhibit anti-breathing.
This suggests that breathing mode transitions might be frequent and the
structure of the BLRs could be more complex than expected in the simple
breathing model.

Meanwhile, it is notable that more than half of the sources in our EVQ sample
have a maximum spectral gap $>$3000 days which is much longer than the
duration of common RM projects. Therefore the absence of significant
breathing/anti-breathing on average in H$\beta$ (and in Mg II, and marginal in
C IV) in our EVQs could suggest possibly different evolution of the BLR (i.e.,
no breathing on average) on longer timescales.

##  5 Conclusions

In this work, we measured all 20,069 available SDSS spectra for the 14,012
EVQs built in Paper I with PyQSOFit. The catalog is presented in §A. From this
sample, we select 1259 EVQs with multi-epoch SDSS spectra after eliminating
spectra with suspicious calibrations, and study the spectral variations in
individual EVQs.

We find clear ”bluer-when-brighter” relation in multi-epoch spectra of EVQs,
consistent with previous works on normal quasars and AGNs. The iBeff of the
broad Mg II and C IV in EVQs are significantly seen. However, no robust iBeff
of H$\beta$ is detected, which could be attributed to stronger host
contamination at longer rest-frame wavelengths. Meanwhile, no systematical
variation of the broad line shape and asymmetry with continuum flux is found.
We do not detect statistically significant broad line breathing of H$\beta$,
Mg II or C IV either, suggesting the BLR evolution on longer timescales
($\sim$ 3000 days in the observed frame) could be different from those on
shorter timescales.

More interestingly, through comparing the iBeff of our EVQs with CLQs selected
in literature, we show that CLQs are more likely a biased (because of its
definition) sub-population of normal quasars with extreme continuum variation,
instead of a distinct population.

The work is supported by National Natural Science Foundation of China (grants
No. 11890693, 12033006, & 12192221), and the Cyrus Chung Ying Tang
Foundations.

##  Appendix A The EVQ Catalog

We publish the EVQs catalog along with all measured spectral quantities in
this paper. The format of the catalog is described in Table 1 and the full
catalog is available at https://doi.org/10.5281/zenodo.8328175.

Table 1: FITS catalog format

Number | Column Name | Format | Unit | Description
---|---|---|---|---
0 | SDSS_NAME | string |  |  Unique identifier from the SDSS DR14 quasar catalog
1 | RA | float32 | deg |  Right ascension (J2000)
2 | DEC | float32 | deg |  Declination (J2000)
3 | Z | float32 |  |  Redshift
4 | MEANMAG_G | float32 | mag |  Weighted mean magnitude of g-band light curve
5 | NG | int32 |  |  Number of observations in g-band light curve
6 | PLATE | int32 |  |  SDSS plate number
7 | MJD | int32 |  |  MJD when spectrum was observed
8 | FIBER | int32 |  |  SDSS fiber ID
9 | SPECPRIMARY | bool |  |  If the spectrum is the primary observation of object
10 | SPECMAG_G | float32 | mag |  Spectrum projected onto g filters
11 | SPECMAG_G_ERR | float32 | mag |  Error in SPECMAG_G
12 | SPECMAG_R | float32 | mag |  Spectrum projected onto r filters
13 | SPECMAG_R_ERR | float32 | mag |  Error in SPECMAG_R
14 | SN_CONTI | float32 |  |  Mean signal-to-noise ratio per pixel of continuum estimated at wavelength around 1350, 3000, and 5100 depending on the spectral coverage
15 | PL_SLOPE | float32 |  |  Slope of AGN power law
16 | PL_SLOPE_ERR | float32 |  |  Error in PL_SLOPE
17 | LOGL5100 | float32 | erg s-1 |  Logarithmic continuum luminosity at rest-frame 5100Å
18 | LOGL5100_ERR | float32 | erg s-1 |  Error in LOGL5100
19 | LOGL3000 | float32 | erg s-1 |  Logarithmic continuum luminosity at rest-frame 3000Å
20 | LOGL3000_ERR | float32 | erg s-1 |  Error in LOGL3000
21 | LOGL1350 | float32 | erg s-1 |  Logarithmic continuum luminosity at rest-frame 1350Å
22 | LOGL1350_ERR | float32 | erg s-1 |  Error in LOGL1350
23 | FE_2240_2650_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Rest-frame flux of the UV Fe II complex within the 2240-2650Å
24 | FE_2240_2650_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in FE_2240_2650_FLUX
25 | FE_2240_2650_EW | float32 | Å |  Rest-frame EW of the UV Fe II complex within the 2240-2650Å
26 | FE_2240_2650_EW_ERR | float32 | Å |  Error in FE_2240_2650_EW
27 | FE_4435_4685_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Rest-frame flux of the optical Fe II complex within the 4435-4685Å
28 | FE_4435_4685_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in FE_4435_4685_FLUX
29 | FE_4435_4685_EW | float32 | Å |  Rest-frame EW of the optical Fe II complex within the 4435_4685Å
30 | FE_4435_4685_EW_ERR | float32 | Å |  Error in FE_4435_4685_EW
31 | HB_BR_PEAK | float32 | Å |  Peak wavelength of Hb broad component
32 | HB_BR_PEAK_ERR | float32 | Å |  Error in HB_BR_PEAK
33 | HB_BR_BISECT | float32 | Å |  Wavelength bisect the area of Hb broad component
34 | HB_BR_BISECT_ERR | float32 | Å |  Error in HB_BR_BISECT
35 | HB_BR_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Flux of Hb broad component
36 | HB_BR_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in HB_BR_FLUX
37 | HB_BR_EW | float32 | Å |  Rest-frame EW of Hb broad component
38 | HB_BR_EW_ERR | float32 | Å |  Error in HB_BR_EW
39 | HB_BR_SIGMA | float32 | km s-1 |  Line dispersion of Hb broad component
40 | HB_BR_SIGMA_ERR | float32 | km s-1 |  Error in HB_BR_SIGMA
41 | HB_BR_FWHM | float32 | km s-1 |  FWHM of Hb broad component
42 | HB_BR_FWHM_ERR | float32 | km s-1 |  Error in HB_BR_FWHM
43 | HB_BR_FWQM | float32 | km s-1 |  FWQM of Hb broad component
44 | HB_BR_FWQM_ERR | float32 | km s-1 |  Error in HB_BR_FWQM
45 | HB_BR_FW10M | float32 | km s-1 |  FW10M of Hb broad component
46 | HB_BR_FW10M_ERR | float32 | km s-1 |  Error in HB_BR_FW10M
47 | HB_BR_Z50 | float32 | km s-1 |  The half maximum center shift of Hb broad component
48 | HB_BR_Z50_ERR | float32 | km s-1 |  Error in HB_BR_Z50
49 | HB_BR_Z25 | float32 | km s-1 |  The quarter maximum center shift of Hb broad component
50 | HB_BR_Z25_ERR | float32 | km s-1 |  Error in HB_BR_Z25
51 | HB_BR_Z10 | float32 | km s-1 |  The 10 percent maximum center shift of Hb broad component
52 | HB_BR_Z10_ERR | float32 | km s-1 |  Error in HB_BR_Z10
53 | HB_NA_PEAK | float32 | Å |  Peak wavelength of Hb narrow component
54 | HB_NA_PEAK_ERR | float32 | Å |  Error in HB_NA_PEAK
55 | HB_NA_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Flux of Hb narrow component
56 | HB_NA_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in HB_NA_FLUX
57 | HB_NA_EW | float32 | Å |  Rest-frame EW of Hb narrow component
58 | HB_NA_EW_ERR | float32 | Å |  Error in HB_NA_EW
59 | HB_NA_SIGMA | float32 | km s-1 |  Line dispersion of Hb narrow component
60 | HB_NA_SIGMA_ERR | float32 | km s-1 |  Error in HB_NA_SIGMA
61 | HB_NA_FWHM | float32 | km s-1 |  FWHM of Hb narrow component
62 | HB_NA_FWHM_ERR | float32 | km s-1 |  Error in HB_NA_FWHM
63 | OIII4959C_PEAK | float32 | Å |  Peak wavelength of OIII4959 core component
64 | OIII4959C_PEAK_ERR | float32 | Å |  Error in OIII4959C_PEAK
65 | OIII4959C_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Flux of OIII4959 core component
66 | OIII4959C_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in OIII4959C_FLUX
67 | OIII4959C_EW | float32 | Å |  Rest-frame EW of OIII4959 core component
68 | OIII4959C_EW_ERR | float32 | Å |  Error in OIII4959C_EW
69 | OIII4959C_SIGMA | float32 | km s-1 |  Line dispersion of OIII4959 core component
70 | OIII4959C_SIGMA_ERR | float32 | km s-1 |  Error in OIII4959C_SIGMA
71 | OIII4959C_FWHM | float32 | km s-1 |  FWHM of OIII4959 core component
72 | OIII4959C_FWHM_ERR | float32 | km s-1 |  Error in OIII4959C_FWHM
73 | OIII4959W_PEAK | float32 | Å |  Peak wavelength of OIII4959 wing component
74 | OIII4959W_PEAK_ERR | float32 | Å |  Error in OIII4959W_PEAK
75 | OIII4959W_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Flux of OIII4959 wing component
76 | OIII4959W_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in OIII4959W_FLUX
77 | OIII4959W_EW | float32 | Å |  Rest-frame EW of OIII4959 wing component
78 | OIII4959W_EW_ERR | float32 | Å |  Error in OIII4959W_EW
79 | OIII4959W_SIGMA | float32 | km s-1 |  Line dispersion of OIII4959 wing component
80 | OIII4959W_SIGMA_ERR | float32 | km s-1 |  Error in OIII4959W_SIGMA
81 | OIII4959W_FWHM | float32 | km s-1 |  FWHM of OIII4959 wing component
82 | OIII4959W_FWHM_ERR | float32 | km s-1 |  Error in OIII4959W_FWHM
83 | OIII5007C_PEAK | float32 | Å |  Peak wavelength of OIII5007 core component
84 | OIII5007C_PEAK_ERR | float32 | Å |  Error in OIII5007C_PEAK
85 | OIII5007C_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Flux of OIII5007 core component
86 | OIII5007C_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in OIII5007C_FLUX
87 | OIII5007C_EW | float32 | Å |  Rest-frame EW of OIII5007 core component
88 | OIII5007C_EW_ERR | float32 | Å |  Error in OIII5007C_EW
89 | OIII5007C_SIGMA | float32 | km s-1 |  Line dispersion of OIII5007 core component
90 | OIII5007C_SIGMA_ERR | float32 | km s-1 |  Error in OIII5007C_SIGMA
91 | OIII5007C_FWHM | float32 | km s-1 |  FWHM of OIII5007 core component
92 | OIII5007C_FWHM_ERR | float32 | km s-1 |  Error in OIII5007C_FWHM
93 | OIII5007W_PEAK | float32 | Å |  Peak wavelength of OIII5007 wing component
94 | OIII5007W_PEAK_ERR | float32 | Å |  Error in OIII5007W_PEAK
95 | OIII5007W_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Flux of OIII5007 wing component
96 | OIII5007W_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in OIII5007W_FLUX
97 | OIII5007W_EW | float32 | Å |  Rest-frame EW of OIII5007 wing component
98 | OIII5007W_EW_ERR | float32 | Å |  Error in OIII5007W_EW
99 | OIII5007W_SIGMA | float32 | km s-1 |  Line dispersion of OIII5007 wing component
100 | OIII5007W_SIGMA_ERR | float32 | km s-1 |  Error in OIII5007W_SIGMA
101 | OIII5007W_FWHM | float32 | km s-1 |  FWHM of OIII5007 wing component
102 | OIII5007W_FWHM_ERR | float32 | km s-1 |  Error in OIII5007W_FWHM
103 | MGII_BR_PEAK | float32 | Å |  Peak wavelength of MgII broad component
104 | MGII_BR_PEAK_ERR | float32 | Å |  Error in MGII_BR_PEAK
105 | MGII_BR_BISECT | float32 | Å |  Wavelength bisect the area of MgII broad component
106 | MGII_BR_BISECT_ERR | float32 | Å |  Error in MGII_BR_BISECT
107 | MGII_BR_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Flux of MgII broad component
108 | MGII_BR_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in MGII_BR_FLUX
109 | MGII_BR_EW | float32 | Å |  Rest-frame EW of MgII broad component
110 | MGII_BR_EW_ERR | float32 | Å |  Error in MGII_BR_EW
111 | MGII_BR_SIGMA | float32 | km s-1 |  Line dispersion of MgII broad component
112 | MGII_BR_SIGMA_ERR | float32 | km s-1 |  Error in MGII_BR_SIGMA
113 | MGII_BR_FWHM | float32 | km s-1 |  FWHM of MgII broad component
114 | MGII_BR_FWHM_ERR | float32 | km s-1 |  Error in MGII_BR_FWHM
115 | MGII_BR_FWQM | float32 | km s-1 |  FWQM of MgII broad component
116 | MGII_BR_FWQM_ERR | float32 | km s-1 |  Error in MGII_BR_FWQM
117 | MGII_BR_FW10M | float32 | km s-1 |  FW10M of MgII broad component
118 | MGII_BR_FW10M_ERR | float32 | km s-1 |  Error in MGII_BR_FW10M
119 | MGII_BR_Z50 | float32 | km s-1 |  The half maximum center shift of MgII broad component
120 | MGII_BR_Z50_ERR | float32 | km s-1 |  Error in MGII_BR_Z50
121 | MGII_BR_Z25 | float32 | km s-1 |  The quarter maximum center shift of MgII broad component
122 | MGII_BR_Z25_ERR | float32 | km s-1 |  Error in MGII_BR_Z25
123 | MGII_BR_Z10 | float32 | km s-1 |  The 10 percent maximum center shift of MgII broad component
124 | MGII_BR_Z10_ERR | float32 | km s-1 |  Error in MGII_BR_Z10
125 | MGII_NA_PEAK | float32 | Å |  Peak wavelength of MgII narrow component
126 | MGII_NA_PEAK_ERR | float32 | Å |  Error in MGII_NA_PEAK
127 | MGII_NA_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Flux of MgII narrow component
128 | MGII_NA_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in MGII_NA_FLUX
129 | MGII_NA_EW | float32 | Å |  Rest-frame EW of MgII narrow component
130 | MGII_NA_EW_ERR | float32 | Å |  Error in MGII_NA_EW
131 | MGII_NA_SIGMA | float32 | km s-1 |  Line dispersion of MgII narrow component
132 | MGII_NA_SIGMA_ERR | float32 | km s-1 |  Error in MGII_NA_SIGMA
133 | MGII_NA_FWHM | float32 | km s-1 |  FWHM of MgII narrow component
134 | MGII_NA_FWHM_ERR | float32 | km s-1 |  Error in MGII_NA_FWHM
135 | CIV_BR_PEAK | float32 | Å |  Peak wavelength of CIV broad component
136 | CIV_BR_PEAK_ERR | float32 | Å |  Error in CIV_BR_PEAK
137 | CIV_BR_BISECT | float32 | Å |  Wavelength bisect the area of CIV broad component
138 | CIV_BR_BISECT_ERR | float32 | Å |  Error in CIV_BR_BISECT
139 | CIV_BR_FLUX | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Flux of CIV broad component
140 | CIV_BR_FLUX_ERR | float32 | $10^{-17}\text{erg cm${}^{-2}$ s${}^{-1}$}$ |  Error in CIV_BR_FLUX
141 | CIV_BR_EW | float32 | Å |  Rest-frame EW of CIV broad component
142 | CIV_BR_EW_ERR | float32 | Å |  Error in CIV_BR_EW
143 | CIV_BR_SIGMA | float32 | km s-1 |  Line dispersion of CIV broad component
144 | CIV_BR_SIGMA_ERR | float32 | km s-1 |  Error in CIV_BR_SIGMA
145 | CIV_BR_FWHM | float32 | km s-1 |  FWHM of CIV broad component
146 | CIV_BR_FWHM_ERR | float32 | km s-1 |  Error in CIV_BR_FWHM
147 | CIV_BR_FWQM | float32 | km s-1 |  FWQM of CIV broad component
148 | CIV_BR_FWQM_ERR | float32 | km s-1 |  Error in CIV_BR_FWQM
149 | CIV_BR_FW10M | float32 | km s-1 |  FW10M of CIV broad component
150 | CIV_BR_FW10M_ERR | float32 | km s-1 |  Error in CIV_BR_FW10M
151 | CIV_BR_Z50 | float32 | km s-1 |  The half maximum center shift of CIV broad component
152 | CIV_BR_Z50_ERR | float32 | km s-1 |  Error in CIV_BR_Z50
153 | CIV_BR_Z25 | float32 | km s-1 |  The quarter maximum center shift of CIV broad component
154 | CIV_BR_Z25_ERR | float32 | km s-1 |  Error in CIV_BR_Z25
155 | CIV_BR_Z10 | float32 | km s-1 |  The 10 percent maximum center shift of CIV broad component
156 | CIV_BR_Z10_ERR | float32 | km s-1 |  Error in CIV_BR_Z10
157 | LOGLBOL_5100 | float32 | erg s-1 |  Logarithmic bolometric luminosity estimated based on LOGL5100
158 | LOGLBOL_5100_ERR | float32 | erg s-1 |  Error in LOGLBOL_5100
159 | LOGLBOL_3000 | float32 | erg s-1 |  Logarithmic bolometric luminosity estimated based on LOGL3000
160 | LOGLBOL_3000_ERR | float32 | erg s-1 |  Error in LOGLBOL_3000
161 | LOGLBOL_1350 | float32 | erg s-1 |  Logarithmic bolometric luminosity estimated based on LOGL1350
162 | LOGLBOL_1350_ERR | float32 | erg s-1 |  Error in LOGLBOL_1350
163 | LOGMBH_HB | float32 | $M_{\odot}$ |  Logarithmic black hole mass estimated based on broad Hb line
164 | LOGMBH_HB_ERR | float32 | $M_{\odot}$ |  Error in LOGMBH_HB
165 | LOGMBH_MGII | float32 | $M_{\odot}$ |  Logarithmic black hole mass estimated based on broad MgII line
166 | LOGMBH_MGII_ERR | float32 | $M_{\odot}$ |  Error in LOGMBH_MGII
167 | LOGMBH_CIV | float32 | $M_{\odot}$ |  Logarithmic black hole mass estimated based on broad CIV line
168 | LOGMBH_CIV_ERR | float32 | $M_{\odot}$ |  Error in LOGMBH_CIV
169 | LOGLBOL | float32 | erg s-1 |  The adopted fiducial bolometric luminosity
170 | LOGLBOL_ERR | float32 | erg s-1 |  Error in LOGLBOL
171 | LOGMBH | float32 | $M_{\odot}$ |  The adopted fiducial black hole mass
172 | LOGMBH_ERR | float32 | $M_{\odot}$ |  Error in LOGMBH
173 | LOGREDD | float32 |  |  Logarithmic Eddington ratio based on fiducial BH mass and bolometric luminosity

Note. — We provide all 20,069 available spectral measurements of 14,012 EVQs
selected in Ren et al. 2022. Repeated observed spectra of a same EVQ will have
the same SDSS_NAME with different spectral info (PLATE,MJD, and FIBER). We
include the SPECPRIMARY flag to indicate if the spectrum is the best
observation of this object. The unmeasurable parameters are set to -999. The
errors are obtained from 100 iterations of Monte Carlo simulation. The
complete table is available on https://doi.org/10.5281/zenodo.8328175

##  Appendix B Evaluate spectral calibration

The fiber-drop (Dawson et al., 2013) is an unavoidable problem in SDSS
spectrophotometry leading to varying degrees of underestimation of spectral
flux and has been widely reported in literature (e.g. Shen et al., 2014; Sun
et al., 2015; Guo et al., 2020). In certain cases the spectrophotometry could
also be overestimated due to contamination to the line of sight by spurious or
unrelated signals. Such effects may produce artificial spectral variability
which is hard to distinguish from the intrinsic spectral variation. This would
be particularly more relevant if one aims to search for rare events (such as
CLQs) out of a huge number of spectra. In this work our EVQs are pre-selected
based on multi-epoch photometric observations. Below we develop a technique to
evaluate the reliability of the SDSS spectrophotometry of the EVQs through
comparing the spectroscopic photometry with the photometric light curve
modeled with damped random walk (DRW). The spectroscopic photometry which
severely deviates from the DRW yielded ranges would be excluded as potentially
unreliable.

###  B.1 Extend the light curves

The light curves we built for EVQ selection, consisted of data from SDSS and
PS1, only cover a span of 1998 to 2014. About half of the SDSS spectra were
conducted by SDSS IV since 2014 (Blanton et al., 2017), significantly beyond
the time span of our photometric light curves. We then introduce the PTF/iPTF
and ZTF observations to extend the photometric light curves to cover the
epochs when SDSS IV spectra were obtained.

The Palomar Transient Factory (PTF) is a fully-automated, wide-field survey
conducted by the Palomar 48-inch Samuel Oschin Schmidt telescope with
12K$\times$8K CCD array during 2009 to 2012 (Rau et al., 2009; Law et al.,
2009). The intermediate Palomar Transient Factory (iPTF) is the successor PTF
ran from 2013 to 2017 with a relatively higher cadence. The Zwicky Transient
Facility (ZTF) is a next-generation optical time-domain survey build upon the
PTF/iPTF. It started from 2018, scanned 3750 $\rm deg^{2}/hour$, and is still
in operation (Masci et al., 2018). The PTF/iPTF and ZTF surveys are shallower
than the SDSS and the PS1 with relatively larger photometry error.
Nevertheless, among the 3042 EVQs with repeated SDSS spectra, only 11 of them
have no detection in either PTF, iPTF or ZTF.

###  B.2 Evaluating the reliability of spectroscopic photometry

We use the celerite (Foreman-Mackey et al., 2017) to model our $g$\- and $r$\-
band light curves and estimate their DRW parameters (timescale $\tau$ and
variability amplitude $\sigma$) with maximum likelihood estimation (MLE)
method in scipy. The mean magnitude $\bar{m}$ of DRW model is assumed to be
the mean of the light curve. The detailed explanation of fitting a light curve
with DRW process is presented in Zu et al. (2011, 2013, 2016).

Next, we utilize the DRW parameters to determine the probability distribution
function of the magnitude at a given epoch. According to the DRW process, the
next epoch signal ($s_{n+1}$ at $t=t_{n+1}$) is only relative to the present
one $s_{n}$, which can be expressed as:

$p(s_{n+1}=x_{n+1}|s_{n}=x_{n})=N(\mu_{n+1},\eta_{n+1}^{2}),$ (B1)

where

$\mu_{n+1}=e^{-\frac{\Delta t}{\tau}}x_{n}+\bar{m}(1-e^{-\frac{\Delta
t}{\tau}}),$ (B2) $\eta_{n+1}^{2}=\sigma^{2}(1-e^{-\frac{2\Delta
t}{\tau}})+err_{n}^{2}+err_{n+1}^{2}.$ (B3)

The $s_{n}$ and $s_{n+1}$ are the magnitude of the couple of adjacent light
curve data at $t_{n}$ and $t_{n+1}$ respectively, $err$ the observational
error, and $\Delta t=t_{n+1}-t_{n}$.

In any gap of the light curve, with the magnitude known at both ends, we can
calculate the probability distribution function (PDF) of the magnitude at any
epoch ($t_{n}<t_{0}<t_{n+1}$) within this duration. According to Bayes
formula, the PDF of the signal at $t_{0}$ given its adjacent data can be
determined as:

$\begin{split}&p(s_{0}=x_{0}|s_{n}=x_{n},~{}s_{n+1}=x_{n+1})\\\
=&\frac{p(s_{0}=x_{0},~{}s_{n+1}=x_{n+1}|s_{n}=x_{n})}{p(s_{n+1}=x_{n+1}|s_{n}=x_{n})}\\\
=&\frac{p(s_{0}=x_{0}|s_{n}=x_{n})p(s_{n+1}=x_{n+1}|s_{0}=x_{0},~{}s_{n}=x_{n})}{p(s_{n+1}=x_{n+1}|s_{n}=x_{n})}\\\
=&\frac{p(s_{0}=x_{0}|s_{n}=x_{n})p(s_{n+1}=x_{n+1}|s_{0}=x_{0})}{p(s_{n+1}=x_{n+1}|s_{n}=x_{n})}.\end{split}$
(B4)

To evaluate the reliability of a spectroscopic photometry, we first set
$t_{0}$ to the date of the spectrum and find its adjacent photometric data in
the corresponding light curve to calculate the PDF of the magnitude at this
epoch. By integrating this PDF, we can give the confidence interval of the
true magnitude of this spectrum. The sketch of the 99.7% confidence interval
is shown in Fig. 8.

Figure 8: An example of $g$-band light curve showing the 99.7% confidence
intervals of the magnitude in the light curve gaps. The observed photometric
magnitudes from different surveys are plotted with distinct colors and marks.
The spectroscopic magnitudes of three SDSS spectra are marked with stars and
vertical dashed lines. As can be seen, the spectroscopic magnitude of the
first spectrum is out of the 99.7% confidence intervals, which is likely due
to the fiber-drop effect.

In this work, we consider the spectroscopic photometry reliable if both its
$g$\- and $r$\- band spectroscopic magnitude fall within the 99.7% confidence
intervals. Under this criterion, 301 spectra are excluded as potentially
unreliable, leaving 1259 sources with repeated observations.

###  B.3 The validity of the approach

Precisely evaluating whether this method can properly rule out unreliable
spectroscopic photometry could be challenging. However, we can provide an
indirect approach to demonstrate its validity. We further select a sample of
low-variability quasars (LVQs) with $|\Delta g|_{\text{max}}<0.2$ mag, where
$|\Delta g|_{\text{max}}$ represents the maximum difference in magnitude
between any two epochs in the light curves constructed by Paper I. Since we do
not expect a sudden change in flux for these low-variability quasars, spectra
with significantly deviating flux from their light curve means in LVQs are
likely caused by calibration issues.

In Fig. 9, we show the normalized distribution of
$g_{\text{spec}}-g_{\text{mean}}$ (the difference between the $g$-band
magnitude from the spectrum and the light curve mean, see §3.1 in Paper I for
details) for our sample and the LVQs. The distribution of
$g_{\text{spec}}-g_{\text{mean}}$ of LVQs (for which we only expect small
deviation from zero) does span a broad range, indicating the calibration
issues are not negligible at large $|g_{\text{spec}}-g_{\text{mean}}|$. For
our EVQs, our screening process does exclude considerable fraction of epochs
at large $|g_{\text{spec}}-g_{\text{mean}}|$ but minor fraction at small
$|g_{\text{spec}}-g_{\text{mean}}|$ (the ratio of blue to yellow line in Fig.
9 could represent the excluded fraction as a function of
$g_{\text{spec}}-g_{\text{mean}}$). At large
$|g_{\text{spec}}-g_{\text{mean}}|$ (e.g. $>$ 0.8), we have excluded
relatively more epochs than the prediction from LVQs (blue line lying above
the green line), suggesting our approach is somehow conservative that some
good spectra could have been excluded (which is however unavoidable).

Figure 9: Normalized distribution of $g_{\text{spec}}-g_{\text{mean}}$ for
multi-spectra EVQs (yellow) and LVQs (green). The blue solid line plots the
excluded portion of the yellow line. A logarithmic y-axis is used to better
visualize the outliers, i.e., with large $|g_{\text{spec}}-g_{\text{mean}}|$.

## References

  * Aretxaga et al. (1999) Aretxaga, I., Joguet, B., Kunth, D., Melnick, J., & Terlevich, R. J. 1999, The Astrophysical Journal, 519, L123, doi: 10.1086/312114
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Barth et al. (2015) Barth, A. J., Bennert, V. N., Canalizo, G., et al. 2015, ApJS, 217, 26, doi: 10.1088/0067-0049/217/2/26
  * Blanton et al. (2017) Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28, doi: 10.3847/1538-3881/AA7567
  * Boggs & Rogers (1989) Boggs, P. T. P., & Rogers, J. J. E. 1989, Nistir 89-4197, 1989, doi: 10.6028/nist.ir.89-4197
  * Cai et al. (2016) Cai, Z.-Y., Wang, J.-X., Gu, W.-M., et al. 2016, ApJ, 826, 7, doi: 10.3847/0004-637X/826/1/7
  * Cai et al. (2020) Cai, Z.-Y., Wang, J.-X., & Sun, M. 2020, ApJ, 892, 63, doi: 10.3847/1538-4357/ab7991
  * Cai et al. (2018) Cai, Z.-Y., Wang, J.-X., Zhu, F.-F., et al. 2018, ApJ, 855, 117, doi: 10.3847/1538-4357/aab091
  * Cohen et al. (1986) Cohen, R. D., Puetter, R. C., Rudy, R. J., Ake, T. B., & Foltz, C. B. 1986, ApJ, 311, 135, doi: 10.1086/164758
  * Dawson et al. (2013) Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10, doi: 10.1088/0004-6256/145/1/10
  * Denney et al. (2008) Denney, K. D., Peterson, B. M., Dietrich, M., Vestergaard, M., & Bentz, M. C. 2008, ApJ, 692, 246, doi: 10.1088/0004-637X/692/1/246
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, The Astronomical Journal, 154, 220, doi: 10.3847/1538-3881/AA9332
  * Gaskell et al. (2021) Gaskell, C. M., Bartel, K., Deffner, J. N., & Xia, I. 2021, MNRAS, 508, 6077, doi: 10.1093/mnras/stab2443
  * Goad et al. (2004) Goad, M. R., Korista, K. T., & Knigge, C. 2004, MNRAS, 352, 277, doi: 10.1111/j.1365-2966.2004.07919.x
  * Goodrich (1995) Goodrich, R. W. 1995, The Astrophysical Journal, 440, 141, doi: 10.1086/175256
  * Green et al. (2022) Green, P. J., Pulgarin-Duque, L., Anderson, S. F., et al. 2022, ApJ, 933, 180, doi: 10.3847/1538-4357/ac743f
  * Grier et al. (2017) Grier, C. J., Trump, J. R., Shen, Y., et al. 2017, ApJ, 851, 21, doi: 10.3847/1538-4357/aa98dc
  * Grier et al. (2019) Grier, C. J., Shen, Y., Horne, K., et al. 2019, ApJ, 887, 38, doi: 10.3847/1538-4357/ab4ea5
  * Guo & Gu (2016) Guo, H., & Gu, M. 2016, ApJ, 822, 26, doi: 10.3847/0004-637x/822/1/26
  * Guo et al. (2018) Guo, H., Shen, Y., & Wang, S. 2018, PyQSOFit: Python code to fit the spectrum of quasars, Astrophysics Source Code Library.  http://ascl.net/1809.008
  * Guo et al. (2019) Guo, H., Sun, M., Liu, X., et al. 2019, ApJ, 883, L44, doi: 10.3847/2041-8213/ab4138
  * Guo et al. (2020) Guo, H., Shen, Y., He, Z., et al. 2020, ApJ, 888, 58, doi: 10.3847/1538-4357/ab5db0
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Homan et al. (2020) Homan, D., MacLeod, C. L., Lawrence, A., Ross, N. P., & Bruce, A. 2020, MNRAS, 496, 309, doi: 10.1093/mnras/staa1467
  * Homayouni et al. (2020) Homayouni, Y., Trump, J. R., Grier, C. J., et al. 2020, ApJ, 901, 55, doi: 10.3847/1538-4357/ababa9
  * Hu et al. (2015) Hu, C., Du, P., Lu, K.-X., et al. 2015, ApJ, 804, 138, doi: 10.1088/0004-637X/804/2/138
  * Kang et al. (2021) Kang, W.-Y., Wang, J.-X., Cai, Z.-y., & Ren, W.-K. 2021, The Astrophysical Journal, 911, 148, doi: 10.3847/1538-4357/abeb69
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895, doi: 10.1088/0004-637X/698/1/895
  * Kong et al. (2006) Kong, M. Z., Wu, X. B., Wang, R., Liu, F. K., & Han, J. L. 2006, A&A, 456, 473, doi: 10.1051/0004-6361:20064871
  * LaMassa et al. (2015) LaMassa, S. M., Cales, S., Moran, E. C., et al. 2015, ApJ, 800, 144, doi: 10.1088/0004-637X/800/2/144
  * Law et al. (2009) Law, N. M., Kulkarni, S. R., Dekany, R. G., et al. 2009, Publications of the Astronomical Society of the Pacific, 121, 1395, doi: 10.1086/648598/XML
  * López-Navas et al. (2022) López-Navas, E., Martínez-Aldama, M. L., Bernal, S., et al. 2022, Mon. Not. R. Astron. Soc. Lett., 513, L57, doi: 10.1093/mnrasl/slac033
  * Lu et al. (2022) Lu, K.-X., Bai, J.-M., Wang, J.-M., et al. 2022, ApJS, 263, 10, doi: 10.3847/1538-4365/ac94d3
  * MacLeod et al. (2010) MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, The Astrophysical Journal, 721, 1014, doi: 10.1088/0004-637X/721/2/1014
  * MacLeod et al. (2016) MacLeod, C. L., Ross, N. P., Lawrence, A., et al. 2016, MNRAS, 457, 389, doi: 10.1093/mnras/stv2997
  * MacLeod et al. (2019) MacLeod, C. L., Green, P. J., Anderson, S. F., et al. 2019, ApJ, 874, 8, doi: 10.3847/1538-4357/ab05e2
  * Masci et al. (2018) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2018, Publications of the Astronomical Society of the Pacific, 131, 018003, doi: 10.1088/1538-3873/AAE8AC
  * pandas development team (2020) pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.3509134
  * Pâris et al. (2018) Pâris, I., Petitjean, P., Aubourg, É., et al. 2018, A&A, 613, A51, doi: 10.1051/0004-6361/201732445
  * Park et al. (2011) Park, D., Woo, J.-H., Treu, T., et al. 2011, ApJ, 747, 30, doi: 10.1088/0004-637X/747/1/30
  * Pogge & Peterson (1992) Pogge, R. W., & Peterson, B. M. 1992, AJ, 103, 1084, doi: 10.1086/116127
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Pu et al. (2006) Pu, X., Bian, W., & Huang, K. 2006, MNRAS, 372, 246, doi: 10.1111/j.1365-2966.2006.10850.x
  * Rakić (2022) Rakić, N. 2022, MNRAS, 516, 1624, doi: 10.1093/mnras/stac2259
  * Rakić et al. (2017) Rakić, N., La Mura, G., Ilić, D., et al. 2017, A&A, 603, A49, doi: 10.1051/0004-6361/201630085
  * Rau et al. (2009) Rau, A., Kulkarni, S. R., Law, N. M., et al. 2009, Publications of the Astronomical Society of the Pacific, 121, 1334, doi: 10.1086/605911/XML
  * Ren et al. (2022) Ren, W., Wang, J., Cai, Z., & Guo, H. 2022, ApJ, 925, 50, doi: 10.3847/1538-4357/ac3828
  * Risaliti et al. (2007) Risaliti, G., Elvis, M., Fabbiano, G., et al. 2007, The Astrophysical Journal, 659, L111, doi: 10.1086/517884
  * Ross et al. (2020) Ross, N. P., Graham, M. J., Calderone, G., et al. 2020, MNRAS, 498, 2339, doi: 10.1093/mnras/staa2415
  * Ross et al. (2018) Ross, N. P., Ford, K. E., Graham, M., et al. 2018, Monthly Notices of the Royal Astronomical Society, 480, 4468, doi: 10.1093/mnras/sty2002
  * Ruan et al. (2019) Ruan, J. J., Anderson, S. F., Eracleous, M., et al. 2019, doi: 10.48550/arXiv.1909.04676
  * Ruan et al. (2016) Ruan, J. J., Anderson, S. F., Cales, S. L., et al. 2016, ApJ, 826, 188, doi: 10.3847/0004-637X/826/2/188
  * Rumbaugh et al. (2018) Rumbaugh, N., Shen, Y., Morganson, E., et al. 2018, ApJ, 854, 160, doi: 10.3847/1538-4357/aaa9b6
  * Sesar et al. (2007) Sesar, B., Ivezić, Ž., Lupton, R. H., et al. 2007, The Astronomical Journal, 134, 2236, doi: 10.1086/521819
  * Shen et al. (2011) Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45, doi: 10.1088/0067-0049/194/2/45
  * Shen et al. (2014) Shen, Y., Brandt, W. N., Dawson, K. S., et al. 2014, ApJS, 216, 4, doi: 10.1088/0067-0049/216/1/4
  * Sheng et al. (2020) Sheng, Z., Wang, T., Jiang, N., et al. 2020, ApJ, 889, 46, doi: 10.3847/1538-4357/ab5af9
  * Sou et al. (2022) Sou, H., Wang, J.-X., Xie, Z.-L., Kang, W.-Y., & Cai, Z.-Y. 2022, MNRAS, 512, 5511, doi: 10.1093/mnras/stac738
  * Sun et al. (2018) Sun, M., Xue, Y., Richards, G. T., et al. 2018, ApJ, 854, 128, doi: 10.3847/1538-4357/aaa890
  * Sun et al. (2015) Sun, M., Trump, J. R., Shen, Y., et al. 2015, ApJ, 811, 42, doi: 10.1088/0004-637X/811/1/42
  * Sun et al. (2020) Sun, M., Xue, Y., Brandt, W. N., et al. 2020, The Astrophysical Journal, 891, 178, doi: 10.3847/1538-4357/ab789e
  * Sun et al. (2014) Sun, Y.-H., Wang, J.-X., Chen, X.-Y., & Zheng, Z.-Y. 2014, ApJ, 792, 54, doi: 10.1088/0004-637X/792/1/54
  * Tholine & Osterbrock (1976) Tholine, J. E., & Osterbrock, D. E. 1976, ApJ, 210, L117, doi: 10.1086/182317
  * Trevese & Vagnetti (2002) Trevese, D., & Vagnetti, F. 2002, ApJ, 564, 624, doi: 10.1086/324541
  * Vanden Berk et al. (2004) Vanden Berk, D. E., Wilhite, B. C., Kron, R. G., et al. 2004, The Astrophysical Journal, 601, 692, doi: 10.1086/380563
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang et al. (2020) Wang, S., Shen, Y., Jiang, L., et al. 2020, ApJ, 903, 51, doi: 10.3847/1538-4357/abb36d
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * Wu & Shen (2022) Wu, Q., & Shen, Y. 2022, Astrophys. J. Suppl., 263, 42, doi: 10.3847/1538-4365/ac9ead
  * Xin et al. (2020) Xin, C., Charisi, M., Haiman, Z., & Schiminovich, D. 2020, MNRAS, 495, 1403, doi: 10.1093/mnras/staa1258
  * Yang et al. (2018) Yang, Q., Wu, X.-B., Fan, X., et al. 2018, ApJ, 862, 109, doi: 10.3847/1538-4357/aaca3a
  * Yang et al. (2023) Yang, Q., Green, P. J., MacLeod, C. L., et al. 2023, The Astrophysical Journal, 953, 61, doi: 10.3847/1538-4357/acdedd
  * Zajaček et al. (2020) Zajaček, M., Czerny, B., Martinez–Aldama, M. L., et al. 2020, ApJ, 896, 146, doi: 10.3847/1538-4357/ab94ae
  * Zhang (2013) Zhang, X. G. 2013, MNRAS, 435, 2141, doi: 10.1093/mnras/stt1432
  * Zhu et al. (2016) Zhu, F.-F., Wang, J.-X., Cai, Z.-Y., & Sun, Y.-H. 2016, ApJ, 832, 75, doi: 10.3847/0004-637X/832/1/75
  * Zu et al. (2016) Zu, Y., Kochanek, C. S., Kozłowski, S., & Peterson, B. M. 2016, ApJ, 819, 122, doi: 10.3847/0004-637X/819/2/122
  * Zu et al. (2013) Zu, Y., Kochanek, C. S., Kozłowski, S., & Udalski, A. 2013, ApJ, 765, 106, doi: 10.1088/0004-637X/765/2/106
  * Zu et al. (2011) Zu, Y., Kochanek, C. S., & Peterson, B. M. 2011, ApJ, 735, 80, doi: 10.1088/0004-637X/735/2/80