A surprising excess of radio emission in extremely stable quasars: a unique clue to jet launching?

Abstract

Quasars are generally divided into jetted radio-loud and non-jetted radio-quiet ones, but why only 10% quasars are radio loud has been puzzling for decades. Other than jet-induced-phenomena, black hole mass, or Eddington ratio, prominent difference between jetted and non-jetted quasars has scarcely been detected. Here we show a unique distinction between them and the mystery of jet launching could be disclosed by a prominent excess of radio emission in extremely stable quasars (ESQs, i.e., type 1 quasars with extremely weak variability in UV/optical over 10 years). Specifically, we find that $>$ 25% of the ESQs are detected by the FIRST/VLASS radio survey, while only $\sim$ 6-8% of the control sample, matched in redshift, luminosity, and Eddington ratio, are radio-detected. The excess of radio detection in ESQs has a significance of 4.4 $\sigma$ (99.9995%), and dominantly occurs at intermediate radio loudness with R $\sim$ 10 – 60. The radio detection fraction of ESQs also tends to increase in the ESQ samples selected with more stringent thresholds. Our results are in contrast to the common view that RL quasars are likely more variable in UV/optical due to jet contribution. New clues/challenge posed by our findings highlight the importance of extensive follow-up observations to probe the nature of jets in ESQs, and theoretical studies on the link between jet launching and ESQs. Moreover, our results makes ESQs, an essential population which has never been explored, unique targets in the burgeoning era of time domain astronomy, like their opposite counterparts of quasars exhibiting extreme variability or changing-look features.

1. Introduction

Based on the relative radio intensity, i.e. radio loudness, which is defined as the ratio of the radio flux density to the optical one (e.g. Kellermann et al. 1989), quasars could be divided into radio-loud (RL) and radio-quiet (RQ) ones, or jetted and non-jetted ones (Padovani et al. 2017; Panessa et al. 2019). The RL quasars, with powerful relativistic and well collimated jets (Blandford et al. 2019), are typically 1000 times brighter in radio than the RQ ones (Miller et al. 1993; Panessa et al. 2019). While it is generally believed that the rotational energies of the black hole (BH) or the inner accretion flow could be extracted to power jets (Blandford & Znajek 1977; Blandford & Payne 1982; Blandford et al. 2019), why the relativistic jets have only been launched in a small population ($\sim$ 10%) of quasars has been a mystery for many decades. This is particularly puzzling considering that both the RL and RQ quasars have rather similar, except in radio, spectral energy distributions (e.g. Elvis et al. 1994; Shang et al. 2011), suggesting both populations are powered by a similar accretion process. Furthermore, Dunlop et al. (2003); Floyd et al. (2004, 2013) found that hosts of both RL and RQ quasars are similarly massive galaxies, and they remain star forming (Floyd et al. 2013; Kellermann et al. 2016), suggesting the host does not dominantly determine the radio loudness of quasars. While observations have revealed that the RL quasars have more massive BHs and smaller Eddington ratios (but with strong overlap in range) compared with the RQ ones (e.g. Laor 2000; Lacy et al. 2001; Ho 2002), searching for observational differences (other than jet induced phenomena) between RL and RQ quasars with matched BH mass and Eddington ratio may yield essential clues to understanding the jet launching mystery.

Variability has been a defining characteristic of AGNs and quasars (e.g. Ulrich et al. 1997). Studying the variability, particularly in UV/optical which is believed to be predominantly driven by magnetic turbulence in the accretion disc, can be of great help in probing the yet unclear underlying physics of the inner accretion process. Plentiful observational studies have reported correlations between the UV/optical variation and several known parameters, including the luminosity, rest frame wavelength (i.e., stronger variability at shorter wavelengths), Eddington ratio, BH mass, redshift (Vanden Berk et al. 2004; Wilhite et al. 2005; Wold et al. 2007; Wilhite et al. 2008; Bauer et al. 2009; Ai et al. 2010; MacLeod et al. 2010; Meusinger et al. 2011; Zuo et al. 2012; Meusinger & Weiss 2013; Kozłowski 2016; Sun et al. 2018), and new parameters including X-ray loudness and equivalent width of emission lines (Kang et al. 2018, 2021). As for the RL quasars, they often exhibit stronger rapid (e.g., intra-night) variability (e.g. Gupta & Joshi 2005) and marginally stronger long-term variability (e.g. Helfand et al. 2001; Vanden Berk et al. 2004), compared with the RQ ones, which could be attributed to the contribution of the UV/optical emission from the jet which could be more variable than the disk component.

In the era of time domain astronomy, great attention of the community has been attracted to AGNs exhibiting violent variability, e.g., changing-look (CL) AGNs and quasars (e.g. Cohen et al. 1986; LaMassa et al. 2015; MacLeod et al. 2016; Sheng et al. 2017; Yang et al. 2018; Sheng et al. 2020; Green et al. 2022; Ricci & Trakhtenbrot 2022), and extremely variable quasars (EVQs, Rumbaugh et al. 2018; MacLeod et al. 2019; Ren et al. 2022), which are a small population of sources showing strong variability likely driven by yet unclear violent changes in the inner accretion disc activity. In contrast, in this work we focus on extremely stable quasars (ESQs), which exhibit rather weak or undetectable long-term (over 10 years) variability in UV/optical. Such quasars, the opposite counterparts of CL quasars and EVQs in the parameter space of variability, have never been studied in literature. We find that ESQs exhibit significant excess of radio emission compared with the normal quasars, providing unique new clues to jet launching in quasars. Throughout this work, cosmological parameters of $H_{0}=70\rm\ km~{}s^{-1}~{}Mpc^{-1}$, $\Omega_{\rm m}=0.3$, and $\Omega_{\Lambda}=0.7$ are adopted.

