{'fraction_non_alphanumeric': 0.04534922211188348, 'fraction_numerical': 0.04592850049652433, 'mean_word_length': 4.263283972125436, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Quasi-periodic Oscillation in the gamma-ray Emission from the Non-blazar Active Galactic Nucleus PKS 0521-36

Long-term gamma-ray variability of a non-blazar Active Galactic Nucleus (AGN) PKS 0521-36 is investigated by using Fermi-LAT pass 8 data covering from 2008 August to 2021 March. The results show that the histogram of the gamma-ray fluxes follows a log-normal distribution. Interestingly, in the analysis of about 5.8-year (from MJD 56317 to 58447) LAT data between two outbursts (occurring during 2012 October and 2019 May respectively), a quasi-periodic oscillation (QPO) with a period of about 1.1 years (about 5 sigma of significance) is found in the Lomb-Scargle Periodogram (LSP), the Weighted Wavelet Z-transform (WWZ) and the REDFIT results. This quasi-periodic signal also appears in the results of Gaussian process modeling the light curve. Therefore, the robustness of the QPO is examined by four different methods. This is the first gamma-ray QPO found in a mildly beamed jet. Our results imply that the gamma-ray outbursts play an important role in the formation of the gamma-ray QPO.


INTRODUCTION
It is generally accepted that Active Galactic Nuclei (AGNs) existing in the center of a small fraction of galaxies are powerful emitters for the entire galaxies (e.g., Padovani et al. 2017). Usually, an AGN is powered by a central super massive black hole (SMBH) and emits variable emission covering from radio to high-energy gamma-ray bands (e.g., Barkov et al. 2010;Madejski & Sikora 2016). Non-thermal GeV-TeV γ-rays from AGNs are produced in their powerful relativistic jets.
The largest class of AGNs detected by Fermi Large Area Telescope (Fermi-LAT) is of blazar type (Ajello et al. 2020). A typical blazar has relativistic jet oriented at very small angle (θ <5 • ) to the line of our sight (Ackermann et al. 2015), and therefore its jet emission is strongly Doppler-boosted (Urry & Padovani 1995). Rapid flares on timescale of a few minutes in blazars have been observed at GeV γ-ray energies (e.g., Ackermann et al. 2016;Meyer et al. 2019) and TeV energies (e.g., Aleksić et al. 2011). This kind of rapid flare suggests an extremely compact γ-ray emission region which likely locates in the broad-line region (BLR) of flat spectrum radio quasars (FSRQs). However, no signs of γ-ray absorption by BLR photons in γ-ray spectrum suggests that γ-ray emission region should be outside of the BLR (e.g., Costamante et al. 2018). To reconcile the observations of fast variability and simultaneous energy spectrum, new models are proposed (see Böttcher 2019;, for recent reviews). Moreover, timing analysis, such as the exploration of periodicity, would provide important information about the nature of these sources. Quasi-periodic variabilities of blazars in γ-ray energies have been investigated, and possible quasi-periodic variabilities in more than 20 blazars have been reported based on analyses of Fermi-LAT data (e.g., Sandrinelli et al. 2014Sandrinelli et al. , 2016Sandrinelli et al. , 2017Ackermann et al. 2015;Prokhorov & Moraghan 2017;Zhang et al. 2017a,c,b;Zhou et al. 2018;Bhatta 2019;Peñil et al. 2020;Zhang et al. 2020). However, the cause for the γ-ray quasi-periodic variabilities remains controversial. Yan et al. (2018) investigated the relation between γ-ray variability amplitude and spectral index for PG 1553+113, and suggested that the γ-ray quasi-periodic variability of PG 1553+113 (Ackermann et al. 2015) could be caused by a periodic modulation in Doppler factor. Several possible explanations have been proposed for the γ-ray quasi-periodic variabilities in blazars, for instance, presence of a binary system of SMBHs (Valtonen et al. 2016;Tavani et al. 2018), modulation in the accretion flow feeding the jet (Gracia et al. 2003), and a periodic lighthouse effect (Camenzind & Krockenberger 1992).
Seventy non-blazar AGNs are included in the Fourth Catalog of AGNs detected by Fermi-LAT (4LAC; Ajello et al. 2020). Among them, there are 41 radio galaxies (RGs), 9 narrow-line seyfert 1 galaxies (NLSy1s), 2 steep-spectrum radio quasars (SSRQs), 5 compact steep spectrum radio sources (CSSs), 1 Seyfert galaxy, and 11 other AGNs. 6 NLSy1s, 7 RGs, 3 CSSs, 1 SSRQ, and 3 other AGNs are found to be variable (Ajello et al. 2020). The jet properties of NLSy1 are similar to those of FSRQs, but with lower powers (e.g., D'Ammando et al. 2013D'Ammando et al. , 2015aD'Ammando et al. , 2016Zhang et al. 2013;Paliya et al. 2019). γ-ray flares on a daily timescale are found in several NLSy1s (D'Ammando 2019). RGs are misaligned AGNs, and therefore the jet Doppler boosting effect is not strong. Intraday variability is found in the Fermi-LAT data of the RG NGC 1275 (Sahakyan et al. 2018;Tanada et al. 2018). All of the non-blazar AGNs are not bright enough to have good sampling in their long-term γ-ray light curves, except for NGC 1275 and PKS 0521-36. The γ-ray variability and spectral evolution of NGC 1275 have been studied in detail (e.g., Baghmanyan et al. 2017;Sahakyan et al. 2018;Tanada et al. 2018). Here, we focus on the γ-ray properties of PKS 0521-36. A combination of γ-ray properties of blazars and non-blazar AGNs can provide insightful information on AGN jet physics.
PKS 0521-36 is a peculiar non-blazar AGN with broad emission lines in the optical and ultraviolet bands and the steep radio spectrum (Scarpa et al. 1995;Keel 1986). D'Ammando et al. (2015b) raised a possibility that PKS 0521-36 is an intermediate source between broad-line radio galaxy and SSRQ. It was previously classified as an N galaxy, then a BL Lac object (D'Ammando et al. 2015b). The latest LAT source catalog (4FGL; Abdollahi et al. 2020) suggests that the existing data cannot determine a clear classification of this source. Very Long Baseline Array (VLBA) observations reveal that the radio structure of its jet is similar to misaligned AGN (D'Ammando et al. 2015b). Recently, Angioni et al. (2019) confirmed that the jet of PKS 0521-36 is not highly beamed.
A study on the long-term variability behavior of PKS 0521-36 in γ-ray energies is lack. The characteristics of γ-ray variability are helpful to understand the nature of this source. In this work, we perform a full time-domain analysis on the latest Fermi-LAT data of PKS 0521-36. The paper is organized as follows. The basic introduction of Fermi-LAT and results of searching for periodicity by different methods are presented in Section 2. Our discussion of the results is given in Section 3. Finally, a summary is presented in Section 4.

Fermi-LAT Likelihood Analysis
Fermi-LAT is a pair conversion γ-ray detector with wide field of view (∼ 2.4 sr), large effective area (> 8000 cm 2 at ∼1 GeV), and broad-energy covering from ∼ 20 MeV to > 300 GeV. Thanks to its all-sky monitoring capabilities, one can get long-coverage and uniform observations of γ-ray sources. A detailed description of LAT can be found in Atwood et al. (2009).
The data used here were collected during ∼12.6 years of Fermi-LAT operation from 2008 August 4 to 2021 March 29 [Modified Julian Day (MJD): 54682 − 59302] in the 0.1-500 GeV energy range. The analysis of data follows the standard criteria for the point-source analysis 1 . We employ the Science Tools package of version v11r05p3 available from Fermi Science Support Center 2 (FSSC), and use the Pass 8 LAT database in which we choose a region of interest (ROI) of 15 • centered at PKS 0521-36. Two background templates, Galactic and extragalactic diffuse emissions, are modeled with the file gll iem v07.fits and iso P8R3 SOURCE V2 v1.txt, respectively. The instrumental response function "P8R3 SOURCE V2" is adopted in data processing. In order to get events with high probability of being photons, we only keep the SOURCE class events through setting "evclass = 128, evtype = 3" in the Fermi tool gtselect. In addition, the events with zenith angle > 90 • are filtered out to avoid contaminating from Earth's limb. Good time intervals (GTIs) are selected by using the criterion of "(DATAQUAL > 0)&&(LATCONFIG == 1)" in the Fermi tool gtmktime. Above settings ensure that we adopt high-quality observations.

γ-ray SED and Light Curves
In the process of analyzing the entire data, an initial binned maximum likelihood analysis (Mattox et al. 1996) is performed. The spectral parameters and the normalization factors of all sources lying within 5 • of PKS 0521-36 as well as the normalization factors of the sources located beyond 5 • but within 10 • from our target are set to be free, while the sources from ROI but beyond 10 • have their parameters fixed to the values in Fermi-LAT fourth source catalog (Abdollahi et al. 2020). In the iterative fittings, the best-fitting results are obtained, and are saved to be a new model file which is the basis of constructing light curves and spectrum energy distribution (SED). The spectral model of LogParabola (LP) form is used here, which is . (1) The global fitting results in an integrated average flux of (1.06 ± 0.02) × 10 −7 photons cm −2 s −1 and a total Test Statistics (TS) value of 18447.2. The TS is defined as 2log(L/L 0 ), where L is the maximum likelihood of the model with a point source at the target position, and L 0 is the maximum likelihood without the source (Mattox et al. 1996). The best-fitting spectral parameters are the spectral slope α = 2.43±0.02 and the curvature parameter β = 0.05±0.01. The average SED of PKS 0521-36 in the whole time period is shown in Figure 1. The unbinned maximum likelihood optimization is used to produce γ-ray light curve of PKS 0521-36. In each time bin, the new model file mentioned above is employed, but the spectral parameters of all sources are freezed. The time bins having T S <25 and the number of the model-predicted photon N pred <10 are excluded . This ensures that the robust detections of the source are adopted in the analysis (Kapanadze et al. 2020;Raiteri et al. 2013). A 30-day binning light curve is produced. The light curve shown in Figure    X-ray bands. The second flaring occurred in 2012 October has been studied by Angioni et al. (2019), and the latest flaring activity occurred in 2019 May.

Flux Distribution
One important property of a time series is its probability density function (PDF). One can explore the mechanism that drives variability by analyzing the flux distribution of a light curve (e.g., Rieger 2019). PDF can be estimated by fitting a histogram of observed data. The distribution of the γ-ray photon fluxes of PKS 0521-36 is presented in a   histogram with 25 bins, and is fitted by log-normal and Gaussian functions in the form of respectively. The parameter σ is the standard deviation and the parameter µ is the mean of the distribution. The flux histogram and best-fitting results are shown in Figure 3. The errors of data points are estimated by the method given by Gehrels (1986). The best-fitting parameters as well as the reduced χ 2 obtained from the fittings with Gaussian and log-normal models are listed in table 1. It is found that the log-normal model is preferred over the Gaussian model.

Searching for Periodicity
The Lomb-Scargle periodogram (LSP;Lomb 1976;Scargle 1982)   astronomy. The LSP and WWZ methods are then applied to the light curve of PKS 0521-36 to search for periodicity. A clear periodic modulation appears from MJD 56317 to 58447 (Figure 4). A significant peak at ∼400 days appears in both LSP and WWZ powers, which indicates a quasi-periodic oscillation (QPO). The period is estimated by fitting the power peak with a Gaussian function. It is 411±33 days given by LSP and 417±39 days given by WWZ. The bootstrap approximation method is used to calculate the false alarm probability (FAP), which is a robust way to estimate FAP   The algorithm given by Emmanoulopoulos et al. (2013) is also used to evaluate the significance of this QPO signal. The properties of observed variation, power spectral density (PSD) and PDF, are well applied in this method. The PSD model of the form P (f ) = Af −α is used, where α is the power-law index. The PDF is in the form of a log-normal distribution. The 4σ and 5σ significance curves are obtained by simulating the 30-day binning light curve 10 5 times, and are shown in Figure 5. The significance of the QPO signal is ≈ 5σ.
Further evidence for this QPO is presented in the folded light curve constructed by performing the likelihood analysis method with a 411-day period. Phase zero is set at MJD 56317. The variation with the phase is clear (see Figure 6). The fitting to the folded light curve in one phase with a constant results in χ 2 red =449.26/14, suggesting a significant variability in the folded light curve.
An additional method, REDFIT (Schulz & Mudelsee 2002), is also employed for the periodicity detection. This method estimates the red-noise spectrum by fitting the data with a first-order autoregressive (AR1) process. It can precisely evaluate the significance of the peaks in the PSD against the red-noise background. We use the program REDFIT3.8e 3 . The bias of Fourier transform for unevenly spaced data is also removed in the program. Two conditions should be satisfied when the program is used for data: (1) the data points of light curve should not be too clustered; and (2) the noise background of light curve can be approximated by an AR1 process. Obviously, our data meets the first condition. Moreover, the result of the non-parametric runs test for our data is rtest = 15 (inside the 98% acceptance region [9,22]), which indicates that the noise background spectrum is consistent with an AR1 model (Schulz & Mudelsee 2002). Therefore, REDFIT can be used to detect periodicity in our data. We set n 50 =2 and select a hanning window to reduce spectral leakage (Schulz & Stattegger 1997). The results show that there is a peak with more than 99% significance level, and the corresponding periodicity is 440.9±78.6 days (Figure 7). This period is consistent with the LSP and the WWZ results. Note that the REDFIT method only provides a maximum of significance of 99%.

Gaussian Process Modeling the γ-ray Light Curve
Besides LSP, WWZ and REDFIT, an alternative approach for the analysis of astronomical variability is the Gaussian process modeling light curves, which treats the variability in the time rather than frequency domain (e.g., Kelly et al. 2009Kelly et al. , 2014Li & Wang 2018). The PSD can be constructed by using modeling results. One popular Gaussian process  model is the continuous time autoregressive moving average (CARMA) model developed by Kelly et al. (2014), which has been applied to γ-ray variability of blazars (e.g., Goyal et al. 2018;Ryan et al. 2019;Yang et al. 2021). Here a newly developed Gaussian process model celerite (Foreman-Mackey et al. 2017), is used to investigate the γ-ray variability of PKS 0521-36.
In this work, a stochastically driven damped simple harmonic oscillator (SHO) is used, and its differential equation is given by where ω 0 , Q, (t) are the frequency of the undamped oscillator, the quality factor of the oscillator, and a stochastic driving force (white noise here), respectively. The PSD in the case is where S 0 is proportional to the power when ω = ω 0 . A specific model for a time series can be built up with a mixture of determined number of SHO terms. The Markov Chain Monte Carlo (MCMC) algorithm (Foreman-Mackey et al. 2013) is used to perform fitting to light curve. From MCMC sampling, we can evaluate values and uncertainties of model parameters. In detail, 5 × 10 4 samples are generated in our analysis, and the first 2 × 10 4 samples are taken as burn-in sampling. Figure 8 shows the fitting results to the 30-day binning light curve of PKS 0521-36 from MJD 56317 to 58447 with the celerite model (SHO×2). For the modeling results, we use Kolmogorov-Smirnov (KS) test to check the consistency between the standardized residuals and a normal distribution with mean zero and standard deviation of unity. The p-values of the test is p=0.23, which means this consistency is credible (Yang et al. 2021). The auto-correlation functions (ACFs) of the residuals and squared residuals in Figure 8 are almost inside the 95% confidence limits, so that the correlation structures have been captured by the model. These results indicate that the γ-ray variability can be described well by the two-SHO model. The posterior probability densities for the model parameters are shown in Figure 9 and the value of these parameters are given in Table 2. A QPO feature with period of ∼ 400 days appears in the PSD constructed from the modeling results ( Figure 10).
There is an aperiodic component and a periodic component in the γ-ray variability during the time range of MJD 56317-58447. The period of the QPO is estimated as 406 ± 10 days which is consistent with the results given by LSP, WWZ and REDFIT. The corresponding quality factor of the oscillator is large (  Flux (ph cm 2 s 1 ) Figure 11. One-day binning light curve (blue points) and 12-hour binning light curve (red points) for the bright flare during 2019 May 9-June 8 (MJD 58612-58642).

Short Timescale Variability
Intraday γ-ray variability of PKS 0521-36 has been found in the LAT data (D'Ammando et al. 2015b;Angioni et al. 2019). A recent bright γ-ray flare occurred in 2019 May 9-June 8 (MJD 58612-58642) has not been investigated. In our work, intraday variability during this flare is observed in the one-day and half-day binning light curves ( Figure 11). Here we use unbinned maximum likelihood method to construct the light curves, and select time bins having TS 4. This flare is not bright enough to use the selection criteria of TS 25 for such a short time bin.
In order to determine variability timescale, the 12-hour binning light curve is fitted by using the following function (Abdo et al. 2010) where F c , F 0 , t 0 , T r and T d are a constant flux, peak flux, peak time, and rise and decay times, respectively. This function can be used to describe one-peak flare. The light curve with multi-peaks is considered as a whole, and is fitted by using the model comprising several components of Equation 6. The best-fitting results are presented in Figure 12, and the parameters are given in Table 3. The residuals are fitted by the constant of zero successfully with χ 2 red = 0.58, suggesting that our fitting to the light curve is acceptable. It is found T r ∼ T d ∼ 7 hours. This is the third short timescale γ-ray flare in PKS 0521-36. However, the loose data selection criteria of TS 4 decreases the significance of the derived timescale.
The SEDs during 2019 May 9 -June 8 (MJD 58612-58642) and 2012 October 12 -November 11 (MJD 56212-56242) are also produced. The spectral parameters are listed in Table 4, together with the parameters for the average SED  Note-(1) time period, (2) peak number, (3)-(7) parameters in Equations 6 , (8) the reduced χ 2 , (3)(4) in unit of 10 −7 photons cm −2 s −1 . The value with * represents that it is fixed in fitting. shown in Figure 1. It is found that during the two flares, the SEDs are the same, which can be described by a power-law form, different from the average LP SED.

DISCUSSION
Our results show that the PDF obtained from the γ-ray light curve of PKS 0521-36 from 2008 August 4 to 2021 March 29 is preferentially log-normally distributed. The log-normal distribution of flux is found in X-ray binaries (e.g., Uttley & McHardy 2001), X-ray NLSy1s (e.g., Uttley et al. 2005), and blazars (e.g., Sinha et al. 2016;Shah et al. 2018;Ait Benkhali et al. 2020;Bhatta & Dhital 2020). Similar as X-ray binaries, the log-normal distribution of γ-ray fluxes may indicate an influence of the accretion disk on the jet (Giebels & Degrange 2009), i.e., the multiplicative behavior in an accretion disc impacting the jet. In this scenario, the variability timescale is dominated by the viscous timescale of the accretion disc (Rieger 2019). However, the tentative timescale of ∼ 7 hours γ-ray variability found in PKS 0521-036 is unlikely achieved in this scenario (Czerny 2006). The cascade-related emission processes may also lead to a log-normal distribution of fluxes (Rieger 2019), for instance, the cascade emission in lepto-hadronic jet model (e.g., Mannheim 1993;Böttcher et al. 2013;Yan & Zhang 2015). In this case, the observed GeV emission is produced by the cascade processes, which is feasible. However, it would require an extremely high jet power (e.g., Petropoulou et al. 2015;Cerruti 2020).
A simple and natural scenario is that the log-normal distribution of variability results from the emitting particles with the log-normal distribution Rieger 2019). This interpretation is associated with the fluctuations in the particle acceleration rate . The average LP SED suggests an average LP electron distribution. Both the log-normal and LP electron distributions could be the result of stochastic acceleration (Tramacere et al. 2011).
Here, we report a γ-ray QPO with a significance of ∼ 5σ in non-blazar AGN PKS 0521-36. The period is ∼1.1 years. This QPO behavior only appears in the temporal range from MJD 56317 to 58447 which is between two outbursts (occurred in 2012 October and 2019 May respectively). This significant QPO is obtained by three different periodicity detection methods (LSP, WWZ and REDFIT). Evidence for this QPO is also found in the Gaussian process modeling the 30-day binning γ-ray light curve. It is confirmed that the jet of PKS 0521-36 is a mildly beamed (Leon et al. 2016;Angioni et al. 2019). Therefore, the QPO of PKS 0521-36 is the first γ-ray QPO found in a mildly beamed jet. The average SED during MJD 56317-58447 is described by a LP model (see Table 4) which is different from the power-law SEDs in the two outbursts.
The detection of a helicoidal motion in the optical jet of PKS 0521-36 was reported (Jiménez-Andrade et al. 2017), which may be associated with helical magnetic field or jet precession. A simultaneous VLBI core flux activity during the γ-ray flare in 2010-2011 was reported by Angioni et al. (2019). This suggests that the γ-ray emission region is located inside the radio core rather in the jet. γ-ray variability on timescale of ∼7 hours is observed before and after the QPO. If the γ-ray photons produced in a single region, the helicoidal motion along the jet at pc scales (Jiménez-Andrade et al. 2017) could not lead to the γ-ray QPO of PKS 0521-36.
It seems that the outburst during 2012 October plays a key role in producing this QPO. It is not clear whether the outburst during 2019 May destroys this QPO or triggers another QPO. The data collected in future several years will answer this question.

SUMMARY
We have carried out a temporal analysis of the γ-ray emissions from an interesting non-blazar AGN PKS 0521-36 with the Fermi-LAT data from 2008 August 4 to 2021 March 29. Our main results are as follows.
(i) The PDF of the γ-ray fluxes is more consistent with a log-normal distribution rather than a Gaussian distribution. It may be due to the log-normal distribution of emitting electrons, which could be the result of stochastic acceleration. Moreover, the total average LP SED also can be produced by the stochastic acceleration.
(ii) Variability on timescale of ∼ 7 hours is found during the outburst in 2019 May. It is the third short timescale γ-ray flare in PKS 0521-36. Such a variability indicates that the γ-ray emission region is compact. However, it should be noted that the timescale may be in low significance.
(iii) A QPO with the period of ∼1.1 years is found at ∼ 5σ confidence level in the time interval (MJD 56317-58447) between two outbursts. Four different methods are used to examine the robustness of the QPO. All the methods give significant evidence for this QPO. This is the first γ-ray QPO found in a mildly beamed jet. The SED during the QPO interval is different from that during the two γ-ray outbursts in 2012 October and 2019 May. It is argued that the two γ-ray outbursts play a key role in the formation of this γ-ray QPO. The γ-ray variability during the QPO interval can be well modeled by two SHO components, an aperiodic component and a periodic component.