Optical and γ-ray variability analysis of BL Lacertae object TXS 1902+556

ABSTRACT

The variability data for the BL Lacertae object TXS 1902+556 in the optical and |$\gamma$|-ray wavebands were obtained from the 0.76-m Katzman Automatic Imaging Telescope and the Fermi Large Area Telescope (Fermi-LAT), covering periods of 14.4 and 14.7 yr, respectively. The variability properties were systematically analysed, with particular emphasis on the first comprehensive investigation of radiation variation in the optical waveband. Four well-established techniques were employed for this purpose: the Lomb–Scargle periodogram, REDFIT program, Jurkevich method, and discrete correlation function (DCF) approach. The optical waveband exhibits quasi-periodic oscillations (QPO) with a time-scale of |$P_{\rm O}=276.8\pm 6.1$| d at a significance level |$3.87\sigma$|⁠, while the |$\gamma$|-ray waveband does not exhibit any significant periodicity. However, it should be noted that the QPO time-scale is consistent with the Sun-gaps in the optical light curve within 2|$\sigma$| uncertainties. The optical QPO behaviour is most likely attributed to the helical motion of the jet driven by the orbital motion in a supermassive black hole binary system. Moreover, we have provided an explanation for the absence of QPO in the |$\gamma$|-ray light curves. Furthermore, utilizing the DCF method, a weak correlation between the variability in the optical and |$\gamma$|-ray wavebands was observed, suggesting that the emission of TXS 1902+556 may be generated through a combination of synchrotron self-Compton (SSC) and external Compton (EC) processes, or a leptonic–hadronic hybrid process.

1 INTRODUCTION

Blazars, which comprise flat spectrum radio quasars (FSRQs) and BL Lacertae objects (BL Lacs), as a special subclass of active galactic nuclei (AGNs) distinguished by their rapid and violent variability across an extensive range of wavebands (Urry & Padovani 1995; Böttcher 2019). Differentiating FSRQs and BL Lacs can be done based on their emission lines: FSRQs display a strong emission line, while BL Lacs exhibit weak or absent ones, characterized by an equivalent width of |${\rm EW} < 5\,\mathring{\rm A}$| (Böttcher 2019). Furthermore, blazars are classified into three categories according to the synchrotron component’s peak frequency |$\nu _{\rm peak}$| in the spectral energy distribution (SED): low energy peak blazar with |$\nu_{\rm peak} < 10^{14}\,{\rm Hz}$|⁠, intermediate energy peak blazar with |$10^{14}\,{\rm Hz} < \nu_{\rm peak} < 10^{15}\,{\rm Hz}$|⁠, and high energy peak blazar with |$\nu_{\rm peak} > 10^{15}\,{\rm Hz}$| (Abdo et al. 2010). Their emission is predominantly characterized by non-thermal radiation from the jet, evincing robust variability (Li et al. 2015; Bhatta 2021).

Variability is an influential tool for studying and understanding the physics of blazars (Li et al. 2022a; Otero-Santos et al. 2023). The analysis of variability can facilitate the investigation of the radiation mechanism and provide constraints on the parameters of the physical model (Lainela et al. 1999; Chandra et al. 2014; Li et al. 2022a; Otero-Santos et al. 2023). Detecting periodicity and correlation between different wavebands in blazars can aid in interpreting the relativistic emission within various theoretical frameworks. Moreover, the significant correlation between different wavebands supports lepton processes rather than hadron processes (Liao et al. 2014). The time-scales of variability in blazars span a wide spectrum, ranging from minutes to years (Otero-Santos et al. 2020), classified into three categories based on time-scales: intra-day variability (IDV), short-term variability (STV), and long-term variability (LTV) (Gupta et al. 2008; Gaur et al. 2012; Li et al. 2015, 2016; Otero-Santos et al. 2020).