2. Selection of extremely stable quasars

We adopt the 10-year-long light curves (MacLeod et al. 2012) for the 9258 spectroscopically confirmed quasars in the Sloan Digital Sky Survey (SDSS) Stripe 82, which has been scanned around 60 times in the $ugriz$ bands (Sesar et al. 2007), to calculate the long-term variability amplitude of quasars. We do not include additional photometric data of them from other time domain surveys, such as Pan-STARRS1 (Flewelling et al. 2020), because the different filter transmissions between surveys would cause extra scatter to the variability measurements. As the SDSS $g$ and $r$ bands have the best photometry for quasars in Stripe 82 (see Fig. 2 of Sun et al. 2014), and the intrinsic variability of quasars is expected to be stronger at shorter wavelength, in this work we primarily adopt $g$ band to select ESQs, and utilize the other bands to secure the selection. We focus on 9146 quasars from MacLeod et al. (2012) which have at least 10 photometric measurements in the $g$-band light curve, after excluding a few unphysical data points that may be present. In addition, to use the up-to-date physical properties, such as BH mass and Eddington ratio, for these quasars, we also drop a few sources for which no counterpart in the SDSS data release 16 quasar (DR16Q) catalogue is found within 2 arcsec or the matched counterpart does not have a valid measurement on the BH mass by Wu & Shen (2022).

The excess variance, $\sigma_{\mathrm{rms}}^{2}$, has widely been utilized to quantify the quasar variability with a single model-independent parameter (e.g. Vaughan et al. 2003). However, as demonstrated in Appendix A for the ESQs, the canonical $\sigma_{\mathrm{rms}}^{2}$ could be significantly biased by a minority of epochs with too large photometric uncertainties. We thus revise the canonical calculation of $\sigma_{\mathrm{rms}}^{2}$ by adding weight to each photometric measurement (see Appendix A).

Figure 1.— Distribution of the $g$-band $\sigma_{\mathrm{rms}}$ of 9146 SDSS quasars in Stripe 82. For sources whose $\sigma_{\mathrm{rms}}^{2}$ are statistically non-detected with S/N($\sigma_{\mathrm{rms}}^{2}$) $<$ 2, upper limits to $\sigma_{\mathrm{rms}}$ are assigned (red). To the left of the vertical solid black line ($\sigma_{\mathrm{rms}}$ = 0.05), there are 317 quasars with smallest $g$-band variability in the sample. For comparison, the vertical dashed line marks (to its right) equal number of sources with the strongest $g$-band variability.

Figure 2.— An illustration of the $g$-band light curves of an EVQ (maximal $\Delta g>1$ mag; top panel), a normal quasar (median panel), and an ESQ (bottom panel). The quasar IDs marked in plot are from MacLeod et al. (2012).

We calculate $\sigma_{\mathrm{rms}}^{2}$ and error($\sigma_{\mathrm{rms}}^{2}$) for each quasar in all SDSS bands. For sources with variability statistically non-detected in one band, i.e., with S/N($\sigma_{\mathrm{rms}}^{2}$) = $\sigma_{\mathrm{rms}}^{2}/\emph{error}(\sigma_{\mathrm{rms}}^{2})<2$, we adopt $2\times$error($\sigma_{\mathrm{rms}}^{2}$) as the upper limit to $\sigma_{\mathrm{rms}}^{2}$. Distribution of the derived $\sigma_{\mathrm{rms}}$ in $g$ band is illustrated in Fig. 1, where we can see that the ESQs are the opposite counterparts of the EVQs in the parameter space of variability amplitude.

We consider a series of thresholds of $\sigma_{\mathrm{rms}}$ to define ESQs. Any source would be identified as ESQs if

1. the $g$-band $\sigma_{\mathrm{rms}}$ is less than the threshold, and

2. the $u$-, $r$-, $i$-, and $z$-band $\sigma_{\mathrm{rms}}$ are all either less than the threshold or non-detected.

Using thresholds of $\sigma_{\mathrm{rms}}$ $\leq$ 0.02, 0.03, 0.04, and 0.05 magnitudes, we identify 3, 25, 53, and 136 ESQs, respectively. In Fig. 2 we illustrate a typical $g$-band light curve of an ESQ, in comparison with a normal quasar and an EVQ. Note the SDSS provides PSF magnitude for quasars which is unbiased and optimal for point-like sources (Stoughton et al. 2002). The host contamination is expected to be weak for luminous quasars, especially for our ESQs which tend to have even higher luminosities compared with all SDSS quasars (see Fig. 3).

3. Excess of radio emission in ESQs

Figure 3.— Normalized distributions of redshift, bolometric luminosity, BH mass, and Eddington ratio of all ESQs (red) selected with threshold of $\sigma_{\mathrm{rms}}\leq 0.05$, compared to those of the control quasars (CQ, blue) and all quasars in the SDSS DR16Q with valid measurements (black). All values are taken from Wu & Shen (2022).

