2. Multi-wavelength observations

3. Periodicity search

We adopted various methodologies in search of a potential periodic signal in the $\gamma$ -ray and optical light curves of blazar ON 246. Figure 1 illustrates the 10-day binned $\gamma$ -ray light curve along with optimal Bayesian Block representation (top panel) and an ASAS-SN optical light curve in the bottom panel.

The utilised methodologies include the Lomb-Scargle periodogram (LSP), Weighted Wavelet Z-Transform (WWZ), and a first-order autoregressive model $\left(AR(1)\right)$ in this QPO investigation. A detailed description and the observed findings from all the utilised methods are given in the following Section 3.1, 3.2.

3.1 Lomb-Scargle Periodogram

The Lomb-Scargle periodogram (LSP) (Lomb Reference Lomb1976; Scargle Reference Scargle1979) is one of the most widely used approaches in the literature to identify any potential periodic signal in a time series. In which a sinusoidal function fits the time series using the least square method. This approach is capable of handling the non-uniform sampling in the time series data efficiently by reducing the impact of noise and gaps and providing a precise measurement of the identified periodicity. In this study, we used the GLSP package to compute the Lomb-Scargle (LS) power. The expression of LS power is given as VanderPlas (2018):

(1) \begin{equation}\begin{split}P_{LS}(f) = \frac{1}{2} \bigg[ &\frac{\left(\sum_{i=1}^{N} x_i \cos(2\pi f (t_i - \tau))\right)^2}{\sum_{i=1}^{N} \cos^2(2\pi f (t_i - \tau))} \\[4pt] &+ \frac{\left(\sum_{i=1}^N x_i \sin(2\pi f (t_i - \tau))\right)^2}{\sum_{i=1}^N \sin^2(2\pi f (t_i - \tau))} \bigg]\end{split}\end{equation}

where $\tau$ is specified for each f to ensure time-shift invariance:

(2) \begin{equation} \tau = \frac{1}{4 \pi f} \tan^{-1} \left( \frac{\sum_{i=1}^N \sin\left( 4 \pi f t_i \right)}{\sum_{i=1}^N \cos\left( 4 \pi f t_i \right)} \right)\end{equation}

where we selected the minimum frequency $\left( f_{\min} \right)$ and maximum frequency $\left( f_{\max} \right)$ in temporal frequency range as 1/T and 1/2 $\Delta T$ , respectively, and here T and $\Delta T$ represent the total observation time frame and the time difference between two consecutive points in the light curve, respectively.

The LSP analysis reveals prominent peaks at frequencies of $\sim$ 0.00134 day⁻¹ ( $746\pm68$ days) in the $\gamma$ -ray LSP (see Figure 2) and $\sim$ 0.00132 day⁻¹ ( $757\pm106$ days) in the optical LSP (see Figure 3). Both peaks have a local significance level exceeding 99.73 $\%$ . The uncertainty on the observed period is estimated by fitting a Gaussian function to the dominant LSP peak, and the obtained half-width and half maxima (HWHM) value is used as an uncertainty on period. The distribution of LS power as a function of frequency is given in Figure 2.

Figure 3. The detected QPO signals in the optical emissions from ON 246. The top panel shows the LSP with a dominant peak at $\sim$ 0.00132 day $^{-1}$ has a local significance level exceeding 99.73 $\%$ . The bottom panels display the wavelet map (bottom left panel) and avg. wavelet power at frequency of $\sim0.00131$ day $^{-1}$ with a significance level greater than 99.73 $\%$ .

3.2 Weighted Wavelet Z-transform

In contrast to the LSP approach, the Weighted Wavelet Z-transform (WWZ) (Foster Reference Foster1996) emerges as a powerful, robust, and widely used method in astronomical studies to identify any potential periodic pattern in irregularly sampled light curves. The WWZ method incorporates wavelet analysis, enhancing the LSP’s capabilities by providing better localisation of periodic signals in both temporal and spectral space. In studying the evolution of a periodic signal over time, this approach emerges as a powerful tool, enabling us to identify and characterise the nature of a periodic signal.

In this study, we adopted the abbreviated Morlet kernel that has the following functional form (Grossmann & Morlet Reference Grossmann and Morlet1984):

(3) \begin{equation} f[\omega (t - \tau)] = \exp[i \omega (t - \tau) - c \omega^2 (t - \tau)^2]\end{equation}

and the corresponding WWZ map is given by,

(4) \begin{equation} W[\omega, \tau \;:\; x(t)] = \omega^{1/2} \int x(t)f^* [\omega(t - \tau)] dt\end{equation}

where $f^*$ is the complex conjugate of the wavelet kernel f, $\omega$ is the frequency, and $\tau$ is the time-shift. This kernel acts as a windowed DFT, where the size of the window is determined by both the parameters $\omega$ and a constant c. The resulting WWZ map offers a notable advantage; it not only identifies dominant periodicities but also provides insights into their duration over time.

Figure 4. Analysis of the light curves, left panel represent the REDFIT curve of $\gamma$ -ray emissions and right panel exhibit the REDFIT curve of optical emissions, using the AR(1) process with the REDFIT software. The red noise-corrected power spectrum (black) is presented alongside theoretical (blue) and average AR(1) (cyan) spectra. The significance levels of 99 $\%$ , 95 $\%$ , and 90 $\%$ are indicated in red, green, and brown, respectively.

In this study, we used publicly available python codeFootnote ^d (Aydin Reference Aydin2017) to generate the WWZ map. The observed power concentration is located around 0.00132 day $^{-1}$ ( $757\pm80$ days) in the WWZ map utilising $\gamma$ -ray emissions (see Figure 2) and at $\sim$ 0.00131 day $^{-1}$ ( $763\pm102$ days) in optical WWZ map (see Figure 3). In both cases, the observed local significance surpasses 99.73 $\%$ . The uncertainty on the period was estimated using the method as described in Section 3.1.

3.3 REDFIT

The light curves of AGNs are typically unevenly sampled, of finite duration, and predominantly influenced by red noise, which arises from stochastic processes occurring in the accretion disk or jet. Red noise spectra are characteristic of autoregressive processes, where current activity is related to past behaviour. The emissions from AGNs are effectively modeled using a first-order autoregressive (AR1) process. To model such behaviour, the software programme REDFIT, developed by Schulz & Mudelsee (Reference Schulz and Mudelsee2002), is specifically designed to analyse the stochastic nature of AGNs dominated by red noise. This software fits the light curve AR(1) process, where the current emission ( $r_t$ ) depends linearly on the previous emission ( $r_{t - 1}$ ) and a random error term ( $\epsilon_t$ ). The AR(1) process is defined as:

(5) \begin{equation} r(t_i)=A_i r(t_{i-1}) + \epsilon(t_i)\end{equation}