Although blazar variability is commonly perceived as inherently stochastic, quasi-periodic oscillations (QPO) have been observed in certain blazar sources across the electromagnetic spectrum on diverse time-scales, spanning from minutes to days, weeks to years, and even decades (e.g. Sillanpaa et al. 1988, 1996a; Lainela et al. 1999; Ostorero, Villata & Raiteri 2004; Xie et al. 2005, 2008; Li et al. 2009, 2015, 2016, 2023; Rani et al. 2010; Sandrinelli, Covino & Treves 2014; Ackermann et al. 2015; Sandrinelli et al. 2016, 2018; Hong, Xiong & Bai 2018; Otero-Santos et al. 2020, 2023; Ren et al. 2021; Tripathi et al. 2021, 2024; Roy et al. 2022a, b). The investigation of QPO in blazars provides a robust approach to delineate the regions of variability and their associated physical mechanisms. Several physical models have been proposed to elucidate the phenomenon of QPO in blazars, including supermassive binary black hole (SMBBH) systems, instabilities in the accretion flow, and geometric models, etc. The occurrence of QPO in blazars is a complex phenomenon. In certain blazars, diverse QPO have been reported within the same blazar and across different wavebands. Additionally, there have been instances where identical QPO time-scales were observed at multiple wavebands. For instance, the optical band observations of OJ 287 reveal a long-term QPO with an approximate periodicity of 12 yr, accompanied by a double-peak structure separated by approximately 1.2 yr (Sillanpaa et al. 1988, 1996a). The long-term QPO can be explained under the framework of SMBBH systerm, and the double-peak structure may be caused by the small black hole passing through the accretion disc of the big on twice per revilution (Sillanpaa et al. 1996a, b). However, Gorbachev et al. (2024) proposed that the 12 yr QPO observed in optical wavelengths can be attributed to variations in the Doppler factor resulting from the presence of a helical structure within the jet. In radio band, a QPO with time-scale 1.12 yr was reported by Hughes, Aller & Aller (1998), which can be understood in terms of a shock-in-jet model. Moreover, Kushwaha et al. (2020) reported the presence of a 314 d QPO in |$\gamma$|-ray emissions from OJ 287. They proposed that both variations, based on accretion discs and their impact on the jet, as well as a purely jet-based origin, offer plausible explanations for the observed QPO behaviour. Furthermore, Pihajoki, Valtonen & Ciprini (2013) confirmed potential QPO with time-scale of a few hours and approximately 50 d in the optical band, while Bhatta et al. (2016) reported a QPO with a period of 400 d. Additionally, Sandrinelli et al. (2016) identified periodicities of around 435 d and approximately 412 d in the NIR-optical and |$\gamma$|-ray light curves, respectively. Moreover, S5 0716+714 also exhibits a diverse range of time-scales in its QPO. In the optical waveband, Hong et al. (2018) identified a QPO with a period of 50 min. Heidt & Wagner (1996) detected a QPO lasting for four days and proposed that this periodicity may originate from the rotation of a hotspot on the accretion disc. Furthermore, Dai et al. (2021) and Yuan et al. (2017) reported QPO behaviours in S5 0716+714 with periods of approximately 6 and 24 d, respectively, suggesting that these variabilities are associated with the jet. A long-term QPO signal about 1060 d in the optical was reported by Yang et al. (2023), who used an SMBBH model to explain it. Raiteri et al. (2003) also reported a QPO behaviour with time-scales of 3.3 yr in the optical band and 5–6 yr in radio bands; Li et al. (2023) found similar signals at different wavelengths: 352 d signal at radio light curve and 960 d signal at |$\gamma$|-ray light curve. Li et al. (2018b) also reported the presence of a long-term QPO in the radio band, with a time-scale of approximately 6.1 yr, suggesting the existence of a SMBBH system. Li et al. (2018a) also investigated the variability of S5 0716+714 in both its 15 GHz radio and |$\gamma$|-ray light curves, revealing QPO time-scales of 266 and 344 d, respectively. Additionally, Chen et al. (2022) analysed the |$\gamma$|-ray light curves of S5 0716+714 and identified a QPO signal with a period of about 31 d. They proposed an explanation for this periodicity involving a spiral motion along the jet caused by a moving blob. The observation of QPO phenomena with different time-scales has been made in numerous blazars, including AO 0235+164, PKS 2155+304, 3C 66A, and 3C 454.3, etc.These observations have prompted various theoretical explanations. Therefore, the investigation into QPO remains an open area for research.

The blazar TXS 1902+556, also known as 4FGL J1903.2+5540 and J1903+5540, is classified as an intermediate energy peak BL Lac object (IBL) and has a redshift of |$z=0.727$| (Shaw et al. 2013; Goldoni et al. 2021). Notably, its spectrum indicates the presence of an Mg ii system at the same redshift (Shaw et al. 2013; Goldoni et al. 2021). The emission variability of TXS 1902+556 has been explored in the literature (e.g. Liodakis et al. 2018; Peñil et al. 2020; Yang et al. 2021; Rueda, Glicenstein & Brun 2022). Quasi-periodic variations with time-scales of approximately 3.8 yr have been reported in the |$\gamma$|-ray waveband by Peñil et al. (2020) at a low significance level |$({>} 2.5\sigma)$|⁠. Additionally, a period of approximately 3.1 yr (1120 d) was identified in |$\gamma$|-ray light curve by Rueda et al. (2022). The multiwavelength cross-correlations of TXS 1902+556 were investigated by Liodakis et al. (2018). They found time lags between optical and radio waveband, optical and |$\gamma$|-ray waveband, |$\gamma$|-ray and radio waveband to be |$66.45\pm 1.02$| d, (⁠|$1.54\sigma$| significance) |$5.03\pm 3.51$| d (⁠|$2.47\sigma$| significance), and |$169.91\pm 10.31$| d (⁠|$1.74\sigma$| significance), respectively.

In this paper, we will investigate the variability properties of TXS 1902+556 in the optical and |$\gamma$|-ray wavebands, along with an examination of the emission mechanisms. Section 2 provides a detailed description of the data reduction and analysis procedures. We then proceed with the periodicity analysis in Section 3, followed by correlation analysis in Section 4. Finally, Section 5 presents the discussion and conclusion.

2 DATA REDUCTION AND ANALYSIS

2.1 The optical waveband variability data

The optical waveband observation data were obtained from the 0.76-m Katzman Automatic Imaging Telescope (KAIT¹), an entirely robotic telescope situated at Lick Observatory atop Mount Hamilton, east of San Jose, California. Since March 2009, KAIT has observed 163 blazars with unfiltered magnitudes. The data of KAIT can be converted to the R-band magnitude taking into account the Galactic extinction of |$A=0.26$| (Li et al. 2003).

The optical data for TXS 1902+556 utilized in this paper spanned a duration of approximately 14.4 yr, ranging from 2009 May 11 to 2023 September 27 and consisted of a total of 526 data points. The light curve of TXS 1902+556 is depicted in the left picture of Fig. 1. To assess the significance of the variability, the variability index is calculated using the formula described in Fan et al. (2002) and Kovalev et al. (2005),

and the uncertainty is given by the following equation (Singh et al. 2018)

where |$F_{\max}$| and |$F_{\rm min}$| are the maximal and minimal flux, |$\sigma _{\max}$| and |$\sigma _{\rm min}$| are the error associated with the maximal and minimal flux, respectively. Singh et al. (2018) suggested that the variability is statistically significant when |$V\ge 3\triangle V$|⁠. In Fig. 1, the unit of the data is magnitude which can be converted to flux density in the unit Jansky (Jy) by the equation,

