Two-dimensional correlation analysis of periodicity in active galactic nuclei time series

The active galactic nuclei (AGN) are among the most powerful sources with an inherent, pronounced and random variation of brightness. The randomness of their time series is so subtle as to blur the border between aperiodic fluctuations and noisy oscillations. This poses challenges to analysing of such time series because neither visual inspection nor pre-exisitng methods can identify well oscillatory signals in them. Thus, there is a need for an objective method for periodicity detection. Here we review our a new data analysis method that combines a two-dimensional correlation (2D) of time series with the powerful methods of Gaussian processes. To demonstrate the utility of this technique, we apply it to two example problems which were not exploited enough: damped rednoised artificial time series mimicking AGN time series and newly published observed time series of changing look AGN (CL AGN) NGC 3516. The method successfully detected periodicities in both types of time series. Identified periodicity of $\sim 4$ yr in NGC 3516 allows us to speculate that if the thermal instability formed in its accretion disc (AD) on a time scale resembling detected periodicity then AD radius could be $\sim 0.0024$ pc.


Introduction
Active galactic nuclei (AGNs) vary on time-scales ranging from minutes and hours to years over the entire electromagnetic spectrum, with no apparent indications of periodicities. Nevertheless, in recent years an increasing number of reports on AGN periodicities have been published [see e.g. [1][2][3][4][5][6][7], suggesting that supermassive binary black holes might be detected through periodicity of their observed time series (i.e. light curves, [see 8, and references therein]).
Establishing a method for recurrent patterns detection in the AGN time series is an important step towards this goal. Many methods have been designed for estimating this periodicity [for an excellent review see 9]. These methods share a number of commonalities, as well as differences. Most of them are some variant of Fourier analysis [10] which has restrictive assumptons: equally spaced observations, the time series is stationary, homoscedastic Gaussian noise with purely periodic signals (i.e. sinusoidal shape). Wavelet analysis does not assume stationarity and is therefore able to detect amplitude and period changes over time. Usually in all Fourier based methods peaks which are indicating periodicity can overlap. However, the Fourier transform, the wavelet transform and related period estimation techniques can not tell about the presence of coordinated or independent changes among signals, as well as about relative directions of signal intensity variations. Our hybrid method based on two-dimensional (2D) correlation analysis were devised to deal with above issues [5].
We aim to further illustrate the performance and application of this 2D hybrid method on synthetic data, where results can be judged carefully and also on observed data, where new insights can be gained. We present computations of 2D correlation maps of damped rednoised artificial time series and newly published long-term monitored time series of a changing look (CL) AGN NGC 3516 [11]. There is some indication that better sampled AGN light curves can be modeled as damped harmonic oscillator perturbed by coloured noise [12]. CL AGNs are objects showing the dramatic variability of the emission line profiles and the change of classification type within very short time interval (from days to years). Periodical variability has been discussed for some well-known CL AGN such as NGC 4151 [see 5,13], NGC 5548 [see 5,14] within the context of supermassive binary black hole candidate and pointed out as possibility for typical CL AGN NGC 2617 [see 15]. It would be interesting to know if CL AGN variability is periodic, since it can be a consequence of tidal disruption events [16] or recoiling supermassive black hole [17].

Data and Method
To demonstrate the utility of this technique, we apply it to two example problems which were not exploited enough.
Interestingly, [12] made estimate that availability of better sampled data in the future will necessitate more sophisticated models for AGN light curves as it is damped harmonic oscillator perturbed by colored noise. Thus, we synthetised artificial damped sinusoid signal corrupted by red noise. Figure 1 (left panels) shows both damped sinusoid and normal sinusoid of period 125 arbitrary units [a.u.]. Normal sinusoid is symmetric with respect to time axis. But by introducing damped motion we brake this symmetry. Although the symmetry can be broken in ways that make it difficult to recognize or reconstruct periodicity, it is there nonetheless. And, of course, the more dramatic and complex the nature and magnitude of the damper, the more complex the task of identifying the original periodicity (symmetry). We perturbed both signals with red noise (see right panels) so that the original signal patterns are no longer able to be recognized. Precise mathematical description of red noise corruption of a signal is given in [5] In a recent study of CL AGN NGC 3516 [11], authors applied Lomb Scargle periodogram in order to detect periodicity in observed light curves of this object. However, potential periodic signals were of small significance. Thus, we applied our method to the long-term monitored continuum, Hα and Hβ fluxes covering 22 yr (from 1996 to 2018). The data are presented and carefully described in [11] (see their Figure 5). Thus, we will not repeat the details here.
Here we recapitulate key aspects of our hybrid method, which is discussed in details in [5]. A conversion of set of light curves to 2D correlation maps is relatively easy, providing rich information about the presence of coordinated or independent signals as well as relative directions of signal variations. Some notable features of our 2D hybrid method are: simplification of complex spectra consisting of many overlapped peaks, enhancement of apparent spectral resolution by spreading peaks over the second dimension, and establishment of direction of changes in signal through correlation coefficients. Some generic properties of correlation map are marked on Figure 2. Our hybrid method produces a contour map of correlation intensity on a period plane defined by two independent period axes corresponding to the two light curves. Peaks on the main diagonal represent the simultaneous change in signals at the same period. Cross peaks located at the off diagonal positions represent the simultaneous change in signals at two different periods. However, we have never observed these cross peaks in objects we analyzed up to now. Positive correlation value means the two periodic signals increase or decrease together in the same direction.

Results and Discussion
Firstly, we present a 2D correlation map (see Figure 3) of a pair of synthetic curves (a sinusoid and damped sinusoid of period 125 [a.u] corrupted with red noise, see right column of Figure 1) to illustrate this approach for establishing the presence of common periodicities. Because of no overlap in this map, the peaks can be readily assigned to specific periodic signals. Peak appearing at the lower left of diagonal positions of the 2D map indicate a clear positive correlation at the period of 122 ± 26[a.u.] with correlation coefficient > 0.8. However, the peak in upper right corner is related to the Nyqist frequency and is not relevant. Clearly the width of waves in damped sinusoid is changing over time and they become wider then waves of sinusoid. As a consequence the signals are not in perfect phase. Due to this, Unexpectedly, the strong relationships are present at periods 1580 ± 743 days (i.e. 4.32 ± 2.04yr [compare to Figure 11 in 11]) and 1385 ± 128 days (i.e. 3.8 ± 0.4yr, see Figure 5). Calculated correlation maps provide an even clearer picture of the time-dependent changes of periodicities in the light curves than Lomb periodogram [compare to Figure 11  These periodicities are similar to each other, indicating that continuum, Hα and Hβ fluctuate with similar periodicity characteristics. It seems that periodicity in continuum, triggers similar variation in Hα and Hβ emission lines.
Perturbed AGN accretion disc can be used as an explanation of continuum flux and corresponding spectral variability [see e.g. 18]. It is believed that such perturbations induce the thermal instability in the disc which time scale is [see 19] where r g = GM BH /c 2 is a gravitational radius and r d is accretion disc dimension.

Conclusion
In the present study, we show a new application of our hybrid method for periodicity detection in AGN time series. We extended the results obtained in [5,7], via the analysis of synthetic sinusoidal and damped sinusoidal signals corrupted with red noise as well as observed continuum, Hα and Hβ light curves of CL AGN NGC 3516. Our hybrid method successfully recovered the period in synthetic time series (i.e. light curves). Further, it detected periods of 1580 ± 743 days in the continuum and Hα as well as 1385 ± 128 days in the the continuum and Hβ of NGC 3516. If the thermal instability had formed in the accretion disc of NGC 3516 on a time scale resembling detected periodicity (∼ 4 yr) we inferred that accretion disc radius is ∼ 0.0024 pc. This agrees well with [20] prediction of accretion disc model for this object and with [11] evidence of its dimension variability.
Both experimental results show the robustness of our method against problems of damped oscillations corrupted with red noise and complex time series of CL AGN NGC 3516. However, in order to validate periodicity detection in NGC 3516 time series, a further decade long-term monitoring is needed.