Figure 4.— The radio (FIRST/VLASS) detection fractions of ESQs (filled or open squares) selected with a series of thresholds of $\sigma_{\mathrm{rms}}$, compared with control samples of quasars (CQs; filled or open circles), matched in redshift, magnitude, and BH mass. Numbers of ESQs and CQs (denominator) and of radio-detected ones (numerator) are marked. Error bars are $1\sigma$ confidence limits based on a combination of Poisson and binomial statistics (Gehrels 1986). See Section 2 and Appendix B for details of selections of the ESQs and control samples, respectively.

As shown in Fig. 3, the selected ESQs span broad ranges of redshift, luminosity, BH mass, and Eddington ratio. To erase potential observational biases, we build control samples of quasars (CQs), with redshift ($z$), $g$-band apparent magnitude, and BH mass ($M_{\mathrm{bh}}$), matched to our ESQs. Note the simultaneous match of redshift, apparent magnitude, and BH mass automatically matches the luminosity and Eddington ratio. See Appendix B for details of the control sample selection.

Figure 5.— Distributions of the radio loudness of the radio-detected (upper: FIRST detected; lower: VLASS detected) ESQs and CQs. Error bars are 1$\sigma$ confidence limits for Poisson statistics. The $P_{\mathrm{KS}}$ is the p-value of the Kolmogorov-Smirnov test between ESQs and CQs.

We match our ESQ and the control samples to the FIRST (1.4 GHz, White et al. 1997) and VLASS epoch 1 (3 GHz, Gordon et al. 2021) catalogs, with a matching radius of 5″ (see Appendix C). The matched results are displayed in Fig. 4. The radio detection fraction of ESQs is remarkably high, $>$ 25%, dependent on the $\sigma_{\mathrm{rms}}$ threshold adopted to identify ESQs. Moreover, the FIRST/VLASS detection fractions of the ESQs are significantly higher than that of the control samples ($\sim$ 6.7% – 8.4%), and the difference gradually decreases towards higher threshold of $\sigma_{\mathrm{rms}}$. This indicates that the more stable an ESQ is, the more likely it would be detected in radio. The clear excess of radio emission in ESQs is further confirmed with $\sigma_{\mathrm{rms}}$ of SDSS quasars measured with Pan-STARRS1 (Flewelling et al. 2020) and Gaia DR3 (Gaia Collaboration et al. 2023) time domain photometry (W. K. Kang et al. 2024 and H. C. Wang et al. 2024, in preparations). In this work we focus on ESQs selected using the SDSS Stripe 82 light curves, which have $\sim$ 10 years time span with $\sim$ 60 visits per band, and are significantly longer/more than those of Pan-STARRS1 and Gaia, thus the variability could be more stringently constrained for ESQs.

We further confirm that the ESQs we selected are type 1 quasars with at least one broad emission line significantly detected, and the dominant fraction of them have UV-to-optical color consistent with their control quasars (see Appendix D). These factors, together with their high luminosity and Eddington ratio (see Fig. 3), indicate their extremely low variability is not dominated by strong host contamination. While it is worthy of further exploring whether redder colors of a small fraction of ESQs are intrinsic or due to obscuration, excluding these redder ESQs does not alter our results (see Appendix D).

Figure 6.— Upper: correlation between radio loudness $R$ and $g$-band $\sigma_{\mathrm{rms}}$ for the radio-detected (left: FIRST detected; right: VLASS detected) quasars in the SDSS Stripe 82. The red and blue lines represent the best-fit bisector linear regression (Isobe et al. 1990) and the corresponding 3$\sigma$ confidence interval, respectively. Shown in the upper-left corners are the bisector regression slope ($s$ with the corresponding 1$\sigma$ error), the Pearson’s rank correlation coefficient ($r$), and the significance level ($P_{\mathrm{null}}$). Numbers of radio-detected quasars are shown in the lower-right corners. Upper limits of $\sigma_{\mathrm{rms}}$ for minor quasars with ${\rm S/N}(\sigma^{2}_{\mathrm{rms}})<2$ are indicated by the leftward arrows. Lower: histogram distributions of the $g$-band $\sigma_{\mathrm{rms}}$ for the radio-detected quasars (RD) in the SDSS Stripe 82, compared with the corresponding control samples of quasars (CQ). Error bars to the histogram are 1$\sigma$ confidence limits for Poisson statistics. Shown in the upper-right corners are numbers of radio-detected quasars with matched control ones (These numbers are different from ones in upper panels because part of radio-detected quasars can not find control quasar). The $P_{\mathrm{KS}}$ is the p-value of the Kolmogorov-Smirnov test between the radio-detected sample and the control sample.

4. Discussion

Figure 7.— Left: similar to the lower left panel of Fig. 6, but here we only plot the $g$-band $\sigma_{\mathrm{rms}}$ distribution for the FIRST detected quasars (RD) with $R<60$ and the corresponding control sample (CQ). We omit the Poisson error bars to avoid confusion. We then plot the expected distribution for a single normal variable quasar using mock light curves, considering the effect of stochastic nature of variation, limited length of light curve, and sparse sampling (see Appendix E). Right: similar to the left panel, we here plot the simulated distribution illustrating the effect of the various lengths of the mock light curves (3 years, 10 years, and 30 years).