here, |$F_{0}$| is the zero-point and m is the data in magnitude. The error in the unit Jy can be given by the following formula:

where |$\triangle m$| is the error of the magnitude m. Then, the maximal flux |$F_{\max}$|⁠, minimal flux |$F_{\rm min}$| and corresponding error |$\sigma _{\max}$| and |$\sigma _{\rm min}$| become

in above formula, |$m_{\max}=16.72$| mag and |$m_{\rm min}=15.47$| mag represent the maximal and minimal magnitudes, and |$\triangle m_{\max}=0.029$| mag and |$\triangle m_{\rm min}=0.015$| mag denote the errors associated with the maximum and minimum magnitudes, respectively. Using equations (1), (2), and (5), the calculated variability index in the optical waveband is |$V_{\rm O}=0.52\pm 0.01$|⁠, implying the significant variability of TXS 1902+556 in the optical waveband.

2.2 High-energy |$\gamma$|-ray: Fermi-LAT data

We analyse the |$\gamma$|-ray light curve in the energy range of 0.3–300 GeV of TXS 1902+556 (4FGL J1903.2+5540), as observed by the Fermi Gamma-ray Space Telescope (Fermi-LAT). We selected the Fermi-LAT Pass 8 Front+Back events (evclass = 128, evtype = 3) data collected from 2008 August 4 to 2023 April 23 (mission elapsed time: 239557417–703904511, ∼14.7 yr), and perform within a |$12^{\circ } \times 12^{\circ }$| square region centred at the position of J1903.2+5540 in the incremental version of the fourth Fermi-LAT source catalogue (Abdollahi et al. 2022). We only considered events within a maximum zenith angle of 90° to minimize contamination from the Earth Limb and applied the expression '|${\rm(DATA\_QUAL > 0)\&\&(LAT\_CONFIG=1)} $|' as well as the instrumental response function '|$\rm{\tt P8R3\, SOURCE\, V3} $|' to select the high-quality data in good time intervals. Subsequently, a model file, encompassing the Galactic diffuse background emission (gll_iem_v07) and isotropic emission (iso_P8R3_SOURCE_V2_v1.txt) models, was created using the user-contributed software make4FGLxml.py. In 4FGL-DR3, the emission of the target source is described by a log-parabolic spectrum model, |${\rm d}N/{\rm d}E = K(E/E_{\rm b})^{-\alpha -\beta \log(E/E_{\rm b})}$|⁠. Employing the binned maximum-likelihood analysis, the best-fit photon spectral parameters |$\alpha$|⁠, |$\beta$|⁠, and |$E_{\rm b}$| were determined to be 1.85 |$\pm$| 0.03, 0.06 |$\pm$| 0.01, and 1.15 |$\pm$| 0.14 GeV, respectively. Based on these parameters, the integrated photon flux was calculated to be |$(9.97\pm 0.21)\times 10^{-9}\,{\rm photons\,cm}^{-2}\,{\rm s}^{-1}$|⁠, with a corresponding test statistic (TS) value of approximately 12 245. We constructed light curves for different bins and found that the 30 d bin light curve with TS |$\geqslant$| 16 (⁠|$4\sigma$|⁠) effectively captured the overall variability characteristics of the source (refer to the right panel of Fig. 1). The maximum and minimum |$\gamma$|-ray fluxes were determined to be |$F_{\max}=(24.1\pm 3.91)\times 10^{-8}\, {\rm photons\, cm}^{-2}\,{\rm s}^{-1}$| and |$F_{\rm min} = (9.02\pm 4.82) \times 10^{-9}\,{\rm {\rm photons\, cm}}^{-2}\,{\rm s}^{-1}$|⁠, respectively. Subsequently, utilizing equations (1) and (2), the |$\gamma$|-ray variability index |$V_{\gamma }=0.93\pm 0.04$| was computed, revealing that the emission variability of TXS 1902+556 is significant at |$\gamma$|-ray waveband.

3 PERIODICITY ANALYSIS

Variability analysis is a potent tool for probing emission properties and provides valuable insights into the physical mechanisms of blazars. Confirmed periodic variability can impose stringent constraints on physical parameters and models (Lainela et al. 1999). In this paper, we aim to analyse the variability periodicity in the optical and |$\gamma$|-ray wavebands of TXS 1902+556 using the Lomb–Scargle periodogram (LSP) method (Lomb 1976; Scargle 1982), the Jurkevich method (Jurkevich 1971), and the REDFIT program (Schulz & Mudelsee 2002) to study the emission variability properties.

