---
abstract: 'We study the long-term variability in the optical monitoring database of Ark 120, a nearby radio-quiet active galactic nucleus (AGN) at a distance of 143 Mpc ($z=0.03271$). We compiled the historical archival photometric and spectroscopic data since 1974 and conducted a new two-year monitoring campaign in 2015-2017, resulting in a total temporal baseline over four decades. The long-term variations in the optical continuum exhibit a wave-like pattern and the H$\beta$ integrated flux series varies with a similar behaviour. The broad H$\beta$ profiles have asymmetric double peaks, which change strongly with time and tend to merge into a single peak during some epochs. The period in the optical continuum determined from various period-search methods is about $20$ years and the estimated false alarm probability with null hypothesis simulations is about $1\times10^{-3}$. [The overall variations of the broad H$\beta$ profiles also follow the same period. However, ]{} the present database only covers two cycles of the [suggested]{} period[, which strongly encourages]{} continued monitoring to track more cycles and confirm the periodicity. Nevertheless, in light of the possible periodicity and the complicated H$\beta$ profile, Ark 120 is one candidate of the nearest radio-quiet AGNs with possible periodic variability and [it is thereby]{} a potential candidate [host for]{} a sub-parsec supermassive black hole binary.'
author:
- 'Yan-Rong Li'
- 'Jian-Min Wang'
- 'Zhi-Xiang Zhang'
- Kai Wang
- 'Ying-Ke Huang'
- 'Kai-Xing Lu'
- Chen Hu
- Pu Du
- Edi Bon
- 'Luis C. Ho'
- 'Jin-Ming Bai'
- 'Wei-Hao Bian'
- 'Ye-Fei Yuan'
- Hartmut Winkler
- 'Eduard K. Denissyuk'
- 'Rashit R. Valiullin'
- Nataša Bon
- 'Luka [Č.]{} Popovi[ć]{}'
title: |
    A Possible $\sim$20-Year Periodicity in Long-term Optical [Photometric and Spectral]{} Variations of the Nearby Radio-Quiet\
    Active Galactic Nucleus Ark 120
---

Introduction
============

[lccccccccc]{} D99b & 1 & 42,392-47,944 & 1974-1990 & 63 & $2.0\times7.0$ & 1.35 & $1.008\pm0.017$ & $\phm{+}2.05\pm0.25$\
P89 & 2 & 44,168-47,414 & 1979-1988 & 74 & $7.0$ & 0.925 & $1.005\pm0.021$ & $\phm{+}1.39\pm0.32$\
P98 & 3 & 47,524-50,388 & 1988-1996 & 141 & $5.0\times7.6$ & 0.91 & 1 & 0\
W96 & 4 & 48,619-48,834 & 1991-1992 & 20 & $5.0\times10.0$& 0.70 & $1.078\pm0.005$ & $\phm{+}0.46\pm0.06$\
D08 & 5 & 48,629-53,445 & 1992-2005 & 88 & $3.0\times11.0$& 1.13 & $1.037\pm0.003$ & $\phm{+}0.09\pm0.03$\
L & 6 & 57,295-57,839 & 2015-2017 & 100 & $2.5\times8.5$ & 0.96 & &\
M79 & 7 & 25,622-33,446 & 1935-1950 & 52 & & & &\
P83 & 8 & 42,392-44,226 & 1974-1979 & 36 & $15.0$ & & $0.687\pm0.087$ & $\phm{+}0.45\pm0.95$\
P89 & 2 & 43,426-47,414 & 1977-1988 & 74 & $7.0$ & & $0.607\pm0.297$ & $\phm{+}5.04\pm0.34$\
D99a1 & 9 & 43,963-50,818 & 1979-1998 & 25 & $14.3$ & & $1.132\pm0.070$ & $-5.04\pm0.90$\
D99a2 & 9 & 43,079-51,186 & 1976-1999 & 50 & $27.5$ & & $1.126\pm0.064$ & $-8.07\pm1.03$\
W92/97 & 10 & 47,534-52,282 & 1989-2002 & 37 & $20.0$ & & $0.719\pm0.017$ & $-0.79\pm0.28$\
R12 & 11 & 52,207-55,644 & 2001-2011 & 54 & & & $0.787\pm0.058$ & $-0.13\pm0.85$\
K14 & 12 & 52,237-54,319 & 2001-2007 & 67 & $8.3$ & & $0.903\pm0.030$ & $-0.34\pm0.40$\
C & 13 & 53,709-56,580 & 2005-2013 & 73 & & & $0.358\pm0.015$ & $\phm{+}6.51\pm0.22$\
H & 14 & 55,122-55,268 & 2009-2010 & 96 & 7.5 & & 1 & $-7.86\pm0.05$\
A & 15 & 55,949-58,345 & 2012-2018 & 476 & & & $0.808\pm0.005$ & $-2.97\pm0.07$ \[tab\_lc\]

Periodic brightness variability in long-term monitoring of active galactic nuclei (AGNs) is widely used to search for supermassive black hole (SMBH) binary candidates in modern time-domain surveys (e.g., @Graham2015a [@Graham2015b; @Liu2015; @Liu2016; @Charisi2016; @Zheng2016]). Although there [are]{} alternative explanations, such a periodicity is generally believed to [reveal]{} the binary’s orbital motion, which modulates accretion processes into the black holes, giving rise to periodic variations in disk emissions. To date systematic searches from large surveys have yielded more than one hundred quasar and AGN candidates with periodic variability, distributed over a wide redshift range up to $z\sim3$ (@Graham2015b [@Charisi2016; @Liu2016]). Due to the limited temporal baselines ($\lesssim$10 years), these detected rest-frame periods are confined to short timescales of several years or hundreds of days. There are also a number of quasars and nearby AGNs that were reported individually to exhibit periodicity based on long-term databases, including (but not limited to) OJ 287 ($z=0.31$; @Valtonen2008), NGC 4151 ($z=0.003$; @Guo2006 [@Bon2012]), and NGC 5548 ($z=0.017$; @Li2016 [@Bon2016]). OJ 287 is long-known to have strong radio emission while NGC 4151 and NGC 5548 [have relatively weak radio emission]{} (@Ho2002).

![image](ark120_lc_org.pdf){width="49.50000%"} ![image](ark120_lc.pdf){width="49.50000%"}

In this paper, we study the variability and analyze the possible periodicity in Ark 120 ($z=0.03271$[^1]), a nearby radio-quiet AGN with a radio-loudness $R\approx0.1$ (@Condon1998 [@Ho2002]). Ark 120 has been intensively observed and investigated in optical/UV (e.g., @Kollatschny1981a [@Kollatschny1981b; @Schulz1981; @Alloin1988; @Marziani1992; @Peterson1998; @Stanic2000; @Popovic2001; @Doroshenko2008; @Kuehn2008]) and X-ray bands (e.g., @Vaughan2004 [@Nardini2016; @Reeves2016; @Gliozzi2017; @Lobban2018]). The [published]{} photometric and spectroscopic monitoring data of Ark 120 date back to 1974. Even earlier (1935-1950) there were historical archival plates collected in the Harvard College Observatory that recorded the photometric images of Ark 120 (@Miller1979). [A recent campaign spectroscopically monitored Ark 120 for about 100 days in 2017 with the 2.3m telescope of Wyoming Infrared Observatory (@Du2018). We carried out further observations between 2015 and 2017 using the 2.4m telescope of Yunnan Observatories.]{} We compile publicly available photometric and spectroscopic data of Ark 120 and obtain a total temporal baseline stretching over four decades. By visual inspection, the long-term variations in [the]{} optical continuum display a wave-like pattern with a possible period of $\sim$20 years. Moreover, the broad H$\beta$ profile of Ark 120 is highly asymmetric and varies strongly with time (@Doroshenko2008 [@Du2018]), a [remarkable]{} analogy to [what was discovered in]{} NGC 5548 (@Li2016 [@Bon2016]).

