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

We systematically search for quasiperiodic oscillatory (QPO) signals on the month timescale among the 1525 sources given in the Fermi Large Area Telescope Light Curve Repository. We find a transient QPO of 31.3$\pm$1.8 days in the gamma-ray band light curve of the TeV blazar S5 0716+714, which has seven cycles (MJD 55918-56137) for the first time by weighted wavelet Z-transform and Lomb-Scargle periodogram methods. Monte Carlo simulations based on the power spectral density and probability distribution function were used to evaluate the confidence level of the QPO, and the result is $\sim 4.1\sigma$. Seasonal autoregressive integrated moving average modeling of the light curve revealed it is a significant physical QPO. The physical models to explain the sporadic QPO of the month-timescale QPOs in blazar were discussed. Our studies indicate that the helical jet model and blob move helically in a curved jet model to properly explain this kind of transient QPO.


INTRODUCTION
Active galactic nuclei (AGN) are energetic astrophysical sources that are powered by the accretion of gas onto the super-massive black hole (SMBH) at its center. It shows unique radiation characteristics, and all band radiation from radio to high-energy gamma-ray has been observed (Padovani et al. 2017). Blazar is a subclass of AGN, whose relativistic jet orientates the observer's line of sight. According to the intensity of the emission line in the spectrum, a blazar is further divided into BL Lacertae (BL Lac) object and flat spectrum radio quasar (FSRQ). The spectrum of BL Lac contains 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 of AGN's QPOs were reported in the gamma-ray band, ranging from days to months and even years of timescales. For the year-like QPOs, the physical mechanism of them was mostly attributed to pulsating accretion flow instability produces precession or helical jets, a homogeneous curved helical jet scenario, 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;Zhang et al. 2017a,b,c;Kushwaha et al. 2020;Zhang et al. 2020Gong et al. 2022). For the month-like QPOs, the physical mechanism for the generation of these QPOs 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 days QPO respectively (Zhou et al. 2018;Gupta et al. 2019;Roy et al. 2022;Sarkar et al. 2021). For more short (the day-like) QPOs, the physical mechanism of them is flux enhancements by the magnetic reconnection events in the magnetic islands inside the jet. This kind of QPO was reported in two-source, i.e., CTA 102 (∼ 7.6 days) and PKS 1510-089 (∼ 3.6 days) 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 years 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 time scales 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 time scales of hours or days, and no reliable QPO signals have been reported on long time scales (year-like oscillation) (Gierliński et al. 2008;Lin et al. 2013;Pan et al. 2016;Zhang et al. , 2018Stella & Vietri 1998). 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 time scale of years (King et al. 2013;Bhatta 2017;Ren et al. 2021a,b). S5 0716+714 is a very famous BL Lac object, which has been monitored by many telescopes for a long time, and its redshift z = 0.31 (Donato et al. 2001). S5 0716+714 is found to have intra-day periodic oscillation, and the mass of the central black hole is estimated to be > 2.5 × 10 6 M , no correlation between the flux level and the time scale of intra-day variability (IDV) was found, so it is speculated that there is more than one radiation mechanism in this source ). There is a significant correlation between the change of gamma-ray flux and the change of position angle in the jet measured by very long baseline interferometry (Rani et al. 2014). From January 19 to February 22, 2015, swift/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 of flux, the braking energy is transferred to higher energy, 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 intra-day QPO evidence has been found in multiple monitoring in different time epochs (Liu et al. 2017;Yuan et al. 2017;Hong et al. 2018;Kaur et al. 2018;Liu et al. 2019Liu et al. , 2021. 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, TESS) 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 time scale 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 days 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 days bin. 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.