where $r(t_i)$ is the flux value at time $t_i$ and $A_i = \exp\left( \left[ \frac{t_{i-1} - t_i}{\tau}\right] \right) \in [0,1]$ , A is the average autocorrelation coefficient computed from mean of the sampling intervals, $\tau$ is the time-scale of autoregressive process, and $\epsilon$ is a Gaussian-distributed random variable with zero mean and variance of unit. The power spectrum corresponding to the AR(1) process is given by

(6) \begin{equation} G_{rr}(f_i) = G_0 \frac{1 - A^2}{1 - 2 A \cos\left( \frac{\pi f_i}{f_{Nyq}} \right) + A^2}\end{equation}

where $G_0$ is the average spectral amplitude, $f_i$ are the frequencies, and $f_{Nyq}$ is representing the Nyquist frequency.

In our study, we used the publicly available REDFITFootnote ^e code to analysis the light curve. In this method, the Nyquist frequency is defined as $f_{Nyq} = H_{fac}/ (2 \Delta t)$ , where the factor $H_{fac}$ is introduced to prevent the noisy high-frequency end of the spectrum from influencing the fit, as described by equation (7). The REDFIT analysis detected prominent peaks at frequencies of $\sim$ 0.00128 day⁻¹ ( $781\pm160$ days) in the $\gamma$ -ray light curve (see the left panel of Figure 4) and $\sim$ 0.00132 day⁻¹ ( $757\pm160$ days) in the optical light curve (see the right panel of Figure 4). The uncertainties in the periods were estimated using the methodology described in Section 3.1.

Figure 5. The figure displays the posterior probability distributions of the DRW model parameters, obtained from the $\gamma$ -ray light curve (left panel) and the ASAS-SN light curve (right panel).

4. Gaussian process modelling

The observed variability in AGN is inherently stochastic. The AGN light curves can be well described by the stochastic processes, also known as Continuous Time Autoregressive Moving Average [CARMA(p, q)] processes (Kelly et al. Reference Kelly, Becker, Sobolewska, Siemiginowska and Uttley2014), defined as the solutions to the stochastic differential equation:

(7) \begin{equation}\begin{split}\frac{d^p y(t)}{dt^p} + \alpha_{p-1}\frac{d^{p-1}y(t)}{dt^{p-1}}+\ldots+\alpha_0 y(t) =\\\beta_q \frac{d^q \epsilon(t)}{dt^q}+\beta_{q-1}\frac{d^{q-1}\epsilon(t)}{dt^{q-1}}+\ldots+\beta_0 \epsilon(t),\end{split}\end{equation}

where y(t) represents a time series, $\epsilon (t)$ is a continuous time white noise process, and $\alpha_*$ and $\beta_*$ are the coefficients of autoregressive (AR) and moving average (MA) models, respectively. Here, p and q are the order parameters of AR and MA models, respectively.

The simplest model is a continuous autoregressive [CAR(1)] model, also known as the Ornstein-Uhlenbeck process. It is a popular red noise model (Kelly, Bechtold, & Siemiginowska Reference Kelly, Bechtold and Siemiginowska2009; Kozłowski et al. Reference Kozłowski2009; MacLeod et al. Reference MacLeod2012; Ruan et al. Reference Ruan2012; Zu et al. Reference Zu, Kochanek and Udalski2013; Moreno et al. Reference Moreno, Vogeley, Richards and Yu2019; Burke et al. Reference Burke2021; Zhang, Yan, & Zhang Reference Zhang, Yan and Zhang2022, Reference Zhang, Yan and Zhang2023; Sharma et al. Reference Sharma, Kamaram, Prince, Khatoon and Bose2024; Zhang, Yang, & Dai Reference Zhang, Yang and Dai2024; Sharma, Prince, & Bose Reference Sharma, Prince and Bose2024), usually referred to as the damped random walk (DRW) model, described by the following differential equation:

(8) \begin{equation} \left[ \frac{d}{dt} + \frac{1}{\tau_{DRW}} \right] y(t) = \sigma_{DRW} \epsilon(t)\end{equation}

where $\tau_{DRW}$ and $\sigma_{DRW}$ are the characteristic damping time-scale and amplitude of the DRW process, respectively. The mathematical form of the covariance function of the DRW model is defined as

(9) \begin{equation} k(t_{nm}) = a \cdot \exp(\!-\!t_{nm} \, c),\end{equation}

where $t_{nm} = | t_n -t_m|$ denotes the time lag between measurements m and n, with $a = 2 \sigma_{DRW}^2$ and $c = \frac{1}{\tau_{DRW}}$ . The power spectral density (PSD) of the DRW model is defined as:

(10) \begin{equation} S(\omega) = \sqrt{\frac{2}{\pi}} \frac{a}{c} \frac{1}{1 + (\frac{\omega}{c})^2}\end{equation}

The DRW PSD has a form of Broken Power Law (BPL), where the broken frequency $f_b$ corresponds to the characteristic damping timescale $\tau_{DRW} = \frac{1}{2\pi f_b}$ .

In the best-fit parameters estimation of the DRW model for both light curves, we employed the Markov chain Monte Carlo (MCMC) algorithm provided by the emceeFootnote ^f package (Foreman-Mackey et al. Reference Foreman-Mackey, Hogg, Lang and Goodman2013). For the modelling, we employed the EzTaoFootnote ^g package, which is built on top of celerite.Footnote ^h In this study, we generated the distributions of the posterior parameters by running 10 000 steps as burn-in and 20 000 as burn-out, as illustrated in Figure 5. The resulting best-fitted $\gamma$ -ray and optical light curves are presented in Figures 6 and 7. The autocorrelation functions of the standardised residuals and their squares indicate that the DRW model effectively captures the characteristic variability in both light curves. The corresponding DRW power spectral densities are displayed in Figure 8, where the shaded regions represent areas of unreliability due to the finite duration and sampling cadence of the light curves.

Figure 6. The celerite modelling with DRW model was performed using the 10-day binned $\gamma$ -ray light curve of blazar ON 246 over 4 000 days from the time stamp MJD 55932. In this figure, the top panel depicts the $\gamma$ -ray flux points with their uncertainties, along with the best-fitting profile of the DRW model in blue, including the 1 $\sigma$ confidence interval. The bottom panels represent the autocorrelation functions (ACFs) of the standardised residuals (bottom left) and the squared of standardised residuals (bottom right), respectively, along with 95 $\%$ confidence intervals of the white noise.

Figure 7. This figure demonstrates the modelling of the ASAS-SN light curve of ON 246 with the DRW model. The top panel shows the best-fitting profile of the DRW modelling along with 1 $\sigma$ confidence interval. The bottom panels represent the autocorrelation functions as described in Figure 6.