As aforementioned in §1, the RL quasars are known to show (marginally) stronger variability in UV/optical compared with the RQ ones, likely due to the jet contribution. Contrarily, the prominent radio excess in ESQs reported here is indeed out of expectation. To determine the radio physical properties of these ESQs and the relationship to the general quasars, we adopt the canonical definition of the radio loudness, $R=f_{\rm 5GHz}/f_{2500{\rm\mathring{A}}}$, where $f_{2500{\rm\mathring{A}}}$ is the rest-frame 2500${\rm\mathring{A}}$ flux density from Wu & Shen (2022), and $f_{\rm 5GHz}$ is the 5 GHz flux density inferred from the observed-frame 1.4 GHz (FIRST) or 3 GHz (VLASS) flux densities assuming $f_{\nu}\propto{\nu}^{-\alpha}$ with a radio spectral index of $\alpha$ = 0.5 (Jiang et al. 2007, see also Appendix C). Fig. 5 plots distributions of the radio loudness for the radio-detected ESQs and CQs. We see the excess of the radio-detected ESQs mainly occurs in the radio intermediate regime ($R\sim$ 10 – 60) which clearly requires jet. The excess is also visible at $R\sim$ 1 – 10 which falls in the radio quiet regime (likely due to weaker jet in those sources). However, the excess shows a sharp cutoff and disappears in the radio loud regime ($R>60$). While we see no statistical difference between ESQs and CQs at $R>60$ (partially due to the small sample size), the lack of excess ESQs at $R>60$ is statistically significant (compared with that at $R<60$, with a p-value of $\sim$ 0.001), probably owing to strong jet contribution to the observed UV/optical variability.

We further explore the connection between radio emission and UV/optical variability in quasars from a different perspective. In Fig. 6 we plot the $g$-band $\sigma_{\mathrm{rms}}$ of the radio-detected quasars in the SDSS Stripe 82. We similarly build the control sample for them with matched redshift, $g$-band apparent magnitude, and $M_{\mathrm{bh}}$. Since we would like to measure $\sigma_{\mathrm{rms}}$ for the control sample as well, the control sample is selected from the SDSS Stripe 82 only. Due to the small sample size, for each radio-detected quasar we select only one control quasar. We remove a small fraction of the radio-detected quasars for which we can not find a control quasar in Stripe 82.

In the upper panels of Fig. 6, the clear correlation between the radio loudness and $g$-band $\sigma_{\mathrm{rms}}$ for the radio-detected quasars could be due to increasing jet contribution to the observed variability with increasing $R$. This may again explain the lack of excess ESQs at much larger $R$ (Fig. 5), as ESQs could only be picked out of jetted quasars with low to intermediate $R$ for which jet contribution to the observed variability amplitude is weak.

As shown in the lower panels of Fig. 6, the radio-detected quasars have on average similar $g$-band $\sigma_{\mathrm{rms}}$ compared with the control sample. For the FIRST- and VLASS-detected quasars and CQs, the p-values of the KS test are 2.2% and 13%, respectively. However, we see a clear excess of the radio-detected quasars at $\sigma_{\mathrm{rms}}$ $<$ 0.05 compared with the control sample (i.e., 57 ESQs vs 24 CQs for the FIRST detection, with a p-value of 0.0001, and 51 ESQs vs 24 CQs for the VLASS detection, with a p-value of 0.0009). This is consistent with the pattern we have shown above, and suggests a link between jet production and ESQs.

To avoid strong jet contribution to the observed $g$-band $\sigma_{\mathrm{rms}}$, in Fig. 7 we plot the FIRST-detected quasars similar to the lower left panel of Fig. 6 but excluding sources with $R>60$. The statistical difference between the FIRST-detected quasars and the corresponding control sample is now evident (with a KS test p-value of 0.004), and the excess of the FIRST-detected sources with $g$-band $\sigma_{\mathrm{rms}}<0.05$ remains significant (i.e., 46 ESQs vs 18 CQs, with a p-value of 0.0002).

Before discussing the physical implication of the radio excess in ESQs, there are two notable questions to be addressed: 1) why a large portion of ESQs are radio non-detected (Fig. 4), and 2) why some normal quasars could also have very small $g$-band $\sigma_{\mathrm{rms}}$ (the control sample in Fig. 6)?

For the first question, we median-stack the radio images of radio non-detected ESQs and CQs (see Liao et al. 2022, for the stacking approach), and find higher median radio flux densities for ESQs than the control samples. Based on the FIRST (VLASS) images, the median radio flux density for ESQs is 92${}^{+24}_{-18}$ $\mathrm{\mu}$Jy (76${}^{+19}_{-20}$ $\mathrm{\mu}$Jy) and larger than 59${}^{+5}_{-1}$ $\mathrm{\mu}$Jy (50${}^{+5}_{-4}$ $\mathrm{\mu}$Jy) for the control sample. Though statistically marginal due to the small sample size of ESQs, this suggests that some of the radio non-detected ESQs also exhibit excess of radio emission. Much deeper radio images are desired to detect those radio fainter ESQs, and address whether all true ESQs have jets.

For the second question, one factor we need to consider is the stochastic nature of the variation, that a normal variable quasar may exhibit large variation amplitude at one time, but much weaker variation at a different time. The sparse sampling and limited length of light curves could also play a role in magnifying such observational effect. We build mock light curves for the normal variable quasars assuming the damped random walk (DRW, Kelly et al. 2009) and non-DRW processes (see Appendix E). The simulations confirm that the normal variable quasars could temporarily appear as “ESQs” due to the stochastic nature of the variation and/or the limited length of light curve, thus contaminate the selection of ESQs (see Fig. 7). The results of our simulations could also explain the drop of the radio-detected fraction of ESQs with increasing the threshold of $\sigma_{\mathrm{rms}}$ (see Fig. 4), as the contamination from the normal variable quasars to ESQs is stronger at larger $\sigma_{\mathrm{rms}}$. Our simulations also show that longer light curves are essential to better distinguish real ESQs from the normal variable quasars (see Fig. 7 and Appendix E).

