A Possible 3 yr Quasi-periodic Oscillation in γ-Ray Emission from the FSRQ S5 1044+71

Blazars, with their jets pointing almost directly to the Earth (Urry & Padovani 1995), are a special subclass of active galactic nuclei (AGNs). They show very extreme observational variability over almost the whole electromagnetic wave band. Because of the abundant optical observations, many variability findings have been claimed in the optical band (e.g., Fan & Lin 2000; Li et al. 2009; Bhatta et al. 2016; Fan et al. 2021). The most compelling sample may be OJ 287, which shows an optical periodic signal with a quasi-periodic cycle of ∼12 yr (Sillanpaa et al. 1985; Kidger et al. 1992; Valtonen et al. 2006). Thanks to the launch of the Large Area Telescope (LAT) on board Fermi in 2008 June (Atwood et al. 2009), long-coverage observations on different timescales (from seconds to years) can be provided by taking advantage of LAT's all-sky monitoring capabilities. For PG 1553+113, a 2.18 yr quasi-periodic cycle in γ-ray was first reported by Ackermann et al. (2015).

According to the optical emission line features, blazars are usually divided into two subclasses: flat-spectrum radio quasars (FSRQs), with strong emission lines, and BL Lac objects (BL Lacs), with weak or even no emission lines. Blazar emission ranges from radio to TeV, which is generally dominated by nonthermal radiation. The spectrum energy distribution (SED) shows two humps, and it is generally accepted that the lower energy hump peak of the typical multiwavelength SED of a blazar is dominated by synchrotron emission. The higher energy hump peak in the MeV–GeV band could be produced by inverse Compton scattering of synchrotron photons (Bloom & Marscher 1996; Finke et al. 2008) and external photons (e.g., from the accretion disk, broad-line region, or dusty torus; see Sikora et al. 1994; Kang et al. 2014). The γ-ray emission of FSRQs is generally produced by the external Compton (EC) mechanism.

Quasi-periodic variability studies could give insights into the physics of blazars and black hole (BH)–jet systems. Quasi-periodic oscillations (QPOs) in blazars occasionally present in optical, X-ray, and radio bands on diverse timescales. Variability is also the typical observation feature of Fermi blazars, and it is usually aperiodic. However, nearly 30 blazars are reported to have possible QPOs based on Fermi-LAT data, with timescales ranging from months to several years (e.g., Ackermann et al. 2015; Sandrinelli et al. 2016b; Prokhorov & Moraghan 2017; Zhang et al. 2017a; Zhou et al. 2018; Bhatta 2019; Peñil et al. 2020, and references therein). A year-like timescale for quasi-periodic variation often appears to occur in Fermi blazars. However, there is still no straightforward model for describing these possible periodicities. The cause of the γ-ray quasi-periodic variabilities still remains controversial. Several explanations have been proposed to explain the QPO γ-ray variabilities in blazars: (i) lighthouse effects in jets (Holgado et al. 2018); (ii) the existence of a binary system of supermassive BHs (SMBHs; Komossa & Zensus 2016); (iii) jet precession or helical structure, with periodic changes of the Doppler factor (Ackermann et al. 2015); and (iv) quasi-periodic injections of plasma into the jet, caused by pulsational accretion flow instabilities (Tavani et al. 2018). Here, a binary BH model, including an accretion model and a dual-jet model, is used for this quasi-periodic variability mechanism.