Figure 8. This figure presents the DRW PSDs obtained from $\gamma$ -ray and ASAS-SN observations, along with their 1 $\sigma$ confidence intervals. The two shaded regions highlight biased regions due to observational limitations, such as finite length and the mean cadence of the light curve. The regions with orange star symbols represent the invalid areas in $\gamma$ -ray PSD, while the blue hatch line regions in PSD indicate limitations imposed by the ASAS-SN light curve’s duration and mean cadence.

5. Probe the significance

AGN emissions exhibit stochastic variability and are well characterised by red noise. The combination of red noise, characteristic nature, and uneven sampling in the light curve can lead to spurious peaks in the periodogram. Therefore, it is crucial to evaluate the significance of any periodic signals detected in the light curves. In the estimation of significance, we employed two different approaches.

Figure 9. The Lomb-Scargle periodograms of the $\gamma$ -ray (left panel) and ASAS-SN (right panel) observations of the blazar ON 246 are shown. The figure presents the LSPs of the original light curves (black) and spectral windows (pink). The significance levels (blue) of the dominant observed peaks in both periodograms are estimated using DRW synthetic light curves. The detected peaks at 0.00105 day $^{-1}$ ( $\sim$ 950 days) and 0.00105 day $^{-1}$ ( $\sim$ 950 days) in $\gamma$ -ray and ASAS-SN observations exceed 3 $\sigma$ and 4 $\sigma$ significance levels, respectively. The shaded regions in both figures represent the invalid areas estimated using the mean cadence and baseline of the light curves. Logarithmic versions of the periodograms are provided in the sub-figures.

The periodogram is usually represented as a power spectral density (PSD) of a form $P(\nu) \sim A \nu^{-\beta}$ , where the temporal frequency is represented by $\nu$ and $\beta \gt 0$ represents the spectral slope. In the first approach, to calculate the statistical significance of the periodic feature, we employed the approach developed by Emmanoulopoulos, McHardy, & Papadakis (Reference Emmanoulopoulos, McHardy and Papadakis2013). The methodology involved modelling the observed PSD using a power-law model. For that, we employed DELightcurveSimulation Footnote ⁱ code, which involved randomising the amplitude and phase of the Fourier components, each mimicking the characteristics of the original data, including the best-fitting power-law model slope and similar properties in terms of flux distribution. We performed a Monte Carlo simulation, generating $5\times10^4$ synthetic light curves for each case. The simulations were based on the best-fitting power-law model slopes ( $\beta$ = 0.48 $\pm$ 0.16 for $\gamma$ -rays light curve and $\beta$ = 1 $\pm$ 0.16 for ASAS-SN light curve) utilising their flux distribution properties, respectively. Each simulated light curve for both cases mimics the underlying properties of the original light curves. After that, we generated the Lomb-Scargle periodogram of each simulated light curve as we did for both original light curves. To estimate the significance level of the dominant peak in original LSPs, we calculated the 84th, 97.5th, 99.85th, and 99.995th percentiles of the 50 000 simulated LSP for each frequency value, which correspond to the 1 $\sigma$ , 2 $\sigma$ , 3 $\sigma$ , and 4 $\sigma$ significance levels. In this first approach, the significance level of the dominant peak in $\gamma$ -ray LSP surpasses the 3 $\sigma$ and just touches the 3 $\sigma$ significance level in ASAS-SN LSP, as shown in Figures 2 and 3.

In addition to modelling the PSDs with a simple power-law, we also applied a bending power-law model, as detailed in Appendix A. The prominent peak observed in the $\gamma$ -ray LSP at $\sim 0.00134 \ \mathrm{day^{-1}}$ exhibits a significance of approximately $3\sigma$ , while, the dominant peak at $\sim 0.00132 \ \mathrm{day^{-1}}$ in the optical LSP surpasses a significance threshold of $3.48\sigma$ .

In the second approach, considering that AGN variability is stochastic and well characterised by a red noise process, we employed the Damped Random Walk (DRW) model to describe both light curves and determine the optimal model parameters, including the amplitude and damping timescale. Using the EzTao package, we simulated 10 000 synthetic light curves with a sampling rate consistent with the real observations. After generating these synthetic light curves, we computed Lomb-Scargle periodograms (LSPs) for each, following the same procedure as for the original $\gamma$ -ray and ASAS-SN light curves. The significance levels were then estimated using the same methodology as described earlier. From this analysis, the peaks at 0.00134 day $^{-1}$ in $\gamma$ -ray LSP, while the peak at 0.00132 day $^{-1}$ in the optical LSP surpasses the 4 $\sigma$ threshold. Additionally, we calculated the spectral window periodogram by constructing a light curve with a total number of time stamps ten times larger than the original within the observed temporal frame. In this synthetic light curve, the time stamps matching the original observations were assigned a value of one, while all others were set to zero. Further, we applied the LSP method to generate the periodogram, as shown in pink in Figure 9. The QPO findings are summarised in Table 1.

Table 1. Summary of the quasi-periodic signal detection using three different methodologies. Column (1) lists the 4FGL name of the blazar ON 246, while Column (2) specifies the observational waveband. Columns (3) and (4) present the QPO frequencies obtained from the Lomb-Scargle periodogram and Weighted Wavelet Z-transform methods, along with their uncertainties. The local significance of each detected QPO signal is provided in parentheses next to the corresponding frequency value. Column (5) presents the estimated local significance level of the QPO in LSP based on DRW-modeled mock light curves, and Column (6) provides the QPO frequency and significance level derived from the REDFIT analysis.

Note: The QPO peaks were fitted using a Gaussian function, and the uncertainties correspond to the half-width at half-maximum (HWHM).

6. Gamma-ray/optical cross correlations

We investigate the possible time lags between the $\gamma$ -ray and optical light curves of ON 246 utilising the interpolated cross-correlation function (ICCF: Peterson et al. Reference Peterson, Wanders, Horne, Collier, Alexander, Kaspi and Maoz1998, Reference Peterson2004), which is one of the commonly used methods in the time-series analysis of AGNs. As we know, AGN light curves generally are not regularly sampled in time but are discretely sampled N times at time stamps $t_i$ , where $t_{i+1} - t_i = \Delta t$ for all values $1 \le i \le N-1$ . ICCF method emerges as a powerful technique to estimate the time-leg between two-time series. The method uses the linear interpolation method to deal with unevenly sampled light curves and calculate the cross-correlation coefficient as a function of the time lag for two light curves:

(11) \begin{equation} F_{CCF}(\tau)=\frac{1}{N} \sum_{i=1}^N \frac{\left[ L(t_i) - \bar{L} \right] \left[ C(t_i - \tau) - \bar{C} \right]}{\sigma_L \sigma_C}\end{equation}

