A 31.3 day Transient Quasiperiodic Oscillation in Gamma-ray Emission from Blazar S5 0716+714

Active galactic nuclei (AGNs) are energetic astrophysical sources that are powered by the accretion of gas onto the supermassive black holes (SMBH) at their centers. They show unique radiation characteristics, and multiband radiation from the radio to high-energy gamma rays has been observed (Padovani et al. 2017). Blazars are a subclass of AGNs, whose relativistic jets orientate along the observer's line of sight. According to the intensity of the emission line in the spectrum, blazars is further divided into BL Lacertae (BL Lac) objects and flat-spectrum radio quasar (FSRQ). The spectra of BL Lac objects contain very weak and narrow emission lines, while FSRQs show wide and strong emission lines (Urry & Padovani 1995; Ghisellini et al. 2011).

Based on the Fermi Gamma-ray Space Telescope long-term survey data, a few AGN quasiperiodic oscillations (QPOs) were reported in the gamma-ray band, ranging from days to months and even years. For the year-like QPOs, the physical mechanism was mostly attributed to pulsating accretion flow instability, which produces precession or helical jets, a homogeneous curved helical jet scenario, and Lense–Thirring precession of the flow; Keplerian binary orbital motion would induce periodic accretion perturbations (Ackermann et al. 2015). Several cases were reported, such as PG 1553+113, PKS 0426-380, PKS 0301-243, PKS 2155-304, PKS 0537-441, PMN J0948+0022, OJ 287, PKS 0601-70, PKS 0521-36, 4FGL J0112.1+2245 (Ackermann et al. 2015; Sandrinelli et al. 2016; Zhang et al. 2017a, 2017c, 2017d, 2017e, 2020, 2021; Kushwaha et al. 2020; Gong et al. 2022). For the month-like QPOs, the physical mechanism is mainly attributed to the helical structure in the jet (Zhou et al. 2018). For example, PKS 2247–131, B2 1520+31, PKS 1510–089, and 3C 454.3 were reported to have 34.5, 71, 92, and 47 day QPOs, respectively (Zhou et al. 2018; Gupta et al. 2019; Sarkar et al. 2021; Roy et al. 2022). For shorter (day-like) QPOs, the physical mechanism is flux enhancements by magnetic reconnection events in the magnetic islands inside the jet. This kind of QPO was reported in two sources, i.e., CTA 102 (∼7.6 days) and PKS 1510–089 (∼3.6 days; Sarkar et al. 2020; Roy et al.2022).

The QPO signals of AGN are also found in optical, X-rays, and radio bands. The most well-known case is BL Lac OJ 287, which exhibits ∼12 yr of QPO in its optical light curve based on more than a century of monitoring (Kidger et al. 1992; Valtonen et al. 2006). In addition, AO 0235+16, AO 0235+164, S5 0716+714, and PKS 2155-304 are reported to have QPO signals in the optical band. It was found that these QPO signals are distributed on the timescales of the hours, days, and years (Raiteri et al. 2001; Liu et al. 2006; Gupta et al. 2009; Lachowicz et al. 2009; Zhang et al. 2014). RE J1034+396, 2XMM J123103.2+110648, 1H 0707-495, and Mrk 766 are reported to have QPO signals in X-rays. Unlike those found in gamma-ray, these X-ray QPO signals oscillate on timescales of hours or days, and no reliable QPO signals have been reported on long timescales (year-like oscillation) (Stella & Vietri 1998; Gierliński et al. 2008; Lin et al. 2013; Pan et al. 2016; Zhang et al. 2017b, 2018). FSRQ J1359+4011, PKS J0805-0111 and BL Lac PKS 0219-164, PKS J2134-0153 are reported to have QPO signals in the radio band, most of which are on the timescale of years (King et al. 2013; Bhatta 2017; Ren et al. 2021a, 2021b).