The LSP method is a well established and widely used traditional method (Peñil et al. 2020; Li et al. 2022a) for assessing variability time-scales in astronomical time series. This method can identify periodicity in time series, irrespective of whether the time series is uniform or uneven. The LSP algorithms (Lomb 1976; Scargle 1982) have been thoroughly documented in numerous sources (e.g. Zheng et al. 2008; Li et al. 2015, 2016, 2022a; Yang et al. 2020; Wang & Jiang 2021 ). Based on the LSP method, clear peaks would be evident in the periodogram at time-scales corresponding to the periods present in a time series. Furthermore, we estimated the period errors using the half width at half-maximum (HWHM) of the peak (Li et al. 2009, 2015, 2022b). The confidence level of the LSP results was determined using Monte Carlo simulations (Yang et al. 2020; Wang & Jiang 2021). Initially, 10 000 artificial light curves were simulated to have the same observational properties as the original data set shown in Fig. 1. The simulations were conducted based on a power-law red noise model (⁠|${\rm PSD} \propto f^{-\beta }$|⁠) using the algorithm proposed by Timmer & König (1995), with power-law slopes of |$\beta _{\rm O}=1.36\pm 0.25$| and |$\beta _{\gamma }=0.82\pm 0.17$| for the optical and |$\gamma$|-ray wavebands, respectively. The estimation of these power-law slopes |$\beta$| was performed using the Power Spectrum Response Method (PSRESP²; Uttley, McHardy & Papadakis (2002)), which also provides the uncertainty associated with the estimated slope and quantifies the goodness of fit through a success fraction measurement. The success fraction quantifies the deviation between the data and the fit for each scanned slope value, enabling identification of the optimal slope that accurately reproduces the derived power spectral density (PSD). To estimate |$\beta$|⁠, we employ a simulation approach based on Timmer & König (1995) generating 1000 light curves with identical observational properties as those in the original data set. The |$\beta$| values were scanned over the range of 0 to 3.0 with an interval of 0.02, yielding a |$\beta$| distribution obtained using the PSRESP method in both optical and |$\gamma$| wavebands as illustrated in Fig. 2. Then, the LSP of the 10 000 artificial light curves was computed, and the confidence level of the LSP results was established using the spectral distribution of the artificial light curves. The 99.7 per cent (3|$\sigma$|⁠) confidence interval is evaluated by considering a range of slopes within the uncertainties. The LSP results and the confidence interval are plotted in the top panels of Fig. 3, illustrating the LSP spectrum and the 99.7 per cent (3|$\sigma$|⁠) confidence interval. The left-hand and right-hand panels of the upper part in Fig. 3 depict the LSP results of the optical and |$\gamma$|-ray wavebands, respectively. From the left-hand panel of the upper part in Fig. 3, one can find that a distinct peak is observed at time-scales of |$276.8\pm 6.1$| d, with a local significance level (⁠|$p_{\rm local}$|⁠) exceeding 3|$\sigma$|⁠. This suggests the presence of potential periodicities of |$276.8\pm 6.1$| d in the optical waveband light curve. Furthermore, the right-hand panel of the upper part in Fig. 3 reveals a peak at a time-scale of |$1180.2\pm 273.8$| d, which is consistent with the 3.1 yr periodicity of |$\gamma$|-ray emissions as reported by Rueda et al. (2022). However, the local significance level of the peak is lower than 3|$\sigma$|⁠.

To confirm the behaviour of the QPO with a time-scale of 276.8 d in the optical wavelength, we divided the light curve into two segments: one spanning from MJD 54963 to 57670 and the other ranging from MJD 57900 to 60215. The light curves of these two segments are illustrated in the top panels of Fig. 4. The duration of these two segments is respectively 2707 and 2315 d, which exceeds the QPO time-scale by more than six times. Subsequently, we analysed the QPO behaviour in both segment light curves using the LSP method and presented the results in Fig. 4 at the bottom panels. The bottom panels of Fig. 4 demonstrate the presence of QPO behaviour in both segment light curves, with confidence levels exceeding 3|$\sigma$|⁠. These findings correspond to time-scales of 279.0 and 275.6 d for the two segment light curves, respectively, which are consistent with the QPO time-scales obtained throughout the entire light curve. Hence, the existence of QPO behaviour with a periodicity of |$276.8\pm 6.1$| d in the optical band appears to be plausible. Moreover, the LSP analysis of the sector1 light curve in Fig. 4 reveals an additional periodicity of 398.0 d, exhibiting a significance level of approximately 3|$\sigma$|⁠. Additionally, the sector2 light curve’s LSP results reveal two additional peaks exceeding the 3|$\sigma$| level at time-scales of 1052.0 and 321.5 d. However, there are two distinct peaks in the LSP significance level curve at time-scales of approximately 320 and 400 d (see Figs 3 and 4), suggesting that these time-scales are attributed to red noise. In addition, the significance level of the time-scale of 321.5 and 398.0 d is lower than 3|$\sigma$| in the overall light curve results. Hence, the time-scales of 398.0 d in the first light curve and 321.5 d in the second light curve do not arise from genuine physical variability but may be attributed to red noise. Additionally, considering the time span of 2315 d for the sector2 light curve, which is merely twice as long as the duration of 1052.2 d, it is not possible to definitively confirm the time-scale of 1052.5 d therefore further observational data is crucial to substantiate it. Similarly, we also partitioned the |$\gamma$|-ray light curve into two intervals (MJD 54683−57353 and MJD 57383−60053) and analysed their variations time-scale using LSP method. The segmented light curves and the corresponding LSP results are presented in Fig. 5, indicating that the significance of all time-scales is below 3|$\sigma$|⁠.

In order to further enhance the credibility and authenticity of the identified periods, the global significance |$p_{\rm global}$| was estimated to account for look-elsewhere effect (Gross & Vitells 2010), which corrects for local significance |$p_{\rm local}$| (Cumming 2004; Zechmeister & Kürster 2009). The global significance |$p_{\rm global}$| is evaluated by employing the subsequent equation,

where N, representing the number of independent frequencies, can be estimated by counting the number of peaks observed in the periodogram (Cumming 2004; Zechmeister & Kürster 2009). Furthermore, the estimation of the number of independent frequencies can be achieved by calculating the separation between peaks in the periodogram |$\delta f$| and comparing it with the frequency range searched f, according to the following formula (Cumming 2004; Zechmeister & Kürster 2009),