where N is the number of data points in the light curves L and C. Each light curve has a corresponding mean value ( $\bar{L}$ and $\bar{C}$ ) and uncertainty ( $\sigma_L$ and $\sigma_C$ ).

The ICCF is evaluated for a time lag ( $\tau$ ) in a range [-1 000, 1 000] days with searching step $\Delta \tau$ , which should be smaller than the median sampling time of the light curves. We adopted $\Delta \tau$ = 7 days and used the public PYTHON version of the ICCF, PYCCF (Sun, Grier, & Peterson Reference Sun, Grier and Peterson2018) in this study. As applying the ICCF to the light curves, a strong peak is obtained, and its corresponding centroid is calculated using the ICCF for the time lags around the peak. We adopted the centroid of the CCF ( $\tau_{cent}$ ) using only time lags with $r\gt 0.8 r_{\max}$ , where $r_{\max}$ is the peak value of the CCF. The 1 $\sigma$ confidence on the time lag is estimated using a model-independent Monte Carlo method. We found a lag of $0.736_{-2.73}^{+3.13}$ day with the $\gamma$ -ray leading the optical emissions. In addition to cross-correlation centroid distribution (CCCD), we also estimated the cross-correlation peak distribution (CCPD; see the right panel of Figure 10). To access the significance, we also calculated the ICCF between the $\gamma$ -ray and DRW mock optical light curves (see Section 5); the observed findings are shown in the left panel of Figure 10. After simulating the 10 000 mock ICCFs between $\gamma$ -ray and optical mock light curves, we calculated the significance level of the ICCF peak observed in correlation analysis between the observed $\gamma$ -ray and ASAS-SN light curves. In the left panel of Figure 10, a red dashed curve represents the 99.9999 $\%$ , corresponding to 4 $\sigma$ , significance level. The observed finding indicates a significant correlation with almost zero-day lag between $\gamma$ -ray and optical emissions, suggesting a common origin of them (Cohen et al. Reference Cohen, Romani, Filippenko, Bradley Cenko, Lott, Zheng and Li2014).

Figure 10. Cross-correlation analysis between $\gamma$ -ray and optical flux variations. The left panel shows ICCF (solid black curve) with a significance level (dash red curve) of 4 $\sigma$ . The cross-correlation centroid distribution (CCCD) in orange and cross-correlation peak distribution (CCPD) in blue are given in the right panel, where vertical dashes in orange and blue represent the median of CCCD and CCPD, respectively.

7. Black hole mass estimation

The mass of the central black hole in an AGN is one of the most crucial parameters for understanding its central engine. Several methods have been developed to estimate black hole mass, including stellar and gas kinematics and reverberation mapping (Gupta, Srivastava, & Wiita Reference Gupta, Srivastava and Wiita2008). The stellar and gas kinematics approach requires high spatial resolution spectroscopic observations to resolve the gravitational influence of the black hole. Reverberation mapping relies on detecting high-ionisation emission lines from regions close to the black hole. In addition to the black hole mass estimation methods mentioned above, the variability timescales can also serve as a useful tool for determining the mass of an AGN’s black hole. The timescales of high-amplitude variations may be linked to the black hole mass, offering valuable insights into the central engine of AGNs. The shortest variability timescale provides a constraint on the size of the emitting region in blazars. In relativistic jets, the emitting region is often simplified as a ‘blob’ with a characteristic size D. This size can be constrained by the relation:

(12) \begin{equation} D \lesssim \frac{\delta \Delta t_{{obs}} c}{1+z}\end{equation}

where $\delta$ is the Doppler factor of the relativistic jet, $\Delta t_{obs}$ is the observed variability timescale, z is the redshift of the blazar, and c represents the speed of light.

Several studies (Hartman Reference Hartman, Richard Miller, Webb and Noble1996; Ghisellini & Madau Reference Ghisellini and Madau1996; Xie et al. Reference Xie, Bai, Zhang and Fan1998; Yang & Fan Reference Yang and Fan2010) have suggested that $\gamma$ -ray emissions originate from a region located a few hundred Schwarzschild radii ( $r_g\equiv GM_* /c^2$ ), where G and $M_*$ are the gravitational constant and black hole mass, respectively. The emission region is typically constrained to a size of $r \lt 200 r_g $ . Yang & Fan (Reference Yang and Fan2010) provided a relation to estimate the lower limit of the black hole mass in terms of $\gamma$ -ray luminosity, redshift, and variability timescale, which is defined as

(13) \begin{equation} \frac{\mathrm{M}}{ \ \ \mathrm{M_{\odot}}} \ge 1\times10^3 \ \frac{\Delta \mathrm{ t_{obs}}}{\mathrm{1+z}} \ \left[\frac{L_{\gamma}\mathrm{(1+z)}}{6.3\times 10^{40} \Delta \mathrm{ t_{obs}}}\right]^{\frac{1}{4+\alpha_{\gamma}}}\end{equation}

The $\gamma$ -ray luminosity is defined as

(14) \begin{equation} L_{\gamma} = 4 \pi d_L^2 (1+z)^{\alpha_{\gamma}-2} \times f\end{equation}

where f is the total $\gamma$ -ray flux calculated using equations (19) or (20), $(1+z)^{\alpha_{\gamma} - 2}$ represents a K-correction factor, and $d_L$ is the luminosity distance value estimated following the expression given as

(15) \begin{equation} d_L = \frac{c}{H_0} \int_1^{1+z} \ \frac{1}{\sqrt{\Omega_M x^3 \ + \ 1 \ - \ \Omega_M }} \ dx\end{equation}