S5 0716+714 is a very famous BL Lac object, which has been monitored by many telescopes for a long time, and has a redshift of z = 0.31 (Donato et al. 2001). S5 0716+714 is found to have intraday periodic oscillation, and the mass of the central black hole is estimated to be > 2.5 × 10⁶ M_⊙. No correlation between the flux level and the timescale of intraday variability (IDV) was found, so it is speculated that there is more than one radiation mechanism in this source (Gupta et al. 2009). There is a significant correlation between the change in the gamma-ray flux and the change in the position angle in the jet measured by very long baseline interferometry (Rani et al. 2014). From 2015 January 19 to February 22, Swift X-Ray Telescope (XRT) monitored the continuous burst status of S5 0716+714. The observed X-ray spectra conform to the breaking-power-law model. With the increase in the flux, the braking energy is transferred to higher energies, which indicates that synchrotron radiation is dominant in the observed X-ray spectra (Wierzcholska & Siejkowski 2015). NuStar observed the hard X-ray of S5 0716+714 and found a clear energy band fracture (8 keV; Wierzcholska & Siejkowski 2016). S5 0716+714 has obvious optical IDV, and multiple intraday QPO evidence has been found in multiple monitorings in different time epochs (Liu et al. 2017, 2019, 2021; Yuan et al. 2017; Hong et al. 2018; Kaur et al. 2018). In the burst stage of S5 0716+714, the gamma-ray spectrum shows an obvious broken feature between 0.93 and 6.90 GeV, and five photons with energy E > 0.1 TeV are found (Geng et al. 2020). The TeV emission from S5 0716+714 originated in a superluminal knot, which means that the high-energy radiation may come from the moving emission region in the helical path upstream of the jet (MAGIC Collaboration et al. 2018). The space telescope (Transiting Exoplanet Survey Satellite) in 2019 December–2020 January observed and analyzed the fast variability of the object with unprecedented high-resolution sampling, revealing the characteristics of IDV (Raiteri et al. 2021). Recent studies have shown that there is a significantly positive or negative correlation between S5 0716+714 radio bands (15, 37, and 230 GHz) and gamma-ray (0.1–200 GeV) fluxes in different time ranges (Kim et al.2022).

Among 1525 sources with a variation index greater than 21.67 given in the Fermi Large Area Telescope Light Curve Repository (LCR ⁵ ), we search for the QPO on a timescale of months. The results show that S5 0716+714 is another blazar with months QPO, in addition to the previously reported PKS 2247–131, B2 1520+31, PKS 1510–089, and 3C 454.3. S5 0716+714 has a ∼31.3 day QPO with a 4.1σ confidence level, based on a Large Area Telescope (Fermi–LAT) light curve [Modified Julian Day (MJD): 55918−56137] with 5 day bins. Fermi–LAT data reduction and results of searching for periodicity by different methods are presented in Section 2. The conclusion and discussion are given in Section 3.

2.1. Fermi–LAT Data

Fermi–LAT, ⁶ the primary instrument of the Fermi satellite mission, is an imaging, wide -field-of-view, high-energy gamma-ray telescope, covering the energy range from below 20 MeV to more than 300 GeV. LAT has a large angular field of view of about 2.3 sr, covers the entire sky every three hours, and has a large effective area of >8000 cm² at ∼1 GeV (Atwood et al. 2009).

For data selection, we selected LAT events in the mission elapsed time (MET) 239557417.0–647393497.0 (∼12.9 yr) from the Fermi PASS 8 database and used the Science Tools package of version v11r05p3 provided by Fermi Science Support Center2 (FSSC ⁷ ). First, we chose the events belonging to the SOURCE class (evclass = 128, evtype = 3) within the energy range of 0.1–300 GeV from a circular region of interest (ROI) having a radius of 15° centered at the source 4FGL J0721.9+7120 (S5 0716+714). We excluded events with zenith angles over 90° to improve the point-spread function and reduce diffuse emissions, and use the standard filter (DATA_QUAL>0) && (LAT_CONFIG == 1) to select a good time interval (GTI) to obtain high-quality data. Afterward, we create the source model XML file containing the spectral shapes of all sources within the radius of the ROI+10° around the source location, including two background templates for Galactic and extragalactic diffuse emissions, modeled using files gll_iem_v07. fit and iso_P8R3_SOURCE_V3_v1. txt, respectively. Note that the normalization of the two diffuse emission components is set as a free parameter in the analysis. Finally, we performed a maximum-likelihood analysis of the input XML spectrum file using the gtlike tool and the instrument response function (IRF) P8R3_SOURCE_V3 to obtain the spectrum of the source. In the 4FGL catalog, the final source spectrum is modeled using a logarithmic parabola (Abdollahi et al. 2020) as follows:

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}=K{\left(\displaystyle \frac{E}{{E}_{b}}\right)}^{-\alpha -\beta \mathrm{log}\left(E/{E}_{b}\right)},\end{eqnarray} \tag{ 1 }$$