Based on the equation (7), the number of independent frequencies, N, equal to 107. The local significance level of the time-scale 276.8 d is approximately 3.87|$\sigma$| (⁠|$p_{\rm local,O}=0.00011$|⁠), which has been estimated using previously described Monte Carlo simulations and power-law slopes of |$\beta _{\rm O}=1.36$|⁠. Consequently, considering equation (6) and |$p_{\rm local,O}=0.00011$|⁠, the global significance is calculated as |$p_{\rm global,O}=0.0117$|⁠, indicating a significant period of 276.8 d with a global significance level of 2.52|$\sigma$|⁠. Similar to the calculation of global confidence in optical bands, we have estimated the global confidence level for a time-scale of 1180.2 d in |$\gamma$|-ray light curves using |$p_{\rm local,\gamma }=0.00953$| (2.59|$\sigma$|⁠) and N = 39. The results suggest a global confidence level of approximately 1.0|$\sigma$| for the time-scale of 1180.2 d in |$\gamma$|-ray light curves. This suggests that there is not significant periodicity variability in the |$\gamma$|-ray light curve, which is consistent with the findings reported by Yang et al. (2021).

To cross-validate the reliability of the LSP method results, we conducted a period analysis of the source using the REDFIT program (Schulz & Mudelsee 2002). This program, based on the first-order autoregressive (AR1) model (Schulz & Mudelsee 2002; Li et al. 2022a), can calculate the bias-corrected spectrum of a time-series and the significance of peak frequencies in the spectrum (Schulz & Mudelsee 2002). The maximum false alarm probability (FAP) level of the peaks provided by the REDFIT program is 99 per cent (2.5|$\sigma$|⁠) (Peñil et al. 2020; Li et al. 2022a). The bias-corrected spectrum, red noise spectrum, and the 90 per cent, 95 per cent, and 99 per cent FAP level lines are plotted in the bottom panels of Fig. 3, revealing three significant peaks at time-scales of 156.8, 276.4, and 1167.1 d for the optical waveband and a prominent peak at a time-scale of 1234.5 d for the |$\gamma$|-ray waveband, each with a higher FAP level than 99 per cent. The time-scale of 276.4 d in optical waveband from the REDFIT program is consistent with the results obtained using the LSP methods.

Moreover, the Jurkevich method (Jurkevich 1971) was also utilized to analyse the time-scales of the light curves plotted in Fig. 1. This method is proficient at detecting variability time-scales from irregularly spaced, non-sinusoidal modulated data samples. By testing a series of trial periods based on expected mean square deviation and phase folding techniques, the Jurkevich method divides the data sample into m groups to calculate the squared deviations |$V_{l}^{2}$| of the lth term (Jurkevich 1971),

where |$x_{i}$| is the individual observation and |${m}_{l}$| is the number of observations in the lth group. Then, the sum of the squared deviations in all groups is |$V^{2}_{m}=\sum _{l=1}^{m} V^{2}_{l}$|⁠. Smaller values of |$V_{m}^{2}$| for a trial period closer to the real one indicate a stronger period. To evaluate the significance of the period, a parameter |$f=\frac{1-V_{m}^{2}}{V_{m}^{2}}$|⁠, which was derived by Kidger, Takalo & Sillanpaa (1992), is utilized. A larger value of f corresponds to a higher significance level. Specifically, f > 0.5 suggests a strong period, while f < 0.25 indicates a weak or spurious one. The results of the Jurkevich method, illustrated in Fig. 6, display five possible periods in optical waveband light curve of Fig. 1 with time-scales of 277, 554, 821, 1223, and 1660 d. The f values associated with these periods are found to be 0.48, 0.64, 0.75, 1.08, and 0.82, respectively, indicating those periods are significant. Notably, an interesting discovery is the presence of an astronomical multiple frequency relation among these time-scales, where the time-scales of 554, 821, 1223, and 1660 d are approximately twice, three times, four times, and six times the time-scale of 277 d, respectively. It is worth noting that the time-scale of 277 d exhibits good agreement with the results obtained using the LSP method and REDFIT program. The right-hand panels of Fig. 6 reveal the presence of five distinct minimum values of |$V^{2}_{m}$| at time-scales of 1221, 1668, 2071, 2447, and 3353 d. However, the corresponding f values for these time-scales are less than 0.1, indicating that they represent spurious time-scales. Therefore, no significant periodicity can be observed in the |$\gamma$|-ray light curves, which is consistent with the results obtained using the LSP method.

4 CORRELATION ANALYSIS

In order to investigate the emission mechanism of TXS 1902+556, we conducted an analysis of the correlation between the optical and |$\gamma$|-ray waveband emissions using the DCF method (Edelson & Krolik 1988). The significant correlation between the variability in the optical and |$\gamma$|-ray wavebands implies that both emissions can be suitably explained by leptonic models, suggesting a common origin for both (Liao et al. 2014). Conversely, a weak or no correlation between lower-energy and higher-energy emissions suggests different origins for the emissions (external Compton, EC; Dermer & Schlickeiser 1993; Błażejowski et al. 2000; Dermer et al. 2009) or that the emissions are caused by hadronic processes (Böttcher 2005; Bonning et al. 2012; Liao et al. 2014). The DCF method is used to examine the correlation and time lag between two variable temporal series, utilizing algorithms provided by Edelson & Krolik (1988). The centroid time lags |$\tau$|⁠, estimated from the bumps closer to zero lag in the DCF profile, are utilized for determining. The significance level of the correlation between the optical and |$\gamma$|-ray waveband emission was estimated by the Monte Carle simulation as the foregoing analysis for the confidence level of LSP method.

