Multi-band optical variability on diverse timescales of blazar 1E 1458.8+2249

Abstract

This study presents an analysis of the optical variability of the blazar 1E 1458.8+2249 on diverse timescales using multi-band observations, including observations in the optical BVRI bands carried out with the T60 and T100 telescopes from 2020 to 2023 and ZTF gri data from 2018 to 2023. On seven nights, we searched for intraday variability using the power-enhanced F-test and the nested ANOVA test, but no significant variability was found. The long-term light curve shows a variability behaviour in the optical BVRI bands with amplitudes of $\sim$100% and in the gri bands with amplitudes of $\sim$120%, including short-term variability of up to $\sim$1.1 mag. Correlation analysis revealed a strong correlation between the optical multi-band emissions without any time lag. From 62 nightly spectral energy distributions, we obtained spectral indices between 0.826 and 1.360, with an average of $1.128\pm0.063$. The relationships of both spectral indices and colour with respect to brightness indicate a mild BWB trend throughout the observation period, both intraday and long-term. We also performed a periodicity search using the weighted wavelet Z-transform and Lomb–Scargle methods. A recurrent optical emission pattern with a quasi-periodicity of $\sim$340 days is detected in the combined V- and R-band light curves. The observational results indicate that the blazar 1E 1458.8+2249 has a complex variability, while emphasising the need for future observations to unravel its underlying mechanisms.

1. Introduction

A subclass of active galactic nuclei (AGNs), Blazars contain a relativistic jet closely aligned with the line of sight to the observer and a supermassive black hole (SMBH) with a mass of $10^{6-10} \,\,\mathrm{M}_{\odot}$ at the centre of the galaxy. They are characterised by rapid emission variations across wavelengths ranging from radio to gamma-ray, high and variable polarisation, apparent superluminal motion, non-thermal continuous emission, and high-energy gamma-ray radiation (e.g. Urry & Padovani 1995; Woo & Urry 2002).

The multiwavelength spectral energy distribution (SED) of a blazar has a typical double hump structure (e.g. Fossati et al. 1998). The low energy bump peak ranges from IR/optical to UV/X-rays and is well explained by the synchrotron mechanism of electrons travelling at relativistic speed entering the magnetic field of the jet. The inverse Compton process satisfactorily reproduces the high-energy peak in the MeV-TeV range (e.g. Abdo et al. 2010; Sikora et al. 2009). Blazars can also be divided into two subclasses: BL Lacertae objects (BL Lacs) and flat-spectrum radio quasars (FSRQs). FSRQs have broad emission lines in their optical spectra, whereas BL Lacs have very weak or no emission lines in their featureless optical spectra (Marcha et al. 1996; Beckmann & Shrader 2012).

The nature of blazars can be studied through their variability, which is generally divided into intraday variability (IDV), short-term variability (STV), and long-term variability (LTV). The IDV, also known as microvariability, refers to flux variations over timescales ranging from tenths of minutes to hours, while the STV ranges from days to months, and the LTV has a longer timescale from months to years (e.g. Xie et al. 2004; Gupta et al. 2008; Agarwal et al. 2022). Various models related to jet and accretion disc have been proposed to explain the variability on diverse timescales, e.g. shocks moving within the jet (Marscher & Travis 1996), gravitational microlensing (Schneider & Weiss 1987), and variations in the Doppler factor due to the spiral motion of the emission plasma (Raiteri et al. 2017). However, many details of the models are still under discussion (Bhatta 2021).

The colour and spectral changes associated with the variation of optical flux in blazars provide valuable information about their emission mechanisms. The bluer-when-brighter (BWB) trend is often observed in BL Lac objects. This trend is often explained by the shock-in-jet model and variations in the Doppler factor (e.g. Ghisellini et al. 1997; Villata et al. 2002; Rani et al. 2010; Agarwal et al. 2019). The redder-when-brighter (RWB) trend is more common in FSRQs and is thought to be caused by additional radiation from the accretion disc (e.g. Villata et al. 2006; Gu & Ai 2011). Another colour trend is achromatic behaviour, which is often attributed to variations in the Doppler factor in the context of the geometrical scenario (Villata et al. 2002). However, there are indications of more complex trends at different timescales (e.g. Gupta et al. 2016; Isler et al. 2017; Agarwal et al. 2021; Negi et al. 2022), and a consistent framework has not yet been established for blazars.