The paper is organized as follows. Section 2 introduces briefly our new monitoring campaign and then presents our construction of a long-term database for Ark 120. Section 3 studies the long-term optical variability and measures the power spectral density (PSD). In Section 4, we perform detailed periodicity analysis for the optical light curve and estimate the associated false alarm probability. In Section 5, we summarize the long-term variations of the broad H$\beta$ profile. A brief discussion on the possible periodicity is given in Section 6 and conclusion is given in Sections 7.

Observations and Historical Database
====================================

New Observations between 2015-2017
----------------------------------

We carried out a two-year observation campaign between 2015 and 2017 using [the]{} Lijiang 2.4m telescope, located in Yunnan Province, China, operated by Yunnan Observatories. [We used the Yunnan Faint Object Spectrograph and Camera (YFOSC) with Grism 14, which covers a wavelength range of 3800-7200 [Å]{} and provides a resolution of 1.8 [Å]{} pixel$^{-1}$.]{} This spectrograph is equipped with a slit long enough to allow us to simultaneously observe a nearby comparison star (@Du2014). We use this comparison star as a reference standard to achieve accurate absolute and relative flux calibrations, whereas the absolute flux of the comparison star is calibrated by observations of spectrophotometric standards during nights of good weather conditions. This flux calibration method is better than the common way based on the flux constancy of the narrow \[\] line as described below. However, we confirm that the narrow \[\] $\lambda5007$ lines of the calibrated spectra are constant with a mean of $0.96\times10^{-13}{\rm~erg~s^{-1}~cm^{-2}}$ and a scatter within 2.5%. The slit was fixed at a projected width of $2.5^{\prime\prime}$ and the spectra were extracted using a uniform window of $8.5^{\prime\prime}$. [The raw spectroscopic data are reduced using standard IRAF version 2.16 routines before absolute flux calibration, which includes bias subtraction, flat-field correction, and wavelength calibration.]{} The details for the data reduction and analysis will be presented in a separate paper [that describes the reverberation mapping analysis (see also @Du2014)]{}. We obtained [100 spectroscopic observations in total]{}, with a typical exposure [time]{} of 30 minutes. [The 5100 [Å]{} continuum flux is measured in a band of 20 [Å]{} wide around 5100 [Å]{} in the rest frame. To be consistent with the other datasets, the H$\beta$ flux is measured as follows: we first subtract the continuum underneath the H$\beta$ line by interpolating between two continuum bands 4760-4790 [Å]{} and 5000-5020 [Å]{}. We then integrate the continuum-subtracted H$\beta$ flux between 4810 [Å]{} and 4910 [Å]{}.]{}

The [previous published]{} spectroscopic observations were only to 2005, so there is no temporal overlap with our observations. As a result, to align the fluxes with the historical data, we need to correct for the differences in the host galaxy contamination. Based on the two-dimensional decomposition of [the]{} [*HST*]{} host galaxy image of Ark 120 by [@Bentz2009], we estimate the relative difference of the host galaxy contribution to 5100 [Å]{} fluxes though the apertures used in our observations and in [@Peterson1998] [to be]{} $\Delta F_{\rm gal}=-(0.9\pm0.3)\times10^{-15}~\rm erg~s^{-1}~cm^{-2}~{\text{\AA}^{-1}}$.

![image](Miller_lc.pdf){width="95.00000%"}