S5 1044+71 is a distant FSRQ (with redshift z = 1.15; Polatidis et al. 1995). In the latest LAT source catalog (4FGL-DR2, for Data Release 2; Ballet et al. 2020), 4FGL 1048.4+7143 is associated with S5 1044+71. It was classified as a low-synchrotron-peaked blazar (for sources with synchrotron peak frequency ${\nu }_{\mathrm{peak}}^{{\rm{S}}}\lt {10}^{14}$ Hz) by the second LAT catalog of AGNs (Ackermann et al. 2011). The LAT observed γ-ray flaring activity from S5 1044+71 in 2014 January (D'Ammando & Orienti 2014). Besides, it was reported that S5 1044+71 showed a marked increase in flux activity in γ-ray in 2016 December, which is about a factor of 16 greater than the average flux reported in the third Fermi-LAT catalog (3FGL; Ojha & Carpen 2017). Since the launch of Fermi in 2008 June, flux flares for S5 1044+71 have also been found in multiwavelength. It showed a near-infrared brightening in 2013 January, which was about 1.2 mag brighter than its previous flux (Carrasco et al. 2013). Its R-band flux was observed to be in a flaring state on 2013 October 25, which was ∼1.5 mag substantially brighter than its usual brightness (Blinov & Kougentakis 2013). Later, it showed a high-radio state from 2014 January to February (Trushkin et al. 2014a, 2014b). A significant optical enhancement was observed in 2017 January, with R = 15.44 ± 0.20 mag (Pursimo et al. 2017), which was associated with the flare state in γ-ray, as noted above.

In this paper, we have performed a detailed time series analysis of the FSRQ S5 1044+71, based on the LAT data for the interval between 2008 August and 2021 July. We present the Fermi data analysis as well as the periodicity searching methods and results in Section 2. The results are discussed in Section 3, with a summary given in Section 4.

2.1. Fermi-LAT Observations and Data Reduction

The LAT scans the whole sky every 3 hr in the energy range from 20 MeV to >300 GeV (Atwood et al. 2009). For the data selection, we chose LAT events from the Fermi Pass 8 database in the time period from 2008 August 4 15:43:36 (UTC) to 2021 July 3 00:00:00 (UTC), in the energy range 0.1–300 GeV. For the target S5 1044+71, a 20° × 20° region centered on its position was selected. Following the recommendations of the LAT team, ⁴ we selected events with zenith angles of less than 90° to prevent possible contamination from the Earth's limb. The analysis tool Fermitools 2.0.8 and the instrument response function P8R3_SOURCE_V2 were used. In addition, the background Galactic and extragalactic diffuse emission were added in the source model, using the spectral model gll_iem_v07.fits and the file iso_P8R3_SOURCE_V2_v1.txt, respectively. The normalizations of the two diffuse emission components were set as free parameters in the analysis.

We constructed light curves binned in 5 day time intervals, by performing the standard binned maximum likelihood analysis. The choice of the 5 day binning provided the shortest time intervals that were long enough for all bins to be detected (the maximum likelihood Test Statistic (TS) values being larger than 9). We also tried 1–30 day bins, but the 5 day bins were the most appropriate bin, since they not only show the details of the flux variation, but also ensure that S5 1044+71 can be detected in all bins. The source model is based on the LAT 10 yr source catalog (Abdollahi et al. 2020; Ballet et al. 2020), and the normalization parameters and spectral indices of the sources within 5° of the target, as well as the sources within the region of interest with variable index ≥ 72.44 (Acero et al. 2015), were set as free parameters. All other parameters were fixed at their catalog values in 4FGL-DR2. We used the original spectral models in 4FGL-DR2 for the sources, and a simple power law for S5 1044+71 in the source model.

Using the Fermi-LAT data, we obtained a 5 day binned light curve (12.9 yr long) for S5 1044+71. The light curve is shown in Figure 1, for which only the flux data points with the maximum likelihood TS values being larger than 9 are plotted. We can clearly see a quasi-periodic variability that nearly begins from MJD 56000 (there are actually data points with TS values larger than 9 starting from MJD 56013). The source is in a quiescent state before the possible oscillation cycle. Hence, the following period analysis only uses the LAT data over an interval of ∼9 yr (from MJD 56013 to 59298).

Zoom In Zoom Out Reset image size

2.2. Searching for γ-Ray Periodicity

Many algorithms have been used to search for variability periodicity. Here, in order to obtain the periodic component with a higher significance level, five different methods are applied to the light curve to search for the γ-ray periodicity, as follows: the Date-compensated Discrete Fourier Transform (DCDFT) method, the Jurkevich (JV) method, the Lomb–Scargle Periodogram (LSP) method, a Fortran 90 program (REDFIT), and the Weighted Wavelet Z-transform (WWZ) method are used in this work. Of these, REDFIT is used to obtain the significance of the signal, and we also make light-curve simulations to obtain the robust significance.

• (i)  
  DCDFT+CLEANest is a superior technique (Ferraz-Mello 1981; Foster 1995), which is especially powerful for unevenly spaced data. We applied it to the light curve, which can be done as described in Foster (1995). The CLEANest algorithm can clean false periodicities so as to remove false peaks. This gives a DCDFT period result of 3.06 ± 0.43 yr, and ∼3.03 yr with the CLEANest method. The period is obtained by fitting the power peak with a Gaussian function. The half-width at half-maximum of the Gaussian fitting at the position of the peak is taken as a measure for the uncertainty of the signal (Kidger et al. 1992).
• (ii)  
  The JV method is based on the expected mean square deviations, ${V}_{m}^{2}$ (Jurkevich 1971). It tests a run of trial periods, T, around which the data are folded, split into m terms. The trial period is expected to be equal to a true one when ${V}_{m}^{2}$ reaches its minimum. Later, Kidger et al. (1992) introduced a fraction reduction of the variance, $f=\tfrac{1-{V}_{m}^{2}}{{V}_{m}^{2}}$. A value of f ≥ 0.5 (meaning ${V}_{m}^{2}\leqslant 0.67$) suggests a very strong periodicity. The higher the f value, the larger the confidence of the period. For the present data, the JV method gives the minimum of ${V}_{m}^{2}$ as 0.506 (meaning f = 0.977) for a trial period of 3.06 ± 0.35 yr. The results from both the DCDFT and JV methods are shown in Figure 2.
• (iii)  
  To obtain the robust significance of the signal, we simulate light curves based on the best-fitting result of the power spectral density (PSD) obtained. The details of the simulation and significance estimation methods are given by Emmanoulopoulos et al. (2013) and Bhatta et al. (2016). Following the procedure, we simulated 10⁶ light curves with the DELCgen program, and evaluated the significance of the signal. The result indicates a significance of ∼3.6σ for the period signal of 3.06 yr, which is shown in Figure 3. We thus conclude that a ∼3.6σ QPO exists in the γ-ray light curve during MJD 56013–59298.
• (iv)  
  The REDFIT program ⁵ (Schulz & Mudelsee 2002), which is based on the LSP method (Lomb 1976; Scargle 1982), is often performed to estimate the red noise levels in light curves of blazars. This program estimates the red noise spectrum by fitting the data with a first-order autoregressive process. It can precisely evaluate the significance of the PSD peaks against the red noise background. When REDFIT is applied to the 5 day binned data, the result is as shown in Figure 4, which shows that the significance of the signal peak is higher than a 99% confidence level. The 2.96 yr peak is within the error range of the results of the previous methods, and the significance of the two smaller peaks (1.52 and 0.81 yr) is relatively lower. Note that the REDFIT program only provides a maximum significance of 99%, which corresponds to ∼ 2.5σ.
• (v)  
  We also use the WWZ method (Foster 1996) to search for QPOs. WWZ is a period extraction algorithm, based on wavelet analysis and vector projection. It is very suitable for the analysis of nonstationary signals, and it has advantages in time–frequency local characteristic analysis. When the method is applied to the present γ-ray data, the corresponding DCDFT and time-averaged WWZ powers are as shown in the right-hand panel of Figure 5, while the 2D plane contour map of the WWZ power spectrum is as shown in the left-hand panel of Figure 5. The result shows a clear peak at ∼3.08 yr, with an uncertainty resulting from a Gaussian fitting of ± 0.36 yr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For clarity, we list all the results obtained by using the five methods in Table 1. We also show the fitting results of the average periodicity using these methods as the dashed red curve in Figure 1.

Table 1. Periodicity Searching Results

Method	Period (yr)
DCDFT	3.06 ± 0.43
JV	3.04 ± 0.35
LSP	3.06 ± 0.44
REDFIT	2.96 ± 0.40
WWZ	3.08 ± 0.36

Download table as:  ASCIITypeset image

We folded the γ-ray light curve using a phase-resolved binned likelihood analysis with a 3.06 yr period, obtained by using the DCDFT method. The folded light curve with phase zero corresponding to MJD 56012.66 is shown in Figure 6, in which 16 phase ranges are set. This result also clearly confirms the signal that the amplitude of the γ-ray flux varies with phase. The folded γ-ray spectral photon index is also given in Figure 7. When comparing the spectral shapes of the different phases, a correlation between flux and photon index is clearly visible, suggesting the tendency of a harder-when-brighter pattern, which is usually seen in blazar flares (e.g., Hayashida et al. 2015; Shukla et al. 2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have carried out a temporal analysis of γ-ray observations of the FSRQ S5 1044+71 by Fermi-LAT from 2008 to 2021. Our results reveal a quasi-periodic variability in γ-ray, with a period cycle of 3.06 ± 0.43 yr at a significance level of 3.6σ. γ-ray QPOs with a significance of ≥3σ have mainly been reported in BL Lacs, especially in high-synchrotron-peaked BL Lacs (HBLs; Sandrinelli et al. 2014, 2016b; Ackermann et al. 2015). It is interesting that quasi-periodic variabilities are also found in other subclasses of blazars. A quasi-period of 3.35 ± 0.68 yr in the γ-ray light curve was reported for the FSRQ PKS 0426–380. (Zhang et al. 2017b). The 3.06 yr quasi-period of S5 1044+71 is very similar to that of PKS 0426–380. Interestingly, although these two FSRQs have longer observed periods (T_obs ∼ 3 yr) than the three HBLs (PKS 2155–304, PKS 0301–243, and PG 1553+113; T_obs ∼ 2 yr), their intrinsic periods (T_sou) are almost the same as for the three HBLs on account of T_obs = T_sou(1 + z). For the FSRQ subclass of blazars, relatively efficient broad-line-region (BLR) emission lines and accretion disk emission are present (D'Ammando et al. 2011). Since S5 1044+71 is an FSRQ, the emission from the accretion disk, the BLR, and the jet will be expected to contribute to the total γ-ray emission from the blazar by the EC mechanism.

Since the launch of LAT, more QPOs in γ-ray have been reported. There have been several analyses of systematic searches for QPOs in Fermi-LAT γ-ray sources based on 3FGL (e.g., Prokhorov & Moraghan 2017; Bhatta & Dhital 2020; Peñil et al. 2020; Zhang et al. 2020b). Together with the previous studies (e.g., Sandrinelli et al. 2014, 2016b; Ackermann et al. 2015; Zhang et al. 2017a; Bhatta 2019, and references therein), there are nearly 30 possible blazar QPOs. Almost all of them have year-long periods. PKS 2247–131 is the first case that exhibits a clear month-like oscillation of 34.5 days. Relatively shorter QPOs have also been detected. Gupta et al. (2019) reported a ∼71 day period for B2 1520+31. Sarkar et al. (2021) reported a dominant period of ∼47 days in γ-ray and optical light curves for 3C 454.3, covering nine cycles over 450 days of observations, which is the highest number of cycles ever detected in a blazar light curve. Bhatta (2019) reported a 330 day subyear timescale γ-ray QPO that persisted for nearly seven cycles in Mrk 501. Such cases are rare, as there are not many γ-ray QPOs have been detected as lasting more than five cycles. Prokhorov & Moraghan (2017) performed a systematic search for QPOs over a period ranging from days to years in Fermi-LAT γ-ray sources. They confirmed three γ-ray blazar QPOs that were claimed previously, including PKS 2155–304, PG 1553+113, and BL Lacertae. In addition, they also found evidence of possible periodic behaviors of four other blazars: S5 0716+71, 4C 01.28, PKS 0805–07, and PKS 2052–47. Zhang et al. (2020b) found periodic signals in 4C 01.28 and S5 0716+71 during the observation period from 2008 August to 2016 December, but the signals disappeared in the interval from 2008 August to 2018 February. This reminds us of the complexity of AGN QPO analysis, and a similar concern is also proposed in Covino et al. (2019).

Peñil et al. (2020) performed a systematic periodicity search study using nine years of LAT data with 10 different techniques, and found 11 AGNs showing periodicity signals at a higher than 4σ significance level from at least four algorithms. The periods of nine of the 11 AGNs had not been reported before. This condition is identified with the high-significance tag, while in Bhatta & Dhital (2020) this criterion corresponds to a significance higher than 99%. Zhang et al. (2020b) also investigated whether there is a relation between γ-ray QPO frequency and BH masses in AGNs, and found no significant correlation. It is important to note that the number of cycles covered by LAT is unavoidably small, with year-long periods, and this may affect the estimates of claimed periodicity and the significance. Furthermore, the methods and the criteria for confirming a high-significance QPO still remain controversial, so it is quite hard to make a definite list of confirmed γ-ray blazar QPOs. From the literature, it is interesting to find that 30 blazars are reported as showing periodic signals, as listed in Table 2. Out of the 30 blazars, 13 are FSRQs and 17 are BL Lacs. Of the 17 BL Lacs, there are 7 HBLs, 5 intermediate-synchrotron-peaked BL Lacs (IBLs), and 5 low-synchrotron-peaked BL Lacs (LBLs), if we adopt the classifications of Abdo et al. (2010) and Fan et al. (2016).

A periodically changing viewing angle causing a varying Doppler factor is a possible interpretation for γ-ray QPOs, which is related to a helical jet (Camenzind & Krockenberger 1992). The presence of a binary SMBH is one way of interpreting the helical structure (Sobacchi et al. 2017). A long-term quasi-periodic variability, as in our case, is well explained by the binary BH model (Sillanpaa et al. 1988; Fan et al. 2002, 2007, 2021; Valtonen et al. 2008). This mechanism can be interpreted by two different models, as described in Qian et al. (2007, 2014). The first model is the accretion model. This could be described as the accretion rate increasing when the secondary BH passes through the primary BH, so as to cause periodic flux flares. The orbital period is the time interval between the two flaring peaks. The second model, the dual-jet model, does not consider the change in the accretion rate. Rather, the change in the observational angle will cause the periodic change of the Doppler factor δ, and finally lead to the periodic flux flares in observations, as shown in Qian et al.(2007).

Accretion model: in order to obtain the intrinsic orbital parameters of a binary system, the observed period is corrected for the cosmological expansion effect from the orbital period using the redshift, T_sou = T_obs/(1 + z). For the binary BH system, Kepler's law of motion gives the following relationship:

$$\begin{eqnarray}&&{T}_{\mathrm{sou}}^{2}=\displaystyle \frac{4{\pi }^{2}{(a+b)}^{3}}{G(M+m)},\end{eqnarray} \tag{ 1 }$$

where T_sou represents the intrinsic orbital period (=1.42 yr in the present source) and a and b represent the major and minor axes, respectively. G represents the universal gravitational constant. M and m represent the masses of the main BH and the secondary BH. Equation (1) can be equivalent to the following formula (Fan et al. 2010; see also Fan et al. 2021):

$$\begin{eqnarray}&&{T}_{\mathrm{sou}}\sim 1.72{M}_{8}^{-1/2}{r}_{16}^{3/2}{\left(1+\displaystyle \frac{m}{M}\right)}^{-1/2}\,\mathrm{yr},\end{eqnarray} \tag{ 2 }$$

where M₈ is the primary BH mass in units of 10⁸ solar mass and r₁₆ = a + b is the orbital radius in units of 10¹⁶ cm, respectively. The period calculated by this model from the light curve is considered to be the orbital period of the binary BH system. Thus, assuming m/M ∼ 0.001, as for OJ 287 in Sillanpaa et al. (1985), and adopting the mass of the BH M = 14.5M₈ (Paliya et al. 2021), we obtain r₁₆ = 2.14, namely a + b = 2.14 × 10¹⁶ cm.

Dual-jet model: from the γ-ray light curve, we can see three distinct flux peaks, as shown in Figure 8. This phenomenon is very similar to the radio and optical periodic flares of 3C 454.3 in Qian et al. (2007) and Fan et al. (2021). Therefore, we also consider the dual-j et model. This model involves two jets coming from two massive BHs that rotate periodically. In this way, the two observation angles (θ) corresponding to the two jets can be obtained (Qian et al. 2007):

$$\begin{eqnarray}&&\cos {\theta }_{1}(t)=\sin {\psi }_{1}\cos (\omega t+{\phi }_{1})\sin i+\cos {\psi }_{1}\cos i,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\cos {\theta }_{2}(t)=\sin {\psi }_{2}\cos (\omega t+{\phi }_{2})\sin i+\cos {\psi }_{2}\cos i.\end{eqnarray} \tag{ 4 }$$

Hence, ϕ₁ and ϕ₂ are the azimuths of the orbital plane. The two angles refer to the same azimuth angle (ϕ₁ = 145°) between the observer and component 1 when t = 0, and the angle (ϕ₂ = 215°) between the observer and component 2. The two components make the angle (ψ₁ and ψ₂) with the orbital normal, and the sight direction of the observer forms the angle (i) with the normal of the orbital plane. Based on this, we can obtain the expression of the Doppler factor (δ₁ and δ₂) and the apparent velocity (β_app1 and β_app2) changing with (t) (Qian et al. 2007):

$$\begin{eqnarray}&&{\delta }_{1}(t)={[{\rm{\Gamma }}(1-{\beta }_{1}\cos {\theta }_{1})]}^{-1},{\delta }_{2}(t)={[{\rm{\Gamma }}(1-{\beta }_{2}\cos {\theta }_{2})]}^{-1},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\beta }_{\mathrm{app}1}^{2}(t)=-{\delta }_{1}^{2}-1+2{\delta }_{1}{\rm{\Gamma }},{\beta }_{\mathrm{app}2}^{2}(t)=-{\delta }_{2}^{2}-1+2{\delta }_{2}{\rm{\Gamma }},\end{eqnarray} \tag{ 6 }$$

where Γ is the Lorentz factor. Then, the change of the Doppler factor (δ₁ and δ₂) eventually leads to the change of the observed flux (S₁ and S₂):

$$\begin{eqnarray}&&{S}_{1}(t)={S}_{b1}{\delta }_{1}^{x},{S}_{2}(t)={S}_{b2}{\delta }_{2}^{x},\end{eqnarray} \tag{ 7 }$$

where S_b1 and S_b2 are normalization constants for fitting the observed light curves. The intrinsic period is T_sou = 1.42 yr (=3.06/(1 + 1.15)). The fitting results of this model are listed in Table 3. We also obtained a quiescent flux of 0.50 × 10⁻⁷ ph cm⁻² s⁻¹ and x = 4. S(t) = S₁(t) + S₂(t) + 0.50 is shown with the solid green line in Figure 8. The three main flux peaks and the three secondary peaks are well fitted by the model.

Zoom In Zoom Out Reset image size

According to the model fitting results, we can further discuss the range of Doppler factors of periodic sources. The range of component 1 is 3.73 < δ₁ < 16.92, and the range of component 2 is 0.84 < δ₂ < 6.51, as shown in Figure 9. For FSRQ cases, a Doppler factor range of 5 < δ < 18 and the mean value δ = 13.16 were obtained from Ghisellini et al. (2014) and Zhang et al. (2020a). Our result is very close to theirs. Meanwhile, Fan et al. (2013) calculated the Doppler factors of 138 Fermi blazars, where the result for S5 1044+71 was δ_γ = 9.33. This result was located in the range of δ₁, corresponding to the time period of the first jet as well as the flux flare state, and it would be interesting to know the description of where the jet came from in our work. In addition, we collected the γ-ray periodic sources reported previously and compared their Doppler factors with S5 1044+71 in column 6 of Table 2. This shows that the Doppler factors of all the periodic sources are roughly in the same range. In addition, we also calculated the apparent velocity (β_app). β_app1 of component 1 is 15.70 < β_app1 < 29.93, and that of component 2 is 7.55 < β_app2 < 20.53, as shown in Figure 10. According to our calculations, S5 1044+71 is expected to be a superluminal source, due to its apparent velocity β_app > 1. Recently, superluminal velocity has been found in the range of 0.53 < β_app < 34.80, in Xiao et al. (2019). Therefore, our conclusion is consistent with the velocity range.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It can be seen from Figure 8 that the fitting of the modeling curve and the observed light curve is good. Therefore, we can not only explain that the periodic flux flare is probably mainly caused by the Doppler boosting effect, but we can also predict the next flux peak, in 2022 September. In Table 4, we give the flare dates of the observations and the model. At the same time, we also obtain the oscillation range of the Doppler factors and the apparent superluminal factors of the two components, so as to understand the oscillations of the two jets caused by the Doppler boosting effect. However, sometimes the observed light curve deviates from the model light curve (e.g., on MJD 56239.5 and MJD 57109.5), which indicates that the change in the accretion rates is likely to be one of the causes, and not just the Doppler boosting effect. Therefore, we propose that the periodic oscillation behavior of S5 1044+71 may be caused by the Doppler boosting effect, with the irregular variation of the accretion rate as a supplement.

In this paper, we have analyzed the Fermi-LAT data of S5 1044+71 from 2008 to 2021. We have used five different methods to search for its γ-ray periodicity, and come to the following conclusions:

• (1)  
  Our results reveal a possible quasi-period of 3.06 ±0.43 yr with a significance of 3.6σ in its γ-ray light curve.
• (2)  
  A binary BH model, including an accretion model and a dual-jet model, is used to explain this year-long possible oscillation, with the results fitting the observation flares. This suggests that the periodic flux flare behavior in γ-ray may be caused by the Doppler boosting effect, supplemented by the irregular variation of the accretion rate.
• (3)  
  We further calculate and discuss the range of the Doppler factor and the apparent velocity for S5 1044+71, which is expected to be a superluminal source and to show a flare in 2022.

Thanks are given to the reviewer for the constructive comments and suggestions. The work is partially supported by the National Natural Science Foundation of China (NSFC U2031201 and NSFC 11733001) and the Guangdong Major Project of Basic and Applied Basic Research (grant No. 2019B030302001).

We also acknowledge the science research grants from the China Manned Space Project with NO. CMS-CSST-2021-A06.