{'fraction_non_alphanumeric': 0.04662756317637715, 'fraction_numerical': 0.03561415763235702, 'mean_word_length': 4.38878003506239, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Local Active Galactic Nuclei with Large Broad-H{\alpha} Variability Reside in Red Galaxies

Inspired by our serendipitous discovery of six AGNs with varying broad-Halpha fluxes over years out of our searching for intermediate-mass black holes (IMBHs), we conduct a systematic investigation of changing-look (CL) and large-variability AGNs. We collect all the CL AGNs at z<0.15 and the reverberation mapped AGNs with strongly variable broad Halpha, and perform careful decomposition fittings to both their images and spectra. We find two observational facts: (1) The host galaxies of local CL and large-variability AGNs, mainly being Seyferts, are in the red (gas-poor) tail of the general Seyfert galaxy population. (2) In contrast, there is a significant trend that their more luminous counterparts namely CL and extremely variable quasars (CLQs and EVQs) are different: CLQs are generally in blue galaxies; in terms of the diagram of SFR and M* local CL Seyfert galaxies are located in the green valley, whereas CLQ hosts are in the star-forming main sequence. We propose explanations for those strongly variable Seyferts and quasars, respectively, under the thought that accretion disks broadly depend on nuclear fueling modes. Local large-variability and CL Seyferts are in nuclear famine mode, where cold-gas clumps can be formed stochastically in the fueling flow, and their episodic infall produces sharp peaks in the accretion-rate curve. CLQs and EVQs are in feast fueling mode, which may account for both their preference to blue galaxies and their variability pattern (high-amplitude tail of the continuous distribution). Lastly, we propose a new thinking: to search for IMBHs by optical variability in red galaxies.


INTRODUCTION
It has been known for about five decades that all active galactic nuclei (AGNs, including Seyfert galaxies and more luminous ones -quasars) are variable at some level in their optical continua and broad emission lines. Generally, the optical variability is of small amplitude, with the rms fluctuation typically being within 10-20% on timescales of weeks to years (e.g., Collier & Peterson 2001;Walsh et al. 2009), and behaves in a stochastic (or chaotic) manner (e.g., Kelly et al. 2009). For a long time, there were only a few exceptions discovered serendipitously, varying with much larger amplitudes (e.g., by a factor of 5; Alloin et al. 1985; see Antonucci 2018 for a historical commentary). Yet this field is being renewed, as a welcome to the so-called time-domain astronomy nowadays. Large or extreme flux changes on timescales of months to years have been observed in 100 AGNs, and the number is still increasing owing to the systematic searches based on large-scale photometric and spectroscopic surveys (e.g., Sheng et al. 2017;Rumbaugh et al. 2018;Yang et al. 2018;MacLeod et al. 2019;Graham et al. 2020). The change can be so dramatic that optical broad emission lines can disappear completely or (re)appear, i.e., the so-called optically changing-look (CL) AGNs (e.g., LaMassa et al. 2015).
It is fair to say that the studies on the large/extreme variability of either the AGN optical continuum (e.g., Rumbaugh et al. 2018) or broad emission lines (e.g., MacLeod et al. 2019;) are just developing. For instance, there are no unified definitions or criteria yet for the terms "large or extreme variability" and "optically changing look" in the literature. The optical CL AGNs (e.g., LaMassa et al. 2015; and see Yang et al. 2018) are often referred to as AGNs that experienced the spectral type transitions between types 1.8-2 (namely having weak or no broad Hβ, or even no broad Hα) and Type 1 (namely having both broad Hβ and Hα). However, the detection of broad-Hβ lines usually depends on the spectral signalto-noise ratio (S/N) level (see §3.3 of MacLeod et al. 2019; also §3.1 below). A consensus has emerges from the recent studies that those optical CL AGNs (including CL quasars) and extremely variable AGNs (mainly being quasars in the literature so far, usually shortened as EVQs) generally reflect variable accretion rate rather than variable obscuration (e.g., LaMassa et al. 2015;Sheng et al. 2017;MacLeod et al. 2019;Jiang & Blaes 2020;Graham et al. 2020). Yet the concrete mechanisms are not clear so far (see, e.g., Lawrence 2018;Dexter & Begelman 2019;Sniegowska et al. 2020;Jiang & Blaes 2020). There is a trend from the latest large-sample studies (Rumbaugh et al. 2018;MacLeod et al. 2019;) that EVQs and optical CL quasars appear to be the high-amplitude tail of the continuous distribution of quasar variability; this variability pattern is different from that of their less luminous counterparts, nearby CL Seyferts such as Fairall 9, Mrk 590 and Mrk 1018 (cf. §5.6 of . Those local Seyferts instead exhibit a secular pattern. For example, according to the available data, Mrk 1018 was a Seyfert 1.9 in the year of 1979 and became Type 1 in 1984; then, the AGN brightness remained constant roughly (as a Seyfert 1), with a small rms (0.15 dex) for at least 10 yr (i.e., 2001-2011 with observations available), and steadily decreased after 2011 (returning to Type 1.9 in 2015) by 3-4 mag until a minimum around October 2016, and again gradually got brightened at a low level (Krumpe et al. 2017;). It appears to keep the slow brightening so far, according to the latest observation reporting it as a faint Type 1.5 (Hutsemékers et al. 2020).
Naturally, one extends the exploration of CL AGNs to their host-galaxy properties. There have been several studies in the literature. For instance, two detailed studies on the nuclear cold gas (Koay et al. 2016) and nuclear warm gas (Raimundo et al. 2019) of the prototypal Mrk 590 reveal complex gas distributions and rich structures in the inner kiloparsec and 100 pc scales. Those nuclear features are normal in local (invariable) Seyfert galaxies and indicate secular processes for AGN fueling. There are also three studies on the host-galaxy properties of small samples. Charlton et al. (2019) analyzed the host galaxies of four faded CL quasars (z 0.2) using Gemini imaging and found that their hosts mostly have major-merger features and reside in the "green valley" between blue star-forming galaxies and dead red galaxies; i.e., the hosts of CL quasars are just "like the majority of [luminous] AGNs", as they concluded. Yu et al. (2020) analyzed in detail the host galaxies of five nearby CL AGNs (z 0.05) in the MaNGA survey based on their spectra, imaging, and IFU datacube. Among their results, there is an interesting finding: in terms of the diagram of global (namely whole-galaxy) star-formation rate (SFR) versus stellar mass (M ⋆ ), the host galaxies of their CL AGNs are located in the star-forming main sequence (SFMS), just like the hosts of their non-CL (i.e., invariable) broad-line AGNs. Most recently, Dodd et al. (2021) compiled a list of 17 CL AGNs from the literature and divided them into two groups: 11 at low redshifts z < 0.15 and 6 at 0.15 < z < 0.25. They presented a surprising result: in the global SFR versus M ⋆ diagram, their z > 0.15 CL AGN sample and their reference AGN sample (namely local Seyfert 2 galaxies) fall in the SFMS, whereas, in stark contrast, their z < 0.15 CL AGN sample is located in the green valley (see their Figure 1). Besides, they found that CL AGNs and highly variable AGNs have high galaxy Sérsic indices and high bulge-to-total light ratios, "implying high stellar density in their cores" (compared with their reference AGN catalog; see Figures 3 and 4 of Dodd et al. 2021). The above results by different teams are not consistent with each other. We believe one major reason is their small sample sizes, which result in their respective biases (e.g., the high fraction of LINERs in the z < 0.15 sample of Dodd et al. (2021), see §3.2 below; see the fourth conclusion in their §6). A second reason lies in the (unwary) use of ready catalogs of physical quantities derived by mass production; for example, SFR values differ significantly between different catalogs, with the difference even being formidable for certain populations (e.g, more than 1 dex systematically; see §8 of Salim et al. 2016), due to their respective systematic errors (e.g., an obvious yet not the most important one: aperture correction). This is particularly the case when a galaxy harbors an AGN and thus its SFR calculation has to account for the AGN contamination (i.e., subtracting the narrow emission lines of the AGN, dust emission heated by the AGN, and to be worse, direct AGN continuum and/or broad emission lines if the AGN is of Type 1).
In this work, we report (1) our discovery of six strongly variable (low-mass) AGNs 1 preferentially in red galaxies, and (2) as an extension of the discovery, our systematic investigation of the host galaxies of low-z CL AGNs and largevariability AGNs, particularly as to any connections between host-galaxy properties and the large/extreme-variability (including CL) phenomenon. 2 The broad-Hα fluxes of the six sources varied by 1.3 − 3.0 times (namely 0.3 − 1.2 mag) on a timescale of years. Their variability was found unexpectedly during our spectroscopic campaign initially planned to search for IMBHs in broad-line AGNs, using the Magellan Echellette Spectrograph (MagE; Marshall et al. 2008) mounted on the 6.5m Magellan Baade telescope. Most surprisingly, among our broad-line AGNs identified by MagE (15 broadline sources in total), those hosted by blue galaxies generally vary little, whereas a significant fraction of those hosted by red galaxies exhibit large broad-Hα variability, i.e., the six sources are predominantly in relatively red galaxies (even redder than the general population of low-z Seyfert 2 galaxies; see §3.3). Their large broad-Hα variability is confirmed by our ensuing spectroscopic observations. In order to avoid the possible bias caused by the small sample size of our broad-line AGNs by MagE, we collect all the low-z CL AGNs (z < 0.15) available in the literature and the strongly variable AGNs observed by reverberation mapping (RM AGNs in short) with available variability measures. Those data secured our conclusion: in terms of all diagnostic diagrams of color versus M ⋆ , stellar absorption-line index Hδ A versus 4000Å break (D 4000 ), and SFR versus M ⋆ , the large-variability and CL AGNs (mostly being Seyferts not LINERs) are predominantly in red galaxies and redder than the general population of low-z Seyfert 2 galaxies.
The structure of the present paper is as follows. In §2, we describe the data, data reductions and analyses for the six MagE AGNs, the low-z CL AGNs, and the strongly variable RM AGNs. In order to make the logical flow of the text friendly to the reader, we put two parts of the data analyses in Appendix A, which are specific to the six variable MagE AGNs and handle their multiband light curves and X-ray spectra. In §3, first we present the results: the genuineness of the large variability of the six MagE sources ( §3.1); the BPT and other AGN properties of the large-variability and CL AGNs ( §3.2); our discovery that low-z CL and largevariability Seyferts generally reside in redder, SFR-deficit host galaxies compared with the control sample of local Seyfert 2 galaxies, while CL quasars at relatively higher redshifts tend to be in the opposite, preferring to blue host galaxies that are located in the star forming main sequence ( §3.3). These discoveries inspired us to think about the idea that accretion disks broadly depend on nuclear fueling modes; in §3.4 we propose theoretical explanations for those strongly variable Seyferts and quasars, respectively. Lastly (in §3.5), we return to IMBH research, which was the initial goal of our observational campaign that branched out in the present study, by proposing a new thinking on opticalvariability selection for IMBHs, i.e., an implication of the discovery of this work. §4 is the summary.

Brief Description of the MagE Sources
The six largely variable (low-mass) AGNs are the following: SDSS J083909.65+072431.5 (hereafter J0839+0724), SDSS J111349.83+000733.9 (hereafter J1113+0007), SDSS J125710.76+272417.6 (hereafter J1257+2724), SDSS J134245.69+243524.0 (hereafter J1342+2435), SDSS J141249.63−030720.9 (hereafter J1412−0307), and SDSS J144242.63+011911.2 (hereafter J1442+0119). They are a part of sources observed by MagE in 2017 March and July, with the initial goal of identifying IMBHs in broad-line AGNs (namely low-mass AGNs). Throughout the whole MagE run of target selection and observing, AGN variability did not come to our mind. From the 2017 observing run, there are 15 broad-line AGNs confirmed in total. The host-galaxy colors g − i (measured from their Petrosian magnitudes) of the 15 sources range from 0.43 to 1.14, with eight sources being red with g − i > 1. Details of the MagE observation, data reduction, and all sources will be presented in another paper (W.-J. Liu et al. 2021, in preparation). In this paper hereafter, we will call the six largely variable (low-mass) AGNs the MagE sample.
J0839+0724 decreased in broad-Hα flux between the two epochs of the SDSS and MagE spectra, while the other five sources increased. We took follow-up spectroscopic observations for the five sources if possible, which confirmed their broad-Hα variability (see §3.1 for details).
Here, we define the "large or strong variability" of the six sources to be a broad-Hα flux change (the ratio of the maximum flux to the minimum flux) greater than 1.3 . This ratio is defined essentially in the same way as the variability measure R max commonly used in the RM literature (e.g., Peterson et al. 2004;Barth et al. 2015). Another common measure of the overall light-curve variability is the so-called fractional variability amplitude F var (see the Appendix §A.1 for details), which is designed to measure the intrinsic variability amplitude (i.e., correcting for the effects of measurement errors) and is thus more robust than R max to noises and outlier values (see §A.1; also Barth et al. 2015). Certainly, the two measures are significantly correlated. R max = 1.3 for Hα or Hβ (i.e., our above threshold value) roughly corresponds to F var = 0.07 (see Table 3 of Barth et al. 2015). According to Barth et al. (2015), F var > 0.1 means strong variability, roughly corresponding to R max 1.55 (see their Table 3). This number is higher than the broad-Hα flux changes of three sources in our MagE sample, which are as follows: J0839+0724 (1.30), J1342+2435 (1.41), and J1412−0307 (1.32). We keep the three as strongly variable sources because of the variabilities of their continua in other wavelength bands (see §3.1), particularly because of the behavior of J1412−0307. Although belonging to the two sources least variable in broad Hα, J1412−0307 has R max > 10 in soft X-ray, 1.75 (0.61mag) in the WISE W2-band light curve (binned every half years), and 1.51 (0.45 mag) even in the optical V -band light curve (binned every 90 days) without removing the dominant starlight component. Mentioned in passing, R max > 1.3 (i.e., 0.3 mag) is larger than the general variability amplitude of low-mass AGNs selected by optical continuum variability (e.g., Baldassare et al. 2018;Martínez-Palomera et al. 2020;Baldassare et al. 2020; ?).

Collected Low-z CL AGNs
Because our MagE sample is small and the emission-line variability of the sources is not extreme, we decide to extend our exploration to the sources in the extreme and dramatic end. Thus, we collect all the CL AGNs at z 0.15 reported in the literature so far. The redshift cut is set to ensure a better decomposition of the AGN and starlight components of the SDSS images (see also Dodd et al. 2021). We obtain 31 CL AGNs, constituting the low-z CL sample. All the sources have SDSS images, and 29 of them (except NGC 2617 and NGC 3516) have SDSS spectra. For the two CL AGNs, we find one optical spectrum publicly available in the 6dF survey (see also Shappee et al. 2014) for NGC 2617 and one nuclear optical spectrum in Kennicutt (1992) for NGC 3516, respectively.

RM AGNs with Large Broad-Hα Variability
We also collected RM AGNs with large/strong variability in Hα or Hβ emission lines, which should belong to the same population as the 6 variable MagE sources. The criterion is this: RM AGNs with measured F var 0.1 (Barth et al. 2015; see our §2.1.1) for Hα or Hβ available in the literature. We obtain 15 RM AGNs finally. They turn out to be at very low redshifts (0.002 < z < 0.04), and from three RM samples (Peterson et al. 2004;Bentz et al. 2009;Barth et al. 2015).

Observations and Data Reductions for the MagE Sample
Brief information on all the spectroscopic observations for the six variable sources in the aforementioned MagE sample is listed in Table 1. Below we describe the instrumental and observing matter of every observation and the corresponding data reduction.
All the six variable sources in the MagE sample have archival SDSS spectra, and J1113+0007 was spectroscopically observed twice. The SDSS is an imaging and spectroscopic survey, using a dedicated 2.5 m telescope to image one-quarter of the sky and to perform follow-up spectroscopic observations. Fibers that feed the SDSS spectrographs have an aperture of 3 ′′ diameter. The nominal total exposure time of the survey spectra is 45 minutes, which typically yields an S/N of 4.5 pixel −1 for sources with a g-band magnitude of 20.2. The spectra are flux-and wavelengthcalibrated by the SDSS pipeline, with 4096 pixels from 3800 to 9200 Å at a resolution R ≡ λ/∆λ ≈ 1800 (i.e., instrumental dispersion σ instr ≈ 70 km s −1 ). • MagE.
We conducted spectroscopic observation with MagE aboard the Magellan Baade telescope on UT 2017 March 25, using a 1 ′′ -slit, which gives an instrumental dispersion σ instr ≈ 26 km s −1 as measured from the arc lamp spectra. The spectral coverage is approximately 3200-10000Å across 15 echellette orders. One-dimensional spectra were extracted and wavelength-calibrated using MASE reduction pipeline (Bochanski et al. 2009). The telluric H 2 O absorption features longward of 6800Å were removed by dividing with the normalized spectrum of a standard star. More details of the data reduction can be found in W. -J Liu et al. (2021, in preparation). • Follow-up by P200/DBSP and Xinglong/BFOSC.
To further check and monitor the broad-Hα variability for the six sources, we carried out several follow-up optical spectroscopic observations, as summarized below. For J1257+2724, J1412−0307, and J1442+0119, spectra were taken using the Double Spectrograph (DBSP) on the Palomar 200 inch Hale telescope (P200) at Palomar Observatory on UT 2019 June 24-25. The instrument is configured with a D55 dichroic, 1. ′′ 5 slit, 600/4900 grating for the blue side and 600/7600 grating for the red side, which gives a spectral resolution of σ instr ≈ 108 km s −1 for the blue side, and σ instr ≈ 70 km s −1 for the red side (i.e., almost the same as the SDSS resolution), and provides a wavelength coverage of 3400-9300Å. The three sources were observed with two separate 500-900 s exposures.
For J1257+2724, we took a spectrum with the Beijing Faint Object Spectrograph and Camera (BFOSC) mounted on the 2.16m telescope at Xinglong Observatory on UT 2020 January 1. The grating G4 and a slit of width 2. ′′ 3 were used. This setup gives a spectral resolution of σ instr ≈ 230 km s −1 , as measured from the sky emission lines and comparison arcs, and provides a wavelength coverage of 3850-8000Å. Two separate 1800s exposures were taken.
The two-dimensional spectral data of the DBSP and BFOSC observations were reduced with the standard routines for longslit spectra in IRAF. The telluric H 2 O absorption features longward of 6800Å were removed by dividing with the normalized spectrum of a standard star.

Spectral fitting and analysis
For the six variable MagE sources and all the 31 lowz CL AGNs, as well as several sources as comparison (e.g., the four CL quasars of Charlton et al. 2019), we fit their optical spectra, in order to measure their host-galaxy properties (such as D 4000 and Hδ A ) and emission lines. Based on the measurements of emission lines, we in turn can derive SFR (from [O II] λ3727) and AGN parameters such as M BH and L/L Edd . The variable MagE sources have multiepoch spectra as described in §2.2.1. Among the CL sources, 29 use the SDSS spectra, NGC 2617 uses the 6dF spectrum, and NGC 3516 uses the spectrum observed by Steward Observatory 2.3 m telescope (see §2.1.2) Regarding the 15 strongly variable RM sources, because their publicly available optical spectra are heterogeneous and incomplete, and generally dominated by the AGN emission, we give up to fit them. We do not need to use the spectrumbased quantities for their host-galaxy properties and instead use the quantities based on imaging decomposition (such as color, magnitude, and derived stellar mass, after correcting for AGN contamination; see Table 4 and §3.3.1).
Regarding the comparison samples for our low-z CL or large-variability AGNs, such as the four CL quasars at z 0.2 in Charlton et al. (2019), the Seyfert 2 sample selected by Dong et al. (2010), and the non-AGN galaxy sample selected by Dong et al. (2012), the aforementioned spectral properties are based on their SDSS spectra; we either fit their SDSS spectra by ourselves (for the CL quasars) or obtain the fitting results from the corresponding researchers (for the latter two control samples).
Prior to fitting, all spectra were de-redshifted to the vacuum rest-frame wavelength and corrected for Galactic extinction using the extinction map of Schlegel et al. (1998) and the reddening curve of Fitzpatrick (1999). We use the same procedures as Dong et al. (2012) for the continuum modeling and emission-line profile fitting; we only provide a brief description here. We model the continuum of each of the spectra using the starlight templates and a single reddened power law to represent the AGN continuum in the spectral fitting. For SDSS, DBSP, and BFOSC spectra, we use the starlight templates from Lu et al. (2006), built from the spectra of simple stellar populations of Bruzual & Charlot (2003). The spectral resolution of the starlight templates used is 3 Å across the wavelength range of 3200-9500Å, corresponding to a median resolving power λ/∆λ ≈ 2000, which is comparable with the spectral resolution of SDSS spectra (69 km s −1 ). The DBSP spectra have identical spectral resolution to the SDSS spectra, and the BFOSC spectra have a much lower one. The resolution of MagE spectra is about twice as good as that of SDSS spectra; a starlight template built from the MagE spectra of four types of stars (B0, A5, K0III, and M4III), which were observed by MagE on the same day with the six MagE AGNs, is adopted in continuum fitting. More detailed information on this will be provided in W.-J. Liu et al.(2021, in preparation). The starlight templates are broadened by convolving with a Gaussian of different widths and shifted to match the stellar velocity dispersion, so that the stellar absorption lines could be well subtracted. A χ 2 minimization is performed iteratively over the whole spectroscopic wavelength range, except for the regions emission lines (e.g., Hα, Hβ, [O III]λ5007, [S II]λλ6716, 6731, [O I]λ6300, and [O II] λ3727). In the fitting, we treat each spectrum of the same AGN independently to find the best-fit model for each spectrum.
Mentioned in passing, we also measured the stellar velocity dispersions from the MagE spectra for the six MagE sources along with the entire MagE sample. Details are in W.-J. Liu et al. (2021, in preparation). Briefly, the fitting method is almost the same as described in Xiao et al. (2011), in which the spectra of velocity template stars are broadened and fit to the galaxy spectra locally in a specific spectral region, e.g., Mg I b (5040-5430Å), or Ca II triplet (8450-8750Å); the spectra of velocity template stars were also observed by MagE on the same day as the MagE AGNs. We do not use the σ ⋆ data in this work and only list them in Table 2 for any possible interest of the reader.
We then fit the continuum-subtracted spectrum ( • Calibrating and scaling multiepoch spectra (MagE sources).
Due to the difference in adopted apertures/slit widths, seeings, and flux calibrations of the spectra taken by different instruments, their nuclei-emission-line flux levels are not exactly the same. Therefore, the emission-line-flux calibrations between spectra of different epochs are needed, and the [O I]λ6300 is the ideal forbidden line to do this.
[O I]λ6300 is mainly produced in the "partially ionized zone" related primarily to AGN radiation, so it is little affected by star formation in galaxies. behavior is also seen in an "Turn-on" AGN SDSS J1115+0544 (Yan et al. 2019), which is explained by the increase of the ionization continuum. The fitting results also indicate that the continua of the six AGNs are dominated by starlight. Their AGN continuum fluxes at 5100Å account for only 5%-22% of the total continuum fluxes.
We also calculate the AGN spectral types for every spectra of the six MagE AGNs. We adopt the definition of spectral subtypes introduced by Winkler (1992, see also Osterbrock 1981 andVéron-Cetty &Véron 2001), which uses a combination of the presence or absence of broad Hα and broad Hβ, and the flux ratio of total Hβ to [O III]. We adopt S/N > 3 as the detection threshold of broad emission lines. The final emission-line parameters of every spectra, as well as the spectral subtypes, are listed in Table 1.
• Calculating spectrum-based quantities for host galaxies.
To quantify the properties of host galaxies, we calculate three spectrum-based quantities: the 4000 Å break D 4000 (Balogh et al. 1999), the Balmer absorption-line index Hδ A (Worthey & Ottaviani 1997), and the SFR, from their SDSS spectra 3 for the six MagE sources and for the low-z CL sources except two. The two exceptional sources have too strong AGN emission in the optical, with the AGN fraction > 30% in the r-band images (see Table 3) and larger in the 3 ′′ aperture SDSS spectra. As for the RM sources, only a fraction have spectra in SDSS, and their SDSS fiber spectra are generally dominated by AGNs (because those sources are very nearby). Thus we do not calculate or use the three quantities for the RM sources in this work. Besides, we also calculate the 3 quantities based on their SDSS spectra for the four CL quasars of Charlton et al. (2019).
D 4000 is a good indicator of stellar age, and its value is small for young stellar populations and is large for old, metal-rich galaxies. It increases monotonically as the luminosity-weighted mean age of the stellar population increases (Kauffmann et al. 2003a). Hδ A is defined to measure the strength of stellar Hδ absorption line and indicates the burstiness of the star forming activity. Strong stellar Hδ absorption (positive large Hδ A values) means that the galaxies experienced a burst of star formation that ended within 0.1-1 Gyr ago. The two indices used in combination can give the mean stellar age and tell the star formation is bursty or secular. We calculate the two indices from the stellar continua decomposed from the SDSS spectra of the MagE and CL samples and listed them in Tables 2 and 3. The SFRs are calculated according to the measured [O II] λ3727 flux from their SDSS spectra using the calibration for the 3 NGC 2617, a source in the low-z CL sample, has no SDSS spectrum, thus we actually use its available 6dF spectrum publicly. Because [O II] λ3727 is not covered by the 6dF spectrum, we cannot give its SFR in Table 3. NGC 3516 is another low-z CL AGN that has no SDSS spectrum and image. Koss et al. (2011) observed its u, g, r, i, z image using the Kitt Peak 2.1 m telescope. In Table 3, we directly use the galactic g, r, i magnitudes, rband PSF fraction, and stellar mass given by Koss et al. (2011), and use the bulge and disk Sérsic indexes decomposed from the Hubble Space Telescope (HST) F814W image by Kim et al. (2017). The D 4000 and Hδ A are measured from its nuclear optical spectrum observed using the Steward Observatory 2.3 m telescope during 1989-1991 (Kennicutt 1992). Because the flux units of the spectrum are counts, so we cannot give its SFR in Table 3. AGN hosts of Zhuang & Ho (2019), which removes the [O II] λ3727 flux from the AGN and accounts for the influence of metal abundance on the SFR. Note that SDSS and BOSS spectra are taken within a fiber aperture of 3 ′′ or 2 ′′ diameter, while the angular sizes of most of the above low-z AGN hosts (see Tables 2 and 3 for their r-band Petrosian radii) are much larger. Thus, the three spectrum-based quantities mainly reflect the host-galaxy properties in the inner regions.

Image Fitting and Analysis
All of the sources in the MagE and low-z CL samples, as well as 10 sources in the RM sample, have archival SDSS images. For those sources except two, we fit their SDSS images in order to remove the AGN emission and measure the host-galaxy properties. The two exceptions are the following: J111803.22+450646.8 in the low-z CL sample, because a bright foreground star sits in its central region; and NGC 4051 in the RM sample, because it is very nearby and the AGN is bright, and thus there are too many saturated pixels in the central region of its SDSS images. As for the five rest RM sources without SDSS images, as well as NGC 4051, fortunately Koss et al. (2011) have taken deeper optical images and performed two-dimensional fitting; we simply adopt the host-galaxy magnitudes and colors (with the AGN contamination removed) as well as the AGN fraction in the r band there.
The SDSS imaging data were collected in a drift-scan mode in five bandpasses (u, g, r, i, and z) on nights of pristine conditions, with a typical seeing of 1. ′′ 5 in the r band. The images have a total exposure time of 54 s per filter. The photometric calibration is accurate to 5%, 3%, 3%, 3%, and 5%, respectively. Figure 2 displays the g-r-i composite images for all the 6 MagE sources (panel a), 30 low-z CL AGNs that have SDSS images (panel b), and the aforementioned 10 sources in the RM sample (panel c). Almost all of the host galaxies are well resolved, except several objects in the low-z CL sample (roughly with Petrosian radius < 4 ′′ in the r rand; see Table 2,3,4): 4 : J0909+4747, J1003+3525, J1132+0357, J1358+4934, J1428+1723, and J1545+2511. It is also obvious from Figure 2 that in most of the images the AGN emission does not swallow up the host galaxies.
We perform a two-dimensional (2D) decomposition of AGNs and host galaxies into their SDSS images in the g, r and i bands of every aforementioned source (44 sources in total) using GALFIT (Peng et al. 2002(Peng et al. , 2010. Prior to the fitting of every image, all foreground or background photometric objects identified by the SDSS pipeline are carefully masked out manually. In the fitting, the sky background is set to be free. The AGNs are represented by a PSF component. The PSF model images are reconstructed from the psfFiled files provided by SDSS for every given field. Concerning the host galaxies, for most sources the galactic morphologies are basically regular (axisymmetric), and we adopt the traditional fitting approach implemented in GALFIT. Yet for about a dozen sources, their images are rather nonaxisymmetric, e.g., with the starlight being dominated by grand-design spirals or bars (e.g., Mrk 1511, J1533+4432), with merger features and so on, or have prominent rings, and we make use of the new machinery implemented in the "new-generation" (NG) GALFIT (namely version 3; Peng et al. 2010); hereafter, we call this approach "NG fitting" because the new methodology is somehow antitraditional in certain respects, as can be sensed in our description below.
In traditional fittings of the regular galaxies, we model the host galaxies with the following three models (with an increasing number of free parameters): a single Sérsic r 1/n function (Sérsic 1968), or one Sérsic plus one exponential, or two Sérsic functions. Here the Sérsic and exponential functions represent the azimuthally averaged, radial SB profiles, which describe the intensity falloff away from the peak (assuming the isophotes being ellipses in shape, i.e., azimuthal symmetry). The exponential function is the special case of a Sérsic function with the index fixed to be n = 1, commonly used to represent galactic disks; yet nowadays researchers realize that the surface-brightness (SB) profiles of (the outer) disks are not necessarily exponential, and Sérsic functions with n < 1 are also common (see also Peng et al. 2010). Following the common practice of this kind of imaging fitting, we begin the fitting of every image with the simplest scheme ("PSF + Sérsic") allowing the scaling factor of the PSF (AGN) and all parameters of the Sérsic (host galaxy) to vary. Then, we try schemes with the other two advanced models for host galaxies, if the fit can be improved significantly in terms of χ 2 (in the spirit of the F-test) and the analysis of residual images; see Dong et al. (2007) for details. In the case of two-component models for the host galaxies, following the common practice hereinafter we interpret the inner Sérsic as the bulge or pseudo-bulge, and the outer one as the disk.
In NG fittings, any galactic component is still represented by the same basic SB profile models used in the traditional fittings, but their azimuthal shapes can be modified from the fundamental ellipse shape by four novel devices: bending, Fourier, coordinate rotation, and truncation modes. Those devices, or called shape operators in the mathematical physics fashion, are carefully devised so that the traditional profile parameters almost keep their original, intuitive meaning intact; in fact, they are merely higher-order modifying functions of the fundamental shape from the perspective of mathematics. In terms of the practical effects, however, those shape operators are magic (particularly when used in combination): for instance, they can metamorphose a Sérsic component into almost any arbitrary shape, thus representing a galactic component/substructure of any realistic morphology (e.g., a spiral with any number of arms, a spiral with a bar, or a tidal tail, let alone a bulge, bar, or their subcomponents; see §7 of Peng et al. 2010). This power naturally results in a new methodology, in contrast to the traditional wisdom of decomposing galactic images. Specifically, for the fitting of the galactic disks we concerned with in this work, the traditional wisdom is this: use a global Sérsic or exponential model to represent the (presumably) main body of the disk, then add localized subcomponents onto it (e.g., a bar). This tradition is actually not a choice, limited by previously available tools (see Peng et al. 2010). And worse, it has reinforced the misimpression that there should be a dominant disk component as the main body and the other (sub)components such as spirals and bars be secondary, which certainly is not true in many galaxies (see §7 of Peng et al. 2010). On the contrary, the NG fitting may divide a galactic "disk" into several (somehow localized) pieces, without a backbone disk component; e.g., the galactic disk of NGC 289 is decomposed into three spiral components each modeled with a morphed Sérsic each by Peng et al. 2010, see their Figure 22; i.e., the sum of the three spiral components is the commonly called "disk" component. Following the advice of Peng et al. (2010), in order to get accurate bulge/disk decomposition, particularly the accurate index values (n) for the inner Sérsic components (see the thorough analyses and remarks in §7 of Peng et al. 2010), for about a dozen galaxies with complex morphologies (such as NGC 2611, Mrk 1511, J1533+4432 and J1605+4526 in the CL and RM samples, as well as two mild ones in the MagE sample -J1113+0007 and J1343+2435; see their demonstrations in Figure 3), we adopt the NG fitting approach. The fitting procedure is the same as described in detail in §7 of Peng et al. (2010). All the inner components (roughly corresponding to bulges or pseudo-bulges) of the galaxies can be well modeled with a Sérsic model (morphed or not), and the outer components (the commonly called "galactic disks") usually need two morphed Sérsic models. To the purpose of this work, we do not list all the fitting results here. Instead, because for their outer components we only need the global/averaged Sérsic indexes (see Tables 2  and 3), again we run a second fitting for every one of those galaxies with two or more Sérsic models for the outer components in their aforementioned NG fittings. In the second fit of an image, we fix the inner component to the best fit of its first NG fitting and use just one (morphed) Sérsic model to fit the outer. Then we adopt the best-fit Sérsic indexes of the second fittings as the final ones for the outer components.
The final g, r and i magnitudes of the host-galaxy starlight (i.e., with AGN light removed), as well the Sérsic indexes of the inner and outer (if present) components, are summarized in Tables 2 and 3. The listed magnitudes have been corrected for Galactic extinction using the dust map of Schlegel et al. (1998) and the reddening curve of Fitzpatrick (1999); we list the magnitudes without k-corrections on purpose to avoid the uncertainty caused by k-correction. The global (namely whole-galaxy) g − i colors of the host galaxies are also listed in Table 2, with the k-corrections applied with the public code kcor provided by Chilingarian et al. (2010). Six galaxies in the CL sample are of too small size (Petrosian radius < 4 ′′ in the r band) to get reliable indexes for their inner Sérsic components (e.g., J1358+4934 particularly, which is small, and to be worse has a bright AGN and a relatively high redshift); for them we mark their inner Sérsic indexes with a superscript ' ? ' in the table. For several additional galaxies that are of relatively small size or somehow with a bright AGN, or at a little bit high redshift, it is not certain whether their outer Sérsic component is present or not, and we denote it with a '?' symbol in the Outer Sérsic Index column. If our best-fit model for a galaxy does not require an outer Sérsic component, we denote with a '-' symbol in that column. Mentioned in passing, J0915+4814 (in the CL sample) has a very small outer Sérsic index of ≈0.1, which actually indicates that there exists an outer ring with a radius of ≈8 ′′ .
Regarding the M ⋆ estimation for the host galaxies, our strategy is as follows. According to Tables 2 and 3, all the 6 MagE sources and most of the low-z CL sources have a small fraction of AGN emission to the total (AGN + starlight) light (< 10% in the r band). Thus, we uniformly estimate M ⋆ for the sources of the two samples based on their NIR magnitudes and do not perform AGN-starlight decomposition into their NIR images. 5 We use the massto-light ratio in the K band (M/L K ), which is relatively insensitive to either dust absorption or stellar population age. The M/L K formula was calibrated by Into & Portinari (2013) and involved the K S magnitude and g − i color, and log M/L Ks = 0.794(g − i) − 0.997 (with a scatter of ±0.1 dex; see their Table 3). All of the above sources were detected in the Two Micron All-Sky Survey (2MASS, Skrutskie et al. 2006), and most are extended sources. we prefer the photometric data in the 2MASS Extended Source Catalog (XSC), for the emission from the whole galaxies; as for some sources not in XSC, we then use their 2MASS Point Source Catalog (PSC) data. Galactic extinction corrections and the k-corrections were performed on the g, i, and K s magnitudes for the M ⋆ estimation. As for the sources in the RM sample, we just use the M ⋆ data available in the two references-Bentz & Manne-Nicholas (2018) and Koss et al. (2011), both with AGN contamination removed; see Table 4 (except one source, IC 4218, for which no ready M ⋆ is available, and we estimate it in the same way as for the MagE and low-z CL sources). For sources listed in Bentz & Manne-Nicholas (2018), we adopt their M ⋆ values estimated using the Into & Portinari (2013) formula, consistent with our above estimation. If a source has M ⋆ in both references, we prefer Bentz & Manne-Nicholas (2018), because their AGN-starlight decomposition was based on NIR images observed with HST. 5 For several sources in the low-z CL sample, the AGN fraction in the r band exceeds 10%, thus our estimation based on NIR magnitudes would overestimate the M⋆ values to some degree (albeit less significantly than in the r band). But this overestimation only weakens our related conclusion in §3.3.

RESULTS AND DISCUSSION 3.1. Large Variability of MagE Sources
• Broad-Hα and continuum variability.
The broad-Hα fluxes of the six variable MagE sources change by 1.3 − 3.0 times (namely 0.3 − 1.2 mag) on a timescale of years, exceeding the common variability amplitudes observed in reverberation mapping experiments (Kelly et al. 2009;Walsh et al. 2009). First of all, we need to verify that the broad-Hα changes of the six sources are real, not due to instrumental effect. Because the spectral resolution of our MagE observations is more than twice that of the SDSS (as well as better observing condition), there are possibilities that the broad-Hα fluxes of the followup MagE spectra are different from the SDSS ones, either increasing (e.g., recovering the broad-Hα wing to a deeper flux density level) or decreasing (e.g., extracting a finer or even recognizing a false profile). Thus we conducted ensuing spectroscopic observations instead with similar (or somehow inferior) spectral resolution and observing conditions to the SDSS, for as many sources as possible by P200/DBSP and Xinglong 2.16m/BFOSC (see §2.2.1). For all the four MagE sources with ensuing lower-resolution spectroscopy, their DBSP or BFOSC observations confirmed the presence of broad-Hα lines, with S/N > 10 (see Table 1); moreover, the measured broad-Hα fluxes are still higher than their SDSS ones, consistent with the increase of their MagE fluxes over the SDSS. Hence their broad-Hα flux changes are secured. For the remaining two sources (with RA smaller than 12 hr), we did not find opportunities to take decent spectra with spectral resolution smaller than the MagE ones. Fortunately for J1113+0007, there are two SDSS spectra taken with a time separation of about one year; the later one gives a larger broad-Hα flux larger by 1.6 times than the former one (albeit with a 2σ significance only, according to the flux difference and errors listed in Table 1); this increasing trend is consistent with the MagE observation taken six additional years later.
The genuineness of their broad-Hα flux changes, as we mentioned in the definition of this large-variability MagE sample ( §2.1.1), is supported by their continuum variability in the optical, MIR, and/or X-ray. Actually, the optical and MIR light curves displayed in Figure A1 are seriously diluted by host-galaxy starlight, and any dramatic variability is not expected. But we are still able to observe significant variability of several sources. First, as listed in Table A1, all of the light curves by ASAS-SN and WISE for the six MagE AGNs satisfy σ 2 rms − err(σ 2 rms ) > 0, indicative of continuum variability (see §A.1). Furthermore, in terms of the threshold for "strong variability" by Barth et al. (2015), F var > 0.1 (see also §2.1.1), three sources (J0839+0724, J1113+0007 and J1412−0307) satisfy it according to their ASAS-SN light curves; this is remarkable because those optical light curves are actually dominated by starlight. Third, according to the peak-to-valley changes measured from the binned light curves (corresponding to R max in Table A1), there are considerable optical or MIR continuum variabilities as follows (only listing those changing by > 0.2 mag peak to valley): J0839+0724, 0.3 mag in the W 2 band; J1113+0007, 0.4 mag in the W 2 band; and the most dramatic one, J1412−0307, which is in fact the delimiter of this MagE sample (i.e., its broad-Hα R max ≈ 1.3), 0.45 mag in the ASAS-SN V band, and 0.61 mag in the W 2 band (see also §2.1.1). Such a large change in the WISE W 2 band reaches the same variability degree as the CL AGNs in the literature (e.g., Sheng et al. 2017). Lastly, among the three sources having two-or multiepoch X-ray observations, the two-epoch fluxes of J1442+0119 change by 2.74 times, and J1257+2724 and J1412−0307 exhibit peak-to-valley changes by > 20 times. We should note that strong/extreme variability in the X-ray usually has different origins from the strong variability or CL events in the optical of the same AGNs, e.g., due to variable obscuration by dust-free gas (see, e.g., Risaliti et al. 2009); even in the cases without variable X-ray obscuration, the connection may be rather complicated between X-ray and optical (CL) variabilities (see, e.g., Ricci et al. 2021). Thus, the observational support from the strong variability in the X-ray of the three MagE sources is not conclusive. But for the three sources at least their X-ray spectra exhibit little gas absorption (see below).
The variability of the MagE sources is not likely due to variable dust obscuration; this is supported by several lines of evidence. The direct evidence comes from the constancy of the broad-line Balmer decrements (Hα b /Hβ b ) of their optical spectra taken at different epochs. As Dong et al. (2008) demonstrated, the intrinsic (unreddened) ratios of Balmer recombination lines from the BLR have a very small scatter in the normal radio-quiet AGN population, and thus the observed broad-line Balmer decrements can be used as an indicator of dust obscuration in the broad emission lines. If the dust obscuration of a source is time variable, then the measured Balmer decrements are variable accordingly; equivalently, the contrapositive holds. J1257+2724, J1342+2435, and J1442+0119 have reliable broad-Hβ measurements in their multiepoch spectra (see Table 1). For J1257+2724, the Hα b /Hβ b ratios of the four spectra (ordered by increasing observing Date) are 3.58±0.63, 3.52±0.18, 3.71±0.17, and 3.79±0.34, respectively. For J1342+2435, the ratios of its three spectra are 3.40±0.51, 3.42±0.36, and 2.93±0.34. For J1442+0119, the ratios are 6.23±1.79, 5.88±0.44, and 6.78±1.01. All the three sources have insignificant variation in Hα b /Hβ b , well within 1σ uncertainty, during those years. Besides, for every object, the consistency of the M BH values estimated from the multiepoch spectra is also against the scenario of variable obscuration (see §3.2 below).
Another line of evidence, for the three sources with Xray observations, comes from their X-ray spectral fitting: the model with free intrinsic absorption is not favored by the data of any spectra (either in high-or low-flux states), and even if the free intrinsic absorption is added into the model, the bestfit N H has a very small impact on the fitting (see §A.2 for the detail).
It is worth noting that the X-ray spectral shape of J1257+2724 is soft (Γ = 2.3) in the high-flux state and gets hard (Γ = 1.1) in the low-flux state. This behavior may be similar to the soft-to-hard state transition at L/L Edd ≈ a few ×10 −2 found in Galactic X-ray binaries (Maccarone et al. 2003); if so, the variability may be due to the change of accretion state. • Spectral type transitions.
As stated in §2.2.2 (see also Table 1), we classify every spectra of the 6 sources into spectral subtypes according to the traditional definition (Osterbrock 1981;Winkler 1992;Véron-Cetty & Véron 2001). The spectral type transitions of the six sources are not as "dramatic" as prototypical CL AGNs; i.e., not between type 1 and type 1.8-2 with the dramatic appearance or disappearance of a strong Hβ component (LaMassa et al. 2015). J0839+0724 (type 1.9), J1113+0007 (type 1.9 or 2), J1342+2435 (type 1.2) and J1442+0119 (type 1.5) do not show significant type transitions. J1257+2724 varied between type 1.5 and 1.2 . J1412−0307 transited from type 1.9 to 1.5 during a timescale of 15 years, and then returned to type 1.9 recently. In its type 1.5 state, the broad-Hβ component is measured with S/N > 10. As described above, J1412−0307 also showed significant continuum variabilities in the soft X-ray, optical, and MIR; in particular, the variability measures of its WISE W 2 light curve are similar to those of typical CL AGNs (see, e.g., Sheng et al. 2017). Hence we regard it as a bona fide CL AGN. 6

AGN Properties
For the six variable MagE AGNs ( §2.1.1) and the 31 low-z CL AGNs ( §2.1.2), we begin with exploring their distribution in the diagnostic diagrams of narrow-line ratios (Figure 4), the so-called BPT diagrams (Baldwin et al. 1981;Kewley et al. 2001Kewley et al. , 2006Kauffmann et al. 2003c), which are a powerful tool to separate Seyfert galaxies, low-ionization nuclear emission-line region sources (LINERs; Heckman 1980), and H II galaxies. Following Kauffmann et al. (2003c) and J1412−0307 are Seyferts; J1342+2435 is a LINER; concerning J1442+0119, its two observations are in the LINER region and the third observation in the Seyfert region, and thus we deem it a 50% probability LINER and 50% probability Seyfert. Thus, among the MagE sources the LINER fraction (LINER/total) is 1.5/6 (25.0%). In the 31 sources of the low-z CL sample, 22 are pure Seyferts, 1 transition object, 7 pure LINERs, and the remaining 1 sits just on the boundary line between Seyferts and LINERs; thus, the LINER fraction is 7.5/31 (24.2%).
All of the 11 low-z CL AGNs (z < 0.15) of Dodd et al. (2021) are included in our low-z CL sample. Their BPT classifications are as follows: four Seyferts, one transition objects, and six LINERs. Thus, the LINER fraction of the low-z CL sample of Dodd et al. (2021) is 6/11 (54.5%), much higher than our low-z CL sample.
We also plot the four CL quasars of Charlton et al. (2019) on the BPT diagrams. It turns out to be two Seyferts and two transition objects.
The sources in our low-z CL sample have M BH estimates in the literature (see the M BH values and references in Table 3), and there are estimations and discussions of their accretion rates (or L/L Edd ) and accretion-rate changes (or state changes) also in the literature. Regarding the sources in the RM sample, there are RM-based M BH measurements and corresponding L/L Edd estimates (see the references in Table 4). Thus here we only describe the estimation for the MagE sources.
We follow the common practice to estimate the virial masses based on the measured luminosity and line width of the broad-Hα emission line. The mass formalism is given by Xiao et al. (2011, their Equation (6)), which is based on Greene & Ho (2005) but incorporates the updated relation between broad-line region size and AGN luminosity calibrated by Bentz et al. (2009). The M BH values of the six MagE AGNs are between 10 6 − 10 7.3 M ⊙ (see Table 1). We find that for any source the derived M BH values from the spectra of different epochs are consistent with each other, with difference 0.3 dex, well within the measurement uncertainty of the viral mass estimation method.
Besides, for the three sources with X-ray measurements ( §A.2), J1257+2724, J1412−0307, and J1442+0119, we also use their L X (2-10 keV) to estimate L bol and L/L Edd . The bolometric correction relation of Marconi et al. (2004) is used, which are ≈ 10 for the L X (2-10 keV) of the three sources. We use the mean value of different-epoch M BH values of each source to calculate their Eddington luminosity (L Edd ). For J1257+2724, the L/L Edd of the high-and lowflux states is 0.023 and 0.012, respectively; for J1412−0307, it is 0.016 and 0.0029. For J1442+0119, there is only one measurement, and its L/L Edd is 0.0072. Those X-ray based L/L Edd estimates are quite consistent with the broad-Hα based ones.
The variable MagE sources stay in the relatively lowaccretion regime, with 0.01 L/L Edd 0.04, i.e., definitely not in or close to the Eddington accretion. We note that a similar low-accretion trend was found in EVQs and CL quasars. Rumbaugh et al. (2018) found that EVQs with ∆g > 1 mag have lower L/L Edd than the control sample of quasars matched in redshift and luminosity. MacLeod et al. (2019) reported that CL quasars with a large optical variability (|∆g| > 1 mag and |∆r| > 0.5 mag) have lower L/L Edd compared with the overall quasar population. Such a trend of large-and extreme-variability AGNs in a relatively low-accretion state may shed light on the underlying accretion-flow physics.

Hosts Redder than General Seyfert 2 Galaxies
Broadly viewed in Figure 2, almost all of the host galaxies of the large-variability AGNs in the three samples are in secular evolution, without ongoing major-merger activity (except two in the RM sample). It is easy to judge that they are not luminous elliptical galaxies; many of them have more or less signatures of galactic disks, such as bars and spirals, while some may be (low-luminosity) S0 or spheroidal galaxies. These characteristics are consistent with the general trend of local AGNs and also the general trend of local galaxies (Kormendy & Ho 2013). Beyond the 0thorder general similarity, however, there appears a perceivable discrepancy from our general impression of the host galaxies of local AGNs: the images in Figure 2 are seldom blue, but slant seriously toward red; this is particularly true when looking at their inner regions.
Based on the derived quantities in §2.3 (magnitudes, colors, inner and/or outer Sérsic indexes, M ⋆ , D 4000 , Hδ A and SFR; corrected for AGN contamination), below we investigate the host-galaxy properties and explore any possible clues of these quantities against or for the large variability of the central AGN. We first utilize the diagnostic tools based on the global (namely whole-galaxy) quantities ( §3.3.1), and then invoke the tools based on the quantities of the inner regions (i.e., quantities derived from the fiber-aperture spectra, properties of the inner Sérsic imaging components, etc.; see §3.3.2).

In Terms of Whole-galaxy Properties
A traditional tool to diagnose galaxies is the diagram of global color versus M ⋆ (e.g., Kauffmann et al. 2003b;Schawinski et al. 2007). Here we use rest-frame g − i to represent the color, just following Charlton et al. (2019) in their study of the host galaxies of CL quasars; this is because the SDSS g, r and i images have similar quality (better than u and z) and the wavelength span between g and i is the largest among all combinations. To represent the general color−M ⋆ parameter space populated by local Seyfert hosts, we use the ∼26000 low-z Seyfert 2 galaxies homogeneously selected from SDSS DR4 by Dong et al. (2010). The g − i values of those Seyfert 2 galaxies are calculated from the SDSS g and i Petrosian magnitudes, accounting for both Galactic extinction and k-correction. The stellar masses of the Seyfert 2 galaxies are basically derived in the same way as described in §2.3. A fraction of Seyfert2 galaxies have no 2MASS K S data, and we then use their WISE W 1 3.4µm luminosity to estimate M ⋆ , with the mass-to-light ratio, log M ⋆ /νL ν (3.4µm) calibrated by Wen et al. (2013).
In addition, we retrieve the data of the 4 CL quasars at z 0.2 (Charlton et al. 2019), and the 52 Seyfert 1 hosts at 0.145 < z < 0.3 in the SDSS Stripe 82 region provided by Bettoni et al. (2015) for comparison. Figure 5 shows the distributions on the g − i versus M ⋆ diagram. Seyfert 2 galaxies, being the reference sample, are plotted as background (gray dots). We plot the contours of the Seyfert 2 distribution, representing 90% (the outermost), 50%, and 20%, respectively, the number fraction of enclosed sources. We also group the Seyfert 2s into nine M ⋆ bins, and plot the median g − i colors of every bins and corresponding standard deviation (black solid circles with ±1σ error bars).
The majority of the large-variability AGNs in any samples are redder than the median-color line of general Seyfert 2s (black dashed line): 3 of the 5 MagE non-CL variable sources (60%), 22 of the 31 z < 0.15 CL AGNs (including 30 z < 0.15 CL AGNs in Table 3 and the MagE CL AGN J1412−0307; 73.3%) 7 , and 10 of the 15 RM sources (66.7%). We perform 2-dimensional Kolmogorov-Smirnov (KS) tests between the variable-AGN samples and the reference Seyfert 2 sample. The p-value (change probability) between all the three variable-AGN samples (51 sources in total) and the Seyfert 2 sample is 0.0021, that between the the above 31 CL sources and the Seyfert 2, 0.0035 . That is, there is a statistically significant difference in the color-M ⋆ diagram between large-variability AGN hosts and general Seyfert 2 galaxies.
Regarding the galactic structure of these host galaxies, according to our 2D imaging fittings it is clear that they are predominantly disk galaxies with a bulge or pseudobulge. Among the 24 galaxies with both reliable inner Sérsic and outer Sérsic indexes in the MagE and CL samples (see Tables 2 and 3), there are merely about four or even fewer spheroidal galaxies. 8 In light of the limited resolution and depth of the SDSS images, as well as the sample incomplete-7 Here we formally take the single MagE CL AGN from the category of MagE large-variability (non-CL) sources and bring it into the CL category. But just as we point out in footnote 2, we actually need not to elaborate such strict separation between the two categories in this study concerning their host-galaxy properties. Anyway, this formal operation involves only one source (J1412−0307), and does not impact the statistics of either the MagE sample or low-z CL sample.
8 If low-luminosity (dwarf) galaxies are well fitted by a single Sérsic with a small index (e.g., n = 1.32 for J1636+4102), then it is somehow arbitrary to classify them as pure disk or spheroidal galaxies. ness induced, we refrain from advanced explorations of the structural and morphological properties of the host galaxies.
Just as we stated in the beginning of this §3.3, although there are a few galaxies in the three samples showing somehow blue colors (say, g − i < 1; see also the galaxies below the line of the median colors of Seyfert 2 galaxies in the g − i versus M ⋆ diagram), their central regions are almost all red. For instance, the two MagE sources with the bluest global color in the sample (J0839+0724 and J1113+0007) have g − i < 1 (see Figure 5); however, if we use the SDSS fiber magnitudes to derive their colors, then their g − i values (after Galactic reddening correction and k-correction) would be 1.17 and 1.28, respectively, and well above the mediancolor line of Seyfert 2 galaxies in Figure 5. As it is not easy to perform AGN-starlight decomposition to the fiber magnitudes, we do not delve into the use of fiber magnitudes but use instead the quantities based on fiber spectra and the properties of the inner Sérsic imaging components in the following.

In Terms of Inner-region Properties
For the six MagE sources and all except two sources in the CL sample, 9 we measured the D 4000 and Hδ A values of the host galaxies on the basis of their SDSS spectra. The SDSS fiber aperture is 1. ′′ 5 in radius, corresponding to 1.47 kpc at the median redshift (z = 0.05) of the 35 MagE and CL sources.
As for the large-variability RM sources, because of the lack of useful spectra as stated in the last part of §2.2.2, we did not use them in the investigation concerning the inner regions of the host galaxies.
First, we employ the traditional diagnostic diagram of D 4000 and Hδ A . Figure 6 shows the distributions of the MagE and low-z CL AGNs, together with non-AGN galaxies and Seyfert 2 galaxies as comparison. The distribution of Seyfert 2 galaxies is peaked at D 4000 ≈ 1.3 and Hδ A ≈ 1.7Å. All but one MagE source are outside the densest 40% contour of Seyfert 2s, with larger D 4000 and/or smaller Hδ A ; this trend is similarly followed by the low-z CL sources. We perform a two-dimensional K-S test to the distributions in the Hδ A versus D 4000 diagram. The p-value (change probability) between the low-z CL plus MagE AGNs and Seyfert 2 galaxies are 0.0007. If we apply to the Seyfert 2 sample the same redshift cut z < 0.15 as the definition of the lowz CL sample, then the p-value is 0.0025. Thus the difference between large-variability AGNs and general Seyfert 2s is significant in terms of the Hδ A versus D 4000 diagram. The difference is more significant than in the color-M ⋆ diagram, just as we expected.
It is easy to see the advantage of using inner-region properties over whole-galaxy properties. For instance, among the five bluest galaxies in the low-z CL sample, 9 Concerning the two CL sources with the AGN fraction in the r band > 30% (see Table 3 and §2.2.2), their SDSS spectra have significant AGN contribution and cannot give accurate D 4000 and Hδ A of the host galaxies. with global color 0.81 g − i 0.95 (see Table 3), four (except J1126+5134) are bluer than the median-color curve of Seyfert 2 galaxies in the diagram of global color and M ⋆ . In contrast, now in the diagram of Hδ A and D 4000 , four (except J1316+3015) have D 4000 > 1.5 indicating the mean stellar age older than 1 Gyr (for the inner-region stars), and four (except J1126+5134) are not above the median-Hδ A curve of Seyfert galaxies. By the way, the latest IFU observation by MUSE/VLT of Mrk 590 demonstrated that most of the stars within the central 10 ′′ (4.9 kpc) have older ages >5 Gyr (Raimundo et al. 2019).
Because those old so-called Lick indexes such as D 4000 and Hδ A do not directly denote the physical properties of galaxies, and because the measurement errors of D 4000 is considerable -particularly considering the small dynamic range of D 4000 for the entire galaxy population, below we employ the relatively new diagnostic tool: the diagram of SFR versus M ⋆ , which is favored by recent studies (e.g., Yu et al. 2020, Dodd et al. 2021).
In Figure  All six MagE sources, in SFR, are below the median SFRs of Seyfert 2 galaxies. All except two of the low-z CL sources are below the median SFR line of Seyfert 2s. Thus, it is unquestionable that large-variability (including CL) AGNs have smaller SFRs than Seyfert 2 galaxies of similar stellar masses. Our large-variability AGNs are mostly Seyfert 1 galaxies (not biased to being LINERs) according to the BPT diagram (see §3.2). On the other hand, Seyfert 1s are the same as Seyfert 2s in terms of the cold-gas content and star formation activity of their host galaxies (e.g., Zou et al. 2019; even for the high-luminosity counterparts namely quasars, Shangguan & Ho 2019). Thus, we conclude that the host galaxies of large-variability and CL Seyferts are in the red tail (i.e., the gas-poor, SFR-deficit tail) of the general Seyfert galaxy population.
We performed various checks and tests to our above conclusion. First of all, because the SFR values we use are based on the calibration of Zhuang & Ho (2019) with AGNemitted [O II] λ3727 flux being removed, we replot the SFR-M ⋆ diagram using SFR values from the [O II] based estimator without accounting for such nonstellar contamination; the result is displayed in Figure 7 (left panel), and clearly the above trend changes little. We tried several other SFR estimators, and the conclusion is the same. We also tested by using the 22 pure Seyferts in the low-z CL AGN sample, and again the conclusion remains intact. Recently, as described in §1, Dodd et al. (2021) presented a surprising result: local Seyfert 2 galaxies fall on the SFMS, whereas their z < 0.15 CL AGNs are on the green valley of the SFR-M ⋆ diagram. Because their z < 0.15 CL AGN sample is small (11 sources in total) and seriously dominated by LINERs (6 LINERs, accounting for 54.5%; see §3.2), we believe that their hostgalaxy difference is mainly between LINERs and Seyferts, and thus we do not make further comparison.
We noticed that there were attempts to the link central concentration of host galaxies to CL variability (e.g., Dodd et al. 2021). We plot such a diagram of inner Sérsic index and M BH (see  2021), we can see that neither our low-z sources nor the CL quasars of Charlton et al. (2019) are different from the reference AGN sample (namely Seyfert 2s) in this respect. 10 That is, the claim is not supported for a connection between high stellar density in the core region (high central concentration) and CL AGN phenomenon. Again, just as we warn in Section 1 against the unconscious use of the SFR values in the ready-made catalogs, there is a similar caveat here: the galactic bulge and disk parameters listed in the ready, mass-produced catalogs (e.g., Simard et al. 2011) were fitted with a general scheme aimed at normal galaxies, and thus special treatments are usually required for specific galaxies, particularly when the galaxies of interest are significantly nonaxisymmetric and/or have a Type 1 AGN (see §2.3). In particular, the abnormally large Sérsic indexes (n > 4) are generally due to poor fitting (e.g., when a galaxy is not well resolved).

Difference between CL quasars and local CL Seyferts.
More interestingly, we note that CL quasars (namely highluminosity AGNs) seem to be different from low-z CL/largevariability Seyferts in host-galaxy properties. In terms of the diagrams of color versus M ⋆ , Hδ A versus D 4000 , and SFR versus M ⋆ (see Figures 5, 6 and 7), the host galaxies of the four CL quasars of Charlton et al. (2019) are relatively young, above the upper SFR boundary of the green valley (consistent with the SFMS); the host galaxies of the six CL AGNs at z > 0.15 of Dodd et al. (2021) hold similar properties. There is a trivial factor contributing to this difference between CL quasars at relatively high redshifts versus CL Seyferts at low redshifts: the aperture effect, wherein the spectroscopic aperture takes more outer-region starlight for high-z AGNs than for low-z AGNs, and the outer regions of galaxies are generally bluer and younger and have more SFRs than their inner regions. But this factor seems unlikely to explain here, because even in terms of whole-galaxy properties (global color, Charlton et al. 2019;global SFR, Dodd et al. 2021) see also §3.3.1) the difference between CL quasars and local CL Seyferts is evident already. It seems that there is no difference in global colors and SFR between CL quasars and normal quasars.
Thus we tend to believe that the host-galaxy difference is real and physical (albeit the evidence comes only from the above 10 CL quasars in total): such a difference corresponds to the difference in variability pattern, namely secular pattern (low-z Seyferts) versus high-amplitude tail (EVQs and CL quasars) as stated in the second paragraph of §1. If this were the case, there should be a deep physical link from nuclear fueling flows to the structure of accretion disks; we will follow this line of thought in §3.4, in detail for the part of the local large-variability (CL) Seyferts and briefly for CL quasars and EVQs.

Dependence of Accretion Disks on Nuclear Fuels
Why do the accretion rates of these large-variability (CL) AGNs, indicated in the continuum and emission-line light curves, change so significantly on timescales of years? This is the current hotly debated theoretical question concerning the physics and structure of accretion disks (e.g., Ross et al. 2018; Dexter & Begelman 2019; Jiang & Blaes 2020), which is beyond the scope of this work. Rather than exploring accretion disks that directly produce the optical continuum emission, here, following the discussions on Mrk 1018 and Mrk 590 (e.g., Denney et al. 2014;Husemann et al. 2016;Raimundo et al. 2019), we would like to discuss broadly from the perspective of nuclear fuel supply. This is inspired by two factors: the preference of large-variability (including CL) Seyferts for red, SFR-deficit (gas-poor) galaxies as discovered in this work and the curve of fueling rate as a function of time produced by numerical simulations of nuclear fuel supply from the ∼1pc scale down toward central massive black holes (e.g., Cuadra et al. 2006Cuadra et al. , 2008Ressler et al. 2018). Also, this is along the line of thought suggested by the host-galaxy difference between CL quasars and CL Seyferts as discussed in §3.3.2: a deep physical link from nuclear fueling flows (or fueling modes) to accretion disks. At the end of this subsection, we discuss to some extent the difference between CL quasars and local CL Seyferts.

Scenario for Local CL/Large-variability Seyferts (Nuclear Famine Fueling): Cold-clump Formation and Episodic Accretion
Low-z Seyferts are generally triggered and fueled by various secular processes (internal or environmental; Kormendy & Kennicutt 2004), and their AGN activity does not depend on the host-galaxy properties on larger scales than, e.g., nuclear star clusters or even smaller structures (Ho 2008;Kormendy & Ho 2013). Kauffmann & Heckman (2009) proposed two distinct regimes of AGN fueling in nearby galaxies. One is the "feast mode" in blue galaxies with huge deposit of cold gas ready on 1 pc scales, which could sustain a steady fueling flow inward. The other one is the so-called "famine mode" in red galaxies with old stellar populations, where the SMBH fuel supply is mainly from slow stellar winds produced by evolved stars (e.g., Davies et al. 2007), or from external cold gas via minor mergers or via accretion of cold-gas clumps in the intergalactic medium particularly when the gas-poor galaxies are in small galaxy groups as realized later (e.g., Davies et al. 2014, their §8.1, §8.2.2 and §8.2.3; Davies et al. 2017). In this context, the 51 large-variability (including CL) AGNs in the 3 samples of the present study should be in the "famine mode," being in old red galaxies.
We must point out that the theoretical scenario we propose here is actually a modified and narrowed version of the "feast vs. famine" notion of Kauffmann & Heckman (2009): while Kauffmann & Heckman (2009) took into account the amount of cold gas in galactic bulges (i.e., on kiloparsec scales), we adopt the inference of King & Pringle (2007) that only the cold gas in the nuclear region (on 1 pc scales) can affect BH accretion. 11 Thus we deliberately use the terms "nuclear feast" and "nuclear famine." In the nuclear famine mode, whether from slow stellar winds or externally from the galactic environment, the information on the cold-gas origin is still retained in the nuclear region (see, e.g., Hobbs et al. 2011); in other words, the BH accretion process does not need larger-scale information than the nuclear cold gas can provide. In the nuclear feast mode, the cold-gas supply on a ≈1 pc scale is much more than required to feed the central BH, and the accretion-disk properties do not care what is the origin of the nuclear cold gas, whether via wet major mergers or secularly from the galactic disks of late-type galaxies (see, e.g., Davies et al. 2014).
However, presently there is no observational census for cold gas on 1 pc or smaller scales of AGNs, and the fueling passage from the 1 pc scale (i.e., AGN tori, if present) down to the outer boundary of accretion disks is beyond the capacity of current observing facilities. It is just now that even ALMA has been launching observations for nearby AGNs and can merely resolve the molecular tori of only a few parsecs (e.g., Combes et al. 2019), let alone the subparsec fuel flows.
Fortunately, we can get insights from the three-dimensional simulations of the fueling passage on these scales in the literature, although those simulations are for the fueling of Sgr A * (as a dim AGN) in the Galactic center where the fuel is fast stellar winds from young stars (Cuadra et al. 2006, 11 Although the smallest spatial scale observationally probed in the present work is still on the bulge or pseudo-bulge scale (the spectra as well as the decomposed inner Sérsic imaging component; §3.3.2), we are convinced by the arguments of King & Pringle (2007) for local Seyferts and believe that the related cold-gas supply is on the 1pc scale. See also Ho (2008) and Kormendy & Ho (2013). 2008Ressler et al. 2018). According to Ressler et al. (2018, see their Fig. 10), the accretion rate measured at 1.5 × 10 −4 pc (namely 370 times the Schwarzschild radius of the Sgr A * BH) can vary by a factor of 9 within tens of years. This is consistent with the simulation of Cuadra et al. (2008, see their Figures 1-3 where the time sampling is every 30 years). Cuadra et al. (2008) clearly demonstrated that while the accretion-rate curve of hot gas is smooth and has only small-amplitude fluctuation, the episodic infall of cold-gas clumps produces sharp peaks in the accretion rate! If the time resolution gets higher (e.g., without averaging of snapshots; see Cuadra et al. 2006), the increasing of accretion rate by cold clumps would get more steep and abrupt, and the duration of the peaks could be shorter (e.g., a few years).
We can imagine that such sharp accretion-rate peaks 12 should be more common in red galaxies, which do not have young stars but AGB and red giant stars to produce slow stellar winds (accordingly, the AGN is not as dim as Sgr A * ); see, e.g., Cuadra et al. (2006) and Shcherbakov et al. (2014) for detailed discussions. Certainly, if we make an analogy of AGN fueling by those simulations of Sgr A * , those sharp peaks actually reflect merely reflect the variability in fueling rate into the outer boundary of accretion disks and will be modulated (or smoothed?) by accretion disks later. At present it is unclear what the signal of the fueling-rate peaks would look like in the optical continuum emission of accretion disks. On the other hand, in fact, the notion of episodic fueling and thus episodic accretion was perceived and formulated from different perspectives (R. Davies 2021, private communication). Initially in the study of cooling flows in galaxy clusters, there are theoretical models, numerical simulations and observations that coldgas clumps (or called blobs, streams, filaments or alike in the literature) can condensate out of the hot intergalactic gas, and several researchers have further argued that some of the cold clumps can sink toward the central BHs of galaxies (e.g., Pizzolato & Soker 2005;Gaspari et al. 2013). In particular, Gaspari et al. (2013) coined the notion of "chaotic cold accretion," in which the large-scale (namely intergalactic and galaxy-scale) cold clouds (clumps) can lose angular momentum and some of them finally "rain" into the central BH via recurrent collisions among clouds and and between clouds and the clumpy AGN torus; the raining is the actual accretion process. They also proposed that chaotic cold accretion seems to be an excellent model to explain AGN variability (see their §7.3). King & Pringle (2006) put forward and derived a "chaotic accretion" mechanism, characteristic of a series of small-scale, randomly oriented accretion events to feed BHs, which was specifically applied to (and analyzed in) the case of nearby Seyferts (King & Pringle 2007); later on, following this line of thought, Hobbs et al. (2011) put forward a concrete scenario emphasizing the role of turbulence. Turbulence is also important in the simulations of subparsec fueling of the aforementioned Cuadra et al. 12 To put more precisely, as described above, the fueling rate peaks. (2006,2008) and Ressler et al. (2018). After their longterm observational studies on nuclear fueling of AGNs, Davies et al. (2014Davies et al. ( , 2017 proposed that there are two modes of inflows feeding low-redshift AGNs, one being quasicontinuous with a plentiful internal supply of gas (e.g., in gas-rich spiral galaxies), the other being stochastic events accreted from the galactic environments (e.g., important for gas-poor galaxies in moderately dense galaxy groups). 13 We can see that, besides our approach on AGN variability, the notion of two (nuclear) fueling modes has been reached from various perspectives in recent years. Again, just as we point out in the second paragraph of this scenario, we would like to stress that for local CL and large-variability AGNs only the nuclear fuel (on 1 pc) scale matters.
A last point we would like to mention: in red galaxies it is probable that even the broad-line region per se is episodic (see Denney et al. 2014). When the aforementioned cold-gas clumps fall into the right radii, they get ionized accordingly by the AGN continuum and become effective in producing certain optical broad lines. This is an interesting picture, and can be modeled by taking the "locally optimally emitting clouds" (LOC; Baldwin et al. 1995) approach, just as, e.g., Korista & Goad (2004) did. Without data from either observations or modelings so far, we however refrain from a full discussion on episodic broad-line regions in this paper.

Speculation for CL Quasars (Nuclear Feast Fueling)
More excitedly, we cannot help thinking about the trend wherein the host galaxies of CL quasars and EVQs appear quite different from that of local CL Seyfert galaxies. Although the number of such quasars with host-galaxy properties analyzed is small (10 sources at most; Charlton et al. 2019, Dodd et al. 2021, there are reasons to believe this trend is real ( §3.3). We think this trend is consistent with the difference between local CL/large-variability Seyferts and EVQs in variability pattern (secular variation versus high-amplitude tail), both differences being suggestive of a deep link between nuclear fuel supply and the structure of accretion disks. We speculate that the "feast mode" of nuclear fueling may account for both the preference of CL quasars and EVQs for blue galaxies and the corresponding accretion-disk structure required by CL quasars and EVQs (see, e.g., Jiang & Blaes 2020). With more than enough cold gas available on ≈1 pc scale (recalling AGN tori as the fuel reservoir), surely the structure of the accretion disks in the feast fueling mode is starkly different from that of the accretion flows for local large-variability (CL) Seyferts described in the above. We have to admit this idea for CL quasars and EVQs is quite speculative, without solid supports from either observations or simulations so far, and thus we defer the exploration for the future.

New Thinking on the Variability Selection for IMBHs
The initial goal of the project branching out from this work is to search for IMBHs. We now return to the IMBH topic with the implication of the finding of this work to the IMBH searching by optical variability. In the flowering time-domain astronomy era, variability selection for lowmass AGNs is promising (see §4.5 of Greene et al. 2020 for a brief review). This field is just developing, and there are only a few such searches based on optical variability so far (Morokuma et al. 2016;Baldassare et al. 2018Baldassare et al. , 2020Guo et al. 2020;Martínez-Palomera et al. 2020).
A general trend in those studies is that the optical continuum variability is generally of low level, less than 0.1 mag, i.e., with the peak-to-valley flux ratio being < 1.1 (see Martínez-Palomera et al. 2020). Note that the "nucleus magnitudes" (or fluxes) reported in Baldassare et al. (2018Baldassare et al. ( , 2020 and Martínez-Palomera et al. (2020) are measured through a small aperture (matching to the seeing in order to collect AGN flux), either by adding together the fluxes of the difference image and template image for each data point (Baldassare et al. 2018(Baldassare et al. , 2020, or by directly performing aperture photometry on the source images (Martínez-Palomera et al. 2020); this is consistent with the data points for the nuclear magnitudes/fluxes displayed in most light curves, such as those used in this work. To clarify possible confusion, we would like to mention in passing the following point. Morokuma et al. (2016) used a different convention to report the data points of their light curves (see the "Flux (sub)" and "Mag (sub)" columns in their Table 1): they adopted the same difference imaging methodology as in Baldassare et al. (2018), but they presented the fluxes and magnitudes measured from the difference images (i.e., with the template subtracted, which is actually more physically motivated). Thus Morokuma et al. (2016) reported a large peak-to-valley magnitude change ∆g = 0.42 (i.e., a flux change by 1.5 times), but this change is with respect to the flux scale of 1 µJy; in fact the peak-to-valley flux change in physical units is tiny, 1.2 µJy, un-surprisingly. Recently, Guo et al. (2020) reported a low-mass AGN at z = 0.823, identified from the Dark Energy Survey (DES) Supernova field by optical variability. Reading from the light curve (their Figure 1), the peak-to-valley magnitude change is really large, ∆g ≈ 0.5 mag (varying around g ≈ 22.0, PSF magnitude measured from the source images).
Besides, in the literature optical-variability-selected lowmass AGNs as a sample were reported to tend to be in galaxies bluer than the general galaxy population (Baldassare et al. 2018), which agreed with the conclusion of Heinis et al. (2016) on the host-galaxy colors of their variability-selected AGNs. Moreover, variability-selected low-mass AGNs of Baldassare et al. (2020) tended to be even bluer (in g − r color) than the low-mass AGNs selected in terms of optical narrow lines by Reines et al. (2013), and mostly resided in the star-forming region of the BPT diagram. All the above reports on variability-selected low-mass AGNs are inconsistent with the discovery of our present work. We are not aware of the reason of the inconsistency at this point.
Instead, we would like to note that the present work reminds us not to ignore red galaxies. It is rewarding to search for low-mass AGNs in red galaxies by optical variability. First of all, as discovered in the present study, low-mass AGNs in red galaxies would have continuum variability of larger amplitude than those in blue galaxies; or in other words, red galaxies have a larger fraction hosting variable low-mass AGNs than blue galaxies. Additionally, in stark contrast to blue galaxies, red galaxies have little starformation dilution of the AGN emission. We could think a little bit further about designing this kind of variability search. Galaxy groups and even clusters would be ideal target fields, where red galaxies such as S0 and spheroidal galaxies are plentiful, and thus the survey efficiency and productivity would be high. Besides, if the above interpretation invoking episodic fueling ( §3.4) is correct, the environmental secular processes red galaxies suffering in galaxy groups and clusters (Kormendy & Kennicutt 2004) can enhance the intermittency in fuel supply (Davies et al. 2017).
4. SUMMARY During our spectroscopic MagE campaign initially planned to search for intermediate-mass black holes (IMBHs) in nearby broad-line AGNs, we unexpectedly found six variable AGNs with relatively small black hole masses (M BH ≈ 10 6−7 M ⊙ ), their broad-Hα fluxes varying by 1.3 − 3.0 times (namely 0.3 − 1.2 mag) on timescale of years. Most surprisingly, among our broad-line AGNs identified by MagE (15 in total), those hosted by blue galaxies generally vary little, whereas a significant fraction of those hosted by red galaxies exhibit large broad-Hα variability; in other words, the six variable sources are predominantly in red galaxies.
The above unexpected "bonus" motivated us to explore AGNs with large variability (including CL AGNs), aiming at their host-galaxy properties (particularly as to any connections between those properties and the CL and largevariability AGN phenomenon) in a systematic way. We collected all the low-z changing-look AGNs available in the literature and RM AGNs with broad-Hα F var > 0.1, and performed careful imaging and spectral fittings. From our observational investigations, we draw the following two main conclusions about the connection between host-galaxy properties and the CL and large-variability AGN phenomenon: • Local CL and large-variability AGNs (mainly Seyferts) reside in redder, more SFR-deficit (presumably gaspoor) galaxies than the control sample of local Seyfert 2 galaxies. That is, the host galaxies of those strongly variable Seyferts are in the red tail (i.e., the gaspoor, SFR-deficit tail) of the general Seyfert galaxy population.
• In contrast, there is a significant trend that their more luminous counterparts, namely CL quasars (CLQs) and extremely variable quasars (EVQs), are different from local CL Seyferts in host-galaxy properties. For instance, in terms of the diagnostic diagram of global color and M ⋆ , the host galaxies of CLQs are generally blue (see also Charlton et al. 2019); in terms of the diagram of SFR and M ⋆ (of the inner regions, Figure 7; see also Dodd et al. 2021), local CL Seyfert galaxies are located in the green valley, whereas CLQ hosts are in the so-called SFMS.
These two discoveries inspired our theoretical thinkings about the physical link between nuclear fuel supply and accretion-disk structure, and about the implication to the field of IMBH research in turn. We proposed a physical scenario for local CL and large-variability Seyferts, a speculation for CLQs and EVQs, and an implication for IMBH searches based on optical variability, as follows: • We presented an explanation for the preference of local CL and large-variability Seyferts to old red host galaxies from the perspective of the nuclear fueling mode, which is a modified and narrowed concept of "famine mode" proposed by Kauffmann & Heckman (2009). The concrete mechanism may be revealed by three-dimensional simulations of the fueling passage from the 1 pc scale down to the outer boundary of accretion disks. In similar existing simulations, coldgas clumps can be formed stochastically in the fueling flow on 1 pc scales. While the accretion-rate curve of hot gas is smooth and has only small-amplitude fluctuation, the episodic infall of the cold-gas clumps produces sharp peaks in the accretion rate (measured at the outer boundary of the accretion disk). We discussed the feasibility of this scenario, namely the time scales of the rising, lasting, and falling of this kind of cold-clump accretion activity ( §3.4).
• We speculated that the "nuclear feast mode" may account for both the preference of CL quasars and EVQs to blue galaxies and their variability pattern (high-amplitude tail of the continuous distribution) that is different from the secular variation of local CL Seyferts. With more than enough cold gas piled up on ≈1 pc scale (say, AGN tori as the reservoir), surely the structure of the accretion disks in the feast fueling mode (see, e.g., Jiang & Blaes 2020) is starkly different from that of accretion flows for local CL Seyferts. We defer the exploration along this line of thought for the future.
• We proposed a new thinking on the design of opticalvariability selection for IMBHs: to launch variability searches in red galaxies. This strategy would be more efficient than usual blind surveys. And it can be regarded as a kind of deliberate debiasing and reminder, because the variability-selected low-mass AGNs so far tend to be in blue galaxies.
We really became excited by the observational discoveries about the connection between host-galaxy properties and the large-variability (including CL) AGN phenomenon, triggered by the unexpected "bonus" (the six variable lowmass AGNs) out of our observing campaign actually aimed at IMBHs; also excited ourselves by the ensuing theoretical thinkings listed above. There appear to be a lot of lines of fruitful work for the future. From an observational standpoint, the direct measurement of cold-gas content of the host galaxies-particularly on 1 pc scale-of such strongly variable Seyferts and quasars is in demand; comparison studies of considerably large samples, i.e., being statistically robust, are necessary. Through numerical experiments, it would be instructive to carry out threedimensional simulations for different nuclear fueling modes, simulating the fueling passage from ≈1 pc scale down toward (the outer boundary of) the accretion disks, or even handling both the nuclear fueling flows and some part of accretion flows simultaneously. And finally, theoretical insights (see Antonucci 2018) are always needed.
We thank the anonymous referee for a helpful report that significantly improved this paper. We also thank Richard Davies for helpful comments, particularly on the two modes of AGN fueling; C. Martin Gaskell and Ski (Robert Antonucci) for valuable comments and discussions, particularly on accretion and on early discoveries and facts long ignored (ahead of their time), Ning Jiang for his guidance on the 2Ddecomposition of SDSS images with GALFIT; and Fuguo Xie for advice and discussions during the course of this work. This work is supported by the Natural Science Foundation of China grants (NSFC 11703079, 11873083) and the "Light of West China" Program of Chinese Academy of Sciences (CAS). Su Yao acknowledges support by an Alexander von Humboldt Foundation Fellowship. D.W.X. and J.W. are supported by the Natural Science Foundation of China grant 11773036.This work has made use of the spectra obtained with the Magellan Baade Telescope/MagE, Hale Telescope/DBSP, and Xinglong 2.16m telescope/BFOSC. The Hale Telescope/DBSP data were obtained through the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories, Chinese Academy of Sciences, and the Special Fund for Astronomy from the Ministry of Finance. Observations obtained with the Hale Telescope at the Palomar Observatory were obtained as part of an agreement between the National Astronomical Observatories, Chinese Academy of Sciences, and the California Institute of Technology. This project is supported by the National Natural Science Foundation of China (NSFC-11421303, NSFC-11603021, and NSFC-11833007). We also acknowledge the support of the staff of the Xinglong 2.16m telescope, which was partially supported by the Open Project Program of the Key Laboratory of Optical Astronomy, National Astronomical Observatories, Chinese Academy of Sciences.

VARIABILITIES OF THE SIX MAGE AGNS
In this Appendix, we present two parts of the data analyses for the six variable MagE sources: analysis of the multiband light curves and X-ray spectra fitting.

A.1. Multiband Light Curves of MagE sources
We construct the multiband light curves for the six MagE sources using all of the photometric data publicly available (see Figure A1).
The optical light curves were constructed using the data from Catalina Real-Time Transient Survey (CRTS) (Drake et al. 2009) and the V -band and g-band magnitudes from ASAS-SN (Shappee et al. 2014;Kochanek et al. 2017). CRTS is one of the largest time-domain optical surveys currently operating, which is performed using unfiltered light and nominally transformed to the V -band zero point. ASAS-SN is a long-term project designed to monitor the entire sky on a rapid cadence to find nearby supernova and other bright transients, providing V -band and g-band photometric data. We first removed the data points with large uncertainties and then binned the data with a bin size of 90 days, which roughly correspond to the natural observing gaps (see Figure A1). In every panel, every original observed data point is plotted together with their ±1σ error bars. For each binned data point, we report the median value within a bin, and the error is the sum of two parts (in quadratic form): (1) the random error of the mean, calculated according to the error propagation formula from the quoted statistical errors of every measured data points in the bin, and (2) the standard error of the mean, i.e., the standard deviation (namely difference) of the measured data points divided by √ N (N being the number of the data points in a bin). The MIR light curves are constructed using the W 1 and W 2 of the WISE (Wright et al. 2010) and NEOWISE-R (Mainzer et al. 2014) surveys. Following Sheng et al. (2017), we removed bad data points with poor image quality ("qi_fact" < 1) or with flagged moon masking ("moon mask" = 1). We binned the data points every half a year. The average values and their errors in every bin are calculated in the same way as the above optical light curves.
The X-ray data are obtained from the literature and archives, by XMM-Newton, Chandra, and ROSAT. Three sources, J1257+2724, J1412−0307, and J1442+0119, have multiple observations over the past decades. For the purpose of variability, we simply use the X-ray fluxes retrieved from 4XMM-DR9 catalog (Webb et al. 2020), Chandra source catalog (CSC, Evans et al. 2010), and ROSAT Catalogs (White et al. 1994and Salvato et al. 2018Panzera et al. 2003Anderson et al. 2007 for J1442+0119). The XMM-Newton and Chandra observations cover the energy range 0.2-10 keV, while ROSAT only covers 0.1-2.4 keV. To build the X-ray light curves, we adopt the 0.2-2 keV X-ray fluxes for XMM-Newton and Chandra observations and convert the ROSAT fluxes into the same energy range using the measured spectral slope of each source (or the mean slope if two or more slopes are measured; see §A.2). The XMM-Newton 0.2-2 keV fluxes are calculated by summing up the fluxes of EPIC band 1 (0.2-0.5 keV), band 2 (0.5-1.0 keV), and band 3(1.0-2.0 keV) in the 4XMM-DR9 catalog. The Chandra 0.2-2.0keV X-ray fluxes are calculated by summing up the fluxes of "U_Flux_Ap", "S_Flux_Ap", and "M_Flux_Ap" in CSC.
To characterize the variability of those multiband light curves, we calculate a set of variability statistics, which (except for R max ) are otherwise not applicable when the data points are few (such as the multiepoch broad-Hα flux data). Besides the maximum variability R max and fractional variability amplitude F var mentioned in §2.1.1, we use a third measure, normalized excess variance (σ 2 rms ), and its error (see, e.g., Martínez-Palomera et al. 2020). Below we give their functional definitions. F var , called fractional variability amplitude historically in the literature (see Edelson et al. 2002), is defined to be the squared root of the excess variance of the total light curve (the excess is presumably the intrinsic, with the variance of random measurement errors subtracted), which is then divided by the mean of the light curve, Here the quantity σ 2 is the variance of the light curve (i.e., of all the observational data points), x is the mean, and δ 2 is the mean of the squared errors (δ i ) associated with the observations (x i ), i.e., where N is the number of the observational data points (fluxes or magnitudes) of the light curve. F var is commonly used in the literature of both X-ray time-series analysis and optical reverberation mapping. Yet it has two limitations: (1) by definition, F var cannot apply when σ 2 < δ 2 (see Table A1); (2) because both σ 2 and δ 2 involve summing over all the data points of the light curve, it does not appear necessary (and even not good in principle) to do the summing twice separately. Thus we introduce also σ 2 rms , called the normalized excess variance, which is similar to F var being overall measures of the intrinsic variability of light curves that correct for measurement errors from photon counting and detector read noise, but is free from the above two limitations. The definition is straightforwardly the calculation of the excess (presumably intrinsic) variance, normalized by the square of the mean (x), as follows: The error of σ 2 rms , due to Poisson noise, is as follows (Martínez-Palomera et al. 2020), Light curves with (σ 2 rms − err(σ 2 rms ) ) > 0 can be regarded to be intrinsically variable (Martínez-Palomera et al. 2020). All the above measures are calculated for every light curves, and listed in Table A1. R max , being the simplest and not accounting for measurement errors, is calculated based on the binned light curves (see the above); the other are are calculated based on the original data points and their associated measurement errors.
A.2. X-ray Spectra of MagE sources As described above, J1257+2724, J1412−0307, and J1442+0119 were observed by XMM-Newton or Chandra. This subsection analyzes their X-ray spectra to check their flux changes and to investigate whether there are any the possible absorption features.
As demonstrated in Figure A1, both J1257+2724 and J1412−0307 have a dozen X-ray observations and have shown violent variability in X-ray during the past 30 yr. For J1257−2724, we select XMM-Newton observation runs of 2006 June 14 and 2010 December 5 to extract the X-ray spectra, which have a high enough S/N and represent the low-and high-flux states, respectively. For J1412−0307, in a similar way, we extract the the XMM-Newton spectrum from the observation on 2008 July 27 as its highflux state. Its lowest-flux states appear in the observations between 2001 and 2004, which however have too few statistics to provide a meaningful constraint on the absorption. Thus we select the observation on 2015 as a test of the absorption in the lowflux state. The third source, J1442+0119, has only one observation run, by Chandra. For the XMM-Newton data, the spectra are extracted preferentially from PN for its large effective area. When PN data are not available, we combine spectra from two MOS CCD arrays and response files. The XMM-Newton spectra are re-binned so that each bin achieves S/N 3. For the Chandra data, the spectra are extracted following standard procedures. The Chandra spectra are rebinned so that each bin contains at least 25 counts.
We started with a single power-law model with hydrogen absorption fixed at the Galactic value for each source, which gives a well fit to all the spectra. Then we added into the model an intrinsic hydrogen absorption component with the column density (N H ) being a free parameter. But this model does not improve the fit in all cases (the five spectra), and the best-fit N H values are all small, as follows. For each source, the best-fitting models give upper limits of N H at 90% confidence level: < 4 × 10 20 and < 2.8 × 10 21 cm −2 for the high-state and low-state spectra of J1257+2727, respectively; < 7 × 10 19 and < 3.6 × 10 20 cm −2 for the high and low states of J1412−0307, respectively; < 1.1 × 10 21 cm −2 for J1442+0119. Thus we suggest that the intrinsic absorption is insignificant.
We also tried adding a thermal component (bbody in XSPEC) to the model for J1257+2727 and J1412−0307. It turns out that only the fitting to the high-state spectrum of J1412−0307 can be improved in a statistical sense, with an decrease ∆χ 2 = 9 after the adding of two more free parameters (i.e., the decrease of the degree of freedom ∆ dof = 2); the F-test is marginally significant, with the p-value being 0.017. This improvement may indicate a possibility that the X-ray variability can be explained by the variation of a thermal component or similar. Yet at present only a single spectrum (the high-flux one) can marginally constrain this additional component, and thus we cannot say anything about the variation of this component. Therefore, we leave the advanced investigation of X-ray properties for a future paper, and in this work adopt the single power-law fitting as the final results. Figure A2 shows the X-ray spectra and their best fits by the single power-law model, as well as the respective residuals. For J1257+2724, the best-fit photon indexes Γ of its high-and low-flux states are 2.3 and 1.1, respectively. Its highest and lowest 2-10 keV luminosities are 4×10 41 erg s −1 and 2.2×10 41 erg s −1 . For J1412−0307, the best-fit photon indexes keep at 1.8 in both its high and low states, and its highest and lowest 2-10 keV luminosities are 2.4×10 42 erg s −1 and 4.5×10 41 erg s −1 . The photon index of J1442+0119 is 1.64, with the 2-10 keV luminosity of 3.92×10 41 erg s −1 .     (2015) and Bentz & Manne-Nicholas (2018). For IC 4218, whose reverberation mapped MBH is not available, we measured its MBH based on its SDSS spectrum in the same way as the MagE sources (see §3.2). Col. (6) Fractional variability amplitude of broad Hβ, characteristic of the overall variability of broad-Hβ light curves, taken from the papers noted in the reference column. Note that here broad Hβ is used because its Fvar is more commonly available than broad Hα in the literature. Cols.(7)-(9) Host-galaxy magnitude, calculated with AGN contamination removed and Galactic extinction corrected. The data with ±1σ errors are from our GALFIT fittings, and the others are taken from Koss et al. (2011). Col.(10) g − i color of the host galaxies, with Galactic extinction and k-corrections performed. Col.(11) Fraction of AGN light to the total in the r-band images, either from our GALFIT fittings or from Koss et al. 2011 (see Cols. (7)- (9)). Col.(12) Stellar mass, estimated with AGN contamination removed. The data marked with a + are taken from Bentz & Manne-Nicholas (2018), those with an * are from Koss et al. (2011), and the remaining one is our estimate in the same way as for the variable MagE sources.   (5) The uncertainty of σ 2 rms due to Poisson noise (see Martínez-Palomera et al. 2020). Col.(6) Light curves with (σ 2 rms − err(σ 2 rms ) > 0) can be regarded to be intrinsically variable (see Martínez-Palomera et al. 2020). Col. (7) Fractional variability amplitude Fvar, calculated from the flux data of light curves (see Barth et al. 2015). Light curves with σ 2 − δ 2 < 0 (see Equation A1) are marked with '−', indicating that Fvar is not applicable. Col.(8) Maximum variability amplitude of light curves. The X-ray data of J1442+0119 are of two epochs only and thus not possible to calculate the statistics in Cols.(4-7).      Table 2), 30 CL AGNs at z < 0.15 (dark red filled diamonds; see Table 3), and the 15 large-variability RM AGNs (brown filled squares; see 4) on the diagram of the global g − i color versus M⋆. Also shown are the 4 CL quasars with z 0.2 of Charlton et al. (2019, light blue open diamonds), and the host galaxies of 52 QSOs from SDSS Stripe 82 analyzed by Bettoni et al. (2015, blue open circles). As a comparison, we plot (by gray dots) about 26,000 low-z Seyfert 2 galaxies selected in terms of BPT diagram by Dong et al. (2010). The contours represent 90% (the outermost), 50%, and 20% of the Seyfert 2 galaxies enclosed, respectively. The black filled dots and dashed line denote the median g − i of Seyfert 2 galaxies in the M⋆ bins (grouped by vertical gray dashed lines); see the text in §3.3.1 for the detail.  (Table 2), 29 low-z CL AGNs (Table 3), as well as the 4 CL quasars at z 0.2 (Table 3), on the diagram of D4000 and HδA. Those AGNs are denoted in the same way as in Figure 5. As comparison, we also plot the aforementioned ≈26,000 Seyfert 2 galaxies (green dots), as well as the ≈32,000 low-z normal galaxies selected by Dong et al. (2012, gray dots). The gray contours represent 99.9% (the outermost), 95%, 68%, and 40% of the normal galaxies enclosed, respectively; the green contours represent 95%, 68%, and 40% of the Seyfert 2 galaxies enclosed, respectively. All of the D4000 and HδA are measured from the SDSS spectra.    (Table 2), 29 low-z CL AGNs (Table 3), as well as the 4 CL quasars at z 0.2 (Table 3), on the diagram of SFR versus M⋆. Those AGNs are denoted in the same way as in Figure 5. As comparison, the aforementioned ≈26,000 Seyfert 2 galaxies (gray dots) are also plotted. The SFR values are derived from the [O II]λ3727 lines in the SDSS spectra, without (left panel) and with (right panel) correcting for the AGN contamination. The contours represent 90% (the outermost), 50%, and 20% of Seyfert 2 galaxies enclosed. THe blue dashed line represents the star forming main sequence, the red dashed-dot line represents the so-called quiescent sequence of old red galaxies, and the region between the two green dotted lines is the green valley; see the text in §3.3.2 for details.  (Table 2), 29 low-z CL AGNs (Table 3), as well as the 4 CL quasars at z 0.2 (Table 3) on the diagram of inner Sérsic index versus MBH. Among low-z CL AGNs, those with r-band Petrosian radii < 4 ′′ (i.e., resulting in unreliable inner-Sérsic indexes) are marked in gray; otherwise, they are in dark red.