Now, it is intriguing to investigate how the prominent radio excess in ESQs reported here could be fitted in the most popular scenarios for the observed broad range of jet production efficiency in quasars and AGNs, i.e., the spin paradigm (e.g. Wilson & Colbert 1995; Sikora et al. 2007). It is commonly believed that the strong magnetic field is a key ingredient in the jet launching (e.g. Blandford & Znajek 1977; Blandford & Payne 1982; Livio et al. 1999; Blandford et al. 2019), and the UV/optical variations in quasars could be driven by the disc turbulence induced by the magneto-rotational instability (e.g. Kelly et al. 2009). As proposed by Cai et al. (2019) presenting a tentative evidence for the more stable inner accretion disc in the RL quasars from the optical color variability study, could the inner discs in the RL quasars have been stabilized by the strong magnetic field (e.g. Begelman & Pringle 2007; Zheng et al. 2011; Li & Begelman 2014; Sadowski 2016)?

Moreover, apart from the prominent radio excess in ESQs, another puzzle is why ESQs are so rare. Only $\sim$ 1.5% of the quasars in the SDSS Stripe 82 (i.e., 135 out of 9146, based on the $\sigma_{\mathrm{rms}}$ threshold of 0.05 mag in all five SDSS bands) have been selected as ESQs, and only $\sim$ 7% of the radio-detected quasars are ESQs (i.e., 34 out of 517 FIRST-detected ones, and 34 out of 499 VLASS-detected ones). Assuming the magnetic field is a key factor on understanding the radio excess in ESQs as well as the rarity, we discuss below possible interpretations that could be further explored from both observational and theoretical aspects.

The rarity of ESQs may suggest that suppressing AGN variability by strong magnetic field in radio quasars could only occur in a minor fraction of sources (see Figs. 6 & 7). One possibility is that the suppressing is only significant with sufficiently strong magnetic field. This seems nicely in line with the theoretical analyses that the magneto-rotational instability (MRI) in the accretion disc could be stabilized beyond a critical magnetic field. Both Pessah & Psaltis (2005) and Das et al. (2018) showed that a toroidal field with $\beta=8\pi P_{\rm gas}/B^{2}<1$ (where $P_{\rm gas}$ is the gas pressure and $B$ the magnetic field) could suppress the MRI, while Bai & Stone (2013) found a critical $\beta$ $<$ 100 for a poloidal field to suppress the MRI. Therefore, we speculate that ESQs could be a minor population with magnetic field above a critical value (or $\beta$ below a critical value) in the inner disc. Note that the jet core-shift has been widely used to measure the magnetic field in jet at pc scale (e.g. Lobanov 1998; Zamaninasab et al. 2014). While a quantitative ratio between magnetic field in the inner disc and that in the pc-scale jet is unclear, it would be helpful to investigate with future core-shift observations whether ESQs have stronger pc-scale magnetic field compared with the control sample.

In the spin diagram, the jet power, $P_{\rm jet}$, depends on both the black hole spin and the poloidal magnetic field ($P_{\rm jet}\sim\Omega_{H}^{2}B_{p}^{2}$, where $\Omega_{H}$ is the angular velocity at the black hole horizon and $B_{p}$ the poloidal magnetic field; e.g., Livio et al. 1999). Under this scheme, the proposed beyond-critical magnetic field could naturally explain the excess of the radio emission in ESQs. Meanwhile, to explain their intermediate radio loudness ($R\sim 10-60$), ESQs should have relatively small black hole spins resulting in small $\Omega_{H}$ thus moderate $P_{\rm jet}$. These ESQs would belong to the dominant radio quiet population if without the beyond-critical magnetic field. Contrarily, some ESQs with both high spin and beyond-critical poloidal magnetic field could be very radio loud. However they could be even rarer compared with the radio intermediate ESQs, and their optical/UV variability could be significantly elevated due to the jet contribution, thus missed by our selection of ESQs. On the other hand, if there are ESQs with strong toroidal field but too weak poloidal field, jet launching would not be expected (Beckwith et al. 2008), which may also explain some of the radio non-detected ESQs.

Overall, our results and the above plausible interpretations make ESQs an essential and unique population under potentially extreme condition, which strongly necessitates extensive theoretical and observational follow-up studies. The quantitative relation between strong magnetic field suppressing and the observed AGN variability is to be established theoretically. Much deeper radio images could draw a panorama of the radio loudness of ESQs. Multi-band radio SEDs and high-resolution radio images could reveal the physical nature of the jets in ESQs, to verify the existence of strong magnetic field. In additional to follow-up radio observations, our work could also makes ESQs unique targets for the time domain astronomy, because the longer duration/higher cadence sampling of time domain observations, the better that stableness of a quasar could be constrained. Studying the relation between ESQs and other multi-band observed parameters of quasars may also yield new hints (see Appendix F, Figs. 15 and 16 for instance). Extending the study (of the connection between disc stableness and jet) to the low mass regime, i.e, stellar mass BHs, is also alluring.

Appendix A A: Weighted excess variance

Figure 8.— Distribution of the simulated $\sigma_{\mathrm{rms}}^{2}$ (zero intrinsic variance is assumed) using different weighted factors.