where α is the spectral index at the break energy (E_b ); we kept E_b fixed during likelihood fitting. The optimizer finds the most suitable spectral parameters but not the location. In other words, the fit tool does not fit the source coordinates. Therefore, test statistic (TS) is added here for each location to calculate the saliency of the point source at that location, TS = −2 $(\mathrm{ln}L0-\mathrm{ln}L)$, where L0 represents the maximum-likelihood value fitted by the background model, and L represents the maximum-likelihood value fitted by adding a test point source to the background model.

Based on the data analysis, we generated flux data for 5 days, excluding very few data points with TS less than 25. The light curve of the flaring state (MJD 55632–56312) is shown in Figure 1 (top panel), and the gray shaded part (MJD 55918–56137) is shown in Figure 1 (bottom panel), where the blue dotted line is the average value of the flux, and the red arrows indicate the peak positions for each cycle.

Zoom In Zoom Out Reset image size

2.2. QPO Analysis and Results

We found a sine-like periodic variability by visual inspection, so we fitted it with a simple sine function, $y={y}_{0}+A\sin (2\pi (t-{t}_{c})/T)$, where T is the period of sine fitting. The fitting results show that the flux changes tend to be sinusoidal (Figure 1, bottom panel, green dotted line). To study the periodicity of the gamma-ray light curve of S5 0716+714, we used three methods to analyze the QPO signal. The first method is the weighted wavelet Z-transform (WWZ; Foster 1996). The WWZ method can decompose the data into the time domain and frequency domain, and convolute the light curve with the kernel related to time and frequency. We used the abbreviated Morlet kernel (Sarkar et al. 2021), and its functional form is as follows:

$$\begin{eqnarray}&&f[\omega (t-\tau )]=\exp [i\omega (t-\tau )-c{\omega }^{2}{\left(t-\tau \right)}^{2}].\end{eqnarray} \tag{ 2 }$$

The WWZ map is then given by,

$$\begin{eqnarray}&&W[\omega ,\tau :x(t)]={\omega }^{\tfrac{1}{2}}\int x(t)f* [\omega (t-\tau )]{dt}.\end{eqnarray} \tag{ 3 }$$

In formula (2) and (3), f^* is the complex conjugate of the Morlet kernel f, ω is the frequency, and τ is the time-shift. The advantage of the WWZ map is that it can search the periodicity in the data by decomposing the signal into the frequency time space at the same time, and detect any main periodicity and its duration time span. In other words, this technique measures the nonstationarity of the data set and indicates the time range of any possible periodic characteristics. For noncontinuous periodic signals (transient periodic signals), WWZ can show a decrease when the periodic signal starts to disappear. The advantage of the WWZ method is that it can detect any major periodicity and its duration span. We obtained the color scale of the WWZ power spectrum of S5 0716+714 in the time interval, MJD 55593−56312 (the flare state), as shown in the left panel of Figure 2. For the light curve during MJD 55918−56137, the WWZ method provides clearer results (see the middle panel of Figure 2). The WWZ power and its time average power strongly indicate the existence of a QPO for about 31.3 ± 1.8 days. The uncertainty is estimated as the half width at half maximum of the Gaussian function fitting the power peak. After that, we also used the Lomb–Scargle periodogram (LSP) method, which is more common in astronomy, to detect the potential quasiperiodic signals in S5 0716+714. The LSP method gives the power of the flux modulation at different frequencies (Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009). For the uniformly sampled light curve, the square of the discrete Fourier transform modulus gives the periodogram. However, for irregular sampling, the LSP method iteratively fits sinusoidal curves with different frequencies to light curves and constructs the periodogram according to the goodness of fit, and can provide accurate frequency and power spectrum intensities. Figure 3 shows the power spectrum. It is obvious that there is a peak around 31.2 ± 1.6 days. This also confirms the signal detection results of the WWZ method. Here, we use false-alarm probability (FAP) to evaluate the confidence level of the peak. $\mathrm{FAP}{({P}_{n})=1-(1-\mathrm{prob}(P\gt {P}_{n}))}^{M}$, where FAP denotes the probability that at least one of the M independent power values in a given frequency band of the white noise periodogram is greater than or equal to the power threshold P_n (Horne & Baliunas 1986; VanderPlas 2018). The green dotted line in Figure 3 shows the FAP = 0.01%, and the confidence level of the peak (31.2 ± 1.6 days) is higher than 4σ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The third method, REDFIT (Schulz & Mudelsee 2002), was also used to search for QPO signals in the light curve. For time-series data in the astronomy field, it is difficult to accurately estimate the red-noise spectrum because the sampling time interval is usually not uniform. REDFIT was proposed to overcome this problem by fitting a first-order autoregressive (AR1) process directly to unevenly spaced time series, avoiding interpolation in the time domain and its inevitable bias. This procedure was also used to test the significance of the time-series flux peaks against the background of red noise in the AR1 process. The emission fluxes of AGN are usually autoregressive (Schulz & Mudelsee 2002; Kushwaha et al. 2020), the current emission flux depends on its immediate previous emission flux, so we can use the AR1 process to model the emission red-noise spectrum of S5 0716+714. The program REDFIT3.8e ⁸ can estimate the spectrum using LSP and Welch overlapped segment averaging (WOSA), with a number of WOSA segments (n50 = 2); selecting a Welch window reduces spectral leakage. As shown in Figure 4, there is a 31.1 ± 1.9 day peak with a confidence level of >99% (it is worth noting that the REDFIT3.8e code provides a maximum significance of 99%). This result is consistent with that of the WWZ and LSP methods. However, the above three methods we use are all similar analysis and not completely independent methods to test the QPO, so we also need to consider the other models to test the physical QPO in this work.