Using the $\Lambda-\mathrm{CDM}$ model with $\Omega_M$ , we calculated $d_L$ , where $H_0$ is a hubble constant ( $H_0=73\,{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$ ).

The $\gamma$ -ray photon number per unit energy can be defined as a power-law distribution

(16) \begin{equation} \frac{dN}{dE} = N_0 E^{-\alpha_{\gamma}}\end{equation}

where $\alpha_{\gamma}$ represents the photon spectral index and $N_0$ is the normalisation constant. $N_0$ can be estimated by integrating equation (16)

(17) \begin{equation} N_0 = N_{(E_L - E_U)} \left( \frac{1}{E_L} - \frac{1}{E_U}\right)\!, \quad \mathrm{if} \ \ \alpha_{\gamma} = 2;\end{equation}

(18) \begin{equation} N_0 = N_{(E_L - E_U)} . \frac{1 - \alpha_{\gamma}}{E_U^{1 - \alpha_{\gamma}} - E_L^{1 - \alpha_{\gamma}}}, \quad \mathrm{if} \ \ \alpha_{\gamma} \neq 2\end{equation}

where $N_{(E_L - E_U)}$ is the integral photons in units of photon cm⁻² s⁻¹ in the energy range $E_L - E_U$ , $E_L$ and $E_U$ represent the lower and upper energy limits respectively. In this study, we set $E_L$ = 100 MeV and $E_U$ = 300 GeV. Thus, the total $\gamma$ -ray flux can be calculated by $f = \int_{E_L}^{E_U} \ EdN$ , which elaborated as

(19) \begin{equation} f = N_{(E_L - E_U)} \left( \frac{1}{E_L} - \frac{1}{E_U}\right)ln \frac{E_U}{E_L}, \quad \mathrm{if} \ \ \alpha_{\gamma} = 2; \end{equation}

(20) \begin{equation} f = N_{(E_L - E_U)} \frac{1 - \alpha_{\gamma}}{2 - \alpha_{\gamma}} \ \frac{\left( E_U^{2-\alpha_{\gamma}} \ - E_L^{2-\alpha_{\gamma}} \right)}{E_U^{1-\alpha_{\gamma}} \ - E_L^{1-\alpha_{\gamma}}}, \quad \mathrm{otherwise} \end{equation}

in units of GeV cm⁻² s $^{-1}$ .

Acharyya et al. (Reference Acharyya2023) reported the spectral parameters of power-law distribution equation (16), including $N_0 = (6.52 \pm 1.12 )\times 10^{-11} \ \mathrm{cm^{-2}} \ \mathrm{s^{-1}} \ \mathrm{MeV^{-1}}$ and $\alpha_{\gamma} = 1.79 \pm 0.14$ . The estimated $\gamma$ -ray luminosity is $\sim 1.18\times 10^{47}\,{\rm erg}\,{\rm s}^{-1}$ and the mass of the black hole to be $M_* \gt 1.42 \times 10^8\,{\rm M}_{\odot}$ using the equations (14) and (13), respectively.

Additionally, Liu & Bai (Reference Liu and Bai2015) proposed a model to determine the upper limit of the black hole mass based on variability studies, which is briefly discussed here. In this model, a blob in the jet-production region initially has a size of $D_0$ and expands to $D_R$ at a distance $R_{jet}$ from the central black hole. As the blob propagates outward along the jet, it expands with an average velocity $\left( \bar{v}_{exp} \right)$ . Since the expansion velocity is non-relativistic ( $\bar{v}_{exp} \lt\lt \bar{v}_{jet}$ , where $\bar{v}_{jet}$ is the relativistic jet velocity), we have the condition $D_0 \lesssim D_R$ . Following equation (12), we can express this relationship as:

(21) \begin{equation} D_0 \lesssim D_R \leqslant \frac{\delta \Delta t_{\mathrm{min}}^{obs}}{1+z} c\end{equation}

where $\Delta t_{\min}^{obs}$ is the observed minimum timescale of variability. Relativistic jets are believed to originate from the inner regions of the accretion disk, close to the central black hole (Blandford & Znajek Reference Blandford and Znajek1977; Blandford & Payne Reference Blandford and Payne1982; Meier, Koide, & Uchida Reference Meier, Koide and Uchida2001). The inner radius of the disk is typically assumed to be near the marginally stable orbit, also known as the innermost stable circular orbit (ISCO). The ISCO radius ( $r_{ms}$ ) is expressed in terms of the gravitational radius ( $r_g$ ) and the dimensionless spin parameter $j = J/J_{\max}$ , where the maximum angular momentum is given by $J_{\max}=GM_*^2 / c$ , as defined in Equation (5) of Liu & Bai (Reference Liu and Bai2015). For a Schwarzschild black hole ( $j=0$ , non-rotating), the ISCO is located at $r_{ms}=6r_g$ , where $r_g=GM_* / c^2$ . The corresponding size of the emitting region is $D_{ms}=12r_g$ .

In the case of Kerr black hole (j = 1, maximally rotating), for prograde orbits, the emitting region is $D_e = 2r_e = 4r_g$ , where $r_e$ represents the equatorial boundary of the ergosphere. These estimated sizes are consistent with findings from general relativistic magnetohydrodynamic (GRMHD) simulations (Meier, Koide, & Uchida Reference Meier, Koide and Uchida2001). Thus, the size of the blob spans from $D_0=4-12 r_g$ for $0\leqslant j \lesssim 1$ . Following equation (13), we have (as given in equations 6a and 6b in Liu & Bai Reference Liu and Bai2015)

(22) \begin{align} M_{{var}}^{{Sch}} &\lesssim 1.695 \times 10^4 \frac{\delta \Delta t_{\mathrm{min}}^{obs}}{1+z} {\rm M}_{\odot}, \quad (j \sim 0) \end{align}

(23) \begin{align} M_{{var}}^{{Kerr}} &\lesssim 5.086 \times 10^4 \frac{\delta \Delta t_{\mathrm{min}}^{obs}}{1+z} {\rm M}_{\odot}, \quad (j \sim 1) \end{align}

where $\Delta t_\mathrm{min}^{obs}$ is in units of seconds.

Acharyya et al. (Reference Acharyya2023) reported the fasted observed variability timescale in $\gamma$ -rays of 6.2 $\pm$ 0.9 hr around MJD 57178. We utilised this timescale to constrain the black hole mass. The calculated black hole masses are $M_{var}^{Sch} \lesssim2.74 \times 10^9\,{\rm M}_{\odot}$ and $M_{var}^{Kerr} \lesssim8.22\times 10^9\,{\rm M}_{\odot}$ . In the calculation, we adopted the Doppler factor value, $\delta=9.6$ (Zhou et al. Reference Zhou, Zheng, Zhu and Kang2021; Acharyya et al. Reference Acharyya2023). The uncertainties in the estimated black hole masses were determined based on the errors in the variability timescale. However, the lack of precise information regarding the exact location of the emission region relative to the central black hole could introduce additional uncertainties, potentially increasing the error bars on the mass estimates.

By combining the estimates of the lower and upper limits of the black hole mass, we obtain a mass range of $(0.142 \lt M_* \lt 8.22)\times 10^9\,{\rm M}_{\odot} $ . Notably, the lower limit derived from the variability study is higher than the estimate reported by Pei et al. (Reference Pei, Fan, Yang, Huang and Li2022).

In recent years, several studies have demonstrated a correlation between variability timescales and black hole mass. Burke et al. (Reference Burke2021) identified a relationship between optical timescales in accretion disks and black hole masses, spanning the entire mass range of SMBH and even extending to stellar-mass BH systems. This correlation has been further explored across different EM bands, including $\gamma$ -rays (Ryan et al. Reference Ryan, Siemiginowska, Sobolewska and Grindlay2019; Zhang, Yan, & Zhang Reference Zhang, Yan and Zhang2022; Zhang, Yan, & Zhang Reference Zhang, Yan and Zhang2023; Sharma et al. Reference Sharma, Kamaram, Prince, Khatoon and Bose2024; Zhang, Yang, & Dai Reference Zhang, Yang and Dai2024; Sharma, Prince, & Bose Reference Sharma, Prince and Bose2024), X-rays (Zhang, Yang, & Dai Reference Zhang, Yang and Dai2024), and sub-millimeter wavelengths (Chen et al. Reference Chen, Bower, Dexter, Markoff, Ridenour, Gurwell, Rao and Wallström2023). Notably, the observed $\gamma$ -ray variability timescales in AGNs have been found to overlap with those in the optical band (Burke et al. Reference Burke2021), suggesting a link between non-thermal jet emissions and thermal emissions from accretion disks. Thus, by utilising the relationship between the rest-frame variability timescale $\left(\tau_{rest}^{Damping}\right)$ and BH mass $\left(M_{BH}\right)$ as established by Burke et al. (Reference Burke2021), $\tau_{rest}^{damping}=107_{-12}^{+12} \left(\frac{M_{BH}}{10^8\,{\rm M}_{\odot}}\right)^{0.38_{-0.04}^{+0.05}} \ \mathrm{days}$ , the mass of the SMBH can be estimated.

In our study, we estimated the mass of the SMBH to be approximately $(7.3 \pm 6) \times 10^9\,{\rm M}_{\odot}$ based on a rest-frame damping timescale of $\sim$ 547 days observed in $\gamma$ -ray emission of the blazar ON 246. While the uncertainty in the SMBH mass estimation is considerably high, the derived mass value is less and close to the upper limit of $M_{var}^{Kerr}$ and higher than the $M_{var}^{Sch}$ , as calculated using the model proposed by Liu & Bai (Reference Liu and Bai2015). Further, the SMBH mass of ON 246, estimated using the optical timescale derived from DRW modelling, is $(1.56 \pm 0.66)\times 10^8\,{\rm M}_{\odot}$ .

The estimated black hole mass of ON 246, derived from variability characteristics in both energy bands, falls within the range of $(0.142 - 8.22) \times 10^9\,{\rm M}_{\odot}$ . The estimated mass also lies in the range (10 $^8$ –10 $^9$ M $_{\odot}$ ) of mass of SMBH for FSRQ as derived in various studies (Ghisellini & Tavecchio Reference Ghisellini and Tavecchio2008; Castignani et al. Reference Castignani, Haardt, Lapi, De Zotti, Celotti and Danese2013; Xiong & Zhang Reference Xiong and Zhang2014; Paliya et al. Reference Paliya, Domnguez, Ajello, Olmo-Garca and Hartmann2021; Zhang et al. Reference Zhang, Yang and Dai2024) using various sample of FSRQs. Using the gamma-ray variability time, Pei et al. (Reference Pei, Fan, Yang, Huang and Li2022) have also estimated the mass of the black hole, but the value is close to $\sim$ 8.08 $\times$ 10 $^{7}$ M $_{\odot}$ , which is much smaller than the value estimated above. This discrepancy is because in Pei et al. (Reference Pei, Fan, Yang, Huang and Li2022), authors have used a fixed variability time of 1 day and the Doppler factor ( $\delta$ ) as a very small value of 0.48.

8. Result and discussion

In this study, we investigated quasi-periodic signals in the $\gamma$ -ray and optical emissions of the blazar ON 246 (4FGL J1230.2+2517) over the period of 11.6 yr, from MJD 55932 to 60081, employing three different methodologies as outlined in Section 3. Our analysis reveals a distinct quasi-periodic signal in the $\gamma$ -ray emissions of the blazar ON 246, with periods of approximately 746 days, 757 days, and 781 days, as identified using the Lomb-Scargle Periodogram (LSP), Wavelet Weighted Z-transform (WWZ), and REDFIT analysis, respectively. The significance of the detected periodicity was assessed through two independent approaches: Monte Carlo simulations and stochastic modelling using the damped random walk (DRW) model. The dominant peaks detected in both LSP and WWZ have a local significance level above 3 $\sigma$ , based on Monte Carlo simulations, and exceed 3 $\sigma$ when evaluated using DRW modelling. The uncertainties in the observed periods were estimated as the half-width at half-maximum (HWHM) by fitting the peak profiles with a Gaussian function. The QPO detected in $\gamma$ -ray emissions is further corroborated by an independent analysis of the optical light curve. Using the same methodologies, we searched for oscillatory signals in the optical data and assessed their significance. A similar periodicities of $\sim$ 757, $\sim$ 763, and $\sim$ 757 days were found in LSP, WWZ, and REDFIT analyses, respectively, with a significance level exceeding 3 $\sigma$ , reinforcing the presence of the detected QPO. Additionally, the phase-folded light curves of Fermi-LAT and ASAS-SN were fitted with sinusoidal functions corresponding to frequencies of $\sim 0.00134 \ \mathrm{day^{-1}}$ and $\sim 0.00132 \ \mathrm{day^{-1}}$ , respectively, as illustrated in Figure 11.

Figure 11. The folded Fermi-LAT and ASAS-SN light curves of ON 246 with a period of 746 and 757 days are shown in the top and bottom panels, respectively. The dashed blue line represents the mean value, and the sine functions (red) with frequencies of 0.00134 and 0.00132 day $^{-1}$ were fitted to the folded $\gamma$ -ray and optical light curves, respectively. Two full period cycles are shown for better clarity.

In addition, we also employed the interpolation cross-correlation function (ICCF) to analyse the correlation between $\gamma$ -ray and optical emissions. As display in the Figure 10, it can be seen that the correlation coefficient is maximum at near zero time-lag, with $\tau_{cent}=0.736_{-2.73}^{+3.13}$ d, have a significance level exceeding 4 $\sigma$ , indicating that the variability features between these two bands are similar in nature and is believed to be originated by the lepton single-zone scenario of blazar emission (Cohen et al. Reference Cohen, Romani, Filippenko, Bradley Cenko, Lott, Zheng and Li2014). Studying correlations across multiple wavebands is crucial for gaining deeper insights into the emission mechanisms of these systems. A strong correlation between low-energy and high-energy emissions can be well explained by the leptonic model. In this scenario, low-energy radiation originates from synchrotron emission, while the same seed photons undergo synchrotron self-Compton (SSC) and inverse Compton scattering to produce high-energy radiation, leading to a significant correlation between the two energy bands (Maraschi, Ghisellini, & Celotti Reference Maraschi, Ghisellini and Celotti1992; Giommi et al. Reference Giommi1999; Tagliaferri et al. Reference Tagliaferri, Ravasio, Ghisellini, Tavecchio, Giommi, Massaro, Nesci and Tosti2003; Zheng et al. Reference Zheng, Zhang, Huang and Kang2013; Liao et al. Reference Liao, Bai, Liu, Weng, Chen and Li2014; Li et al. Reference Li, Jiang, Guo, Chen and Yi2016).

Conversely, when high-energy emission is generated through the external Compton process, which involves seed photons originating outside the jet (Malmrose et al. Reference Malmrose, Marscher, Jorstad, Nikutta and Elitzur2011; Liao et al. Reference Liao, Bai, Liu, Weng, Chen and Li2014), the correlation between low-energy and high-energy emissions tends to weaken or become insignificant. Our findings reveal a strong correlation between the optical and $\gamma$ -ray emissions of the blazar ON 246, supporting the leptonic model’s predictions. Acharyya et al. (Reference Acharyya2023) observed a strong correlation between MeV and GeV emissions with a peak at zero time-lag, suggesting the same origin, and also found the positive time-lags with radio and optical emissions.

In this study, we employed three different methods to estimate the black hole (BH) mass of the blazar ON 246 based on its variability characteristics. First, we determined the lower and upper limits of the BH mass using the models proposed by Yang & Fan (Reference Yang and Fan2010) and Liu & Bai (Reference Liu and Bai2015), respectively. The minimum variability timescale observed in the $\gamma$ -ray band yielded a lower mass limit of $M_* \gt 0.142 \times 10^9\,{\rm M}_{\odot}$ . For the upper limit, we obtained $M_* \lt 2.74 \times 10^9\,{\rm M}_{\odot}$ for a Schwarzschild black hole, and $M_* \lt 8.22 \times 10^9\,{\rm M}_{\odot}$ for a maximally rotating Kerr black hole. In a third approach, we utilised the damping timescales of $\gamma$ -ray and optical light curves, applying an empirical correlation between the rest-frame damping timescale ( $\tau_{{rest}}^{\textrm{Damping}}$ ) and the BH mass, as established by Burke et al. (Reference Burke2021). This method yielded BH mass estimates of $(7.3 \pm 6) \times 10^9\,{\rm M}_{\odot}$ from $\gamma$ -ray variability and $(1.56 \pm 0.66) \times 10^8\,{\rm M}_{\odot}$ from optical variability. Combining these results, we constrain the BH mass of ON 246 to lie within the range of approximately $(0.142 - 8.22) \times 10^9\,{\rm M}_{\odot}$ .

8.1 Potential mechanisms for QPO

Several physical models have been proposed in the literature to explain the phenomenon of quasi-periodic oscillations (QPOs) in blazars. These include scenarios involving supermassive binary black hole (SMBBH) systems, precession or helical motion of relativistic jets (geometric effects), and instabilities in the accretion flow within the disk. A more detailed discussion on each of these scenarios is given below.

The supermassive binary black hole (SMBBH) scenario provides a plausible explanation for long-term flux modulations observed in blazars (Sillanpaa et al. Reference Sillanpaa, Haarala, Valtonen, Sundelius and Byrd1988; Xie et al. Reference Xie, Yi, Li, Zhou and Chen2008; Valtonen et al. Reference Valtonen2008; Villforth et al. Reference Villforth2010; Graham et al. Reference Graham2015; Wang, Yin, & Xiang Reference Wang, Yin and Xiang2017; Gong et al. Reference Gong, Zhou, Yuan, Zhang, Yi and Fang2022b, 2024). Several sources have been identified as potential candidates exhibiting long-term periodic flux modulations across multiple wavebands. Notable examples include a $\sim$ 12-yr QPO in OJ 287 (Sillanpaa et al. Reference Sillanpaa, Haarala, Valtonen, Sundelius and Byrd1988; Valtonen et al. Reference Valtonen2008; Villforth et al. Reference Villforth2010; Sandrinelli et al. Reference Sandrinelli, Covino, Dotti and Treves2016), a $2.18 \pm 0.08$ -yr periodicity in PG 1335 + 113 (Ackermann et al. Reference Ackermann2015), a $1.84 \pm 0.1$ -yr period in PKS 1510-089 (Xie et al. Reference Xie, Yi, Li, Zhou and Chen2008), and a 3-yr periodicity in 3C 66A (Otero-Santos et al. Reference Otero-Santos2020). Several other sources are also considered potential candidates exhibiting long-term periodic flux variations consistent with the SMBBH scenario. These periodic modulations in blazars are often attributed to the orbital motion of the secondary black hole within the binary system.

Rieger (Reference Rieger2004) investigated the potential geometrical origins of periodicity in blazar-type sources, proposing that periodic variations in emission can result from three possible mechanisms, including orbital motion in a binary black hole system, jet precession, and intrinsic jet rotation. The precessional-driven ballistic motion is unlikely to produce observable periods shorter than several decades, While the orbital motion in a close SMBBH system produces a period of $P_{obs}\gtrsim$ 10 days and Newtonian-driven precession in a close SMBBH can be a possible mechanisms of periodicity $P_{obs} \gtrsim 1$ yr. Therefore, the observed periodic flux modulations in emission from these sources canbe reasonably explained by orbital-driven (nonballistic) helical motion in a close SMBBH system.

In our study, the observed oscillation period is $P_{obs}$ = 746 days, which is significantly shorter than the actual physical period $P_d$ due to the light-travel time effect (Rieger Reference Rieger2004). The period is related by the equation, $P_d = \frac{P_{obs} \Gamma^2}{1+z}$ , where $\Gamma$ is the bulk Lorentz factor. To estimate $\Gamma$ , we used a relation from (Megan & Padovani Reference Megan and Padovani1995) and Li et al. (Reference Li, Qin, Gong, Liu, Guo, Gao, Yi and Liu2024), $\Gamma \le \frac{1}{2} \left(\delta + \frac{1}{\delta}\right)$ , which is basically gives the lower limit to the $\Gamma$ at a given value of $\delta$ . This condition is valid only for $\delta \gt1$ . Using a Doppler factor of $\delta$ = 9.6 (Zhou et al. Reference Zhou, Zheng, Zhu and Kang2021), we estimate the lower limit of the Lorentz factor as $\Gamma \geq 4.85$ . Consequently, the intrinsic period of the blazar ON 246 is $\sim$ 36.28 yr. Based on this corrected QPO period, the estimated mass of the primary black hole is

(24) \begin{equation} M \simeq P_{\mathrm{corrected, yr}}^{8/5} \ R^{3/5} \ 10^6\,{\rm M}_{\odot}\end{equation}

where $P_{\mathrm{corrected}}$ is the real physical period in units of years, R represents the mass ratio of the primary black hole to the secondary black hole, $R=\frac{M}{m}$ . Previous studies (Roland et al. Reference Roland, Britzen, Caproni, Fromm, Glück and Zensus2013; Yang et al. Reference Yang, Yan, Zhang, Dai and Zhang2021; Gong et al. Reference Gong, Zhou, Yuan, Zhang, Yi and Fang2022a) have suggested that the R lies within the ranges of 4–10.5, 1–100, and 10–100, respectively. In our analysis, we adopt R within the range of 1–100 to calculate the mass of the primary black hole. Based on this, the estimated mass is found to be in the range of $ M \sim 3.12\times 10^8 - 5 \times 10^9\,{\rm M}_{\odot} $ . The estimated BH mass range is comparable with the derived mass range in Section 7, indicating that the observed year-like QPO is likely caused by non-ballistic helical motion driven by the orbital dynamics of a close SMBBH system.

Additionally, instabilities in the accretion flow within the disk may also contribute to the flux modulations in blazars. In this scenario, oscillations in the innermost regions of the accretion disk or Kelvin-Helmholtz instabilities could lead to quasi-periodic plasma injections into the jet, as a result producing quasi-periodic oscillations in the jet emissions (Gupta, Srivastava, & Wiita Reference Gupta, Srivastava and Wiita2008; Wang et al. Reference Wang, An, Baan and Lu2014; Bhatta et al. Reference Bhatta2016; Sandrinelli et al. Reference Sandrinelli, Covino, Dotti and Treves2016; Tavani et al. Reference Tavani, Cavaliere, Munar-Adrover and Argan2018). The mass of the SMBH can be estimated using the following equation

(25) \begin{equation} M=\frac{3.23 \ \times \ 10^4 \ \delta \ P_{obs} }{\left(r^{3/2}+a\right) \left(1+z\right)}\,{\rm M}_{\odot}\end{equation}

where $P_{obs}$ is in units of seconds, $\delta$ is the Doppler factor, parameter r is the radius of this source zone in untis of $GM/c^2$ , parameter a is the spin parameter of SMBH, and z is redshift. In our study, using equation (17), the estimated mass of the SMBH is $1.02 \times 10^{12}\,{\rm M}_{\odot}$ for a Schwarzschild black hole with r = 6 and a = 0, while for a maximal Kerr black hole with r = 1.2 and a = 0.9982, the estimated mass of SMBH is $6.52 \times 10^{12}\,{\rm M}_{\odot}$ . The estimated masses of SMBH significantly exceed the black hole mass calculated in Section 7. Therefore, the observed flux modulations in both bands are unlikely to be associated with this scenario.

As previously discussed, blazar emission is primarily jet-dominated, and variability in jet emission can also arise due to geometric effects. One such scenario involves a plasma blob moving along a helical trajectory within the jet, causing observed emission variations (Villata & Raiteri Reference Villata and Raiteri1999; Rieger Reference Rieger2004; Li et al. Reference Li, Xie, Chen, Dai, Lei, Yi and Ren2009; Li et al. Reference Li, Chen, Yi, Jiang, Chen, Lü and Li2015; Mohan & Mangalam Reference Mohan and Mangalam2015; Sobacchi, Sormani, & Stamerra Reference Sobacchi, Sormani and Stamerra2017; Zhou et al. Reference Zhou, Wang, Chen, Wiita, Vadakkumthani, Morrell, Zhang and Zhang2018; Otero-Santos et al. Reference Otero-Santos2020; Gong et al. Reference Gong, Zhou, Yuan, Zhang, Yi and Fang2022b; Gong et al. Reference Gong, Tian, Zhou, Yi and Fang2023; Prince et al. Reference Prince2023; Sharma et al. Reference Sharma, Kamaram, Prince, Khatoon and Bose2024). In this geometrical model, emission is enhanced due to relativistic beaming, and as the blob follows a helical path within the jet, the viewing angle changes over time. Such helical jet trajectories can result from jet bending, wiggling motion, or the presence of helical magnetic fields within the jet, leading to periodic emission patterns. Consequently, the blob’s motion causes a continuous variation in the viewing angle relative to the observer’s line of sight, which can be characterised as

(26) \begin{equation} \mathrm{cos} \ \theta_{{obs}}(t)= \mathrm{sin} \phi \ \mathrm{sin} \psi \ \mathrm{cos} (2\pi t / P_{{obs}}) \ + \ \mathrm{cos} \phi \ \mathrm{cos} \psi\end{equation}

where $\phi$ represents the pitch angle between the blob motion and the jet axis, $\psi$ is the viewing angle or inclination angle between the observer’s line of sight and jet axis, and $P_{{obs}}$ represents the observed period of oscillations in the light curve. Due to the helical motion of emitting region in the jet, the Doppler factor undergoes temporal variations, defined as

(27) \begin{equation} \delta = \frac{1}{\Gamma \left(1 \ -\ \beta \ \mathrm{cos} \theta_{{obs}}(t)\right)}\end{equation}

where $\Gamma = \frac{1}{\sqrt{1 - \beta^2}}$ represents the bulk Lorentz factor and $\beta = v_{\mathrm{jet}/c}$ . Following this, the observed flux is defined as $F_{\nu} \propto \delta^3$ . The relationship between the observed period and rest frame period of blob can be defined by the following expression:

(28) \begin{equation} P_{{rest}}=\frac{P_{{obs}}}{1 \ - \ \beta \ {\cos} \psi \ \mathrm{cos} \phi}\end{equation}

The blazer ON 246 is BL Lac type, therefore, we considered typical values of $\phi=2^{\circ}, \ \psi = 5^{\circ}$ (Sobacchi, Sormani, & Stamerra Reference Sobacchi, Sormani and Stamerra2017; Zhou et al. Reference Zhou, Wang, Chen, Wiita, Vadakkumthani, Morrell, Zhang and Zhang2018; Sarkar et al. Reference Sarkar2019; Prince et al. Reference Prince2023), and $\Gamma=4.85$ as estimated above. The period of oscillations in the rest frame to be $P_{{rest}} \sim$ 78.85 yr.

As discussed earlier, Rieger (Reference Rieger2004) outlined three potential mechanisms, one of which suggests that if the observed periodicity exceeds 1 yr, it can be reasonably attributed to the helical motion of the jet driven by the orbital motion of an SMBBH system. Using the expression in equation (16) and periodicity in rest frame, we estimate the mass of the primary black hole to be $M \sim 1.083 \times 10^9, \ 4.314 \times 10^9, \ 1.71 \times 10^{10}\,{\rm M}_{\odot}$ for $R=1,10$ , and 100, respectively. The estimated primary black hole mass in SMBBH for $R=1,10$ is consistent with the derived mass value utilising different approaches as described in Section 7. Therefore, our comprehensive analysis suggests that the observed periodicity is most likely attributed to the jet being driven by the orbital motion of the SMBBH system.