We calculate the revised $\sigma_{\mathrm{rms}}^{2}$ for a quasar light curve with $N$ epochs as

$$\sigma_{\mathrm{rms}}^{2}=\frac{1}{\sum\sigma^{w}_{i}}\sum\sigma^{w}_{i}[\frac{N}{N-1}(X_{i}-\bar{X})^{2}-\sigma_{i}^{2}],$$ (A1)

where $X_{i}$ the observed magnitude at $i$ epoch, $\sigma_{i}$ is the photometric uncertainty at $i$ epoch, and $\bar{X}$ is the weighted average magnitude in which the weight is $\sigma^{2}_{i}$. We consider three different weighted factors $\sigma_{i}^{w}$, with $w=$ 0, -2, and -4, respectively. When $w$ = 0, the equation returns to the canonical form (Vaughan et al. 2003). The smaller (more negative) $w$ is, the smaller the influence of the observed epochs with large photometric uncertainties to the variability amplitude.

Utilizing the SDSS $g$ band as an example, we perform Monte Carlo simulations to select the optimal weighted factor for the selection of ESQs with very weak intrinsic variability. Suppose a quasar with zero intrinsic variability and $g$ = 19.5 has been observed $n$ times, we simulate its magnitude at each epoch as:

$$g_{i}=g+Gau(0,\sigma_{i})$$ (A2)

where $\sigma_{i}$ refers to the photometric uncertainty. We adopt $n$ = 60 which is the average observation number for each source in the SDSS Stripe 82 (MacLeod et al. 2012). To mimic the true observational effects, $\sigma_{i}$ is randomly selected from the real data (out of all photometric measurements of the SDSS Stripe 82 quasars with observed $g=19.5\pm 0.1$ at any epoch). Distribution of the simulated $\sigma_{\mathrm{rms}}^{2}$ is displayed in Fig. 8. Clearly, the canonical form of $\sigma_{\mathrm{rms}}^{2}$ ($w$ = 0) yields the largest scatter, while $w$ = -4 behaves slightly better than $w$ = -2. Simulations assuming various $g$ band magnitudes yield similar results. Thus in this work, we adopt $w$ = -4 to calculate $\sigma_{\mathrm{\mathrm{rms}}}^{2}$.

We also perform similar Monte Carlo simulations to derive the statistical uncertainty of a measured $\sigma_{\mathrm{rms}}^{2}$, i.e., error($\sigma_{\mathrm{rms}}^{2}$), assuming a quasar has zero intrinsic variability and estimating the expected standard deviation of $\sigma_{\mathrm{rms}}^{2}$ (such as the scatter of the mock $\sigma_{\mathrm{rms}}^{2}$ plotted in Fig. 8) simply due to photometric errors. Note this approach is only valid for sources with extremely weak variability, but sufficient for this study.

Appendix B B: build the control samples

For each ESQ, we randomly select 10 control quasars out of the SDSS DR 16 quasar catalog (Lyke et al. 2020) with $\Delta z/z$ $\leq$ 10%, $\Delta g$ $\leq$ 0.2 mag, and $\Delta\log M_{\mathrm{bh}}$ $\leq$ 0.4 dex. The reason we select 10 (instead of one) CQs for each ESQ is to reduce Poisson noise of the control sample. Note while the VLASS fully covers the SDSS survey, the FIRST does not. We only select CQs within the FIRST footprint, and exclude one ESQ which is out of the FIRST footprint from Fig. 4 when presenting the FIRST detection fraction. We further drop one more ESQ which has less than 5 CQs available. For the two other ESQs with at least 5 but less than 10 CQs, their available CQs are repeatedly and randomly selected to make up 10 CQs.

We further note the FIRST sensitivity is slightly deeper in a small region along the equatorial strip (RA = 21.3 to 3.3 hr, Dec = -1 to 1 deg, with a typical detection threshold of 0.75 mJy, instead of 1.0 mJy for the rest dominant FIRST survey area). The deeper strip is right within the SDSS Stripe 82 field. While comparing the FIRST detection fraction of ESQs with the control samples, we need to correct the effect of the non-uniform depth of the FIRST survey. To do so, for any ESQ within the deeper strip, we select its CQs from the DR16Q but with $g$-band flux density brighter by a factor of 1/0.75 to compensate the effect (slightly brighter CQs are needed to be detected in slightly shallower FIRST image). We also request the control quasar to have BH mass larger by a factor of 1/0.75 to ensure the match of Eddington ratio, as the RL fraction could be dependent on the Eddington ratio. After this correction, the FIRST detection fraction of the CQs only slightly increases (by around 1%) and has a negligible effect to the results presented in this study. Note the results we present in Fig. 4 are already after this correction. Furthermore, the FIRST detection fraction of our ESQs in and out the deeper FIRST strip exhibits no significant difference, also indicating the effect of the non-uniform FIRST depth is negligible.

Appendix C C: Cross-match the SDSS DR16Q with the radio catalogs

Figure 9.— The observed offset distribution of cross-matching the SDSS DR16Q catalog with the FIRST (upper panel) and VLASS catalogs (lower panel). We adopt 5″ as the cross-matching threshold (blue dashed lines). The red lines show the expected number of chance alignments with background sources as a function of offset between SDSS and radio sources.

Figure 10.— Radio flux densities from the VLASS epoch 1 (version 3) vs epoch 2 (version 2). Both of them have been converted to the rest-frame 5 GHz assuming a spectral index of $\alpha$ = 0.5. Numbers of the VLASS-detected sources are marked. For a few sources detected in only one epoch, the upper limits are given for the other epoch.