Zoom In Zoom Out Reset image size

Here, we model the light curve using the autoregressive integrated moving average (ARIMA) model, which has been widely used in disciplines or fields such as time-series analysis, signal processing, and econometrics. This is a flexible method consisting of three parts, in which the autoregressive (AR) process can quantify the coefficient of the dependence of the current value on the most recent past value, the integrated (I) process is used to reduce the trend, and the moving average (MA) process can quantify the coefficient of the current value's dependence on the system's recent random shocks (Scargle1981; Vecchia 1985; Feigelson et al. 2018). It is defined in statistics as

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{d}F({t}_{i})=\sum _{j=1}^{p}{\theta }_{j}{{\rm{\Delta }}}^{d}F({t}_{i-j})+\sum _{j=1}^{q}{\phi }_{j}\epsilon ({t}_{i-j})+\epsilon ({t}_{i}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\mathrm{or},\left(1-\sum _{j=1}^{p}{\theta }_{j}{L}^{j}\right){{\rm{\Delta }}}^{d}F({t}_{i})=\left(1-\sum _{j=1}^{p}{\phi }_{j}{L}^{j}\right)\epsilon ({t}_{i}),\end{eqnarray} \tag{ 5 }$$

where p, d, and q are the AR order, difference operator, and MA order, respectively. The seasonal autoregressive integrated moving average model or SARIMA (p, d, q) × (P, D, Q)_s is a model based on ARIMA (p, d, q) evolved methods of identifying physical periodicity in light curves (Adhikari & Agrawal 2013; Sarkar et al. 2020), defined in statistics as

$$\begin{eqnarray}&&\begin{array}{l}\left(1-\sum _{j=1}^{p}{\theta }_{j}{L}^{j}\right)\left(1-\sum _{j=1}^{p}{{\rm{\Theta }}}_{j}{L}^{{sj}}\right){{\rm{\Delta }}}^{d}{{\rm{\Delta }}}_{s}^{D}F({t}_{i})\\ \quad =\,\left(1-\sum _{j=1}^{p}{\phi }_{j}{L}^{j}\right)\left(1-\sum _{j=1}^{p}{{\rm{\Phi }}}_{j}{L}^{{sj}}\right)\epsilon ({t}_{i})+A(t),\end{array}\end{eqnarray} \tag{ 6 }$$

where P, D, Q, and s are the seasonal AR order, difference operator, MA order, and period parameter, respectively. Note that the detailed derivation process for Equations (4)–(6) is shown in the Appendix at the end of this paper. Next, we fit the light curves using the ARIMA and SARIMA models and compare their goodness of fit using the Akaike information criterion (AIC; Akaike 1974). AIC = −2 $\mathrm{ln}L+2k$, where L is the likelihood function, and k is the number of free parameters in the model. AIC rewards models that fit the data better, while penalizing models that use more parameters. Therefore, we build the parameter space to search for the model with the smallest AIC

$$\begin{eqnarray*}&&\varphi =\left\{\begin{array}{l}p,q\in \ [0,8]\\ d,D\in \ [0,1]\\ P,Q\in \ [0,5]\\ s\in \ [0,16]\times 5\,\mathrm{days}.\end{array}\right.\end{eqnarray*}$$

The fitting results are shown in Figure 5. When using the seasonal model SARIMA(4, 0, 1) × (6, 1, 2)_{30 days}, the AIC has the global smallest value of 28.0 (Figure 5(b)), while the AIC value of the nonperiodic best-fit model ARIMA(4, 1, 3) is 489.55 (Figure 5(a)), which also means that the periodic model can better fit the S5 0716+714 light curve. Figure 5(c) presents the AIC values considering different periods, and we find that the best AIC occurs at the 30 ± 2.5 day position. The amplitude of the uncertainty is estimated to be half of the light-curve time bin (2.5 days).

Zoom In Zoom Out Reset image size

2.3. Significance Estimation

As described above, we used several different methods and techniques to analyze possible QPO signals in the light curve of blazar S5 0716+714 and obtained consistent results. However, the statistical characteristics of a large number of light curves in AGN show red-noise-like behavior related to the frequency (for example, the red noise is produced by stochastic flares), which may produce false QPOs (Vaughan et al. 2016). In addition, the sampling instability of the light curve (observation time, weather, seasonal variation, diurnal variation, etc.) will also produce pseudo-QPO signals. Therefore, the estimation of the confidence level is very important in QPO search. We modeled the red noise with a first-order autoregressive process (e.g., REDFIT). The first-order autoregressive process is a pure power-law model. The pure power laws or smooth bending power laws could reasonably model the red-noise power spectral density (PSD), and the latter is closer to the red-noise PSD. Therefore, we further adopt a smooth bending power law to model the PSD. The Monte Carlo (MC) procedure includes three steps: model the original light-curve PSD and probability distribution function (PDF), finally based on the parameters of the PSD and PDF models of the simulated light curves in the LSP and WWZ analyses (Emmanoulopoulos et al. 2013). To estimate the underlying PSD, we model the light-curve red-noise PSD by smooth bending power laws plus a Poisson noise (González-Martín & Vaughan 2012; Max-Moerbeck et al. 2014). Poisson noise is considered here because the observed light curve is a product of the counting detector process, so the observations are subject to Poisson noise, which is imprinted in the corresponding PSD as a constant component. The model form of PSD is

$$\begin{eqnarray}&&p(f:\gamma ,c)=\displaystyle \frac{{{Af}}^{-{\alpha }_{\mathrm{low}}}}{1+{\left(f/{f}_{\mathrm{bend}}\right)}^{{\alpha }_{\mathrm{high}}-{\alpha }_{\mathrm{low}}}}+C,\end{eqnarray} \tag{ 7 }$$

where the model parameters A, α_low, α_high, f_bend, and C are the normalization, low-frequency slope, high-frequency slope, bend frequency, and Poisson noise, respectively. Here we obtain the optimal parameters of the PSD model through maximum-likelihood analysis; i.e., first we calculated the joint probability density function (likelihood function) of obtaining the ensemble of periodogram estimates for the given PSD model, and then we maximize the log-likelihood function probability to get the PSD model optimal parameters. For a more detailed analysis process see Emmanoulopoulos et al. (2013). The optimal likelihood parameter values of the PSD model are A (0.014), α_low (1.001), α_high (4.801), f_bend (0.056), and C (9.638 × 10⁻⁵), respectively. The light curves of AGNs generally exhibit a "burst"-like behavior, and the PDF is distributed according to a right-heavy-tailed distribution. This means that the probability of occurrence of high flux values is greater than that expected from a Gaussian distribution, so the PDF is modeled as a superposition of the two distributions of the two active states (low and flare state). This is modeled using a gamma distribution and a lognormal distribution (Emmanoulopoulos et al. 2013),

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{PDF}(x:\eta ) & = & {\omega }_{1}\displaystyle \frac{{x}^{{\alpha }_{s}-1}{e}^{-x}}{{\rm{\Gamma }}({\alpha }_{s})}+{\omega }_{2}\displaystyle \frac{1}{s(x-\mu )\sqrt{2\pi }}\\ & & \times \,\exp \left(-\displaystyle \frac{{\mathrm{log}}^{2}(\tfrac{x-\mu }{\sigma })}{2s}\right),\end{array}\end{eqnarray} \tag{ 8 }$$

where α_s is the shape parameter of gamma function fitting, and the obtained value is 4.993; $\mathrm{log}\mu$ and σ are the log mean and standard deviations of the lognormal distribution, respectively. The maximum-likelihood analysis obtains the parameter values as 2.886 and 0.401, respectively.

For the PSD- and PDF-based parameter values, we generated 10⁵ artificial light curves using the Python code "emmanoulopoulos." ⁹ (To account for the "red-noise leakage" effect, the generated artificial light curves are much longer (e.g., 100 times) than the original light-curve data set, where a subset with the desired length is randomly selected.) Finally, according to the mean and standard deviations of the distribution of the artificial light-curve PSD at each frequency, the confidence level of the QPO was estimated. The blue and red dashed lines in Figure 2 and the right panel of Figure 3 represent the confidence levels of 3σ (99.73%) and 4σ (99.994%), respectively. The confidence level of the QPO signal detected by the WWZ and LSP methods exceeds 4σ. Therefore, the 31.3 days QPO signal in S5 0716+714 is highly significant.

Through the above analysis, we get that the quasiperiodic signal in S5 0716+714 is significant. We use the ARIMA and SARIMA models to fit the light curve and find that the SARIMA model has an AIC value much better than the best ARIMA model (Figure 5(b)), and this fact proves that the SARIMA model with strict periodic modulation component can better model the light curve. Figure 5(c) shows the AIC values given by the SARIMA model in different seasons (s values). We give 99.994% confidence levels based on the mean and standard deviations of the AIC corresponding to the SARIMA model in different seasons. It can be seen from Figure 5(c) that the best-fitting model of the time-varying curve for s = 30 days has a significance of over 4σ. So we have reason to believe that this QPO signal originates from the physical periodic process. In the next chapter, we will analyze and discuss the physical mechanism of QPO.

In this work, we searched for month-like QPOs in the Fermi Large Area Telescope Light Curve Repository (LCR) catalog of 1525 sources (all of which have variability indexes >21.67). The search results showed a possible QPO in the TeV blazar S5 0716+714. We analyzed the Fermi–LAT gamma-ray (0.1–300 GeV) light curve of this source and found that a QPO of 31.3 ± 1.8 days of this TeV blazar in MJD 55918−56137 with a FAP of less than 0.01%. The probability that the detection due to a false alarm in a sample of 1525 sources is 1 − (1 − FAP)¹⁵²⁵ ∼ 10% which turns out to be roughly 10%. In astrophysical terms, this 10% false alarm is acceptable. The in-depth analysis found that the variations can be modeled as a ∼30 days strictly periodic component superposed on a stochastic autoregressive variability. Together these can be considered QPO behavior. The presence of both components is established independently using wavelet analysis (with the weighted wavelet Z-transform), Fourier analysis (with the Lomb–Scargle periodogram), and autoregressive analysis (comparing the ARIMA and SARIMA models). Compared with the gamma-ray QPO signal previously reported in AGNs, this QPO is the first month-like QPO with seven cycles. The QPO only occurred during the outbreak state in 2011 June (MJD 55918−56137), so this QPO was transient. Most of the QPOs that have been reported are transient signals (Zhang & Wang 2021). We used the simulation method (Emmanoulopoulos et al. 2013) to produce 10⁵ light curves to estimate the significance of the QPO and find the confidence level of ∼4.1σ. However, this kind of month-like transient QPO perhaps originate from the central activity of the blazar or from the influence of background photons. The month-like QPO timescale is very close to the moon's sidereal period of 27.32 days, so it may be modulated by the moon's gamma rays. However, the QPO of S5 0716+714 is transient, while the gamma-ray modulation of the moon is persistent, so the possibility of the influence of the moon can be ruled out. Second, the stochastic processes can produce pseudo-periodic variations, but these 'phantom periodicities' typically show less than three cycles. When the period reaches five cycles, it is relatively straightforward to distinguish between potential physical cycles and stochastic processes (Vaughan et al. 2016). While the QPO behavior in S5 0716+714 appears in MJD 55918−56137 for up to seven cycles, so this QPO signal is likely a physical periodic process. Additionally, the SARIMA model has an AIC value much better than the best ARIMA model, demonstrating that a strictly periodic component must be present in addition to any QPO-like behavior of the autoregressive component. Below, we will discuss based on the fact that the QPO originated from the center's activities.

The possible physical mechanism behind QPO in AGN has been widely discussed, and many models have been proposed. One possible explanation could be a binary SMBH system, in which the periodicity can be induced by the secondary BH by piercing the accretion disk of the primary BH during the orbit motion (Valtonen et al. 2008; Villforth et al. 2010; Ackermann et al. 2015; Sandrinelli et al. 2016). The second is the jet precession model, where a secondary black hole in a noncoplanar orbit generates periodic oscillations due to tidal-induced rapid precession in the inner region of the primary accretion disk (Romero et al. 2000; Liska et al. 2018). The third is the Lense–Thirring precession of accretion disks (Stella & Vietri 1998). The above three models are often used to explain the long-term quasi periods (persistent QPOs with at least year-long periods). However, S5 0716+714 is famous for its rapid and violent intraday variability (Gupta et al. 2009; Wierzcholska & Siejkowski 2015; Liu et al. 2017, 2019, 2021; Yuan et al. 2017; Hong et al. 2018), and the QPO found in this source shows a month-like oscillation, so these models cannot provide appropriate explanation of this short timescale QPO. The models need to be further explored.

As far as we know, the gamma-ray emission in the high-energy bands is attributed to a shock in the blazar jet. The variability of the gamma-ray flux of the blazar can be explained in the scenario of a shock wave propagating along a helical path in the jet (Marscher et al. 2008; Larionov et al. 2013; MAGIC Collaboration et al. 2018). One possible explanation for the 31.3 day QPO is that it comes from a region with enhanced emission, moving helically within the jet (Valtonen et al. 2008; Mohan & Mangalam 2015; Sarkar et al. 2021). When different parts of the helical jet have different angles with the line of sight, the change in the relativistic beam effect will lead to significant flux changes even if there is no internal change in the emission of the jet. So, the blob in the jet moves toward us, and the observation angle of helical motion essentially changes periodically, resulting in the QPO of the flux (Villata & Raiteri 1999). Due to the postulated helical motion of the blob, the viewing angle of the blob with respect to our line of sight (θ_obs) changes periodically with time as (Sobacchi et al. 2017; Roy et al. 2022),

$$\begin{eqnarray}&&\cos {\theta }_{\mathrm{obs}}=\sin \phi \sin \psi \cos (2\pi t/{P}_{\mathrm{obs}})+\cos \phi \cos \psi ,\end{eqnarray} \tag{ 9 }$$

where ϕ, ψ, and P_obs represent the pitch angle of the helical path, the angle of the jet axis concerning our line of sight, and the observed periodicity in the light curve, respectively. The Doppler factor (δ) varies with the viewing angle as $\delta =1/[1-\beta \cos {\theta }_{\mathrm{obs}}]$, where ${\rm{\Gamma }}=1/{(1-{\beta }^{2})}^{1/2}$ is the bulk Lorentz factor of the blob motion with β = ν_jet/c. The value of ϕ usually fluctuates between 0° and 5°; here we choose ϕ = 2°, ψ = 5°, and Γ = 8.5 according to the typical value of the blazar (Sobacchi et al. 2017; Zhou et al. 2018). Due to the relativistic effect and Doppler enhancement effect, the observed period P_obs is much smaller than the actual physical period P_rf. We can estimate the physical period P_rf by the following formula,

$$\begin{eqnarray}&&{P}_{\mathrm{rf}}=\displaystyle \frac{{P}_{\mathrm{obs}}}{1-\beta \,\cos \,\phi \,\cos \,\psi }.\end{eqnarray} \tag{ 10 }$$

Based on the known parameters, P_rf was estimated as 7.57 yr; we can also estimate the distance the blob moves during one period, $D=c\beta {P}_{\mathrm{rf}}\cos \phi \approx 2.09\,\mathrm{pc}$. The total projection distance of seven cycles can be estimated as ${D}_{p}\,=7D\sin \psi \approx 1.27\,\mathrm{pc}$. Parsec-scale helical jets have been detected in several blazars (Bahcall et al. 1995; Vicente et al. 1996; Tateyama et al. 1998). The helical structures of jets have been supported by optical polarization observations (Marscher et al. 2008), but their origin remains unresolved. We can assume that the jet structure is generated under the action of a helical magnetic field (Vlahakis & Königl 2004).

In a modified helical jet model, the blob moving helically in a curved jet (Sarkar et al. 2021; Roy et al. 2022). According to this model, the inclination angle ψ of the jet's axis to the line of sight is no longer constant, but a function of time ψ(t) (Sarkar et al. 2021). In our QPO timescale, the bending of the helical jet will not have a too large angle, so the impact on the value of QPO can be ignored, but the flux value of the QPO signal will slightly increase or decrease due to different inclination angles ψ. This provides a possible explanation for the obvious flux increase in the seventh cycle of the QPO signal.

If there are simultaneous multiband observations, it will provide important evidence for confirming the QPO. Unfortunately, S5 0716+714 did not have available data in lower-energy bands in MJD 55918−56137. The appearance of QPOs in the gamma-ray light curves of blazars is a rare event. It is even more difficult to verify it in multiple bands. So, in the future, long-term and high-temporal-resolution observations of blazars are very important to identify QPOs.

We thank the anonymous referee for useful and constructive comments. This study has made use of the Fermi–LAT data, obtained from the Fermi Science Support Center, provided by NASA's Goddard Space Flight Center (GSFC). The Fermi LAT Collaboration acknowledges the generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. The work is supported by the National Natural Science Foundation of China (grants 11863007, 12063007).

The Fermi–LAT data used in this article are available in the LAT data server at: https://fermi.gsfc.nasa.gov/ssc/data/access/.

The Fermi–LAT data analysis software is available at: https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/.

The Fermi Large Area Telescope Light Curve Repository (LCR): https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository/.

The weighted wavelet Z-transform(WWZ) method: https://github.com/eaydin/WWZ/.

The program REDFIT3.8e: https://www.marum.de/Prof.-Dr.-michael-schulz/Michael-Schulz-Software.html/.

The light-curve simulation method: https://github.com/samconnolly/DELightcurveSimulation/, https://github.com/lena-lin/emmanoulopoulos.

Software for evaluating SARIMA models: https://www.statsmodels.org/v0.13.0/examples/notebooks/generated/statespace_sarimax_stata.html.

As we all know, the light curve of a blazar is a time-series array, which we define as F(t). First, we construct an autoregressive (AR) model for the blazar's light curve,

$$\begin{eqnarray}&&F({t}_{i})=\sum _{j=1}^{p}{\theta }_{j}F({t}_{i-j})+\epsilon ({t}_{i}),\end{eqnarray} \tag{ A1 }$$

where t_i is the time, (t_i ) is the fluctuation or error of the model fitting, p is the time lag (the order of AR) of the past flux affecting the current flux, and θ_j is the coefficient of the AR model. Second, we build a moving average (MA) model for the light curve,

$$\begin{eqnarray}&&F({t}_{i})=\sum _{j=1}^{q}{\phi }_{j}\epsilon ({t}_{i-j})+\epsilon ({t}_{i}),\end{eqnarray} \tag{ A2 }$$

where q is the order of the MA, and ϕ_j is the coefficient of the MA. Next, to consider the nonstationary time series, we compute the d-order continuous difference (integrated) on the light curve,

$$\begin{eqnarray}&&{{\rm{\Delta }}}^{d}F({t}_{i})={\left(1-L\right)}^{d}F({t}_{i}),\end{eqnarray} \tag{ A3 }$$

where L is the lag operator, L^k F(t_i ) = F(t_i−k ), and Δ^d is the continuous difference operator, Δ^d = (1 − L)^d . Finally, AR, I, and MA are combined to form the autoregressive integrated moving average (ARIMA) model,

$$\begin{eqnarray}&&\left(1-\sum _{j=1}^{p}{\theta }_{j}{L}^{j}\right){{\rm{\Delta }}}^{d}F({t}_{i})=\left(1-\sum _{j=1}^{p}{\phi }_{j}{L}^{j}\right)\epsilon ({t}_{i}),\end{eqnarray} \tag{ A4 }$$

where p, q, and d are the orders of the AR, MA, and differencing, respectively, and the lag operator is defined as L^k (t_i ) = (t_i−k ). For light curves with strict periodic modulation, the SARIMA model is introduced with the periodic modulation seasonal (S) component, or SARIMA (p, d, q) ×(P, D, Q)_s ,

$$\begin{eqnarray}&&\begin{array}{l}\left(1-\sum _{j=1}^{p}{\theta }_{j}{L}^{j}\right)\left(1-\sum _{j=1}^{p}{{\rm{\Theta }}}_{j}{L}^{{sj}}\right){{\rm{\Delta }}}^{d}{{\rm{\Delta }}}_{s}^{D}F({t}_{i})\\ \quad =\,\left(1-\sum _{j=1}^{p}{\phi }_{j}{L}^{j}\right)\left(1-\sum _{j=1}^{p}{{\rm{\Phi }}}_{j}{L}^{{sj}}\right)\epsilon ({t}_{i})+A(t),\end{array}\end{eqnarray} \tag{ A5 }$$

where s is the component responsible for seasonality; P and Q are the orders of the seasonal AR and MA models, respectively; Θ_j and Φ_j are the coefficients of the seasonal AR and MA models, respectively; and ${{\rm{\Delta }}}_{s}^{D}$ is the seasonal difference operator, ${{\rm{\Delta }}}_{s}^{D}F{({t}_{i})=(1-{L}^{s})}^{D}F({t}_{i})$. For a more detailed derivation process of the ARIMA and SARIMA models, please refer to Scargle (1981), Adhikari & Agrawal (2013), and Feigelson et al. (2018).