The DCF result, shown in Fig. 7, includes the curves for significance levels 95, 99, and 99.7 per cent. The plot suggests a weak correlation between the variabilities in the optical and |$\gamma$|-ray wavebands, with a significance level surpassing 95 per cent but falling short of 99.7 per cent. This finding is consistent with previous reports in the literature (Liodakis et al. 2018). The results suggest that the emission of TXS 1902+556 cannot be adequately interpreted within the framework of the pure SSC model.

5 DISCUSSION AND CONCLUSIONS

The variability data of the blazar TXS 1902+556 in the optical and |$\gamma$|-ray bands were compiled in this study. We obtained long-term optical waveband variability data for TXS 1902+556 from the KAIT data base, spanning approximately 14.4 yr from 2009 May 11 to 2023 September 27. The |$\gamma$|-ray variability data covering 14.7 yr, from 2008 August 4 to 2023 April 23, was obtained from the Fermi-LAT data base. The emission variability properties of the optical and |$\gamma$|-ray wavebands were systematically studied using four different reliable methods, with particular emphasis on the systematic study of radiation variation in the optical band for the first time. The optical waveband light curve consistently revealed a significant quasi-periodicity with time-scales of |$P_{\rm O}=276.8\pm 6.1$| d using three distinct analysis methods, and the local significance level (⁠|$p_{\rm local}$|⁠) exceeded the |$3\sigma$|⁠. However, the |$\gamma$|-ray waveband did not exhibit significant periodicity, which does not support the findings reported by Peñil et al. (2020) and Rueda et al. (2022). Additionally, periodic Sun-gaps in the light curves may potentially contribute to the observed QPO phenomenon. To investigate whether the observed QPO behaviour in the optical waveband, characterized by time-scales of |$P_{\rm O}=276.8\pm 6.1$| d, can be attributed to specific intervals in the light curve, we performed frequency count analysis on intervals exceeding a duration of 150 d. The frequency count result is presented in Fig. 8, illustrating a Sun-gap of approximately |$233.5\pm 16.8$| d in the optical waveband light curves shown in Fig. 1. It is worth noting that this Sun-gap, with a time-scale of |$233.5\pm 16.8$| d, falls outside the 3|$\sigma$| interval of the QPO time-scale |$P_{\rm O}=276.8\pm 6.1$| d (ranging from 258.5 to 295.1 d). However, upon closer examination using standard statistical arguments, it is evident that these two quantities are consistent with each other within 2|$\sigma$| uncertainties. As demonstrated by the relationship |$233.5 + 2\sigma _{1}$| being approximately equal to |$276.8 - 2\sigma _{2}$|⁠, where |$\sigma _{1}= 16.8$| and |$\sigma _{2}= 6.1$|⁠. Therefore, the QPO time-scale is consistent with the Sun-gaps in the optical light curve within 2|$\sigma$| uncertainties.

Furthermore, we employed the DCF method to investigate the correlation between variability in the optical and |$\gamma$|-ray wavebands. Our findings revealed a weak correlation between the two wavebands, with a level of significance exceeding |$2\sigma$|⁠, which contradicts the results reported by Liodakis et al. (2018). Investigating the correlations in variability is pivotal for advancing our understanding of the emission mechanisms exhibited by blazars. The leptonic model provides a robust explanation for the significant correlation observed between optical and |$\gamma$|-ray variations (Giommi et al. 1999; Tagliaferri et al. 2003; Zheng et al. 2013; Liao et al. 2014; Zheng, Kang & Li 2014). In contrast, the hadronic processes do not generate such a significant correlation (Liao et al. 2014). According to the leptonic model, low-energy radiation is generated by synchrotron radiation, while high-energy radiation is produced through inverse Compton (IC) scattering (Tagliaferri et al. 2003; Li et al. 2016). The process is classified into SSC and EC processes, based on the origin of the seed photons. The SSC process involves source photons derived from synchrotron photons emitted at low energies, resulting in a significant correlation between the variability in low and high energy emissions (Maraschi, Ghisellini & Celotti 1992; Liao et al. 2014). On the other hand, the EC process involves seed photons originating from radiation fields external to the jet, such as the accretion disc, hot dust, broad-line region, and so on (Malmrose et al. 2011; Liao et al. 2014). Our results indicate a weak correlation between variabilities in the optical and |$\gamma$|-ray wavebands for TXS 1902+556. This finding implies that the emission process in the |$\gamma$|-ray waveband is more complex than a single synchrotron self-Compton (SSC) process, suggesting that the emission of TXS 1902+556 may arise from a combination of SSC and EC processes or leptonic–hadronic hybrid processes similar to those observed in other sources such as 1ES 1959+650 (Li et al. 2022b) and S5 0716+714 (Liao et al. 2014).