Figure 11.— The FIRST cutout images (60″$\times$60″) of the 34 FIRST-detected ESQs. In one of these cutouts images, i.e., the fourth column from the left and the third row from the bottom, there seems to be two radio counterparts. However only one radio source (in the center of the cutout, $\simeq 0.1^{\prime\prime}$ offset from SDSS position) is detected in the FIRST catalog, and the other point-like source to the slightly left is likely due to jet structure.

Figure 12.— The VLASS epoch 1 cutout images (11″$\times$11″) of the 34 VLASS-detected ESQs.

Figure 13.— Left: Rest-frame 5GHz flux density of ESQs and CQs (detected by both FIRST and VLASS) derived from the observed frame 3 GHz VLASS flux density ($y$-axis) versus that derived from the 1.4 GHz FIRST flux density ($x$-axis), both assuming a spectral index of $\alpha$ = 0.5. Middle and Right: Rest-frame 5GHz flux density of ESQs and CQs derived through interpolating FIRST and VLASS flux densities ($y$-axis) versus that from FIRST and VLASS ($x$-axis of middle right panel, respectively, both assuming spectral index of $\alpha$ = 0.5).

We cross-match the SDSS DR16Q with both the FIRST and VLASS catalogs, using a matching radius of 5$\arcsec$. Following Reines et al. (2020), in Fig. 9 we plot the offset distribution from cross-matching the SDSS DR16Q catalog with the radio catalogs. As indicated by the red lines, the expected fraction of chance alignments with background sources for the radio-detected quasars below 5″ offset is 4.6% (1082 out of 23458 FIRST detected ones) and 4.3% (1104 out of 25782 VLASS detected ones), respectively. Our Fig. 9 is similar to Fig. 2 of Reines et al. (2020).

For the VLASS, we adopt the epoch 1 Quick Look catalog (Gordon et al. 2021). Utilizing the epoch 2 catalog does not alter our results. Meanwhile, comparing the epoch 1 and epoch 2 VLASS flux densities of our ESQs or the control sample reveals no systematic variation trend (see Fig. 10).

In Fig. 11 (and Fig. 12) we illustrate the FIRST (VLASS) images of the ESQs with FIRST (VLASS) detection. We also compare the rest frame 5 GHz flux densities derived from the FIRST or VLASS assuming a spectral index of 0.5, and from interpolating the FIRST and VLASS flux densities (Fig. 13). We find that these approaches yield generally consistent rest frame flux densities. Fig. 13 also indicates that the observed frame 1.4 – 3 GHz spectral slopes of ESQs and their control sample exhibit no systematic difference.

Appendix D D: the color of ESQs vs the control samples

Figure 14.— Upper: the $u$-$z$ color of our ESQs (selected with a threshold of $\sigma_{\mathrm{rms}}\leq 0.05$) versus that of the control sample. Red circles are ESQs while green squares and diamonds are the FIRST and VLASS detected ESQs, respectively. The statistical errors to the colors due to photometric uncertainties are generally too small to be displayed. For each ESQ we have 10 CQs, and we also plot the mean color (and standard deviation) of the 10 CQs for comparison. Middle: the histogram distribution of the color difference between ESQs and the control quasars. Error bars are 1$\sigma$ confidence limits for Poisson statistics. The $P_{\mathrm{KS}}$ shown in the upper-left corner are the p-values of the KS test between ESQs and FIRST detected ESQs, and between ESQs and VLASS detected ESQs, respectively. Lower: the radio detection fraction of ESQs versus the color difference. Error bars are 1$\sigma$ confidence limits based on a combination of Poisson and binomial statistics.

In Fig. 14 we plot the $u$-$z$ color of ESQs compared with the control samples. While a considerable fraction of ESQs do exhibit redder $u$-$z$ color, the radio detection fraction in ESQs appears independent to the color.

Appendix E E: contamination of normal quasars to ESQs

We perform simulations to assess the contamination to ESQs by the normal variable quasars appearing as “ESQs” as a result of effects of the stochastic nature of variation (i.e., the DRW and non-DRW processes with distinct damping timescales), the limited length of light curve (i.e., 3 years, 10 years, and 30 years), and the sparse sampling. In order to compare with the $g$-band $\sigma_{\rm rms}$ distributions of the FIRST detected quasars with $R<60$ and the corresponding control sample, we implement their real $g$-band samplings in our simulations. Both of them have 303 quasars.

For the DRW process, light curves are simulated using the procedure of Kelly et al. (2009) with three typical damping timescales ($\tau_{\rm DRW}=$ 100 days, 300 days, and 1000 days) and a long-term variation amplitude ($\sigma_{\rm l}=0.1$ mag) for the normal variable quasars. For the non-DRW process, light curves are simulated in terms of the algorithm of Timmer & Koenig (1995) assuming a broken power-law power spectral density (PSD). The non-DRW PSD breaks at a frequency of $(2\pi\tau_{\rm non-DRW})^{-1}$ with $\tau_{\rm non-DRW}=1000$ days, and has a low-frequency slope $\alpha_{\rm l}=-1.3$ suggested by Guo et al. (2017a) and a high-frequency slope $\alpha_{\rm h}=-2$. All light curves are generated in steps of 0.1 day. The length of the DRW light curves are set to 30 years, while $\simeq 90$ years for the non-DRW light curves, whose long-term variation amplitudes are $\sigma_{\rm l}=0.1$ mag as well. A longer length for the non-DRW process is to fully take the diverse power at low frequencies into account; and a length longer than $\simeq 90$ years for the non-DRW process does not alter the results implied by $\simeq 90$-year light curves.