Light Curves of 5100 [Å]{} Continuum and H$\beta$ Fluxes {#sec_data}
--------------------------------------------------------

Table \[tab\_lc\] summarizes the photometric and spectroscopic datasets [published]{} to date for Ark 120 in the literature. These datasets [use]{} various apertures and different flux standards of \[\] $\lambda5007$ for absolute spectral calibrations. Intercalibration is needed to correct for these differences.

For spectroscopic data, we select the dataset of [@Peterson1998] as a reference and scale the other datasets [to]{} a common absolute flux of \[\] $\lambda5007$, $F([\text{\ion{O}{3}}])=0.91\times10^{-13}{\rm~erg~s^{-1}~cm^{-2}}$ as the reference dataset used[^2]. To account for the different aperture sizes, we apply a scale factor $\varphi$ and a flux modulation $G$ to the 5100 [Å]{} flux densities and H$\beta$ fluxes of each dataset with respect to the reference [using]{} (e.g., @Peterson1995 [@Li2014]) $$F_{\lambda}(5100\text{\AA}) = \varphi F_{\lambda}(5100\text{\AA})_{\rm obs} + G,$$ and $$F({\rm H\beta}) = \varphi F({\rm H\beta})_{\rm obs}.$$ The values of $\varphi$ and $G$ are determined by comparing the closely spaced measurements of the two datasets. The interval for comparison is adopted to be 50 days for the datasets D08 and W96 in Table \[tab\_lc\], which yields a good temporal overlap. For the [datasets D99b and P89]{} there is a short interval (2 years) of temporal overlap with the other datasets. Meanwhile, the sampling of these two datasets is poor. We thus increase the interval for comparison up to 100 days to obtain a sufficient number of nearly [contemporary]{} observations. Since we concentrate on long-term variations, the choice of the interval has no influence on our analysis. It is impossible to intercalibrate the H$\beta$ fluxes of our new observations with the other datasets, because there is no temporal overlap. We present the light curve of the H$\beta$ fluxes for the purpose of illustrating the variation behaviour of [the]{} H$\beta$ line. We do not use it for subsequent periodicity analysis.

For $V-$band photometric data, we first convert the magnitudes into flux densities by adopting the zero point $F_\lambda(V=0)=3.92\times10^{-9}\,{\rm erg~s^{-1}~cm^{-2}}\textrm{\AA}^{-1}$ (@Johnson1966). We then scale the flux densities to align with the 5100 [Å]{} flux densities by again applying a scale factor ($\varphi$) and flux adjustment ($G$) to bring them [to]{} a common flux scale with the light curve of [the]{} 5100 [Å]{} continuum [using]{} $$F_{\lambda}(5100~\text{\AA}) = \varphi F_{\lambda}(V)_{\rm obs} + G.$$ The interval for comparison is adopted to be 50 days for all the datasets. Table \[tab\_lc\] also lists the values of $\varphi$ and $G$ for all the spectroscopic and photometric datasets. In the right panels of Figure \[fig\_lc\], we plot the intercalibrated light curves of $V-$band and 5100 [Å]{} flux densities, and H$\beta$ integrated fluxes. The original light curves without intercalibration are plotted in the left panels of Figure \[fig\_lc\] for comparison. Appendix A tabulates a portion of the merged light curve of [the]{} 5100 [Å]{} and $V$-band fluxes, while the entire light curve data is available in the online journal.

The Dataset of [@Miller1979] {#sec_miller}
----------------------------

[@Miller1979] compiled the $B$-band data of Ark 120 over 1935-1950 using the archival plate collection of the Harvard College Observatory. Albeit with poor cadence and low measurement accuracy, this dataset still provides precious historical variation information of Ark 120. Indeed, [@Miller1979] detected significant flares around 1937-1939. We digitalize Figure 1 of [@Miller1979] and obtain the $B$-band magnitude data. It is impossible to intercalibrate this dataset with the other datasets described above. Despite this, a useful comparison is possible by converting the $B$-band magnitudes into fluxes and then simply adjusting the fluxes to make the mean and standard deviation to be at the same scales as those of the other datasets. [*We do not include this dataset for subsequent analysis*]{}, but take it as an auxiliary dataset to illustrate the long-term variations of Ark 120. Figure \[fig\_miller\] plots the light curve data of [@Miller1979] together with the merged light curve of [the]{} $V-$band and 5100 [Å]{} continuum.

![(Top) Power spectral density $P(f)$ and (bottom) $f\times P(f)$ of Ark 120 in the optical and X-ray band (0.3-10 keV). The normalization of the X-ray PSD is adjusted to match the optical PSD. Blue points represent the binned PSD and [the]{} grey line represents the raw PSD. [Slate blue]{} points represent the X-ray PSD from [@Lobban2018]. Vertical dashed lines correspond to the period of 20.5 years.[]{data-label="fig_psd"}](fig_psdall.pdf){width="48.00000%"}

Long-term Optical Variations {#sec_var}
----------------------------

As shown in Figure \[fig\_lc\], all the light curves of [the]{} $V-$band, 5100 [Å]{} continuum and H$\beta$ integrated fluxes oscillate with a wave-like pattern. In the merged light curve plotted in Figure \[fig\_miller\], there [are clearly]{} two major crests around 1985 and 2005 and two major troughs around 1995 and 2015. Meanwhile, [one prominent flare appears]{} around the bottom of each major trough. Both flares rise rapidly in flux by more than 50% within 200 days. The latest flare (2016-2018) then undergoes an even more sudden drop within the same timescale. All the separations between peaks/toughs and between the locations of the flares are roughly 20 years[^3]. If we extrapolate such a repeated pattern to the period of [@Miller1979]’s dataset, we find that the time (1937-1939) of the significant flare reported in [@Miller1979] remarkably coincides with the anticipated time of recurring flares with a period of $\sim20$ years (see Figure \[fig\_miller\]). Unfortunately, the data from [@Miller1979] are very noisy and scattering, therefore, we are unable to be sure whether the variations follow the wave-like pattern seen in the light curve from 1974 to 2018. On [a]{} short timescale of months, Ark 120 also undergoes strong variations, which are commonly believed to originate from local instabilities of the accretion disk (e.g., @Czerny2008 [@Kelly2009]).

[cccccccccc]{} SPL & $-2.2\pm0.2$ & $1.6\pm0.1$ & & & & & & 0 & $8\times10^{-4}$\
DRW & $\phm{+}4.3\pm0.4$ & & & $-3.6\pm0.2$ & & & & $-2.37\pm1.00$ & $5\times10^{-4}$\
BPL & $\phm{+}3.3\pm1.3$ & $2.0\pm0.6$ & $1.6\pm1.1$ & $-3.3\pm0.8$ & & & & $\phm{-}0.31\pm0.82$ & $1\times10^{-3}$\
SPL+Periodic & $-2.1\pm0.3$ & $1.6\pm0.1$ & & & $-1.2\pm1.3$ & $-3.87\pm0.30$ & $-4.5\pm0.8$ & $-3.34\pm1.14$ &\
DRW+Periodic & $\phm{+}3.4\pm0.5$ & & & $-3.2\pm0.3$ & $-1.2\pm1.1$ & $-3.88\pm0.23$ & $-4.7\pm0.8$ & $-1.75\pm0.97$ &\
BPL+Periodic & $\phm{+}2.5\pm1.2$ & $2.3\pm0.8$ & $1.7\pm0.9$ & $-2.8\pm0.8$ & $-1.5\pm1.4$ & $-3.80\pm0.36$ & $-4.5\pm0.9$ & $-2.78\pm0.76$ &

Power Spectral Density Analysis
===============================

In Figure \[fig\_psd\], we [show]{} the PSD of the merged light curve. We resample the light curve into an even time grid, subtract the mean flux to eliminate the zero-frequency power, and then calculate the modulus squared of the discrete Fourier transform at each sampled frequency $f_j$ (@Uttley2002) $$|F_N(f_j)|^2 = \left[\sum_{i=1}^{N} x_i\cos(2\pi f_j t_i) \right]^2 + \left[ \sum_{i=1}^{N} x_i\sin(2\pi f_j t_i) \right]^2,$$ where $x_i$ is the flux at time $t_i$, $f_j=j/N\Delta T$ with $j=1, ..., N/2$, $\Delta T$ is the sampling interval, and $N$ is the number of points in the data. The PSD is defined with an adopted normalization as $$P(f_j) = \frac{2\Delta T}{N}|F_N(f_j)|^2.$$ To compare with the X-ray PSD, we also bin the logarithm of the PSD over a logarithmically spaced frequency grid, with a grid width of $\log(1.3)$ and a minimum of two power measurements per bin (@Uttley2002). The uncertainties of the binned powers are set to be the Poisson noise power.

[@Lobban2018] measured the PSD of Ark 120 in the X-ray band (0.3-10 keV) using $\sim 420$ ks [*XMM- Newton*]{} monitoring data and additional [*Swift*]{}, [*RXTE*]{}, and [*NuSTAR*]{} observation data (see their Figure 13). We adjust the normalization of the X-ray PSD to align with the optical PSD and superpose it upon the optical PSD in Figure \[fig\_psd\]. Figure \[fig\_psd\] illustrates that both the optical and X-ray PSDs can be universally described by a single power law with a slope of $\sim -1.6$, [despite some weak features in the]{} $f\times P(f)$ plot. The similar shapes of the optical and X-ray PSDs potentially imply that the optical and X-ray variabilities may have the same driving mechanism. A plausible scenario is reprocessing of X-ray emission into the UV/optical band (@Guilbert1988).

We apply three widely used PSD models to fit the optical light curve and compare the relative merits of these three models using the framework [`RECON`]{} developed in [@Li2018c]. This framework parameterizes [the]{} AGN time series in [the]{} frequency domain and transforms it back to [the]{} time domain to fit the data, and infers the model parameters using a Markov-chain Monte Carlo (MCMC) algorithm (see @Li2018c for details). Compared with methods that directly fit PSDs (e.g., @Vaughan2016), such an approach can cope with irregularly sampled light curves and also naturally take into account the sampling effects, such as spectral leakage and aliasing. In addition, the framework calculates the Bayesian evidence and hence the Bayes factor that allows us [to]{} perform model comparison. The shortcoming of [this]{} approach is that the (inverse) Fourier transform and MCMC sampling are highly computationally expensive for large numbers of data points. To save computation time, we bin the light curve every 10 days, reducing the number of points from 1473 to 471. Such a binning manipulation degrades the maximum frequency that can be [obtained]{} from the data. Since we mainly concentrate on variations with time scales of years, this manipulation does not influence our main analysis results. [It is worth stressing that in fitting the light curve, we assume independent, Gaussian measurement noise. ]{} Here, [the]{} Bayes factor is defined to be the ratio of the posterior probabilities (@Sivia2006). For two models, say $M_1$ and $M_2$, with equal prior weight, the Bayes factor is simplified to be the ratio of the corresponding Bayesian evidence, $$K = \frac{P(M_2|D)}{P(M_1|D} = \frac{P(D|M_2)}{P(D|M_1)},$$ where $P(M_{1, 2}|D)$ is the posterior probability, $P(D|M_{1, 2})$ is called the Bayesian evidence, which has been marginalized over the whole model parameters, and $D$ represent the observed data. A large value of $K>1$ means that model $M_2$ is [preferable]{} to model $M_1$. A widely used [criterion for quantifying the degree of preference is]{} given by [@Kass1995]. The three PSD adopted models are 1) the single power-law (SPL) model with a form of, $$P_{\rm SPL}(f) = A f^{-\alpha},
\label{equ_spl}$$ 2) the damped random walk (DRW) model with a form of $$P_{\rm DRW}(f) = \frac{A}{1 + (f/f_{\rm b})^2},
\label{equ_drw}$$ and 3) the bending power-law (BPL) model with a form of $$P_{\rm BPL}(f) = A \left\{\begin{array}{cc}
          (f/f_{\rm b})^{-\alpha} & {\rm for}~f > f_{\rm b},\\
          (f/f_{\rm b})^{-\beta} & {\rm otherwise}.
         \end{array}\right.
\label{equ_bpl}$$ Here $A$, $\alpha$, $\beta$, and $f_{\rm b}$ are free parameters with a restriction $\alpha > \beta$. Throughout the calculations, the prior probabilities for $A$ and $f_{\rm b}$ are set to be logarithmic and the value of $f_{\rm b}$ is limited to the frequency range of the data; the prior probabilities for $\alpha$ and $\beta$ are set to be uniform over a range (1, 5) and ($-3$, 5) respectively.

Table \[tab\_par\] summarizes the inferred parameter values for the above three PSD models. The MCMC algorithm adopted in the calculations can not calculate the uncertainty of the obtained Bayesian evidence in a single running. We [ran]{} the algorithm 10 times and set the best estimate of the evidence [based on]{} the mean and the associated uncertainty [based on]{} the standard deviation. The Bayes factor of the DRW and BPL models relative to the SPL model are $\log K=-2.37\pm1.00$ and $0.31\pm0.82$ respectively. This indicates that all the three PSD models give similarly good fits to the observed data.

![image](lsp.pdf){width="100.00000%"}

Periodicity Analysis and False-Alarm Probability
================================================

As described above, visual inspection reveals an interval of 20 years between peaks/troughs, potentially implying that there exists a $\sim20$-year periodicity in the light curve. Below we employ several methods to [attempt to confirm]{} the periodicity and estimate the associated false-alarm probability.

Periodicity {#sec_per}
-----------

We use three methods to [identify possible]{} periodicity in the light curve as follows.

-   [*Lomb-Scargle (LS) periodogram*]{}[^4]. The top right panel of Figure \[fig\_period\] shows the LS periodogram of the merged light curve, which peaks at 20.5$\pm$0.1 years. Here the uncertainty is determined from the uncertainty of the peak frequency given by the formula in @Horne1986 [see Equation (14) therein]. The LS periodogram essentially assumes a sinusoidal model for the data (@Scargle1983 [@VanderPlas2018]).

-   [*Phase dispersion minimization (PDM) analysis*]{}[[^5]]{}. The PDM method is widely used to search for non-sinusoidal periodicity by minimizing the dispersion of the phase-folded data sets (@Stellingwerf1978). The bottom right panel of Figure \[fig\_period\] plots the PDM periodogram, which is minimized at $20.5$ years, consistent with the period from the LS periodogram.

-   [*Sinusoidal fit*]{}. The sampling rate of the light curve is highly irregular and some portions have quite dense cadences. Besides, the light curve fluctuates strongly on time scales of weeks/months. To alleviate these two effects, we bin the light curve every 300 days, resulting in a total of 45 data points over a baseline of 43 years. The uncertainties of points in the newly binned light curve are assigned to be the standard deviations of each bin. The sinusoidal fitting to the binned light curve yields a period of $19.7\pm0.5$ years and sinusoidal amplitude of $(2.32\pm0.20)\times10^{-15}\,\rm erg~s^{-1}cm^{-2}\text{\AA}^{-1}$. The obtained period slightly depends on the adopted bin width, with the value changing between 19.5 and 20.0 years for the bin width at a range of 100-600 days. Throughout the calculations, we use 300 days as the default bin width. In the left panel of Figure \[fig\_period\], we show the binned light curve and the best sinusoidal fit. The original, unbinned light curve is also superposed for the sake of illustration.

    [We note that the LS periodogram is mathematically equivalent to a sinusoidal fit (e.g., @Scargle1983 [@VanderPlas2018]). Here we apply the sinusoidal fit to binned light curve data, for the purpose of alleviating the influences of the highly irregular sampling rate and strong fluctuations on short time scales. By contrast, we calculate the LS periodogram using the unbinned light curve data.]{}

In a nutshell, the above three methods generally [yield]{} periods of $\sim20$ years. We use the period $20.5$ years as the fiducial period in the [analysis that follows]{}.

![image](fig_psd_model.pdf){width="90.00000%"}

Comparison between Periodic and Aperiodic PSD Models
----------------------------------------------------

In Section \[sec\_var\], we apply three aperiodic PSD models to fit the light curve of Ark 120. We also construct periodic models by including a periodic component into the above three PSD models (@Li2018c [@Li2019]). The periodic PSD models have a composite form of $$P(f) = P_{\rm aperiodic}(f) + \frac{A_{\rm p}}{\sqrt{2\pi} \omega_{\rm p}} \exp\left[-\frac{(f-f_{\rm p})^2}{2\omega_{\rm p}^2} \right],$$ where $P_{\rm aperiodic}(f)$ is the aperiodic component that can be $P_{\rm SPL}(f)$, $P_{\rm DRW}(f)$, or $P_{\rm BPL}$ in Equations (\[equ\_spl\])-(\[equ\_bpl\]), and the second term [on]{} the right-hand side is the periodic component, which is parameterized to be a Gaussian with free parameters $A_{\rm p}$, $f_{\rm p}$, and $\omega_{\rm p}$. The corresponding period of the light curve is $1/f_{\rm b}$. Since the values of all these three parameters can span a wide range of magnitudes, their prior probabilities are set to be logarithmic. The values of $f_{\rm p}$ are additionally limited by the frequency range of the data.

We tabulate the inferred parameter values for the above periodic PSD models in Table \[tab\_par\]. Remarkably, the center of the Gaussian for the three models corresponds to a period of $\log(P/{\rm years})=1.31, 1.32, 1.25$ respectively, with a typical uncertainty of 0.30 dex. This is in agreement with the period $\log (P/{\rm year}) = 1.3$ from the LS periodogram. In Figure \[fig\_psd\_model\], we show the best recovered PSDs for both aperiodic and periodic models. The Bayes factors for the three periodic models with respect to the (aperiodic) SPL model are $\log K=-3.34\pm1.00$, $-1.75\pm0.97$ and $-2.78\pm0.76$ respectively. If adopting $\log K<-2.18$ as the criterion for rejecting a model (@Kass1995), the obtained Bayes factors imply that for the current observation data, the periodic models are not preferable [to]{} the aperiodic models.

False-alarm Probability of the Periodicity
------------------------------------------

[@Vaughan2016] showed that red-noise stochastic processes can easily produce false positive periodicity in few-cycle light curves of AGNs. The current light curve data of Ark 120 only has roughly two cycles of the period. Thus, it is important to appropriately calibrate the false alarm probability (e.g., @Koen1990 [@Barth2018; @VanderPlas2018]). False alarm probability essentially measures the probability that an aperiodic variation process (the null hypothesis) can produce spurious periodicity as strong or even stronger than that detected in the observed data. We perform null hypothesis simulations as follows. For a given aperiodic PSD model, we run `RECON` to obtain a posterior sample of model parameters, from which we randomly draw a subsample to generate a set of mock data with exactly the same sampling pattern as the observed data. To simulate the observation noises, independent Gaussian [noise with standard deviation]{} equal to the measurement errors are added into the generated [artificial]{} data. We then apply the three methods described in Section \[sec\_per\] on the [artificial]{} data and analyze the periodicity. For sinusoidal fitting, we again bin the light curve every 300 days. We identify a “periodic candidate” if the following criteria are satisfied:

1.  the LS period of the [artificial]{} data is equal to or smaller than that of the observed data, to ensure that the [artificial]{} data covers as many cycles of the period as the observed data;

2.  the peak power of the LS periodogram is larger than or equal to that of the observed data ($P_{\rm LS}\geqslant0.68$) and the minimum of the PDM periodogram is smaller than or equal to that of the observed data ($P_{\rm PDM}\leqslant0.41$);

3.  a sinusoid improves the goodness of fit compared to a constant with $\Delta \chi^2$ greater than or equal to that of the binned observed data ($\Delta \chi^2 \geqslant 150$; see Section 4.1);

4.  the LS period agrees with the periods from PDM analysis and sinusoidal fitting at a level of $\pm 1.0$ years.

The choice of the second criterion is for the purpose of selecting light curves with strong periodic variation amplitudes; the third criterion is to largely ensure roughly equal fluxes of the crests/troughs as seen in the light curve of Ark 120; and the fourth criterion aims to reduce possible false detections in the LS periodogram (@Graham2013).

The obtained false alarm probabilities are roughly $(5-10)\times10^{-4}$ (equivalent to a significance of $\sim3.3\sigma$) for the above three PSD models. Combining with the Bayes factors calculated in Section 4.2, the values of the false alarm probabilities suggest that the significance of the periodicity in Ark 120 is marginal. Continued monitoring of Ark 120 is therefore [necessary]{} to confirm the periodicity.

![Dependence of the estimated false alarm probability on the monitoring time length. [The]{} vertical dashed line represents the baseline of the observed data. Beyond the data baseline, the false alarm probability is calculated using [artificial]{} light curves for Ark 120. Note that the false alarm probabilities depend on the configurations of [artificial]{} data, which are described in Section 4.4.[]{data-label="fig_fap"}](fig_fap.pdf){width="47.00000%"}

Requirement of Temporal Baseline for the Periodicity Detection
--------------------------------------------------------------

The preceding analysis shows that the [evidence for the periodicity is inconclusive]{} for the current data. A question then arises: how long do we need to monitor Ark 120 to ascertain that the periodicity is significant at a given confidence level? To address this question, we make the following assumptions: 1) the long-term variations are sinusoidal with the parameters given by the sinusoidal fitting in Section 4.1; and 2) the variations at short time scales follow a single power-law PSD with the parameters $\log A=-2.2\pm0.2$ and $\alpha=1.6\pm0.1$ from Section 3. The single power law is set to flatten to a constant below $1.0\times10^{-3}~{\rm day}^{-1}$ to confine the generated variations to short time scales. We generate a [artificial]{} light curve by summing up the above two variation components. Following the typical configurations of the observed data, the cadence of the [artificial]{} light curve is set to be 10 days apart and the seasonal gap is set to last three months. A Gaussian noise with a standard deviation of $0.37\times10^{-15}~\rm erg~s^{-1}~cm^{-2}~\text{\AA}^{-1}$ is added into the [artificial]{} light curve. One quarter of points are discarded to mimic bad weather or [instrumental]{} problems. We then attach the [artificial]{} light curve to the observed data and use this new [artificial]{} data to calculate the false alarm probability as in Section 4.3. Appendix B illustrates three examples of [artificial]{} light curves for Ark 120. In Figure \[fig\_fap\], we show the dependence of the estimated false alarm probability on the monitoring time length for three aperiodic PSD models. The current baseline of the observed database is 43 years, corresponding to a $\sim3.3\sigma$ confidence level (see Section 4.3). Under the observation configurations described above, an additional monitoring period of $\sim$30-60 years (about 1.5-3 cycles) is required to reach a $5\sigma$ confidence level (equivalent to a false alarm probability of $6\times10^{-7}$). We stress that the adopted observation configurations are quite conservative and therefore the obtained monitoring length is just a conservative estimate. Higher data quality (e.g., higher signal to noise ratio and smaller seasonal gap) may reduce the required monitoring length.

![Two examples of the extraction of the broad H$\beta$ profile. The top panel shows the mean spectrum of our Lijiang observations between 2016 and 2017. The bottom panel shows the historical spectrum from November 22, 1976 digitalized from [@Carpriotti1982], which has narrow wavelength coverage.[]{data-label="fig_decomp"}](spec_decomp.pdf){width="48.00000%"}

![(Left) H$\beta$ profiles from 1976 to 2017. (Right) Folded H$\beta$ profile series with phases for a complete cycle using a period [of]{} 20.5 years. Black lines represent the overall mean profiles of each phase. Blue and red lines represent profiles from the first and second cycles of the series respectively. Vertical lines show the wavelength center of [the]{} H$\beta$ line.[]{data-label="fig_prof"}](prof_all.pdf){width="45.00000%"}

H$\beta$ Profile Variations
===========================

We also compile the historical H$\beta$ profile of Ark 120. The H$\beta$ profiles are mainly from [@Stanic2000] and [@Doroshenko2008], which presented 10 averaged H$\beta$ profiles over 1977-1990 (see their Figure 3-11) and 21 H$\beta$ profiles over 1992-2005 (see their Figure 5), respectively. In addition, [@Carpriotti1982] [displayed]{} 6 spectra of Ark 120 over 1976-1981 (see their Figure 3), [@Korista1992] presented 12 spectra over 1981-1989 (see their Figures 1-3), [@Peterson1998] [show]{} one averaged spectrum [from]{} around 1993 (see their Figure 1), [and [@Afanasiev2019] presented one spectrum from 2014 (see their Figure 5)]{}. We digitalize these figures to obtain the profile data. Since we only concentrate on the shape of H$\beta$ profile, the absolute fluxes are no longer important. [Among these references, [@Stanic2000], [@Doroshenko2008], and [@Afanasiev2019] had already extracted the broad H$\beta$ profiles. We directly use their extracted H$\beta$ profiles. For the other spectra, we isolate]{} the broad H$\beta$ profiles following the procedure in [@Stanic2000]. The narrow H$\beta$ line and emission are subtracted using a simple spectral decomposition (see @Stanic2000 for a detail). [The bottom panel of Figure \[fig\_decomp\] shows an example for decomposing the historical spectrum from November 22, 1976 digitalized from [@Carpriotti1982].]{}

![image](fig_prof_period.pdf){width="95.00000%"}

![(Top) H$\beta$ line centroid, (middle) H$\beta$ line dispersion, and (bottom) continuum fluxes of Ark 120. The unit of the continuum flux densities is $10^{-15}~\rm erg~s^{-1}~cm^{-2}~\text{\AA}^{-1}$. The H$\beta$ line centroid and dispersion are calculated over a wavelength range (4741-4981) Å.[]{data-label="fig_hb"}](fig_hb.pdf){width="48.00000%"}

[[@Haas2011] displayed one spectrum from 2009 and [@Marziani1992] presented four spectra from between 1989 and 1990, with both these studies providing us with their electronic data. We queried the NASA/IPAC Extragalactic Database and found one spectrum from the 6dF galaxy survey (@Jones2009) and one spectrum from the updated Zwicky catalog (@Falco1999). We supply two mean H$\beta$ profiles from our two-year monitoring. These spectra are of sufficiently good quality to perform a detailed spectral decomposition. The top panel shows an example of spectral decomposition for the mean spectrum between 2016 and 2017 from our Lijiang observations. The decomposition includes a power law for the continuum, two Gaussians for the narrow \[\] doublet and H$\beta$ line, three Gaussians for the broad emission lines (H$\beta$ and H$\gamma$), an model shelf (@Kovacevi2010), and a host galaxy starlight model based on the ELODIE library (@Prugniel2001 [@Prugniel2007a; @Prugniel2007b]). Some other emission lines, such as , , etc., are also included when necessary (modelled by one or two Gaussians). The flux ratio between \[\] $\lambda$5007 and $\lambda$4959 is fixed to be 3. A Legendre polynomial of the first order is multiplied to the sum of the above components to account for the errors in the Galactic extinction, flux calibration, or any other cause that affects the shape of the continuum. The fitting is implemented using the ULySS code[^6] (@Koleva2009). More details of the spectral decomposition were described in [@Bon2014; @Bon2016].]{}

Furthermore, [@Denissyuk2015] [collected]{} spectroscopic monitoring data of Ark 120 over 1976-2013 using a 70-cm telescope. We use 16 [of their]{} spectra between 2005 and 2011. The spectra are relatively noisy and we decompose the broad H$\beta$ line simply by subtracting the underlying continuum linearly interpolated between 4748 [Å]{} and 5020 [Å]{}. The narrow H$\beta$ and \[\] lines are not subtracted.

[In total, we obtain 76 profiles that cover a time span of 40.1 years from 1976 to 2017 and show them in Figure \[fig\_prof\]. We bear in mind that our compilation is incomplete and does not include those studies that investigated spectral variations of Ark 120 but did not display the spectra.]{} The profiles are highly asymmetric and have double peaks (e.g., around 1981 and 1989), which seem to shift apart and merge again in some epochs. Moreover, the amplitudes of the red and blue peaks vary independently. These strong variations of H$\beta$ profiles potentially reflect that the region emitting the broad H$\beta$ line ([the]{} so called broad-line region) undergoes notable changes. Such behavior had also been found in the H$\beta$ profiles of NGC 5548 (@Li2016 [@Bon2016]). [To determine whether there exist periodic variations in the H$\beta$ profiles, we first normalize the H$\beta$ profiles to a uniform integrated flux. We divide the H$\beta$ wavelength region between 4741 and 4981 [Å]{} into 12 bins and obtain the “flux” in each wavelength bin. We then calculate the “multiband” LS periodogram that combines the variation information from all the bins using the algorithm developed by [@VanderPlas2015]. The right panel of Figure \[fig\_prof\_period\] plots the “multiband” LS periodogram, which peaks around 20 years. This is remarkably consistent with the period found in the continuum light curve. For the sake of comparison, the left panel of Figure \[fig\_prof\_period\] shows the standard LS periodogram in each wavelength bin. The periodograms show diverse patterns among wavelength bins, meaning that it is impossible to detect a period in individual bins under the current data quality. ]{} We fold the profiles into 10 uniformly spaced phases using a period of 20.5 years in the right panel of Figure \[fig\_prof\]. The phase bin width is $\sim2$ years and each phase bin has on average six profiles. It is more evident that in the folded profile series the red and blue peaks shift with phases significantly and their amplitudes vary independently. For the sake of comparison, we plot the H$\beta$ profiles from the first cycle (1976-1996) and second cycle (1996-2017) in blue and red respectively, and the overall mean profile in black. The profiles from the two cycles are generally matched in terms of the separation and amplitude of the double peaks, [further supporting]{} that H$\beta$ profiles may change with a similar period as the continuum and H$\beta$ fluxes.

We calculate the H$\beta$ line centroid and dispersion over a wavelength range (4741-4981) [Å]{} from the above spectra and show their changes with time in Figure \[fig\_hb\]. Despite large scattering, both the H$\beta$ line [centroid]{} and dispersion generally exhibit a wave-like pattern as in the continuum variations. It seems that the centroid and line dispersion increase with the continuum flux density. This is difficult to explain if assuming that the broad-line region is virialized and its size obeys the tight scaling correlation between sizes of broad-line regions and AGN luminosities found in reverberation mapping observations (e.g. @Kaspi2000 [@Bentz2013; @Du2016]). According to this scaling correlation, the size of the broad-line region expands as the AGN luminosity increases, leading to the line width decreasing if the broad-line region is virialized. However, we stress that the present spectra data are compiled from various observations with diverse spectral resolutions, wavelength calibrations, and signal-to-noise ratios. Some of the digitalized spectra only [cover]{} the wavelength region of H$\beta$ line so that it is impossible to recalibrate all the spectra in a self-consistent way. Therefore, more high-quality spectroscopic observations are needed to [test the results]{} in this section. [In Appendix A, we present a portion of the online table for the H$\beta$ line centroid and dispersion data. Again, because it is difficult to reliably estimate the measurement noise for the digitalized spectra, we do not report the uncertainties for the calculated H$\beta$ line centroids and dispersions.]{}

Discussions on the Periodicity
==============================

Regarding the periodicity in AGN variations, there are indeed a variety of theoretic models/interpretations, including the SMBH scenario (see the discussions in @Li2016, @Lu2016, and @Bon2016). Some of the interpretations can be directly excluded in Ark 120. The precessing jet model is implausible as Ark 120 is radio-quiet (@Condon1998 [@Ho2002]) so that the optical continuum is most likely dominated by the disk emission; The other interpretations may be tested with the aid of high-quality reverberation mapping observations of broad emission lines (e.g., @Shen2010 [@Wang2018]). In particular, for the SMBH binary scenario, there are two factors that lead to distinct reverberation mapping signatures: first, the geometry and dynamics of the broad-line region(s) are different compared to [those]{} surrounding a single black hole, considering the complicated gravitational potential jointly governed by both black holes (e.g., @Sepinsky2007 [@Popovic2012]); second, if both black holes are active, there are doublet ionizing sources that illuminate the broad-line region(s), giving rise to distinguishable patterns in velocity-delay maps of broad emission lines (@Wang2018). [Interestingly, a recent reverberation mapping campaign by [@Du2018] found a complicated velocity-delay map for the H$\beta$ broad-line region of Ark 120, which shows a general decreasing trend from the blue ($-3000~{\rm km~s^{-1}}$) to the red ($4000~{\rm km~s^{-1}}$) wing but with a local peak around $1000$-$2000~{\rm km~s^{-1}}$. A detailed investigation is needed to uncover the geometry and dynamics underlying such a velocity delay map.]{}

Ark 120 has a stellar velocity dispersion of $\sigma_\star=192\pm8~\rm km~s^{-1}$ (@Woo2013), resulting in a black hole mass estimate of $M_\bullet=(2.6\pm0.2)\times10^8M_\odot$ from the $M_\bullet-\sigma_\star$ relation for classical bulges (@Kormendy2013). By assuming that an SMBH binary resides at the center of Ark 120, the semi-major axis of the binary’s orbit is $a_\bullet = 27.0$ light-days (0.02 pc or 33 $\mu$as) if using 20.5 years (in the observed frame) as the orbital period. Such a separation can be [potentially]{} spatially resolved either by the Event Horizon Telescope with an angular resolution of $\sim20~\mu$as (@Fish2016) if each black hole produces radio emission, or by [the]{} [*Gaia*]{} satellite with an angular resolution down to $9~\mu$as (@Dorazio2018). On the other hand, the characteristic strain amplitude of the gravitational wave emission is [of]{} the order of $h_s\sim10^{-17}$, which is in the sensitivity [range]{} of next-generation pulsar timing arrays (@Janssen2015). The previous H$\beta$ reverberation mapping observations of Ark 120 show a typical size of 30-40 light-days for the H$\beta$ broad-line region (@Peterson1998 [@Du2018]), indicating that the binary orbit is slightly smaller than or comparable with the H$\beta$ broad-line region. In such cases, the interaction between the binary black holes and the broad-line regions is therefore important to shape the broad-line profiles.

Conclusion
==========

We compile the historical archival photometric and spectroscopic data of Ark 120 over four decades. The long-term variations of both the optical continuum and H$\beta$ integrated fluxes exhibit a wave-like pattern. We analyze the periodicity using various methods and generally obtain a period of $\sim20$ years. The estimated false alarm probability is about $1\times10^{-3}$. The comparison between aperiodic and periodic models based on the Bayes factors suggests that the periodicity is [inconclusive using]{} the current data. Continued monitoring of Ark 120 is [needed]{} to track more cycles to eliminate the false positive rate and confirm the periodicity. The broad H$\beta$ line shows double peaks, which vary strongly with time. [Using the “multiband” LS periodogram developed by [@VanderPlas2015], we find that the overall H$\beta$ profile also varies with the same suggested period in the continuum.]{} These observations make Ark 120 to be [one of]{} the nearest AGNs with possible periodic variability, a remarkable analogy to NGC 5548 (@Li2016 [@Bon2016]).

Although the [evidence for periodicity is inconclusive using the current data, its abundant spectroscopic monitoring data still make Ark 120]{} a good laboratory for studying the origin of asymmetric, rapidly varying broad emission lines in general and evolution of SMBH [binaries, particularly]{} if the possible periodic variations stem from the binary’s orbital motion.

a portion of the data and the entire tables are available in a machine-readable form in the online journal. [Note that most of the compiled spectra are digitalized from published figures with insufficient wavelength coverage to estimate the measurement noise, so we do not report the uncertainties of the calculated H$\beta$ line centroids and dispersions.]{}

[cccc]{} 42392.0000 & 11.68 & 0.84 & P83\
42392.4680 & 8.97 & 0.72 & D99b\
42396.0000 & 10.78 & 0.77 & P83\
42415.0000 & 9.96 & 0.72 & P83\
42476.0000 & 9.12 & 0.67 & P83\
42685.0000 & 6.85 & 0.59 & P83\
...

[ccc]{} 43104.9553 & 4866.66 & 32.27\
43494.3051 & 4873.59 & 35.52\
43879.0498 & 4869.62 & 44.00\
44168.3578 & 4865.34 & 41.01\
44228.8177 & 4872.53 & 38.48\
44545.9451 & 4865.19 & 43.10\
...

Examples of [Artificial]{} Light Curves
=======================================

Figure \[fig\_artificial\] shows three examples of [artificial]{} light curves for Ark 120, which are used to determine the required baseline to detect the periodicity at a given confidence level. The procedure for generating [artificial]{} light curves is detailed in Section 4.4.

![image](fig_artificial.pdf){width="90.00000%"}

[99]{} Alloin, D., Boisson, C., & Pelat, D. 1988, , 200, 17 Astropy Collaboration, Price-Whelan, A. M., Sip[ő]{}cz, B. M., et al. 2018, , 156, 123 Barth, A. J., & Stern, D. 2018, , 859, 10 Bentz, M. C., Denney, K. D., Grier, C. J., et al. 2013, , 767, 149 Bentz, M. C., Peterson, B. M., Netzer, H., Pogge, R. W., & Vestergaard, M. 2009, , 697, 160 Bogdanovi[ć]{} T., Eracleous M., Sigurdsson S., 2009, NewAR, 53, 113 Bogdanovi[ć]{}, T., Smith, B. D., Sigurdsson, S., & Eracleous, M. 2008, , 174, 455 Bon, E., Jovanovi[ć]{}, P., Marziani, P., et al. 2012, , 759, 118 Bon, E., Zucker, S., Netzer, H., et al. 2016, , 225, 29 Bon, N., Popovi[ć]{}, L. [Č]{}., & Bon, E. 2014, Advances in Space Research, 54, 1389 Boroson, T. A., & Green, R. F. 1992, , 80, 109 Capriotti, E. R., Foltz, C. B., & Peterson, B. M. 1982, , 261, 35 Charisi, M., Bartos, I., Haiman, Z., et al. 2016, , 463, 2145 Condon, J. J., Yin, Q. F., Thuan, T. X., & Boller, T. 1998, , 116, 2682 Czerny, B., Siemiginowska, A., Janiuk, A., & Gupta, A. C. 2008, , 386, 1557 Denissyuk, E. K., Valiullin, R. R., & Gaisina, V. N. 2015, Astronomy Reports, 59, 123 D’Orazio, D. J., & Loeb, A. 2018, arXiv:1808.09974 Doroshenko, V. T., Sergeev, S. G., & Pronik, V. I. 2008, Astronomy Reports, 52, 442 Doroshenko, V. T., & Lyuty, V. M. 1999a, Astronomy Letters, 25, 771 Doroshenko, V. T., Sergeev, S. G., Pronik, V. I., & Chuvaev, K. K. 1999b, Astronomy Letters, 25, 569 Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, , 696, 870 Du, P., Brotherton, M. S., Wang, K. et al. 2018, ApJ, 869, 142 Du, P., Hu, C., Lu, K.-X., et al. 2014, , 782, 45 Du, P., Lu, K.-X., Zhang, Z.-X., et al. 2016, , 825, 126 Falco, E. E., Kurtz, M. J., Geller, M. J., et al. 1999, , 111, 438 Fish, V., Akiyama, K., Bouman, K., et al. 2016, Galaxies, 4, 54 Gliozzi, M., Papadakis, I. E., Grupe, D., Brinkmann, W. P., & R[ä]{}th, C. 2017, , 464, 3955 Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015a, , 453, 1562 Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015b, , 518, 74 Graham, M. J., Drake, A. J., Djorgovski, S. G., et al. 2013, , 434, 3423 Guilbert, P. W., & Rees, M. J. 1988, , 233, 475 Guo, D., Tao, J., & Qian, B. 2006, , 58, 503 Haas, M., Chini, R., Ramolla, M., et al. 2011, , 535, A73 Ho, L. C. 2002, , 564, 120 Horne, J. H., & Baliunas, S. L. 1986, , 302, 757 Janssen, G. H., et al. 2015, in Proc. Advancing Astrophysics with the Square Kilometre Array (AASKA14), 37 Johnson, H. L. 1966, , 4, 193 Jones, D. H., Read, M. A., Saunders, W., et al. 2009, , 399, 683 Kaspi, S., Smith, P. S., Netzer, H., et al. 2000, , 533, 631 Kass, R. E., & Raftery, A. E. 1995, Journal of the American Statistical Association, 90, 791 Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, , 698, 895 Kelly, B. C., Becker, A. C., Sobolewska, M., Siemiginowska, A., & Uttley, P. 2014, , 788, 33 Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, , 129, 104502 Koleva, M., Prugniel, P., Bouchard, A., & Wu, Y. 2009, , 501, 1269 Kollatschny, W., Fricke, K. J., Schleicher, H., & Yorke, H. W. 1981a, , 102, L23 Kollatschny, W., Schleicher, H., Fricke, K. J., & Yorke, H. W. 1981b, , 104, 198 Kormendy, J., & Ho, L. C. 2013, , 51, 511 Korista, K. T. 1992, , 79, 285 Koshida, S., Minezaki, T., Yoshii, Y., et al. 2014, , 788, 159 Kova[č]{}evi[ć]{}, J., Popovi[ć]{}, L. [Č]{}., & Dimitrijevi[ć]{}, M. S. 2010, , 189, 15 Kuehn, C. A., Baldwin, J. A., Peterson, B. M., & Korista, K. T. 2008, , 673, 69 Li, Y.-R. 2018a, CyPDM: A fast package to apply the phase disperion minimization algorithm, v0.1.1, Zenodo, doi:10.5281/zenodo.2527265 Li, Y.-R. 2018b, RECON: Measuring spectral power and reconstructing time series in active galactic nuclei, v0.1.0, Zenodo, doi:10.5281/zenodo.2528155 Li, Y.-R., & Wang, J.-M. 2018, , 476, L55 Li, Y.-R., & Wang, J.-M. 2019, , resubmitted Li, Y.-R., Wang, J.-M., Ho, L. C., et al. 2016, , 822, 4 Li, Y.-R., Wang, J.-M., Ho, L. C., Du, P., & Bai, J.-M. 2013, , 779, 110 Li, Y.-R., Wang, J.-M., Hu, C., Du, P., & Bai, J.-M. 2014, , 786, L6 Lindegren, L., Hern[á]{}ndez, J., Bombrun, A., et al. 2018, , 616, A2 Liu, T., Gezari, S., Burgett, W., et al. 2016, , 833, 6 Liu, T., Gezari, S., Heinis, S., et al. 2015, , 803, L16 Lobban, A. P., Porquet, D., Reeves, J. N., et al. 2018, , 474, 3237 Lu, K.-X., Li, Y.-R., Bi, S.-L., & Wang, J.-M. 2016, , 459, L124 MacLeod, C. L., Ivezi[ć]{}, [Ž]{}., Kochanek, C. S., et al. 2010, , 721, 1014 Miller, H. R. 1979, , 35, 387 Miranda, R., Mu[ñ]{}oz, D. J., & Lai, D. 2017, , 466, 1170 Nardini, E., Porquet, D., Reeves, J. N., et al. 2016, , 832, 45 Peterson, B. M., Denney, K. D., De Rosa, G., et al. 2013, , 779, 109 Peterson, B. M., Korista, K. T., & Wagner, R. M. 1989, , 98, 500 Peterson, B. M., Pogge, R. W., Wanders, I., Smith, S. M., & Romanishin, W. 1995, , 107, 579 Peterson, B. M., Wagner, R. M., Crenshaw, D. M., et al. 1983, , 88, 926 Peterson, B. M., Wanders, I., Bertram, R., et al. 1998, , 501, 82 Prugniel, P., Koleva, M., Ocvirk, P., Le Borgne, D., & Soubiran, C. 2007a, Stellar Populations as Building Blocks of Galaxies, 241, 68 Prugniel, P., & Soubiran, C. 2001, , 369, 1048 Prugniel, P., Soubiran, C., Koleva, M., & Le Borgne, D. 2007b, VizieR Online Data Catalog, 3251, Popovi[ć]{}, L. [Č]{}, Stani[ć]{}, N., Kubi[č]{}ela, A., & Bon, E. 2001, , 367, 780 Popovi[ć]{}, L. [Č]{}. 2012, , 56, 74 Reeves, J. N., Porquet, D., Braito, V., et al. 2016, , 828, 98 Roberts, C. A., & Rumstay, K. R. 2012, Journal of the Southeastern Association for Research in Astronomy, 6, 47 Sepinsky, J. F., Willems, B., & Kalogera, V. 2007, , 660, 1624 Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, , 788, 48 Shen, Y., & Loeb, A. 2010, , 725, 249 Sivia, D., & Skilling, D. 2006, Data Analysis: A Bayesian Tutorial (New York: Oxford Univ. Press) Stanic, N., Popovic, L. C., Kubicela, A., & Bon, E. 2000, Serbian Astronomical Journal, 162, 7 Stellingwerf, R. F. 1978, , 224, 953 Uttley, P., McHardy, I. M., & Papadakis, I. E. 2002, , 332, 231 Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, , 452, 851 Vanderplas, J. T. 2015, gatspy: General tools for Astronomical Time Series in Python, v0.1.1, Zenodo, doi:10.5281/zenodo.14833 VanderPlas, J. T., & Ivezi[ć]{}, [Ž]{}. 2015, , 812, 18 Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, , 461, 3145 Vaughan, S., Fabian, A. C., Ballantyne, D. R., et al. 2004, , 351, 193 Wang, J.-M., Songsheng, Y.-Y., Li, Y.-R., & Zhe, Y. 2018, , 862, 171 Winge, C., Peterson, B. M., Pastoriza, M. G., & Storchi-Bergmann, T. 1996, , 469, 648 Winkler, H. 1997, , 292, 273 Winkler, H., Glass, I. S., van Wyk, F., et al. 1992, , 257, 659 Woo, J.-H., Schulze, A., Park, D., et al. 2013, , 772, 49 Zheng, Z.-Y., Butler, N. R., Shen, Y., et al. 2016, , 827, 56

[^1]: The redshift is taken from the NASA/IPAC Extragalactic Database, corresponding to a luminosity distance of 143 Mpc based on the [*WMAP-5*]{} cosmological parameters $H_0=70.5~\rm km~s^{-1}~Mpc^{-1}$, $\Omega_m=0.27$, and $\Omega_\Lambda=0.73$.

[^2]: It is possible that the narrow \[\] line has a long-term secular variation as in NGC 5548 (@Peterson2013). We do not include such a variation in our intercalibration as there were no reliable absolute flux measurements over the whole period of the datasets.

[^3]: We note that some numerical simulations of supermassive black hole binary systems produce a similar variation pattern (crests/troughs and flares) in mass accretion rates onto the binary black holes (e.g, @Bogdanovic2008 [@Bogdanovic2009; @Miranda2017]).

[^4]: We use the module `LombScargle` in the Python package `astropy` (@Astropy2018) to calculate [the]{} LS periodogram. The argument `normalization=’standard’` is switched on to implement the “least-square normalization” of the periodogram. This leads the power of the periodogram to lie between 0 and 1.

[^5]: We use the Python module `CyPDM` to implement PDM analysis (@Li2018a).

[^6]: The ULySS code is available at <http://ulyss.univ-lyon1.fr>.