The QPO behaviour of blazars can be elucidated by various theoretical models, including the SMBBH system (Sillanpaa et al. 1988; Xie et al. 2005, 2008; Graham et al. 2015; Wang, Yin & Xiang 2017; Gong et al. 2022), instabilities in accretion flows (Karouzos et al. 2012; McKinney, Tchekhovskoy & Blandford 2012; Piner & Edwards 2014), and geometrical effects (Rieger 2004; Li et al. 2015; Zhang et al. 2017; Gong et al. 2022, 2023). The inter-day QPO behaviour can be explained within the framework of the shock-in-jet process and the instabilities in accretion flows (e.g. Mangalam & Wiita 1993; Hong et al. 2018; Li et al. 2024). Moreover, the presence of a SMBBH system and geometric effects can provide a plausible explanation for both short-term and long-term QPO (e.g. Sillanpaa et al. 1988, 1996a; Xie et al. 2005; Li et al. 2009; Ackermann et al. 2015; Otero-Santos et al. 2023). The overall radiation of the BL Lac object is primarily dominated by the boosted jet emission, while its main power source originates from accretion (Kushwaha et al. 2020). The periodic variability in emissions observed in blazars can be attributed to the orbital motion of the secondary black hole within the SMBBH system, providing a plausible explanation for various quasi-periodic candidates. These candidates encompass an approximately 12-yr period in OJ 287 (Sillanpaa et al. 1988, 1996a), a period of |$1.84\pm 0.1$| yr in PKS 1510-089 (Xie et al. 2005, 2008), a roughly 3-yr cycle in 3C 66A (Otero-Santos et al. 2020), and a period of |$2.18\pm 0.08$| yr in PG 1553+113 (Ackermann et al. 2015), etc. Several scenarios have been proposed to elucidate QPO, including periodic disturbances in accretion, the double jet model, gravitational lensing effects, and so on (Sillanpaa et al. 1996a; Villata & Raiteri 1999; Otero-Santos et al. 2023). Under this framework, if the primary black hole is surrounded by an accretion disc, the orbital motion induces quasi-periodic tidal perturbations in the accretion disc, leading to periodic gas transfer towards the primary black hole and resulting in variability in emission periodicity. Furthermore, the accretion rate also exhibits fluctuations during the traversal of the primary’s accretion disc by the secondary black hole, leading to recurrent flares (Lehto & Valtonen 1996). Consequently, QPO can be observed due to the orbital motion of the secondary black hole around its primary counterpart (Sillanpaa et al. 1996a). This scenario can also plausibly explain the double-peak structure observed during the outburst phase. The double peak structure can also be explained by the double jet model in the SMBBH system, where each black hole possesses a relativistic jet (Huang et al. 2021). The substantial contribution of each jet is crucial for the overall impact, and the QPO can be attributed to changes in the orientation of both jets. Additionally, the gravitational lensing model offers a plausible interpretation for the quasi-periodic emission (Sillanpaa et al. 1996a; Roy et al. 2022a). In this model, the emission is amplified as it traverses a massive object. Within the SMBBH system, the observed radiation flux is enhanced when the secondary black hole is situated between the primary black hole and the observer. As a result of the orbital motion of the secondary black hole around the primary, the emission would exhibit QPO. Moreover, the emission variability induced by the gravitational lensing effect is expected to possess remarkable symmetry. This model can also offer an explanation for the observed double-peaked structure of emission variability (Roy et al. 2022a).

Assuming that the |$P=276.8\pm 6.1$| d period is attributed to the orbital motion of the SMBBH system, considering the |$P_{\rm rf}=P/(1+z)$| relation with |$z=0.727$|⁠, the periodicity in the rest reference frame amounts to approximately 0.44 yr. The mass of the primary black hole in the SMBBH system can be estimated utilizing the subsequent formula (Begelman, Blandford & Rees 1980; Ostorero et al. 2004; Gong et al. 2022; Li et al. 2022a),

where |$P_{\rm rf,yr}$| is the period in the rest reference frame in years, while R corresponds to the mass ratio between the primary and secondary black holes. Considering |$R=10$| and |$P_{\rm rf}=0.44$| yr, the estimated mass of the primary black hole is approximately |$M\sim 1.06\times 10^{6}\,{\rm M}_{\odot }$|⁠, significantly lower than the reported black hole masses of |$2.366\times 10^{8}\,{\rm M}_{\odot}$| as documented by Pei et al. (2022). Therefore, it may not be appropriate to assume that the variability in periodicity originates solely from the orbital motion of the SMBBH system for TXS 1902+556.

Moreover, the QPO of blazars may be associated with the orbital time-scale of a hotspot, a blob, a flare, or other oscillatory phenomena occurring within the innermost region of the rotating accretion disc (Gupta, Srivastava & Wiita 2009; Gupta et al. 2019). Under the model, the plasma, originating from oscillating accretion discs, will quasi-periodically inject into jets and induce periodic changes in radiant flux (Gupta et al. 2009, 2019; Gong et al. 2022). The mass of supermassive black hole (SMBH) can be estimated by the following formula (Gupta et al. 2009, 2019; Gong et al. 2022)

where |$P_{s}$|⁠, |$\delta$|⁠, r, a, and z denote the observed period in seconds, the Doppler factor, the radius of this source zone in units of |$GM/c^{2}$|⁠, the SMBH spin parameter and the redshift, respectively. The observed period, Doppler factor, and redshift for TXS 1902+556 are |$P=276.8$| d, |$\delta =5.11$| (Fan et al. 2014), and |$z=0.727$|⁠, respectively. The radius r and the spin parameter a for a Schwarzschild black hole are given by |$r=6$| and |$a=0$|⁠, respectively, while for a maximal Kerr black hole they are determined as |$r=1.2$| and |$a=0.9982$| (Gupta et al. 2009, 2019; Gong et al. 2022). Consequently, the mass of the SMBH is approximately |${\sim}1.56\times 10^{11}\,{\rm M}_{\odot }$| for the Schwarzschild case, and |${\sim}1.01\times 10^{12}\,{\rm M}_{\odot }$| for the Kerr scenario. The estimated mass of SMBH, as calculated by equation (10), significantly surpasses the reported black hole masses of |$2.366\times 10^{8}\,{\rm M}_{\odot }$| documented by Pei et al. (2022). Therefore, it may not be appropriate to assume that the QPO originates from perturbations on the innermost stable circular orbit for TXS 1902+556.