For each quasar, we generate 1000 light curves, which are coupled with its real $g$-band sampling. Specifically, each simulated light curve is linearly interpolated at the observed epochs of a quasar, and the interpolated magnitudes are fluctuated by random Gaussian deviates scaled to the observed photometric uncertainties of the quasar. The length of the observed light curves for quasars is 10-year long. The 10-year real observation is repeated 3 times to mimic a 30-year observation, while the last 3-year of the 10-year real observation mocks a 3-year observation.

In Fig. 7 we compare the $g$-band $\sigma_{\rm rms}$ distributions for both the RD quasars with $R<60$ and the corresponding CQ with those simulated for a single normal variable quasar, considering different damping timescales of the DRW process, the effect of the non-DRW process, and the sampling lengths. For each of these effects, the simulated “absolute” distributions of the $g$-band $\sigma_{\rm rms}$ also depend on the adopted $\sigma_{\rm l}$. However, adopting different input $\sigma_{\rm l}$ would just horizontally shift the output distribution of the $g$-band $\sigma_{\rm rms}$ without changing its shape. The similarity of the “relative” distributions for various input $\sigma_{\rm l}$ is confirmed by simulations using a series of $\sigma_{\rm l}$, ranging from 0.01 to 0.4, well blanketing the observed values. Therefore, we adopt the typical $\sigma_{\rm l}=0.1$ for reference and normalized (and shifted horizontally) the resultant distributions of the $g$-band $\sigma_{\rm rms}$ to around the peak of the observed CQ distribution (Fig. 7).

As expected, longer damping timescale, the non-DRW process rather than the DRW one, and shorter sampling length can all increase the probability of a normal variable quasar temporarily appearing as an “ESQ”. As shown in Fig. 7, there are 46 RD quasars and 18 CQs with $g$-band $\sigma_{\rm rms}<0.05$. Assuming the most broad distribution simulated by the non-DRW process with $\tau_{\rm non-DRW}=1000$ days and $\alpha_{\rm l}=-1.3$, $\sim 13\%$ of the RD quasars and $\sim 35\%$ of the CQs are likely normal variable quasars temporarily appearing as “ESQs”. Note that the fractions of contamination by the normal variable quasars are hard to assess as the exact variation model (i.e., DRW vs non-DRW), the model parameters (e.g., the damping timescale and the low frequency PSD slope) and the variance from source to source (currently omitted in the simulations) are yet poorly constrained (e.g. Guo et al. 2017b). Even longer damping timescale (which is very likely as the measured timescale could have been under-estimated due to the limited length of light curves, e.g., Hu et al. 2023) would yield higher contamination by the normal variable quasars (see the left panel of Fig. 7). Meanwhile, longer light curves appear essential to better distinguish real ESQs from the normal variable quasars (see the right panel of Fig. 7).

Appendix F F: Other factors

Figure 15.— Distribution between the equivalent width ($\rm EW$) of broad Mg II line and the $g$-band $\sigma_{\mathrm{rms}}$. Blue lines and grey points are results from Kang et al. (2021). The red circles are our ESQs. The green squares and diamonds are the FIRST and VLASS detected ESQs, respectively.

Figure 16.— Distribution between the X-ray loudness of 0.5 - 12 keV (from the 4XMM-DR13 catalog; Webb et al. 2020) and the $g$-band $\sigma_{\mathrm{rms}}$. The blue lines and grey points are results from Kang et al. (2018). The red points are our ESQs. The green squares and diamonds are the FIRST and VLASS detected ESQs, respectively.

In the manuscript we have controlled the effects of redshift, luminosity, BH mass, and Eddington ratio while comparing ESQs with the control sample. Here we present attempts to explore the effects of other potential factors.

In Kang et al. (2021), we have found that the UV/optical variability amplitude of quasars positively correlates with the equivalent widths (EWs) of C IV, Mg II, and [O III]5007 emission lines. In Fig. 15, we plot our ESQs which have EWs of broad Mg II line from the SDSS DR16Q catalog (Wu & Shen 2022), compared with quasars from Fig. 4 of Kang et al. (2021). Here the linear regression is obtained taking $g$-band $\sigma_{\mathrm{rms}}$ as the independent variable, used to illustrate the expected $EW_{\mathrm{MgII}}$ at given $g$-band $\sigma_{\mathrm{rms}}$. We see that the ESQs are consistent with the correlation that was fit to the Kang et al. (2021) sample, suggesting that compared with the normal variable quasars the ESQs are not distinct in their emission line strength.

The UV/optical variability of quasars was also found be correlate with the X-ray loudness of quasars (Kang et al. 2018). We cross-match our ESQs with the XMM-DR13 catalog (Webb et al. 2020), and compare our ESQs with the normal variable quasars from Kang et al. (2018) in the X-ray loudness versus $g$-band $\sigma_{\mathrm{rms}}$ (Fig. 16). Similarly we find that the ESQs are consistent with the correlation that was fit to the Kang et al. (2018) sample. Meanwhile, the radio-detected ESQs tend to be X-ray louder (statistically yet insignificant due to the small number of ESQs with X-ray coverage).