Across the electromagnetic spectrum, from radio to gamma rays, periodic or quasi-periodic oscillations (QPOs) of flux variations have been reported in some blazars (e.g. Sillanpaa et al. 1988; Jorstad et al. 2001; Raiteri et al. 2001, 2003; Gupta, Srivastava, & Wiita 2009; Villata et al. 2009; Wang et al. 2014; Sandrinelli et al. 2017; Roy et al. 2022b; Chang et al. 2023; Kishore, Gupta, & Wiita 2023; Tripathi et al. 2024). Reported timescales vary from tens of minutes to days, weeks, and even decades. Several mechanisms have been proposed to explain year-long periodicities. These include the presence of a binary SMBH system, a plasma blob moving helically within the jet, and jet precession. Although the nature of the day-like QPOs is still uncertain, it is thought that they originate in the innermost central emission region, i.e. the accretion disc and the black hole.

1E 1458.8+2249 is a relatively bright BL Lac object with $z=0.235$ (Heidt & Wagner 1998). However, there are only a few studies of this object including flux variation and SED since the first epoch (Forman et al. 1978; Stocke et al. 1990). The study by Massaro et al. (2003) analysed optical and X-ray observations of the blazar from 1994 to 2001. The blazar had an R-brightness between 15.5 and 16.5 mag between 1994 and 1998. In January 2000, a flare with a maximum of 14.75 mag was detected in the optical R band. On 8 February 2001, during the second flare, the brightest magnitude observed so far was $R=14.62$ mag. The spectral energy distribution studies showed that the SED classification of 1E 1458.8+2249 is a high energy peaked BL Lac object (HBL; e.g. Fiorucci, Ciprini, & Tosti 2004; Fan et al. 2023). This indicates that the emission of 1E 1458.8 + 2249 is of synchrotron origin. Negi et al. (2022) studied the colour variations of 897 blazars using multi-band light curves obtained with the Zwicky Transient Facility (ZTF). They reported that 1E 1458.8 + 2249 exhibits BWB behaviour. This strong BWB trend ( $r\sim0.8$ ) between $V-I$ colour and V magnitude was also reported by Chang et al. (2023) in their intraday (ID) observations. The authors detected that one of the IDV light curves exhibited a drastic outburst with a variability amplitude of 0.97 mag ( $\%96.5$ ) on 26 May 2011, while the other IDV on 14 April 2020 showed a variability amplitude of only 0.1 mag ( $\%8.7$ ).

In this study, we present multi-colour observations of the blazar 1E 1458.8+2249 to investigate optical flux variations, correlations between different optical bands, and colour behaviour over intraday to long-term timescales. By combining ZTF and our observations, we analyse the quasi-periodic variability behaviour of the blazar.

2. Data

Optical observations were performed in the BVRI bands using a 1.0 m RC telescope (T100) and a 60 cm RC robotic telescope (T60) at the Tubitak National Observatory (TUG) from May 2020 to February 2023. The T100 telescope is equipped with a cryogenic cooled CCD (model SI 1100 Cryo, size 4 096 $\times$ 4 037 px). The T60 is equipped with a thermally cooled CCD (model Andor iKon-L, size 2 048 $\times$ 2 048 px). The exposure times were determined based on the photometric band and the brightness of the source, with a range of 20 to 180 s. A total of 1,555 BVRI frames were acquired over 169 nights, including seven nightly follow-up observations. Standard data reduction and calibration were performed on all CCD frames, i.e. bias subtraction, twilight flat-fielding, cosmic ray removal, and aperture photometry. In order to evaluate the instrumental magnitude of the blazar, we employed star A as the reference, as its position, brightness, and colour are closest to that of the blazar, and stars C1, C2, and C3 as comparisons, ^a as provided by Fiorucci, Tosti, & Rizzi (1998). It is important to note that only instrumental magnitudes were used for the B-band data since there are no B-band magnitudes for the standard stars. Therefore, the B-band light curves were not included in the analyses that require calibrated magnitudes. We provide our observational data, including instrumental magnitudes, exposure times, and dates, in the online supplementary data.

In addition to our observational data, we collected light curves in the gri bands from the ZTF database, ^b spanning from March 2018 (MJD 58194) and July 2023. In order to maintain the quality and goodness of the ZTF light curve, we selected 1,180 ZTF data points with catflags = 0 and chi < 4.

3. Intraday, short-term, and long-term variability

To understand the characteristics of the emission regions and the underlying radiative processes, it is essential to study the flux variability of blazars at different frequencies and timescales, from the shortest to the longest. Our quasi-simultaneous optical observations provide the opportunity to study the variability of the blazar from intraday to long-term scales. In order to obtain a meaningful IDV, we have made seven nightly observations of at least one hour duration, i.e. six of them in the optical R-band and one of them in BVRI. Fig. 1 shows the light curves for these seven observation nights.

Table 1. The test results of the IDV of the blazar 1E 1458.8+2249.

Figure 1. Intraday light curves of the blazar 1E 1458.8+2249 in the optical VRI bands. The blue, gray, and dark red circles represent the V-, R-, and I-band magnitudes, respectively, with a labelled offset. The dates of the observations are given at the top of each plot.

To search for variability in the intraday light curves of the blazar 1E 1458.8+2249, two widely used statistical tests were employed: the power-enhanced F-test and the nested ANOVA (analysis of variance) test (e.g. Gaur et al. 2015; Pandey et al. 2019; Agarwal et al. 2022; Weaver et al. 2020). As noted by de Diego et al. (2015), other commonly used statistical tests to study optical variability, such as the C-test, F-test, and chi-square test, have several limitations that can lead to less reliable results. To improve the robustness of the analyses, it is crucial to include differential light curves from several comparison stars.

The power-enhanced F-test is defined as the comparison of the variance of the blazar’s light curve with the combined variance of multiple comparison stars, as $F_{enh} = \frac{s_{bl}^2}{s_c^2}$ (de Diego 2014). The differential variance of the blazar is represented by $s_{bl}^2$ , while $s_c^2$ is the combined variance of the comparison stars. The degrees of freedom (DOF) are formulated as $u_{bl} = N-1$ and $u_c = k(N-1)$ . The critical value ( $F_{critical}$ ) for the variability test is calculated at a 99% confidence level ( $\alpha = 0.01$ ). If $F_{enh}$ is greater than or equal to $F_{critical}$ , the blazar variability hypothesis is valid. In our research, we used differential light curves of the reference star (the field star A) and three comparison stars to estimate $F_{\text{enh}}$ values and compare them with the critical value ( $F_c$ ) for each observation day.

Nested ANOVA, which is a modified version of ANOVA, is used as a statistical tool to assess mean variations between groups in the differential light curves of the blazar generated by using multiple reference stars. An advantage of nested ANOVA is that it is independent of a specific comparison star, allowing all available stars to be used as reference stars in the analysis. In our study, four standard stars (A, C1, C2, and C3) were used to construct differential LCs of the blazar. These differential LCs were divided into groups of five data points each. Following the methodology in Montgomery (2012) and de Diego et al. (2015), the mean square of the groups ( $MS_G$ ) and the nested observations within the groups ( $MS_{O_{G}}$ ) were calculated. The resulting ratio, $F = MS_G/MS_{O_{G}}$ , follows an F distribution with $(a-1)$ and $a(b-1)$ degrees of freedom in the numerator and denominator, respectively. The variability of the blazar’s light curve was tested for a significance level of $\alpha = 0.01$ .

The results of the $F_{\text{enh}}$ -tests and the nested ANOVA tests are listed in Table 1. We have assumed that the light curve is variable (V) if the F-statistic equals or exceeds the critical value ( $F_c$ ) in both tests; otherwise, it is non-variable (NV). Although the daily monitoring time ranges from 1 to 5 hours, no variability on minute timescales was found in the seven intraday LCs.

Using all the observational data, we have constructed the long-term optical light curve for both the ZTF gri and VRI bands, as shown in Fig. 2. Table 2 gives a brief overview of the LTV light curves, consisting of the minimum, maximum, and mean magnitudes. In addition, we have calculated the variability amplitude (A) as defined by Heidt & Wagner (1996). It is defined as $A = \sqrt{(A_{\text{max}} - A_{\text{min}})^{2} - 2\sigma^{2}}$ , where $A_{\text{max}}$ and $A_{\text{min}}$ are the maximum and minimum magnitudes of the LTV light curve, respectively, and $\sigma$ is the mean error. The long-term light curve shows variability with amplitudes ranging from 100% to 133% depending on the band. During our observations, we recorded the brightest magnitude as $R=15.109$ mag on MJD 59453, while the blazar was at its faintest with $R=16.136$ mag on MJD 59625. The mean R-band magnitude during the monitoring period was $15.716\pm0.052$ , and the variability amplitude was calculated as 102.397%. ZTF observations from 2018 to 2023 show the brightest magnitude at $r = 15.320$ and the faintest magnitudes at $r = 16.427$ , giving an average magnitude of $15.829$ and a variability amplitude of 122%. During these 5 yr of observation, no flare was detected, and the brightness of the blazar didn’t reach the maximum brightness ( $R\sim 14.6$ ) which as in February 2001. The short-term light curve shows a continuous brightening from 2021 April to September (over $\sim$ 150 days) by $\sim$ 1.1 mag in the ZTF gri bands. The LTV light curves, especially between MJD 8900 and 9300, show fluctuations on shorter timescales of the order of days with a smaller amplitude of up to 0.8 mag.

Table 2. Results of the LTV analysis of the blazar.

Figure 2. Long-term light curves of the blazar 1E 1458.8+2249 in the optical gri bands and VRI bands. Vertical lines represent the date of the intraday observation.

4. The multi-band colour/spectral behaviour and the correlation analysis

Quasi-thermal emission from the accretion disc and non-thermal synchrotron emission from the relativistic jet drive the rapid flux and spectral variations in blazars. Understanding the spectral variations of blazars is crucial for identifying the different components that influence the observed flux. To examine the daily multi-band optical spectral energy distributions (SEDs) of the blazar 1E 1458.8+2249, we used nearly simultaneous observations consisting of at least the four VRI- and gri-bands. The calibrated magnitudes of the blazar were corrected for the Galactic extinction using the $A_\lambda$ values taken from the NASA/IPAC Extragalactic Database. ^c . Using the zero-points in Bessell, Castelli, & Plez (1998) for the VRI-bands and in the filter profile service of the Spanish Virtual Observatory (SVO) ^d for the ZTF gri-bands, these extinction-corrected magnitudes were then converted into corresponding fluxes. The optical spectral energy distributions (SEDs) were computed for the 62 nights, as shown in Fig. 3. The SEDs were found to follow a simple power law relationship of the form $F_{\nu} \propto \nu^{-\alpha}$ , where $\alpha$ represents the optical spectral index. The spectral indices for each night were determined by fitting a linear model of $\log(F_{\nu}) = -\alpha~\log(\nu)+ C$ . The analysis reveals that the spectral indices range from 0.826 to 1.360, with an average of $1.128 \pm0.063$ . The standard deviation is calculated to be 0.11. Fig. 4 displays the correlation between the spectral index and V-band brightness. Linear regression analysis shows that the correlation has a slope of $0.333\pm0.048$ (see Table 4), i.e. a mild BWB trend, but there is no significant trend in $\alpha$ over time.

Figure 3. The intranight SEDs of the blazar 1E 1458.8+2249 in optical multi-bands.

Figure 4. Variation of spectral index with V-mag during the entire monitoring period. The colour of the dots indicates the time of observation.

Table 3. Correlation between the colour indices and V-band magnitude of multi-band IDV data.

Table 4. The results of the linear regression between the spectral index ( $\alpha$ )/colours and the corresponding magnitude on the long timescales.

Since the colour variability of the blazars reflects their spectral variability, we performed a correlation analysis between brightness and colours on both intraday and long-term timescales. We have plotted the $V-I$ and $V-R$ colour indices against the V magnitude shown in Fig. 5, and the results of the linear regressions are listed in Table 3. We found mild correlations ( $r\sim0.5$ and $p<0.05$ ) for both the $V-R$ and $V-I$ colour indices, indicating a non-significant BWB trend taking into account the large colour errors with an average of 0.079 mag.

Figure 5. The intraday colour-magnitude plots of the blazar 1E 1458.8+2249 for the date 15 Feb 2023. The silver and dark red dots represent the $V-R$ and $V-I$ colour indices. Black lines represent linear fits.

To investigate the long-term colour variations of the blazar, we analysed our VRI and ZTF gri data sets separately. We computed colour indices by coupling data from the same night. Therefore, we plot the colour indices with larger frequency bases, such as the $V - I$ and $V-R$ colours, against the V-band magnitude in Fig. 6. We have also plotted the $g-i$ and $g-r$ colours against the g-band magnitude for the ZTF dataset. We have listed the results of linear regression on the colour-magnitude relations in Table 4. We found a strong correlation between ZTF $g-i$ colour and g band magnitude with a correlation coefficient of $r=+0.84$ , while moderate correlations ( $r=0.66$ ) were found with the ZTF $g-r$ colour index. However, the colour relations of our data set have weaker correlations with $r\sim0.5$ .

Figure 6. colour-magnitude diagram of the blazar 1E 1458.8+2249 on long timescale. Black lines represent the linear fits.

Understanding the correlations between different spectral bands is crucial as any time lag identified implies spatial distinction between the corresponding emission regions. These results provide valuable insights into the physical processes taking place within the blazar, allowing a nuanced understanding of its complex dynamics and geometry. The Discrete Correlation Function (DCF) is a powerful statistical tool for investigating unevenly distributed time series of multi-band light curves (Agarwal & Gupta 2015, and references therein). The calculation of the DCF involves the use of the unbinned DCF (UDCF), expressed as:

(1) \begin{equation}UDCF_{ij}(\tau) = \frac{(a_i - \bar{a}) (b_j - \bar{b})}{\sqrt{(\sigma_{a^2} - e_{a^2})(\sigma_{b^2} - e_{b^2})}}\end{equation}

Here, $\bar{a}$ and $\bar{b}$ are the mean values of the respective time-series datasets, while $\sigma_{a,b}$ and $e_{a,b}$ are their standard deviations and errors. The time delay between two data points is denoted by $\Delta t_{ij} = (t_{bj} - t_{ai})$ . The DCF is then found by taking average of the UDCF values over the interval $\tau - \frac{\Delta\tau}{2} \leq \tau_{ij} \leq \tau + \frac{\Delta\tau}{2}$ , using the methodology suggested by Edelson & Krolik (1988):

(2) \begin{equation}DCF(\tau) = \frac{\sum_{k=1}^{m} \text{UDCF}_k}{M}\end{equation}

Here, ‘M’ is the number of pairwise time lag values within the specified $\tau$ interval.

In the examination of the long-term light curves, weighted mean magnitudes and mean Julian Dates (JD) were computed from nightly binned observations. The DCF analyses (see Fig. 7) were applied to each pair combination of the nightly binned ZTF gri and optical VRI long-term light curves. Strong correlations were found between all multi-band pairs using a time binning value of 2 days, and no time lag was detected.

Figure 7. Cross-correlation analysis of light curves in the ZTF gri and optical VRI bands using Discrete Correlation Function for the entire monitoring period.

5. Optical quasi-periodic oscillations

The availability of continuous observational data determines the accuracy of the quasi-periodic analysis. Thus, to fill gaps in the optical light curve and extend it, we transformed ZTF magnitudes to the V- and R-bands using the transformation equations ^e (Lupton et al. 2005). The resulting combined light curves, including both the transformed ZTF magnitudes and those obtained in our study, are shown in Fig. 8. It is worth noting that Zhang (2022) was able to identify quasi-periodicity in the quasar SDSS J1321+033055, even without taking into account the difference in magnitude between the CSS (Catalina Sky Survey) light curve and the ZTF light curves.

Figure 8. The V- and R-band light curve combined from the calibrated ZTF data and our observations.

To identify any underlying periodicity in the combined light curve of 1E 1458.8+2249, we employed two widely used methods: the Lomb–Scargle (LS) periodogram (Lomb 1976; Scargle 1982), and the Weighted Wavelet Z-transform (WWZ; Foster 1996).

The LS periodogram is a Fourier transform for detecting periodicity in unevenly distributed time series data. It uses likelihood minimisation to integrate sinusoidal components into the time series data. Utilising the Python Astropy package, ^f the LS periodogram provided a detailed analysis of quasi-periodic signals within the corresponding time series. The WWZ is an improved version of the wavelet transform that is useful for detecting periodic or pseudo-periodic signals, even in uneven time series, and for assessing their persistence throughout the observation. The signal is simultaneously decomposed into both the frequency and time domains. This method uses wavelet functions to fit the observations rather than sinusoidal components. We perform a thorough WWZ analysis ^g of the combined light curves over the entire observation period.

It is important to note that blazar light curves can exhibit both periodic and stochastic variability. The identification of QPOs is attributed to systematic dynamics within the disc or jet. It is critical to make reliable estimates of the significance of any possible detection that could be misinterpreted as periodic peaks due to noise in the power spectrum. Therefore, we simulated 10,000 light curves corresponding to the power spectral density of the blazar light curve using the same methodology ^h as described in O’Neill et al. (2022). By performing identical period analyses using WWZ and LS on the simulated light curves, as we did on the blazar light curve, we obtained significance levels of 95% and 99% for each frequency.

The results of the WWZ and LS analyses for the combined R-band light curve are shown in Fig. 9, including the WWZ power map corresponding to the periodicity and time, the periodogram of the time-averaged WWZ power, and the LS periodogram. The WWZ analysis revealed a dominant signal at a periodicity of $337^{+133}_{-67}$ days in the WWZ periodogram, spanning the entire observation time in the WWZ power map and exceeding the significant level at 99%. On the other hand, the LS periodogram exhibited two peak signals at $301\pm28$ days and $390\pm60$ days, with a combined signal at $390^{+55}_{-112}$ days. These signals exceeded the 95% significance level, but did not reach the 99% threshold. Similar results were obtained for the combined V-band light curve as well as for the ZTF g- and r-band light curves alone. Our results suggest a possible recurrent behaviour in the optical emission of the blazar with a periodicity of about 340 days.

Figure 9. (a) 2D plane contour plot of the WWZ power over time and frequency domain, (b) Time averaged WWZ power periodogram, and (c) the Lomb–Scargle periodogram of the combined R band light curve. The dashed brown and red lines represent the 95% and 99% significance levels, respectively.

6. Discussion and conclusions

The most common explanation for the flux variations of the BL Lac type blazars during an outburst is the jet-in-shock scenario (e.g. Marscher & Gear 1985; Wagner & Witzel 1995; Marscher et al. 2008). In this model, variability arises from a shock or blob propagating through relativistic plasma jets and interacting with the emitting regions. In addition, helical structures, precession or other geometric effects within the jets can change the Doppler factor of the emitting region resulting in flux variability (e.g. Villata & Raiteri 1999; Raiteri et al. 2017). Hotspots or other emission regions on the accretion disc may contaminate the observed optical emission in the blazar during the low state due to the weakened contribution of the jets to the thermal emission (e.g. Chakrabarti & Wiita 1993; Mangalam & Wiita 1993; Marscher & Gear 1985). These mechanisms provide plausible explanations for the observed variations on timescales ranging from hours to months and years.

We present the results of our multi-band observations from 2020 to 2023 and ZTF observations from 2018 to 2023, studying the intraday, short-term, and LTV of the blazar 1E 1458.8+2249 in optical bands. Statistical tests on the seven intraday light curves indicate that the blazar has no significant variation on minute timescales. Combining our IDV results with those obtained in Chang et al. (2023), the duty cycle (DC), which is defined as the fraction of nights showing variability out of the total number of nights (Romero, Cellone, & Combi 1999), is calculated to be about 10%. It’s likely that HBLs, such as 1E 1458.8+2249, have lower DCs and smaller variability amplitudes than low- and intermediate-energy peaked blazars due to the prevention of the magnetic field (Heidt & Wagner 1998; Romero et al. 1999; Gopal-Krishna et al. 2011). Romero (1995) demonstrated that axial magnetic fields inhibit the formation of Kelvin-Helmholtz instabilities in jets spanning from sub-parsec to parsec scales. This prevention is effective if the magnetic fields exceed the critical value $B_c$ , which is defined by the equation $B_c=[4\pi nm_ec^2(\Gamma^2-1)]^{1/2} \Gamma^{-1}$ . Here, n represents the local electron density, $m_e$ refers to the electron rest mass, and $\Gamma$ denotes the jet’s bulk Lorentz factor. Higher magnetic fields ( $B>B_c$ ) in HBLs would suppress the development of small-scale structures, hence decreasing the occurrence of microvariability in the optical light curves. Nevertheless, the DC fraction is below the typical value of $\sim$ 30–50 observed in HBL objects (e.g. Romero et al. 2002; Gopal-Krishna et al. 2011).

On the other hand, the extreme variation of 0.97 mag reported by Chang et al. (2023), in an I-band IDV LC with an average magnitude of 16.43 is comparable to both the STV and even the LTV of the blazar. All other intraday light curves were taken in brighter states of blazar 1E 1458.8+2249, with an average magnitude between 15 and 16 mag. The lack of significant intraday IDV detected in our study and in the literature suggests that the probability of detecting variability increasing with observation time alone may not sufficiently explain the absence of significant IDV. This makes it reasonable to consider alternative instability models for the IDV.

The blazar 1E 1458.8+2249 exhibits moderate STV and LTV, with an amplitude of up to 1.1 mag, and no flare is detected during the five years of observation. The multi-band LTV LCs show a strong correlation without any time lag for a 2-day time binning. This suggests that the optical emission region is likely to be co-spatial or that any spatial differences are too small to be detected with the optical datasets in this study. To understand the spectral variation of the source, we computed the spectral indices of the SEDs in the 62 nightly LCs ranging from 0.826 to 1.360 with an average of $1.128\pm0.063$ which are consistent with those obtained for HBLs within the range of 0.5–1.8 (e.g. Gaur et al. 2012; Meng et al. 2018). However, it is expected that the spectral index should be $\leq$ 0.80 based on simple pure synchrotron emission models and the spectral index of optically thin synchrotron-emitting plasma (Urry & Padovani 1995). This suggests that variations in the contribution of different emission components to the optical emission of the blazar 1E 1458.8+224, such as thermal emission from the accretion disc or host galaxy, or non-thermal emission from different regions of the relativistic jets, are likely responsible for both the deviation in the spectral indices and the larger mean spectral index. Analysis of the variation of the spectral indices with V-band magnitudes over the entire observation period shows a mild detectable BWB trend, consistent with the colour behaviour obtained from VRI-band data set for both IDV and LTV. However, ZTF spectral colours indicate stronger BWB trends. The shock-in-jet model often cited in BL Lac studies can explain the mild BWB trend for both IDV and LTV in our study.

Our results suggest a possible recurrent behaviour in the optical emission of the blazar 1E 1458.8+2249 with a periodicity of about 340 days. The jet itself is most likely driving the observed QPO in blazars, where it dominates the thermal emission and aligns closely with the observer’s line of sight (< ${10}^\circ$ Urry & Padovani 1995). It has been shown that if the blazar is part of a binary SMBH system it could also induce the jet precession effect caused by either closely orbiting binary black holes or warped accretion discs (e.g. Valtonen et al. 2008; Graham et al. 2015; Bhatta et al. 2016). These mechanisms are mostly expected to generate QPOs at year-like timescales. If the angle between the jet axis and the observer’s line of sight is too small and the jet has a high Lorentz factor, the jet precession effect could lead to shorter periodicities, such as the observed QPO of 340 days (e.g. Sandrinelli et al. 2016; Bhatta et al. 2016; Zhang et al. 2017). It is very likely that helical structures are common in blazars (Villata & Raiteri 1999; Rieger 2004). The Doppler boosting effect due to helical or non-ballistic motions of relativistic blobs or shocks within blazar jets could be responsible for the QPO in blazar LCs. (e.g. Camenzind & Krockenberger 1992; Mohan & Mangalam 2015). For the simplest leptonic one-zone model, due to the postulated helical motion of the blob, the viewing angle of the blob with respect to the line of sight ( $\phi_{obs}$ ) changes periodically with time, and the Doppler factor ( $\delta$ ) varies with viewing angle as $\delta=1/[\Gamma(1-\beta\cos{\theta(t)})]$ (Sobacchi, Sormani, & Stamerra 2017; Zhou et al. 2018). Here, $\Gamma=1/\sqrt{1-\beta^2}$ is the bulk Lorentz factor of the blob motion with $\beta=\nu_{jet}/c$ . Given this scenario, the periodicity in the rest frame of the blob is given by $P_{rf}={P_{obs}}/{(1-\beta\cos{\psi}\cos{\phi})}$ . For typical values of the pitch angle of the helical path $\phi=$ ${2}^\circ$ , the angle of the jet axis with respect to the line of sight $\psi=$ ${2}^\circ$ , and $\Gamma=4.6$ (Wu et al. 2007), the rest frame periodicity is calculated as $\sim$ 37 yr for an observed periodicity of $P_{obs}=\sim340$ day. During one period, the blob traverses a distance of $D = c\beta P_{rf}\cos{\phi} \simeq 11$ pc. In order to detect a statistically significant QPO in the present case, at least 4-5 cycles are required, thus the blob would traverse $\sim$ 55 pc throughout the observation time. Under these assumptions, this length scale indicates a very highly curved jet structure to drive the QPO. In addition, a slow increasing trend in the first 5 cycles and then a continuous power trend close to zero can be seen in Fig. 9, but a faster attenuation in the QPO amplitudes after peak is expected (Roy et al. 2022a). Thus, this scenario is less likely for the QPO obtained from the dataset we have. The periodicity driving mechanisms associated with the orbital motions of hot spots, spiral shocks, or other non-axisymmetric phenomena around the innermost stable circular orbit. (e.g. Mangalam & Wiita 1993; McKinney, Tchekhovskoy, & Blandford 2012; Huang et al. 2013; Gupta et al. 2019) are most likely to produce much shorter periodicities than we have found.

Further multi-wavelength observations on diverse timescales, complemented by theoretical modelling of the underlying mechanisms of blazar variability, are needed to improve our understanding of the complex variability of the blazar 1E 1458.8+2249.