Fermi-LAT Data
The Fermi-LAT, the primary instrument of the Fermi satellite mission, is an imaging, wide field-of-view (FoV), 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 and covers the entire sky every three hours and a large effective area of > 8000 cm 2 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 years) time from the fermi PASS 8 database and used the Science Tools package of version v11r05p3 provided by Fermi Science Support Center2 (FSSC). Firstly, 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 pointspread 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.f it and iso P 8R3 SOU RCE V 3 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, using the instrument response function (IRF) P 8R3 SOU RCE V 3 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: 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, T S = −2(lnL0 − lnL), 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 T S 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.

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π(t − t c )/T ), where T is the period of the 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: (2) The WWZ map is then given by, In formula (2) and (3), f * is the complex conjugate of the Morlet kernel f , ω is the frequency, and τ is the timeshift. The advantage of WWZ map is that it can search the periodicity in the data by decomposing the signal into 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 non continuous 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 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). WWZ power and its time average power strongly indicate the existence of QPO for about 31.3±1.8 days. The uncertainty is estimated as the half-width at half the maximum value 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 quasi-periodic signals in S5 0716+714. LSP method gives the power of 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, LSP method iteratively fits sinusoidal curves with different frequencies to light curves and constructs periodogram according to the goodness of fit, and can provide accurate frequency and power spectrum intensity. 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. F AP (P n ) = 1 − (1 − prob(P > P n )) M , F AP 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 F AP = 0.01%, and the confidence level of the peak (31.2±1.6 days) is higher than 4σ. 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 firstorder 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 days 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. 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) 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 (Scargle 1981;Vecchia 1985;Feigelson et al. 2018). It is defined in statistics as or, where p, d, and q are AR order, difference operator, and MA order respectively. The seasonal (S) autoregressive (AR) Integrated (I) moving average (MA) 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 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 section at the end of this paper. Next, we fit the light curves using the ARIMA and SARIMA models respectively, and compare their goodness of fit using the Akaike Information Criterion (AIC, (Akaike 1974)). AIC = −2lnL + 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.
The fitting results are shown in Figure 5. When using the seasonal model SARIMA(4,0,1)×(6,1,2) 30days , the AIC has the global smallest value of 28.0 (Figure 5b), while the AIC value of the nonperiodic best-fit model ARIMA(4,1,3) is 489.55 (Figure 5a), which also means that the periodic model can better fit S5 0716+714 light curve. Figure 5c presents the AIC values considering different periods, we find that the best AIC occurs at the 30 ± 2.5 days position. The amplitude of the uncertainty is estimated to be half of the light curve time bin (2.5 days). SARIMA (4,0,1)×(P,1,Q) s=30days

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 frequency (for example, the red noise is produced by stochastic flares), which may produce the false QPO (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 confidence level is very important in the 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 include three-step: model the original light curve PSD and Probability Distribution Function (PDF), finally based on the parameter of PSD and PDF model to 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 where model parameters A, α low , α high , f bend and C are 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. firstly we calculated the joint probability density function (likelihood function) of obtaining the ensemble of periodogram estimates for the given PSD model, and then 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 −5 ), respectively. The light curve in AGN generally exhibits 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), where α s is the shape parameter of gamma function fitting, and the obtained value is 4.993, logµ and σ are the log mean and standard deviation of the log-normal distribution, respectively. The maximum likelihood analysis obtains the parameter values as 2.886 and 0.401, respectively. PSD and PDF based parameter values, we generated 10 5 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 dataset, where a subset with the desired length is randomly selected). Finally, according to the mean and standard deviation of the distribution of the artificial light curve PSD at each frequency, the confidence level of the QPO was estimated. The blue and red dash 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-day QPO signal in S5 0716+714 is highly significant.
Through the above analysis, we get that the quasi-periodic signal in S5 0716+714 is significant. We use ARIMA and SARIMA models to fit the light curve respectively, we find that the SARIMA model has AIC value much better than the best ARIMA model (Figure 5b), and this fact proves that the SARIMA model with strict periodic modulation component can better model the light curve. Figure 5c 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 deviation of the AIC corresponding to the SARIMA model in different seasons. It can be seen from Figure 5c 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 quasi-periodic oscillations (QPOs) in the Fermi Large Area Telescope Light Curve Repository (LCR) catalog of 1525 sources (all of which have variability index > 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-300GeV) 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 − F AP ) 1525 ∼ 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-day strictly periodic component superposed on a stochastic autoregressive variability. Together these can be considered Quasi-Periodic Oscillator (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 ARIMA and SARIMA models). Compared with the gamma-ray QPO signal previously reported in AGN, this QPO is the first month-like QPO with 7 cycles. The QPO only occurred during the outbreak state in June 2011 , 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 5 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 time scale is very close to the moon's sidereal period of 27.32 days, so it may be modulated by the moon's gamma-rays. But, 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. Secondly, the stochastic processes can produce pseudo periodic variations, but these 'phantom periodicities' typically show less than 3 cycles. When the period reaches 5 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 7 cycles, so this QPO signal is likely 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 supermassive black hole (SMBH) system, periodicity can be induced by the secondary BH by piercing the accretion disc of the primary BH during the orbit motion (Valtonen et al. 2008;Villforth et al. 2010;Ackermann et al. 2015;. The second, jet precession model, a secondary black hole in a non-coplanar 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, 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 intra day variability Wierzcholska & Siejkowski 2015;Liu et al. 2017;Yuan et al. 2017;Hong et al. 2018;Liu et al. 2019Liu et al. , 2021, and the QPO found in this source shows a month-like oscillation, so these models can not provide appropriate explanation of this short time scale QPO. The more details of models need to be 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 gamma-ray flux of the blazar can be explained in the scenario of a shockwave 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 of relativistic beam effect will lead to significant flux change even if there is no internal change in the emission of the jet. So, the blob in the jet moves towards us, 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), cosθ obs = sinφsinψcos(2πt/P obs ) + cosφcosψ, 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 δ = 1/[1 − βcosθ obs ], where Γ = 1/(1 − β 2 ) 1/2 is the bulk Lorentz factor of the blob motion with β = ν jet /c. The value of φ usually fluctuates between 0 • -5 • , here we choose φ = 2 • , ψ = 5 • , and Γ = 8.5 according the typical value of 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, Based on the known parameters, P rf was estimate as 7.57 years, we can also estimate the distance the blob moves in a period, D = cβP rf cosφ ≈ 2.09pc. The total projection distance of 7 cycles can be estimated as D p = 7Dsinψ ≈ 1.27pc. The parsec scale helical jets has been detected in several blazars (Bahcall et al. 1995;Vicente et al. 1996;Tateyama et al. 1998). The helical structure in the jet have been supported by optical polarization observations (Marscher et al. 2008), but its origin has remained unresolved. We can assume that the jet structure is generated under the action of a helical magnetic field (Vlahakis & Königl 2004).
In 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 are no longer constant, but a function of time ψ(t) (Sarkar et al. 2021). In our QPO time scale, the bending of the helical jet will not have 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 angle ψ. This provides a possible explanation for the obvious flux increase in the 7th-cycle of QPO signal.
If there are simultaneous multi-band 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, the long-term and high temporal resolution observations of blazar is very important to identify QPOs of blazar.

ACKNOWLEDGMENTS
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 generous on-going 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).