As previously mentioned, the QPO of blazars can be reasonably interpreted within the framework of the geometrical model and has been extensively discussed in numerous scholarly investigations (e.g. Villata & Raiteri 1999; Rieger 2004; Li et al. 2009, 2015; Mohan & Mangalam 2015; Otero-Santos et al. 2020; Gong et al. 2022, 2023; Gao et al. 2023). The geometrical model proposes that the observed emission flux is enhanced due to beaming effects, in accordance with a power-law relationship of |$F_{\nu }\propto \delta ^{3}$|⁠. Therefore, the variability in emissions is associated with changes in the Doppler factor |$\delta$|⁠, which is a function of the viewing angle |$\theta$|⁠,

where |$\beta =v/c$| and |$\Gamma =(1-\beta ^{2})^{-1}$| represent the normalized velocity and bulk Lorentz factor, respectively. In the geometrical scenario, the emission blob moves along helical trajectories within the jet, causing the viewing angle to change as it moves. The helical trajectories may be attributed to the bending jet, wiggling jet, and the helical magnetic fields in the jet, leading to periodic variations in the viewing angle (Li et al. 2015). Thus, the motion of emitting material along helical trajectories would lead to observable periodic emissions. The relationship between the viewing angle |$\theta$| and the observed period |$P_{\rm obs}$| is described by the following formula:

where |$\phi$| represents the pitch angle between the blob motion and the jet axis, and |$\psi$| represents the inclination angle of the jet to the line of sight. The periodicity in the rest frame of the blob is estimated by the following equation:

For blazar, the typical values of |$\psi$| and |$\phi$| is |$\psi =2^{\circ }$| and |$\phi =5^{\circ }$|⁠, respectively (Banerjee et al. 2023). The value of |$\beta$| can be estimated from the Doppler factor |$\delta$| and bulk Lorentz factor |$\Gamma$|⁠. Given a value of |$\delta$|⁠, a lower limit to the Lorentz factor is given by the following expression (Urry & Padovani 1995),

This expression is valid only for |$\delta > 1$|⁠. For TXS 1902+556, a lower limit of |$\delta =5.11$| was reported by Fan et al. (2014). Then, the lower limit of the Lorentz factor is |$\gamma \ge 2.65$|⁠. Therefore, the value of |$\beta$| is |$\beta =\sqrt{1-\frac{1}{\gamma ^{2}}}\ge 0.93$|⁠. The rest frame periodicity becomes |$P_{\rm rf}\sim 10.24$| yr.

Regarding the periodic variability induced by the helical trajectories of the jet, three potential driving mechanisms were proposed by Rieger (2004): intrinsic rotation of the jet, orbital motion in a SMBBH system, and jet precession. The time-scales of |$P_{\rm obs}$| associated with these driving mechanisms are ≤10 d, ≥10 d, and ≤1 yr, respectively. The period of |$276.8\pm 6.1$| d observed in the light curves of TXS 1902+556 can be reasonably explained by the helical motion of the jet driven by orbital motion in SMBBH system. By utilizing equation (9), the estimated mass of the primary black hole is approximately |$M\sim 1.65\times 10^{8}\,{\rm M}_{\odot }$|⁠, given |$R=10$| and |$P_{\rm rf}=10.24$| yr, which aligns with the reported black hole masses of |$2.366\times 10^{8}\,{\rm M}_{\odot }$| by Pei et al. (2022). Therefore, the QPO of TXS 1902+556 with a time-scale of 276.8 d is most likely attributed to the helical motion of the jet driven by the orbital motion in SMBBH system.

The QPO behaviour identified in our study is limited to the optical band for several reasons. First, the low flux and large relative errors associated with |$\gamma$|-ray measurements may account for this limitation. It is likely that potential QPO in |$\gamma$|-ray are obscured by measurement noise. Furthermore, |$\gamma$|-ray photons travelling over cosmological distances experience attenuation by the extragalactic background light (EBL), leading to the loss of QPO information. The redshift of TXS 1902+556 is |$z=0.727$|⁠, resulting in significant absorption of the EBL, attenuated by |${\rm e}^{-\tau (E_{\gamma})}$|⁠, where |$\tau$| is the optical depth. For Fermi-LAT observations at this redshift, |$E_{\gamma} > 140\, \text{GeV}$| corresponds to |$\tau > 1$| (Finke, Razzaque & Dermer 2010). In addition, the |$\gamma$|-ray emission from BL Lac objects includes an additional spectral component beyond the standard SSC emission. While the QPO in the optical band may originate from a SMBBH system, the secondary emission, such as that from the hadronic model (e.g. proton–photon interactions (Dermer, Murase & Takami 2012; Li et al. 2022b)), could be less affected by helical motion or independent of the primary emission region. This could result in a weakened QPO signature in the observed |$\gamma$|-ray emission, even in the presence of an SMBBH system. Our analysis also shows a weak correlation between the optical and |$\gamma$|-ray bands, suggesting the presence of an additional |$\gamma$|-ray component for TXS 1902+556.

ACKNOWLEDGEMENTS

We thank the anonymous referee for the valuable comments and suggestions. This study utilizes data obtained from KAIT, the authors gratefully acknowledge Weikang Zheng for providing the KAIT data. Additionally, this research employs data observed by the Fermi project, the authors gratefully acknowledge the Fermi team for their invaluable support. This work was supported by the National Natural Science Foundation of China (12063005, 12063006, 12373018), the Young and Middle-aged Discipline and Technology Leaders Reserve Talents in Yunnan Province (202405AC350114), the Program for Innovative Research Team (in Science and Technology) in University of Yunnan Province (IRTSTYN). The authors gratefully acknowledge the computing support provided by the JRT Science Data Center at Yuxi Normal University and the author (QLH) gratefully acknowledges the financial supports from the Hundred Talents Program of Yuxi (grant 2019).

DATA AVAILABILITY

The data utilized in this research were provided by the KAIT and Fermi project and will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES