Searching for quasar candidates with periodic variations from the Zwicky Transient Facility: results and implications

We conduct a systematic search for quasars with periodic variations from the archival photometric data of the Zwicky Transient Facility by cross-matching with the quasar catalogs of the Sloan Digital Sky Survey and V{\'e}ron-Cetty \&V{\'e}ron. We first select out 184 primitive periodic candidates using the generalized Lomb-Scargle periodogram and auto-correlation function and then estimate their statistical significance of periodicity based on two red-noise models, i.e., damped random walk (DRW) and single power-law (SPL) models. As such, we finally identify 106 (DRW) and 86 (SPL) candidates with the most significant periodic variations out of 143,700 quasars. We further compare DRW and SPL models using Bayes factors, which indicate a relative preference of the SPL model for our primitive sample. We thus adopt the candidates identified with SPL as the final sample and summarize its basic properties. We extend the light curves of the selected candidates by supplying other archival survey data to verify their periodicity. However, only three candidates (with 6-8 cycles of periods) meet the selection criteria. This result clearly implies that, instead of being strictly periodic, the variability must be quasi-periodic or caused by stochastic red-noise. This exerts a challenge to the existing search approaches and calls for developing new effective methods.


INTRODUCTION
Periodic variability of quasars has attracted great attention over the past decade because of its possible connection with supermassive black hole binaries on the theoretic side (e.g., Artymowicz & Lubow 1996;D'Orazio et al. 2015) and significant advances towards modern time-domain surveys on the observational side.There had been more than one hundred periodic quasar candidates reported from systematic searches over large surveys, such as from Catalina Real-time Transient Survey (CRTS, Graham et al. 2015a), Palomar Transient Factory (PTF, Charisi et al. 2016), Pan-STARRS1 Medium Deep Survey (PS1 MDS, Liu et al. 2019), and from a combination of the Dark Energy Survey (DES) and Sloan Digital Sky Survey (SDSS; Chen et al. 2020).Several individual quasars were serendipitously found to exhibit (quasi-)periodic variations using long-term archival data (e.g., Valtonen et al. 2008;Li et al. 2016aLi et al. , 2019;;Zhang et al. 2020;O'Neill et al. 2022;Zhang et al. 2021).Those candidates constitute a precious guiding sample for investigating the physical origins of the periodicity and exclusively identifying supermassive black hole binaries.However, we bear in mind that the periodicity might be subject to false positives caused by red-noise stochastic variability of quasars (e.g., Vaughan et al. 2016).A sophisticated estimation of the statistical significance of the periodicity in the presence of red-noise variations is highly required.Meanwhile, a systematic search over new emerging time-domain surveys with improved sampling rates and photometry would provide more valuable candidate samples to testify to the periodicity.
The Zwicky Transient Facility (ZTF) is a new time-domain survey dedicated to a systematic exploration of the northern optical transient sky with a 47 square degree field of view camera equipped on the 48inch Samuel Oschin Telescope at Palomar Observatory (Bellm et al. 2019;Graham et al. 2019).It started operation in 2018 and scans the entire northern visible sky with a typical cadence of three days.Such a moderate sampling cadence makes ZTF well-suitable for exploring quasar periodicity.In this work, we conduct a systematic search for periodic quasar candidates based on ZTF photometric data by crossmatching with the quasar catalogues of the SDSS and Véron-Cetty & Véron (2010).
The paper is organized as follows.In Section 2, we describe the search methodology and estimation of the statistical significance of the periodicity.In Section 3, we briefly summarize the basic properties of our selected periodic quasar candidates and compare our sample with those reported in previous works.In Section 4, we compile other available survey photometric data to extend the temporal baselines to verify the periodicity, and then present a brief discussion on possible explanations for the periodicity and implications of our periodic search.Finally, we summarize the main results in Section 5. Throughout this work, we use a cosmology with Ω M = 0.3, Ω Λ = 0.7, and  0 = 70 km s −1 Mpc −1 .

Initial Quasar Catalogues and ZTF Light Curve Data
We initiate with the quasar catalogues of SDSS DR14 and Véron-Cetty & Véron (2010).The former contains 526,356 spectroscopically identified quasars (Pâris et al. 2018) and the latter contains 168,941 quasars compiled from all known quasars at that time.To discard faint objects, we respectively set a magnitude limit of  ≤ 20 for SDSS DR14 catalogue and  ≤ 20 for Véron-Cetty & Véron (2010) catalogue.These limits are slightly brighter than the ZTF depth of ∼ 20.5 mag (Graham et al. 2019).We then perform a crossmatch within a radius of 5 arcsecs between the two catalogues and remove duplicated objects using the package esutil 1 .We combine the left objects in the two catalogues and obtain a list of 223,061 quasars.
Next, we retrieve the above initial quasar list from the ZTF database using a searching radius of 3 arcsecs.We use the 11th ZTF public data release2 .There are three custom filters, -, -and -band.The photometry was reduced with the Infrared Processing and Analysis Centre pipeline (Masci et al. 2019).The -band generally has the most data points, therefore, unless stated otherwise, we use the band light curves for the following analysis.
In some light curves there apparently appear a few outlier points.We use a fifteenth-order polynomial to fit the light curves and remove those points outside the 3 deviation from the best fits.The polynomial order is adopted somewhat arbitrarily, but generally, it should be high enough to capture all major variation patterns in the light curves.Since we only focus on long-timescale variability, we bin the light curves within an interval of 20 days.The magnitudes and errors of the binned points are assigned the mean and standard error of the points in each 20-day bin, respectively.This binning operation helps to alleviate the possible biases arising from severely uneven sampling and strong short-timescale fluctuations in some light curves.We further discard those objects with less than 40 binned points and are finally left with 143,700 quasars.Below we use the binned light curves for preliminary periodicity searches in Section 2.2.

Identification Methodology of Periodicity
In this section, we design the following two steps to identify raw periodic quasar candidates.In the next section, we will select out the final candidate sample according to the estimated significance of the periodicity.
First, we employ the generalized Lomb-Scargle periodogram (GLSP; Zechmeister & Kürster 2009) implemented in the pyastronomy3 package.The Lomb-Scargle periodogram (LSP) was proposed by Lomb (1976) and Scargle (1982) and had become a widely used method to detect periodic signals in unevenly sampled data.It is mathematically equivalent to a sinusoidal fit with a form of where  is the amplitude,  is the period, and  is the phase.Compared to the LSP, the GLSP additionally includes an offset term in the sinusoidal fit.The advantage of this generalization is that it provides a more accurate period prediction and a better determination of the spectral intensity (Zechmeister & Kürster 2009).We calculate the GLSP over a period range between 50 days and  span , where  span is the time span of the light curve.We assign the best period  GLS corresponding to the largest periodogram peak and then perform a sinusoidal fit with the period fixed to  GLS .As such, we obtain the amplitude  0 of the best-fit sinusoidal function and the residuals by subtracting the best-fit sinusoidal function from the light curve.Following Horne & Baliunas (1986), we define a signal-tonoise (S/N) ratio as  =  2 0 /2 2  , where   is the standard deviation of the residuals.This ratio measures the power arising from the periodic signal relative to the power from the noise.We require the periodic candidates to satisfy the criterion: 1)  > 4.0 so as to select out the statistically significant peak; and 2)  cycle =  span / GLS > 1.5 to ensure at least 1.5 cycles of periodicity in the light curve.
Second, we adopt the auto-correlation function (ACF) to search for periodicity by combining it with the GLSP.ACF describes the degree of auto-correlation of the light curve at different time lags and has also been widely used for detecting periodic signals (e.g., McQuillan et al. 2013;Graham et al. 2015a;Chen et al. 2020).Theoretically, the ACF of periodically driven stochastic systems is expected to follow an exponentially decaying cosine function (Jung 1993) where  is the time lag and  ACF is the period.We use the method of the interpolated cross-correlation function (ICCF; Gaskell & Peterson 1987;White & Peterson 1994) to calculate the ACF and fit it with the exponentially decaying cosine function.We require that 1) the period difference between the GLSP and ACF methods is within 10 percent, namely, |1 −  ACF / GLS | < 0.1; and 2) the best-fit decay rate  < 10 −3 day −1 , which ensures the ACF decaying no more than 1/ over the temporal baseline (about 1500 days) of the light curve (Graham et al. 2015a).
We obtain 184 periodic candidates that satisfy the above selection criterion (see Tabel 1 for details).In Fig. 1, for the sake of illustration, we show the original and binned light curves, along with the GLSP and ACF of the binned light curve for two selected examples.We note that Vaughan et al. (2016) suggested searching over light curves with at least three-period cycles to distinguish between true candidates for periodicity versus light curves that can be fully described via standard stochastic red-noise processes.For these 184 candidates, the present temporal baselines of ZTF data are not long enough to ensure this criterion (the range of period cycles is from 1.5 to 2.8).In Section 4.2, we will supplement these baselines with earlier archival data from other surveys to extend the temporal baselines for selected candidates.

Estimating Significance of the Periodicity
It is known that red-noise-like variability is commonly seen in quasars/active galactic nuclei (AGNs)    sparse and uneven sampling and large photometric uncertainties (e.g., Vaughan et al. 2016;Krishnan et al. 2021).Considering that we only limit 1.5 cycles of periodicity, it is necessary to appropriately estimate the statistical significance of the periodicity.Below we estimate the false-alarm probability (FAP) and significance based on two popular red-noise models, namely, the damped random walk (DRW) model and the single power-law (SPL) model.Here, the SPL model means the power spectral density (PSD) of the light curve follows a SPL.

The False-alarm Probability with the DRW Model
We adopt the DRW process (Kelly et al. 2009) as the null hypothesis.
The DRW process is widely used to describe AGN stochastic variability (e.g., Kelly et al. 2009;Zu et al. 2011;Li et al. 2013Li et al. , 2016b;;Lu et al. 2019;Zhang et al. 2022) and has a covariance function in the form of where   and   are two times,   is the variation amplitude at long timescale, and   is the characteristic damping timescale.For each periodic candidate, we fit the original unbinned light curve by exploring the posterior probability function where  is the long-term mean of the light curve, (, , ) is the prior probability function, and L is the likelihood function, defined as where   represents the observed light curve.The covariance matrix  is given by where   is the measurement uncertainty at the observation time   ,    is the Kronecker delta, and    is given by Equation (3).We set a uniform prior for the logarithm of the parameters   and   and a uniform prior for .We apply the Markov Chain Monte Carlo (MCMC) sampler emcee5 (Foreman-Mackey et al. 2013) to construct the posterior samples of the parameters.We then randomly draw a set of   and   from their posterior samples and generate a simulated light curve with exactly the same cadence as the observed data.The parameter  is a nuisance for the present purpose as we shift the generated mock light curve to match the observed mean magnitude.We further add Gaussian noises with a zero mean and standard deviation equal to the photometric uncertainties to mimic the measurement errors.Again, we bin the mock light curve every 20 days and apply exactly the same identification methods as in Section 2.2.We repeat the above procedures 100,000 times and calculate the FAP6 as FAP =  p / tot , where  p is the number of periodic candidates in the mock light curves and  tot is the total number of simulated light curves.The FAP of the 184 candidates ranges from 1.0 × 10 −5 to 7.8 × 10 −3 (see the right panel of Fig. 2).We mark a significant periodicity if (1 − FAP) ≥ 99.73%, namely, at least a significance of 3.

The Bayesian Information Criterion of Periodicity with the DRW Model
We have estimated the FAP with the null hypothesis that the light curves follow the DRW process.As an alternative approach, we perform a model selection to test if the data favor an additional periodic signal on top of the background DRW variability.To this end, we adopt the Bayesian information criterion (BIC) defined as (Schwarz 1978) where L is the likelihood function,  is the number of free model parameters, and  is the number of data points.As in Section 2.2, we simply adopt a sinusoidal model for the periodic signal and compare such a DRW+sinusoidal model with the pure DRW model.The pure DRW model has three parameters, namely, the red-noise amplitude (  ), characteristic timescale (  ), and mean magnitude ().The DRW+sinusoidal model has three additional parameters, namely, the period (), phase (), and amplitude () of the periodic signal (see Equation 1).Based on Equation ( 5), the likelihood function for the DRW+sinusoidal model is written as where   represents the sinusoidal signal.We again employ the MCMC sampling package emcee to obtain posterior samples of the parameters, from which the best-estimated values and uncertainties are assigned by the means and standard deviations, respectively.
In the left panel of Fig. 2, we compare the best estimated periods with those determined by the GLSP in Section 2.2.We can find the general consistency within uncertainties, despite a few discrepant points that might be caused by the inclusion of the DRW process.The period uncertainties of both approaches generally increase with the period because a longer period means fewer cycles of periodic 2016) arising from searches over a broad period range by counting for all possible false positives within the whole period range.variations, resulting in fewer constraints on the period.We calculate the BIC difference between the above two models as A negative ΔBIC indicates a more preference over the DRW + sinusoidal model.As shown in the right panel of Fig. 2, the obtained ΔBIC are overall negative, meaning that the DRW+sinusoidal model is relatively more preferable.The right panel of Fig. 2 also plots the relation between ΔBIC and the FAP estimated in the preceding section.There appears a generic positive correlation as expected if regardless of the scattering.It is noteworthy that the sinusoidal function exhibits a simple and well-defined structure, whereas a sinusoidalshape light curve stemming from red-noise may exhibit slight deviations from the ideal sinusoidal form.Consequently, even if the fake periodicity is caused by red-noises, there is a possibility of favoring the DRW + sinusoidal model in our fitting.Therefore, it is crucial to carefully choose an appropriate BIC value in order to select the preferable models.To be conservative, in addition to the criterion (1-FAP) ≥ 99.73%, we also require ΔBIC ≤ −107 to select out the best periodic candidates.In the end, we censor 57 less significant quasars from the 184 candidates identified in Section 2.2 and obtain a sample of 127 periodic quasar candidates.

The False-alarm Probability with the SPL Model
In the preceding two sections, we assume that the AGN stochastic variability is characterized by the DRW process.In this section, we estimate FAPs of the 184 candidates selected in Section 2.2 using the SPL model.Due to irregular sampling and seasonal gaps, the calculated PSD from the Fourier transform of the light curves may be distorted from the true spectra (Uttley et al. 2002; see a simple test in Appendix A), making it difficult to obtain reliable model parameters by directly fitting the PSDs.To overcome this limitation, we utilize the framework The model parameters are determined using the MCMC technique so that their posterior samples can be thereby obtained.Here, the prior probabilities for the amplitudes9 and slopes of SPL are independently assigned as natural logarithmic and uniformly distributed, covering a range of (-15, 6) and (1, 5) respectively.After determining SPL parameters with RECON, we calculate the FAP using the same procedure as in Section 2.3.1.We generate mock light curves with the method described by Timmer & Koenig (1995).Fig. 3 compares FAPs derived from the DRW and SPL models.The obtained SPL slopes (listed in Table 2) are generally steeper than the DRW slope (=2) at high frequency (short timescale < ).Overall, the SPL model tends to yield relatively higher FAPs than the DRW model.This is not surprising as the SPL model results in higher rednoise powers at low frequencies (long timescale), which can easily produce more spurious periodic signals.Using the FAPs based on the SPL model, the number of AGNs that satisfy (1-FAP) > 99.73% decreases from 158 (DRW) to 120 (SPL).

Estimation of Significance with the GLS Periodogram
We also adopt the procedure proposed by Vaughan (2005Vaughan ( , 2010) ) to quantify the significance of the periodicity, which applies in the frequency domain and involves calculating the PSD of the light curve.However, as mentioned above, due to irregular sampling and seasonal gaps of the light curves, the PSDs obtained through direct Fourier transform are usually seriously distorted (see details in Appendix A).Therefore, we slightly modify the procedure outlined in  Vaughan (2005Vaughan ( , 2010)).Specifically, instead of using PSDs, we employ the GLSP that applies to irregular sampling (see also Chen et al. 2020).By generating a series of mock light curves from DRW and SPL models, we can construct the respective background levels of red-noises in the GLSP.For each candidate, model parameters are randomly drawn from the corresponding posterior samples obtained by using the RECON framework.We then compare the peak power in the GLSP of the observed light curve with the background red-noise power levels to constrain the significance of the periodicity.In Fig. 1, we plot 1 (68.27%), 2 (95.45%), and 3 (99.73%) significance levels of DRW and SPL red-noises for two exemplary candidates.For the DRW model as the null hypothesis, we identify 145 candidates with the peak powers in the GLSPs reaching the 3 significance level.By combining with the selection criterion outlined in Section 2.3.1 and 2.3.2,we obtain 106 high-significance candidates.For the SPL model, we find 86 periodic quasar candidates that have peak powers in the GLSPs reaching the 3 significance level (out of 138 objects) and (1-FAP) > 99.73%.We summarize these numbers of candidates with different selection criteria in Table 3.

Comparison between DRW and SPL models
The RECON framework can calculate the Bayesian evidence and therefore allows us to compare different red-noise models based on the Bayes factor.For this purpose, we also run RECON with the DRW model on the light curves of the selected candidates.We confirm that the obtained DRW parameters are consistent within uncertainties with those obtained through maximizing the likelihood in Equation 5 in Section 2.3.1.We define the Bayes factor  as the ratio of the Bayesian evidence calculated using the SPL and DRW models.A value of  > 1 means the SPL model is preferable to the DRW model.A widely used criterion for quantifying a strong preference in Bayesian model selection provided by Kass & Raftery (1995) is log  > 1.3.Fig. 4 shows the distribution of log  for the 184 periodic candidates with a median value of log  ≈ 0.94 +1.8  −1.1 .This is a marginal value, indicating that the two models exhibit comparable goodnessof-fit to the observed data.Considering that most of the candidates (151/184) have log  > 0 (i.e., the SPL model is relatively more favorable) and the SPL model gives a more conservative significance estimation 10 , we prefer the sample of 86 periodic candidates with high statistical significance identified using the SPL model.Below we by default use this sample as the fiducial periodic sample and proceed with a summary of its properties and further discussions.

Basic Properties of the Periodic Quasar Sample
In appendix B, we detail how to estimate the 5100 Å luminosity, black hole mass, and accretion rate of our sample from the SDSS archival spectra and tabulate the basic properties of our selected candidates in Tabel 1.
Fig. 5 illustrates the distributions of the redshifts, -band magnitudes, and periods from the GLSPs.Fig. 6 illustrates the distributions of the black hole mass, 5100 Å luminosity, and accretion rate (see 10 It is worth noting that the best-fit parameter   of DRW may be underestimated due to the limited time spans of ZTF light curves, leading to an overestimation of significance of the periodicity. Appendix B for the definition).The redshift distribution peaks around  ∼ 0.5−1.8,generally coincident with the redshift distribution of the whole quasar sample.The magnitude distribution has a narrow range of about 17-20 mag and peaks around 19 mag, corresponding to a black hole mass of ∼ 10 9  ⊙ and luminosity of ∼ 10 45.5 erg s −1 .The period ranges between 500 and 950 days and its distribution peaks around 850 days, but has a hard upper limit.This arises from the temporal baseline (∼1500 days) that restricts the period to be less than ∼1000 days in the observed frame.
Fig. 6 also plots the relations between the rest-frame period (from the GLSPs) and black hole mass, 5100 Å luminosity, and accretion rate.There appear superficial anti-correlations in these relations, especially for the black hole mass and luminosity.A simple inspection indicates that these anti-correlations are caused by the redshift effect as follows.Both the observed periods and magnitudes are concentrated on a narrow range.As a result, the rest-frame period  ∝ (1 + ) −1 and the luminosity  5100 ∝  2  increases with (1 + ), where   is the luminosity distance.The black hole mass can be deemed to be proportional to  5100 if regardless of the accretion rate.These factors lead to the spurious anti-correlations between the rest-frame period and black hole mass/luminosity.This can be further verified by a partial correlation analysis (Kendall & Stuart 1979).The correlation coefficient between  and  excluding the dependence on the third parameter  is defined as

Variability Amplitudes
To investigate the relation of intrinsic variability amplitudes versus luminosity, black hole mass, redshift, and rest wavelength, we measure the variability amplitudes by using the expression (Vanden Berk et al. 2004;Liu et al. 2019) where  0 is the amplitude of the best-fit sinusoidal function (see Section 2.2),  2 =  =1  2  /,   is the magnitude uncertainties of the -th point, and  is the number of observations.Fig. 7a-c plot the relation of log  in -band with log  5100 , log  • and log(1 + ).There show anti-correlations in these relations, similar to previous studies on normal AGNs (e.g., Kelly et al. 2009).Regarding the relation between log  and redshift, however, several  previous studies found a positive correlation (e.g., Vanden Berk et al. 2004).It is unknown whether such a difference is intrinsic or caused by different selected quasar samples.
In Fig. 7d, we show the relation between log  and log( • / rest ).This plot is motivated as follows.The close binary supermassive black holes (SMBHs) have been proposed as a possible origination of periodicity in AGNs (e.g., Bon et al. 2012;Farris et al. 2014;Graham et al. 2015b;D'Orazio et al. 2015;Li et al. 2016aLi et al. , 2019;;Liu et al. 2019;Liao et al. 2021).D'Orazio et al. (2015) used the Doppler boosting of the emissions from the accretion disk surrounding the secondary black hole to explain periodic variations.In this scenario, the velocity of the secondary black hole in a circular orbit is where  =  2 / 1 ≤ 1 is the mass ratio,  • =  1 +  2 ,  1 and  2 are the mass of the primary and secondary black hole, respectively, and  is the orbital period.To the first-order approximation, the variability due to the Doppler boosting is where Δ  is the flux at a specific frequency ,   is the spectral in-dex,  is the phase angle of the orbit, and  is the inclination angle.For a large sample, we expect a relation  ∝ Δ log   ∝ log( • / rest ).However, as shown in Fig. 7d, log  has an anti-correlation with log( • / rest ).We notice that there are correlations of the variation amplitude, black hole mass, and period with redshift, which might cause a superficial anti-correlation between log  and log( • / rest ).Therefore, we perform a partial correlation analysis excluding the dependence on redshift as in Section 3.2.We find a partial correlation coefficient of  = −0.35,which is marginal and possibly indicates that the anti-correlation is not intrinsic.
As mentioned above, the ZTF has three custom filters .Their effective wavelengths are  eff () = 4722.74Å,  eff () = 6339.61Å,  eff () = 7886.13Å, respectively.In the above analysis, we by default use the -band data.We can also calculate the variability amplitudes for the data of the other two bands.As such, we can obtain the relation between the variability amplitude with the rest-frame wavelength.We use the light curves with more than 40 binned points for band and more than 20 binned points for -band (because the -band generally has poor sampling).Fig. 7e illustrates the dependence of the variability amplitudes on rest-frame wavelength, which is calculated by dividing the effective wavelength by the corresponding (1 + ).The variability amplitude of our sample has no obvious declining  trend with rest-frame wavelength.This is opposite to the results of Vanden Berk et al. (2004) and Liu et al. (2019), which, however, found the variability amplitude decreases with rest-frame wavelength as  ∝ exp (−/988Å).A possible cause of this discrepancy may be due to the anti-correlation between variability amplitudes and redshift in our sample (see Fig. 7c).
With future development of photometric sky surveys, such as the ongoing ZTF and the upcoming Vera C. Rubin Observatory's Legacy Survey of Space and Time (Ivezić et al. 2019), a larger sample of periodic quasar candidates with much longer temporal baselines and greater diversity in period and redshift is highly worthwhile to testify the above correlations.

The Asymmetry of Broad Line Profiles
To quantity the asymmetry of the broad emission lines, we adopt the dimensionless asymmetry parameter (Brotherton 1996;Du et al. 2018) where   (3/4) and   (1/4) are the wavelengths at the 3/4 and 1/4 of the peak height, and Δ(1/2) is full with at half maximum (FWHM).
The line profile has a blue asymmetry for  > 0 and a red asymmetry for  < 0. We calculate the asymmetry parameters of the prominent H and C iv profiles obtained by the spectral decomposition described in Appendix B. For the prominent Mg ii line, we find that it is sufficient to use one single Gaussian in the spectral decomposition so that the resulting asymmetry parameters  = 0.In Fig. 8, we plot the distributions of  for the H and C iv lines of our sample.Most of the H profiles are roughly symmetric ( ∼ 0) while the C iv profiles are either red or blue asymmetric.
In the scenario of binary SMBHs, the broad emission lines are expected to be asymmetric due to the orbital motion of the black holes that cause velocity shifts to the lines emitted from the broadline regions (BLRs) surrounding each black hole (e.g., Shen & Loeb 2010; Eracleous et al. 2012;Bon et al. 2012Bon et al. , 2016;;Li et al. 2016aLi et al. , 2019)).However, if the orbital separation is smaller than the BLR size, the binary shares a mutual circumbinary BLR and the influences of the orbital motion might be minimized.The generated broad emission line profiles will be similar to those of normal AGNs.We estimate the orbital separations of our sample using Kepler's law where  g is the gravitational radius.We estimate the H BLR sizes using the well-established size-luminosity relation (see Appendix B).
As shown in Fig. 9, the orbital separations are overall smaller than the BLR sizes.
On the other hand, we note that asymmetric emission lines can also be produced from the BLR of a normal AGN with complicated geometry and dynamics, such as inflows and/or outflows (e.g., Denney et al. 2009;Grier et al. 2013;Fausnaugh et al. 2018;De Rosa et al. 2018;Bao et al. 2022;U et al. 2022;Lu et al. 2022;Chen et al. 2023), spiral arms (e.g., Gilbert et al. 1999;Storchi-Bergmann et al. 2003, 2017;Wang et al. 2022;Du & Wang 2023), and partial dust obscuration (e.g., Gaskell & Harrington 2018;Panda & Śniegowska 2022).In this sense, only using one-epoch line asymmetry, we cannot distinguish binary black hole systems from single black holes.However, long-term spectroscopic monitoring of those periodic candidates can yield time variations of the line asymmetry, shedding more light on testifying binary black hole systems.

Comparison with the Previous Works
With systematic searches over survey data, previously there were 111 objects identified from CRTS (Graham et al. 2015a), 50 from PTF (Charisi et al. 2016), 26 from PS1 MDS (Liu et al. 2019), and 4 from the combination of the DES and SDSS (Chen et al. 2020).Fig. 10 compiles the black hole mass of the previous samples together with our sample and plots the relationship between the black hole mass with redshift.The overall redshift distributions of the black hole mass of those samples are similar.
The CRTS sample of Graham et al. (2015a) was based on the combination of the Million Quasars (MQ) v3.7 catalogue and SDSS DR12 quasar catalog.Charisi et al. (2016) used the Half MQ catalog (a subsample of the MQ Catalogue v4.4) and SDSS DR12 quasar catalog as the input quasar catalogs.As a comparison, we use SDSS DR14 quasar catalog and the Véron-Cetty & Véron catalog.Therefore, we expect that there are significant overlaps among these three input catalogs.However, we cross-match our periodic sample with the samples of Graham et al. (2015a) and Charisi et al. (2016), and do not find matching candidates.The possible reasons that there are no matching candidates among the three samples are as follows: 1) the temporal baseline is crucial and different temporal baselines may lead to different quasars identified.For example, most of the candidates from Graham et al. (2015a) have observed periods larger than 1400 days, the candidates from Charisi et al. (2016) have observed periods less than 500 days, and our sample has observed periods in the range from 500 to 950 days.2) The period and variation amplitude might change with time (e.g., Jiang et al. 2022).3) The periodicity might only appear at some epochs and disappear at other epochs (e.g., O'Neill et al. 2022).The latter two behaviours can be explained once the variability is not strictly periodic, but quasi-periodic in the sense that the corresponding PSD shows a broad peak on top of a red-noise-like shape (see Section 4.2.2 for more discussions).

The Look-elsewhere Effect from Searching over an AGN Sample
In Section 2.3.1, we have estimated the FAP for an individual AGN, which includes the look-elsewhere effect arising from searches over a broad period range.Another type of look-elsewhere effect is the false positive when globally searching over an AGN sample due to random fluctuations.To estimate this kind of "global FAP", we perform the following simulations.We generate 143,700 mock light curves with the same method in Appendix A by assuming that the AGN variability is purely stochastic and modeled by SPL.The SPL slopes are randomly assigned from a range of 1.5-3.5, based on the slope distribution of 184 candidates (see Figure 3).We then search over these mock light curves to identify "fake" periodic candidates using exactly the same procedures employed for the identification of the 86 candidates.We select out a total of 51 candidates, implying a global FAP ∼ 3.6 × 10 −4 .This result indicates that in a statistical sense, most of our 86 periodic candidates might be normal AGNs and their spurious periodic variations are caused by red-noise.
In light of the temporal limitations of ZTF data, additional data will be helpful to reduce the global FAP and verify/falsify the periodicity.In next section, we extend the temporal baselines of ZTF light curves by compiling other time-domain survey data.

Photometric Data Intercalibration
By incorporating archival survey data from CRTS, PTF, PS1, and ATLAS (short for Asteroid Terrestrial-impact Last Alert System), we extend the light curves of the 86 periodic candidates with high statistical significance to further verify the periodicity.Here, we utilize -band data of PTF and PS1 and -band data of ATLAS.Given that these photometric data were obtained with different instruments, filters, and reduction methods, it is imperative to perform intercalibration between them.For the data obtained from ATLAS, the monitoring period temporally overlapped with that of ZTF.We perform .The distribution of the redshift and black hole mass for the periodic candidates found in this work and previous works (Graham et al. 2015a;Charisi et al. 2016;Liu et al. 2019;Chen et al. 2020).intercalibration using the Bayesian package PyCALI11 (Li et al. 2014) by adopting the ZTF light curve as a reference.As for the data from CRTS, PTF and PS1, we first adopt the -band light curve from PS1 as the reference and align other light curves with this reference.We then convolve the SDSS spectrum with the PS1/ZTF filter transmissions to obtain synthetic magnitudes and calculate the magnitude difference between the two photometric systems.Finally, we add this magnitude difference to the aligned data of CRTS, PTF, and PS1 to keep consistent with the ZTF data.In Appendix C, we show the intercalibrated long-term light curves of the 86 periodic candidates.

Verifying the Periodicity
With the extended light curves, we can verify the periodicity using the same methods described in Section 2.212 .Unfortunately, we find that no candidate satisfies the criterion of  > 4.0.This indicates that the light curves of the periodic candidates may not rigorously adhere to the expected sinusoidal form due to contamination of red-noises (see below), the influence of poor photometric uncertainties, or other unknown factors.With this consideration, we reidentify the periodic quasars by discarding the criterion of  > 4.0 but still requiring the peak power in the GLSP exceeding 99.73% (in terms of using the SPL model as the null hypothesis), 1 −  ACF,Extend / GLS,Extend < 0.1, and  < 10 −3 day −1 .Meanwhile, we imposed an additional criterion that the difference between the GLSP periods from extended and ZTF light curves should be within 10 percent.
Finally, three candidates (i.e., SDSS J083111, SDSS J155304 and tion 2.2 because of two main reasons.Firstly, The ATLAS database provides on-the-fly reduced photometry and acquiring 143,700 objects is unrealistic presently.Secondly, accurately intercalibrating fluxes for all archived data poses a challenge for some objects due to the lack of SDSS spectra and/or there being no time overlap between CRTS, PTF, PS1, ATLAS, and ZTF. SDSS J135830; see Figure 11) with 6-8 cycles of periods are identified using the extended data.In particular, the candidates SDSS J083111 and SDSS J155304 reach a significance level of 99.99%.Notably, it seems that the periodic variation in the light curve of SDSS J083111 started to appear around JD 2,454,000 (see the upper panel of Figure 11), after a rapid declining trend that apparently deviates from the expected periodic variability.Such behavior is reminiscent of the candidate PKS 2123-021 reported by O'Neill et al. ( 2022), which exhibits periodic variability only at certain epochs (see Fig. 1 therein).
For the remaining 83 periodic candidates that do not satisfy the criteria, there are two possibilities.First, the periodic variations might be spurious and purely originate from red-noises as discussed in Section 4.1.However, given the sparse sampling and/or large photometric uncertainties of the extended data, more monitoring data are required to finally confirm this point.Alternatively, as mentioned above, the observed behaviours reflect that the variability is quasi-periodic instead of strictly periodic.Apparently, the significance level of the periodicity under the present criteria will be diminished if the period varies over time (e.g., Jiang et al. 2022) or is only present during certain epochs (e.g., O'Neill et al. 2022).This makes it challenging to identify true periodic quasars and estimate the significance of the periodicity.
For the sake of illustration, in Fig. 12 we artificially generate mock light curves using PSDs composed of a SPL with a slope of  = 2 and a Lorentzian function.As can be seen, in the cases of mild Lorentzian widths, the generated light curves show sinusoidal patterns at some epochs but seemingly stochastic patterns elsewhere.Also, it appears that the amplitudes of the sinusoidal patterns vary with time.As expected, the generated periodic signal becomes more visible when the Lorentzian width decreases.Besides the Lorentzian width, the underlying SPL also introduces extra stochasticity to the (quasi-)periodic signals.

Possible Explanations of Periodic Variability
In this section, we present a brief discussion on the explanations of the periodicity based on the three candidates identified with extended light curves by assuming that the periodicity is real (see also Graham et al. 2015a andDotti et al. 2023).

Doppler Boosting of Supermassive Black Hole Binaries
This scenario was proposed by D' Orazio et al. (2015), in which the periodicity arises from the Doppler-boosted emissions of the accretion disk surrounding the secondary black hole.In Section 3.3, we have tested this scenario by checking on the correlation between variability amplitudes and log( • / rest ) for the 86 candidates.We do not find a significant, expected positive correlation by excluding the selection bias.
Another testable point regarding the Doppler boosting scenario is that the resulting relative flux and thus magnitude variations are expected to be wavelength-independent and simultaneous (i.e., no time delay) among different bands.To test this point, we measure time delays among the -, -and -bands ZTF light curves using Table 4. Time delays with the ICCF method and phase differences of the best sinusoidal fits for the -, -and -bands light curves of three candidates.The -band light curves of SDSS J083111 and SDSS J135830 have too few data points to obtain reliable time delays and phase differences.

ID
Δ  Δ  (day) (day) (day) (day) SDSS J083111 −8.1 ± 10.7 -−13.7 ± 4.7 -SDSS J155304 13.5 ± 5.9 35.3 ± 10.9 11.2 ± 4.0 39.0 ± 9.6 SDSS J135830 −8.2 ± 12.5 -0.6 ± 8.5 two methods.The first method relies on calculating the ICCF, from which the time delay is determined as the ICCF centroid above 80% of the peak value.We denote   and   as time delays of the or -bands light curves relative to the -band one, respectively.To reduce the influence of short-timescale variability (ascribed to the accretion disk emissions), we smooth the light curves using a median filter of 7 points prior to ICCF analysis.An alternative method is calculating phase differences among the best sinusoidal fits (with the same fixed period) for the three-band light curves.We define Δ  =(  −   )/2 and Δ  =(  −   )/2, where   ,   and   are the best-fit phases, and  is the period.We list time delays and phase differences (in the observed frame) in Table 4.The obtained time delays and phase differences are consistent with each other within uncertainties.For the candidate SDSS J155304, there are significant time delays between -and /-bands, suggesting that the periodicity might not be caused by the Doppler boosting.For the other two candidates, the time delays and phase differences have large uncertainties and therefore no meaningful constraints can be made.

Periodic Accretion of Supermassive Black Hole Binaries
The orbital motion of supermassive black hole binaries can produce periodic modulation to mass accretion onto each black hole, giving rise to periodic emissions of the accretion disks.Previous hydrodynamic simulations of circumbinary accretion discs predicted that the periodic pattern would exhibit a sawtooth shape rather than a sinusoidal form (e.g., Farris et al. 2014).The quality of the present light curves does not allow us to reliably distinguish the shapes of periodic variations.Nonetheless, in the case of SDSS J155304, the detection of significant inter-band time delays can be explained under this scenario.

Other Explanations
The observed optical flux variations in quasars may be a result of the combination of thermal emissions from the accretion disk and nonthermal emissions from the Doppler-boosted radio jet.The periodic variations arise from the jet precession, which might be driven by the orbital motion of a binary SMBH system or an internally rotating jet flow of a single black hole system.However, after cross-matching with the sample of radio-loud AGNs from Best et al. (2005), we find that none of the three candidates were included in this sample.This rules out the precessing jet scenario.Alternatively, in addition to radio jets, periodic variations could also be attributed to the precession of a warped accretion disk surrounding a single black hole, which results in the obscuration of the emission region (Martin et al. 2007).However, as discussed by Graham et al. (2015b) and Chen et al. ( 2020), for a 10 8  ⊙ SMBH, the precessing timescale is on the order of several tens of years.This is much longer than the periods (∼2.5 years) of the three periodic quasar candidates, whose SMBH masses are in a range of 10 8.3−9.0  ⊙ .

Caveats in the Methodology of Periodicity Identification
There are two caveats in the method we adopted to identify candidates in Section 2.2.
Firstly, we initially select the candidates by using the best-fit sinusoidal fit with the period fixed to  GLS , which corresponds to the largest periodogram peak.In this step, we assume that the periodic signal is solely affected by white noises, namely, the noise power is constant across frequency.However, it is important to note that red-noise power decreases with frequency, as shown in the bottom left panels of Figure 1.Consequently, our approach may result in the exclusion of some periodic candidates with short periods.A possible method to surmount this issue is to pre-select the period corresponding to the periodogram peaks that exhibit the highest significance in light of the red-noise background (the peak power in GLSP reaches 3 significance level at least), and then perform the same subsequent procedures as in Section 2.2 such as sinusoidal function fitting.The red-noise background can be constructed by assuming the PSD follows some specific shapes (e.g., either DRW or SPL; see details in Section 2.3.4).In Section 3.1, the comparison between DRW and SPL models reveals that the latter is more favorable for our quasar sample.However, due to the irregular sampling and seasonal gaps, the SPL parameters cannot be directly recovered by fitting the PSD (see Appendix A).Although the RECON framework offers an alternative method to obtain reliable PSD parameters, it is not feasible to construct the red-noise background for each of 143,700 quasars because of the prohibitive computational time required for running RECON.
Here, we carry out a simple test to obtain a quantitative estimation of the number of candidates with shorter periods that may be missed when using the method in Section 2.2.We adopt the DRW process instead of the SPL model and construct the red-noise background for each object in our parent sample.The subsequent steps and selection criteria fully follow those in Section 2.2.We finally identify 125 candidates covering a period range of 400-900 days, in which 90 objects overlap with our sample of 184 candidates identified in Section 2.2.These results imply that there are no serious missing identifications of short-period candidates in Section 2.2.
Secondly, the adopted criteria in the present work are inclined to select out candidates with sinusoidal variations.Specifically, the limit on S/N ratio  > 4 is sensitive to the amplitude of the sinusoidal variations.In addition, for the ACF selection script, the exponentially decaying cosine function in Equation (2) indeed corresponds to a noise-modulated sinusoidal signal with a form of (Jung 1993) where  is the phase and () is the Wiener noise with a zero mean and correlation function of ⟨( 1 )( 2 )⟩ = (2)min( 1 ,  2 ).Here,  −1 represents the coherence time of the phase fluctuations in the sinusoidal signal.We require  < 10 −3 day −1 , meaning that the coherence time is as long as 1000 days, comparable to the time baselines of ZTF data.As a result, this selection criterion also tends to identify sinusoidal variations, although not as sensitive as the S/N limit.

Future Improvements for the Searching Method of Periodicity
Considering that realistic variations might deviate from a sinusoidal form (e.g., Farris et al. 2014) and also be quasi-periodic with a broad distribution of the period (see Fig. 12), it is necessary to devise new effective identification methods in future works.A possible starting point might be resorting to PSD fitting.However, since AGN light curves usually have highly irregular sampling and significant inhomogeneity in observing instruments and data reduction, the traditional direct PSD fitting as in X-ray analysis (e.g., Uttley et al. 2002) is no longer applicable.Li & Wang (2018) proposed a forward approach to recover any given shapes of PSDs from irregularly sampled light curves with measurement noises.This approach can be easily extended to the case of a quasi-periodic PSD on top of a red-noise PSD.The Bayes factor can be employed to distinguish a quasi-periodic PSD from a pure red-noise PSD.
We note that since a quasi-periodic PSD has distributed frequencies (and thereby periods), reliable identification of such signals would be challenging and require longer temporal baselines and better sampling rates.We expect that as the ZTF data continuously accumulates and other time-domain surveys come into operation in the near future (such as the Rubin-Large Synoptic Survey Telescope; Ivezić et al. 2019), there will be a sizable sample of statistically significant quasi-periodic quasar candidates finally identified.

SUMMARY
Based on the quasar catalogs of SDSS DR14 and Véron-Cetty & Véron (2010), we systematically search for quasars with periodic variations over the ZTF archival light curve data.As a first step, we employ the GLSP and ACF methods to make a primitive selection.We restrict that there are at least 1.5 cycles of periods and sufficiently strong periodic signals in the light curves (see Section 2.2).Also, we require the difference between the obtained periods from the GLSPs and ACF method to be smaller than 10%.An initial sample of 184 periodic quasars candidates is obtained.As a second step, we evaluate the significance of the periodicity by calculating the FAP and comparing the peak power with the background red-noise power levels in the GLSP.We adopt the null hypothesis that quasar stochastic variabilities arise from red-noises modeled by either DRW or SPL.We set a limitation of (1 − FAP) ≥ 99.73% and require the highest GLSP peak exceeding the 3 level.Additionally, for the case of DRW red-noises, we can directly fit the light curves and calculate the BIC to compare between the pure stochastic DRW model and the periodic (sinusoidal + DRW) model.To select out the sinusoidal + DRW model, we adopt a criterion of ΔBIC ≤ −10.Finally, we identify 106 (DRW) and 86 (SPL) significant periodic candidates out of a total of 143,700 quasars, respectively.The Bayes factors derived from RECON demonstrate that for our initial 184 candidates, the SPL model is relatively preferable over the DRW model (see Section 3.1).Therefore, we choose the sample of 86 periodic candidates selected using the SPL model as our final candidate sample.
We perform Monte Carlo simulations to estimate the global lookelsewhere effect, which accounts for false positives resulting from a global search over an AGN sample due to random fluctuations (see Section 4.1).The obtained "global FAP" is ∼ 3.6 × 10 −4 , meaning that there will be about 51 false positives in a sample of 143,700 quasars.A such FAP is comparable to the probability of discovering the above 86 periodic candidates, indicating that, statistically, most of the periodic candidates in our sample might be normal AGNs with their periodic variability caused by red-noise.
Notably, compared to the previous searches made by Graham et al. (2015a); Charisi et al. (2016), their parent quasar catalogs are significantly overlapped with our catalog, however, we find that the identified periodic quasar candidates are not matching.In addition, we extend the ZTF light curves of the 86 candidates by supplying other earlier survey data (including CRTS, PTF, PS1, and ATLAS) to verify the periodicity.With extended light curves, only three quasars with 6-8 cycles of periods satisfy the selection criteria (see Section 4.2.2).These results, together with the high "global FAP" mentioned above, impose an implication that the variations of the candidates are most likely quasi-periodic rather than strictly periodic or that they are purely caused by stochastic red-noise.In this case, the selection criteria and methods used in both this work and previous works (e.g., Graham et al. 2015a;Charisi et al. 2016;Liu et al. 2019;Chen et al. 2020) might be no longer effective.Indeed, those criteria are inclined to identify candidates with sinusoidal variations.In summary, more time-domain data of all identified candidates with longer temporal baselines are necessary to confirm their periodicity and rule out false positives due to red-noise.Also, new effective methods are highly warranted to identify generic quasi-periodic quasar candidates.

Figure 1 .
Figure 1.Two examples of our selected periodic quasar candidates.For each candidate, the upper panel shows the light curves from ZTF DR11.The orange points show the original -band light curve and the black points show the binned light curve with an interval of 20 days.The blue dash line shows the best sinusoidal fit.The two bottom left panels show the GLSP and the significance levels in terms of DRW and SPL red-noise models (see Section 2.3.4).The bottom right panel shows the ACF of the light curve (the black solid line) and the best fit of an exponentially decaying cosine function (the blue dashed line; see Section 2.2).

Figure 2 .
Figure 2. Left panel: a comparison between the periods from the DRW+sinusoidal model fitting and the periods from the GLSP.Right panel: a comparison between the FAP and ΔBIC.The horizontal and vertical dotted lines correspond to (1-FAP) = 99.73% and ΔBIC = −10, respectively.In each panel, the blue points highlight those candidates with (1 − FAP) ≥ 99.73% and ΔBIC ≤ −10.

Figure 3 .
Figure 3.A comparison of FAPs between the DRW and SPL models.The points are color-coded by the best-fit slopes of the SPL from RECON.Both the horizontal and vertical dashed lines correspond to (1-FAP) = 99.73%.

Figure 4 .
Figure 4.The distribution of Bayes factors log , defined as the evidence ratio of the SPL to DRW models (i.e.,  =  SPL / DRW , where  SPL and  DRW are the Bayesian evidence).The vertical red solid line shows log  = 0.

Figure 5 .
Figure5.The distributions of the redshift, magnitude, and observed period (from the GLSP) of the periodic quasar sample.Note that the vertical axis of the rightmost panel is in a logarithm scale.

Figure 6 .
Figure6.The rest-frame period (from the GLSPs) versus black hole mass, 5100 Å luminosity, and accretion rate of our periodic quasar sample.Note that the superficial anti-correlations are caused by the redshift effect.The denotation  shows the partial correlation coefficient excluding the dependence on redshift (see the text).

Figure 7 .
Figure7.Variability amplitudes versus 5100 Å luminosity, black hole mass, redshift, rest-frame wavelength and the ratio of black hole mass to period.The denotation  marks Spearman's rank correlation coefficient.In panels (a)-(d), the blue points are the variability amplitude in -band.In panel (e), the blue, orange and green points are the variability amplitudes in -, -and -band, respectively.The dashed line is from VandenBerk et al. (2004).

Figure 8 .
Figure 8.The distributions of  for the H and C iv lines of our periodic quasar sample.

Figure 9 .
Figure 9.Comparison of the orbital separations versus the H BLR sizes by assuming that the periodicity arises from the orbital motion of supermassive black hole binaries and the BLR size follows the size-luminosity relation.
Figure10.The distribution of the redshift and black hole mass for the periodic candidates found in this work and previous works(Graham et al. 2015a;Charisi et al. 2016;Liu et al. 2019;Chen et al. 2020).

Figure 11 .
Figure11.Three periodic quasar candidates selected from extended light curves.For each candidate, the upper panels show the synthetic light curves from the CRTS (blue), PS1 (brown), PTF (red), ATLAS (light blue), and ZTF (orange) archival survey data.The blue dash line shows the best sinusoidal fit.The bottom left panel shows the GLSP and the significance levels in terms of the SPL red-noise model.The bottom right panel shows the ACF of the extended light curve (the black solid line) and the best-fit using an exponentially decaying cosine function (the blue dashed line).

Figure A1 .
Figure A1.Recovered SPL slopes by directly fitting PSDs from evenly/unevenly sampled light curves versus the input slopes.
4can produce spurious fewcycle periodic signals in AGN light curves, especially for those with

Table 1 .
Properties of the initial 184 periodic quasar candidates identified in this work.The full version is available in a machine-readable form online.

Table 2 .
Statistical significance estimation of the periodicity.The full version is available in a machine-readable form online.   /day) log(   /mag) log FAP ΔBIC >3 in GLSP   log FAP >3 in GLSP  log  HD  HS Whether the highest peaks of GLSPs reach 3 significance levels.Y-yes and N-No.Whether the candidates satisfy our whole selection criteria for high significance when using the DRW (HD) and SPL (HS) models as the null hypothesis.Y-yes and N-No.RECON 8 (Li & Wang 2018; see Appendix A for a brief description of RECON and validity tests.) to determine the parameters of the SPL model.RECON uses a stochastic complex series to parameterize AGN variability in the frequency domain and transform the series back to the time domain to fit the observed light curve.RECON can cope with irregularly sampled light curves and also can recover any forms of PSDs.

Table 3 .
The selection criterion used to identify periodic quasars candidates with high significance and the respective counts derived using the DRW and SPL models as the null hypothesis.The ΔBIC is calculated only for the DRW model.