Study of variability in long-term multiwavelength optical lightcurves of blazar AO 0235+164

Blazars belong to the radio-loud (RL) class of active galactic nuclei (AGNs). This extremely variable class is the union of BL Lacertae objects (BL Lacs) and flat spectrum radio quasars (FSRQs). Blazars host a large-scale relativistic jet of plasma pointing very close to the observer’s line of sight (Urry & Padovani, 1995). The jet is launched from the very near vicinity of the supermassive black hole (SMBH) of mass 10⁶ – 10¹⁰ M_⊙ at the center of the AGN (e.g., Woo & Urry, 2002). Blazars are characterized by highly variable emission throughout the whole electromagnetic (EM) spectrum, from radio to $\gamma$-rays, and their spectral energy distributions (SEDs) are characterized by two broad humps (Fossati et al., 1998). Blazars display high and variable polarization from radio to optical bands, and emit predominately non-thermal emission in the entire EM spectrum. The low-energy hump is ascribed to synchrotron radiation from relativistic leptons, whereas the high-energy hump arises from inverse Compton (IC) processes and sometimes from hadronic processes (e.g., Marscher, 1983; Mücke et al., 2003; Romero et al., 2017, and references therein).

Blazars display flux variability on diverse timescales ranging from a few minutes to several years. Blazar variability has often been divided into three categories, depending on the cadence of the observations: (i) microvariability (Miller et al., 1989), or intraday variability (IDV) (Wagner & Witzel, 1995), or intra-night variability (INV) (Sagar et al., 2004), focusing on the variability over a day or less; (ii) short-term variability (STV), focusing on variability over days to weeks, (iii) and long-term variability (LTV), focusing on timescales of months to years (e.g. Gupta et al., 2004).

The BL Lac object AO 0235+164 is at redshift $z=$ 0.94 (Cohen et al., 1987). Optical spectroscopic and photometric observations of the object have discovered two foreground-absorbing systems at $z=$ 0.524 and $z=$ 0.851 (Cohen et al., 1987; Nilsson et al., 1996; Raiteri et al., 2007). The flux of the source can be both absorbed and contaminated by these foreground systems, and the stars in them may act as gravitational micro-lenses that could contribute to the observed variability. Abraham et al. (1993) did deep CFHT imaging of AO 0235+164 and reported that the source is weakly amplified by macrolensing / microlensing by stars in the foreground.

AO 0235+164 has been extensively observed in the past from radio to $\gamma$-ray bands either in individual EM bands or quasi-simultaneously in multiple EM bands and has shown variations in all those bands on diverse timescales (e.g., Madejski et al., 1996; Rabbette et al., 1996; Takalo et al., 1998; Qian et al., 2000; Webb et al., 2000; Romero et al., 2000; Raiteri et al., 2006, 2008; Hagen-Thorn et al., 2008; Gupta et al., 2008; Agudo et al., 2011; Ackermann et al., 2012; Fan et al., 2017; Kutkin et al., 2018; Wang & Jiang, 2020, and references therein). It is one of the blazars which has displayed very high and variable optical/NIR polarization up to $\sim$45 percent (e.g., Impey et al., 1982; Stickel et al., 1993; Fan & Lin, 1999; Cellone et al., 2007; Ikejiri et al., 2011; Itoh et al., 2016, and references therein). In the Hamburg quasar monitoring program (HQM) this source was observed in the optical R band during 1988–1993, during which a 2.36$\pm$0.25 magnitude variation was detected; a particularly strong brightening in the source of $\sim$1.6 magnitude was reported during February 20–22, 1989 (Schramm et al., 1994). In six nights of optical $B$ and $V$ bands observations during 21–27 September 1992, the blazar was found in an unusually bright state and IDV was detected in both $B$ and $V$ bands (Rabbette et al., 1996). On another occasion, 6 nights of quasi-simultaneous $V$ and $R$ band observations in November 1999, revealed IDV with an amplitude of $\sim$100 percent over timescales of a day, while 0.5 magnitude changes were reported in both bands on a single night (Romero et al., 2000). In multicolor optical/NIR photometric ($BV\!RIJHK$) and $R$-band optical polarimetric observations of AO  0235+164 during its 2006 December outburst, variability on IDV timescales was detected, with increasing minimum timescale of variability from optical to NIR wavelengths; such variations were even detected in the optical polarization (Hagen-Thorn et al., 2008). In three nights of optical observations of the blazar in January – March 2007, IDV and STV were detected (Gupta et al., 2008).

In quasi-simultaneous optical ($V$ and $R$ bands) and radio (22 GHz) observations of AO 0235+164 during 1993–1996, the variability in optical bands showed amplitudes up to 1.5 magnitudes on STV timescales; although the radio variability is less dramatic, in general, it followed the optical behavior (Takalo et al., 1998). For the 1997 AO 0235+164 outburst, quasi-simultaneous multi-wavelength (MW) (radio, optical, NIR, and X-ray) observations were carried out. It was found that the source varied nearly simultaneously over 6 decades in frequency during the outburst and this result was explained in terms of a microlensing event (Webb et al., 2000).

An analysis of this source’s variability over $\sim$25 years led to the suggestion of a $\sim$5.7 years quasi-periodicity of the main radio and optical flares (Raiteri et al., 2001); however, the putative next outburst, predicted to peak around February–March 2004, did not occur, and a new analysis of the optical light curves on a longer time span revealed a characteristic variability timescale of $\sim$8 years, which was also present in the radio data (Raiteri et al., 2006). Recently, optical $R$ band photometric data taken during 1982–2019 showed 5 cycles of double-peaked periodicity of $\sim$8.13 years with a secondary peak following the primary one by $\sim$(1.5–2.0) years (Roy et al., 2022). In another MW campaign from radio to UV bands in 2006–2007, a huge NIR-optical-UV outburst with brightness increase of $\sim$5 magnitudes during February 19 – 21, 2007 was detected (Raiteri et al., 2008). During a major outburst seen in 2009, changes in radio, optical, X-ray, and $\gamma$-ray bands were found to be strongly associated (Agudo et al., 2011). In another simultaneous MW observing campaign of this blazar between 2008 September and 2009 February, $\gamma$-ray activity was found to be well correlated with a series of NIR/optical flares, accompanied by an increase in the optical degree of polarization; the X-ray light curve showed a different 20-day high state of an unusually soft spectrum which did not match the extrapolation of the optical/UV synchrotron spectrum (Ackermann et al., 2012).

AO 0235+164 is one of the sources that often used to be called OVV (optically violently variable). There are several such objects, like 3C 279, 3C 454.3, 4C 29.45, CTA 102, BL Lacertae, etc. Long-term achromaticity and zero lags have widely been found for these sources (Bonning et al., 2012; Zhang et al., 2021; Fan et al., 2006; Raiteri et al., 2017; Guo et al., 2015). AO 0235+164 is peculiar because it is commonly considered a BL Lac, one of the furthest known, but it shares properties with FSRQs. It is also a complex source because its light is contaminated by the southern AGN, ELISA, and absorbed by an intervening galaxy. This paper has undertaken a detailed analysis of the source’s optical brightness and spectral variability over a very long time span ($\sim$5 decades) as well as an investigation of its central engine. Our aim is to shed light on the long and short-term behavior of an emblematic BL Lac object through a detailed analysis of what is likely the most massive data set ever assembled for an object of this kind. The paper is organized as follows. In section 2, we provide descriptions of the observations of AO 0235+164. The section 3 gives our data analysis methods and results. We present a discussion and conclusions in section 4 and section 5, respectively.

Most of the optical $UBVRI$ observations of AO  0235+164 we have employed in this work are taken from The Whole Earth Blazar Telescope¹¹1https://www.oato.inaf.it/blazars/webt (WEBT) (Villata et al., 2002; Raiteri et al., 2017) which is an international collaboration of optical, near-infrared, and radio observers. WEBT has organized several monitoring campaigns on the blazar AO  0235+164, with the participation of many tens of observers and telescopes all around the world. Later, this source was studied by the WEBT and by its GLAST-AGILE Support Program (GASP) (Villata et al., 2008, 2009), which was started in 2007 to record quasi-simultaneous data of various blazars observed by the AGILE and Fermi (formerly GLAST) satellites. WEBT/GASP data on AO  0235+164 were published in Raiteri et al. (2001, 2005, 2006, 2008) and Ackermann et al. (2012). Raiteri et al. (2005) prescribed ways to remove the contribution of the southern galaxy ELISA from the observed optical flux densities and estimated the amount of absorption towards the source in excess of that from our Galaxy in X-ray, ultraviolet, optical, and near-infrared bands.

The WEBT and GASP data were calibrated following a common prescription, i.e., with the same photometry for the same reference stars. For calibration of the AO  0235+164 observations, the adopted photometric sequence includes stars 1, 2, and 3 from Smith et al. (1985). To build a reliable lightcurve for further analysis, clear outliers were removed and minor systematic offsets between various datasets were corrected.

AO  0235+164 was also observed with the 2.2 m telescope of Calar Alto Astronomical Observatory (CAHA, Spain) in November – December 2005, using the CAFOS instrument in imaging polarimetry mode, and photometric data were obtained by adding up the ordinary and extraordinary fluxes from each individual image (Cellone et al., 2007). Photometric data were also obtained with the 2.15 m telescope at Complejo Astronómico El Leoncito (CASLEO, Argentina) along several runs in November 1999, December 2000, August 2004, and January 2005. Results from these data were published in Romero et al. (1999, 2000, 2002) and in two papers by the WEBT collaboration focused on this blazar (Raiteri et al., 2005, 2006). Data from a more recent (December 2019) observing run with the same telescope were used in Roy et al. (2022). Magnitude calibration to the standard system was done using our own photometry of Landolt’s (2009) fields as well as standard stars in the field of AO  0235+164 (Smith et al., 1985; González-Pérez et al., 2001).

We also collected the publicly available optical $R$ and $V$-band data of AO 0235+164, taken at Steward Observatory²²2http://james.as.arizona.edu/~psmith/Fermi/DATA/Rphotdata.html, University of Arizona. These measurements employed the 2.3 m Bok and 1.54 m Kuiper telescopes between 4 October 2008 and 12 February 2018, using the SPOL CCD Imaging/Spectropolarimeter attached to those two telescopes. Details about the instrument, observation, and data analysis are given in Smith et al. (2009). In addition, we included the optical-$BVR$ data from the Small and Moderate Aperture Research Telescope System (SMARTS) public archive³³3http://www.astro.yale.edu/smarts/glast/home.php#. The SMARTS consortium is part of the Cerro Tololo Inter-American Observatory (CTIO), Chile, and has been observing Fermi-Large Area Telescope (LAT)-monitored blazars in the optical $B$, $V$, $R$ and NIR $J$ and $K$ bands. Details about the SMARTS instruments, observations, and data analysis procedures are given in Bonning et al. (2012). These standard magnitudes observed by CASLEO, CAHA, SMARTS, and the Steward observatory were further corrected for the southern galaxy ELISA following Raiteri et al. (2005). We also added other $R$-band optical photometric data from the literature (Takalo et al., 1998; Hagen-Thorn et al., 2008).

We combined all the optical $U$, $B$, $V$, $R$, $I$ band data to plot the long term (1974–2020) MW lightcurves of blazar AO 0235+164 (Figure 1). We removed the observations with errors of more than 0.1 magnitudes and studied long-term and intraday variability, color variation, spectral properties, and inter-band correlations.

In this work, we have presented a detailed temporal and spectral study of the highly variable emission from the blazar AO 0235+164 observed at multiple optical wavebands ($U\!BV\!RI$) from October 1975 to December 2019. The lightcurves have highly uneven data sampling due to gaps in observation seasons and non-uniform observation campaigns. Although $U$-band data are quite sparsely sampled the $BV\!RI$ observations have denser sampling when the source was highly active. Multiple long-term studies suggested that AO 0235+164 shows $\sim$2-year long flaring episodes with multiple sub-flares after intervals of $\sim$8 years (Raiteri et al., 2006; Fan et al., 2017; Roy et al., 2022). Figure 1 shows a difference of about six magnitudes between the quiescent and outburst states in all optical wavebands, corresponding to an energy flux variation of more than two orders of magnitude (Figure 6). The long-term variability amplitudes at all five wavebands are quite similar (Table 1). Also, we found a strong correlation with zero time-lag between the $U\!BV\!I$ observations and the $R$-band data (Figure 2 and Figure 3), which implies a common radiative process at a single emission zone is responsible for the bulk of the emission at the optical wavebands.

Sometimes during the quiescent states of powerful blazars, the disk thermal emission component becomes visible as a big blue bump on top of the synchrotron emission component from the jet in the optical-UV wavebands (e.g., Roy et al., 2021). As the disk emission is bluer than the jet synchrotron emission, an increase in the jet activity during low flux states displays a redder-when-brighter trend. The enhanced jet activity is observed when the charged particles inside the jet get accelerated to higher energies, and then radiate faster. Thus, the jet synchrotron component tends to get bluer with the increase in flux. If the jet emission completely outshines the disk emission, we expect to see a bluer-when-brighter trend (e.g., Isler et al., 2017). The flux increment can also be attributed to the increase in the jet Doppler factor (e.g., Papadakis et al., 2007), which blueshifts the spectrum and produces a bluer-when-brighter trend because of the convexity of the spectrum. Such a trend is seen in the $(B-I)$ vs $R$ magnitude diagram (Figure 4b) and indicates the domination of non-thermal jet emission over the thermal emission component of the accretion disk during both flaring and quiescent states. From the convex shapes of the optical $BV\!R$ SEDs during states ranging from quiescent to flaring (see the accompanying SED video and Figure 6), we may infer that the effect of the disk thermal emission is not significant in optical wavebands even during the low flux states.

This can be explained in terms of the nature of disk thermal emission given the disk luminosity and the central black hole mass computed in subsection 3.3. The primary, and most precise, black hole mass estimation methods are based on stellar and gas kinematics and reverberation mapping (e.g. Vestergaard, 2004). These methods need high spatial resolution spectroscopy data from the host galaxy and/or higher ionization emission lines and are not applicable to most BL Lacertae objects (BL Lacs). But in BL Lacs, if the weak emission lines are present, we can use the empirical methods (Kong et al., 2006) for BH mass estimation. The most common methods used for BH mass estimation for BL Lacs are the shortest variability timescales and periods of QPOs (Gupta et al., 2012). Since BL Lacs are highly variable objects, any BH mass estimation may be treated as an upper limit, and there are possibilities of detection of a shorter variability timescale or shorter QPO period. We obtained a log-scale BH mass of 7.90$\pm$0.25 in solar mass unit. The Steward observatory spectrum we used in our analysis had a narrower Mg II emission line (FWHM=3151 km s^-1) than those of Raiteri et al. (2007) and Paliya et al. (2021), thus resulting in a lower mass estimate. We considered a multi-temperature blackbody type accretion disk model, where the temperature at any portion of the disk is a function of the disk luminosity and the central black hole mass, to compute the thermal emission component. In Figure 10 we plotted the thermal component along with the optical-UV SED during the lowest activity state of AO 0235+164 observed on JD 2452169. It is evident that, as the thermal emission peaks at far UV frequencies ($\sim$3.5$\times 10^{15}$ Hz) in the observer’s frame of reference, the jet emission always dominates in $BV\!RI$ wavebands. We do not see any significant trend in the variation of the $(V-R)$ spectral index ($\alpha_{VR}$) (Figure 7). The sudden rise of the $U$-band flux in Figure 10 is an indicator of a probable UV-soft X-ray bump as discussed in Raiteri et al. (2005, 2006). According to these studies, the source of the bump is either an additional synchrotron component coming from a separate emission region in the jet or the emission of a continuous inhomogeneous jet is suppressed in near UV region due to a discontinuity in opacity or misalignment of that particular emission region. Ackermann et al. (2012) mentioned that the whole optical-UV spectrum is produced by a single synchrotron emitting zone as the shape of the bump does not change with luminosity. They attributed the UV spectral hardening to an artifact due to the overestimation of extinction by Junkkarinen et al. (2004).
For the detection of any statistically significant intraday variability in 33 lightcurves of AO 0235+164 observed at CASLEO/CAHA, we employed different statistical tests widely used in AGN variability studies. The reliability of each of these tests has been disputed (e.g. de Diego et al., 2015; Zibecchi et al., 2017), so we here employed a comparative approach that could allow us to circumvent the limitations affecting any individual test. In the first place, we used the scaled $C$-criterion and the $F$-test. The first compares the dispersion of the blazar lightcurve to the dispersion of a field star (control star), while the latter does so with the variances. According to Zibecchi et al. (2017) and Zibecchi et al. (2020), the F-test has a tendency to classify noisy non-variable curves as a variable (i.e., give false positives), while the $C$-criterion tends to give false negatives. Even though the $C$-criterion (Romero et al., 1999) cannot be considered as an actual statistical test, it may still be a useful parameter to detect variability with high significance. The $F$-test, on the other hand, does not always work as expected, because it is particularly sensitive to non-Gaussian errors (“red noise”), which are usually an issue when analyzing blazars DLCs.

We also used the power-enhanced F-test and the nested ANOVA test, which involve multiple field stars. It is expected that the power-enhanced F-test may also suffer from the same drawback of detecting false variability as the (original) $F$-test. In the nested ANOVA test, in turn, data grouping may lead to false results if data within a time span larger than the (unknown) variability timescale are grouped. Comparing the results of Table 6 and Table 7, while considering the tendencies of giving false results by the respective tests, we can confirm that the source was significantly variable in 4 out of 13 $V$-band lightcurves, and 9 out of 20 $R$-band lightcurves. The source seems to be probably variable in 3 $V$-band and 4 $R$-band lightcurves, and non-variable in the rest. On 1999 November 5, the combination of $C$-criterion and $F$-test indicates non-variability but the combination of power-enhanced $F$-test and nested ANOVA detects variability in the $R$-band lightcurve. The results in the $V$-band lightcurve on that day are exactly the opposite. Similar situations were observed also on 2001 November 9 and 2001 November 12. A visual inspection of the DLCs of these nights reveals that the blazar DLCs were classified as non-variable when either the control star DLC had higher variability (1999 November 5) or the measurement errors of the blazar DLCs were higher due to its low-flux state (2001 November 9 and 12). Higher measurement errors lead to a lower chance of significant variability detection. These strange results may be an example of the drawbacks of the applied methods when trying to recover low-amplitude variations from DLCs affected by non-Gaussian noise (part of the observations on that night were taken at air mass $>2$ and under non-photometric conditions). Otherwise, the combined results of different methods seem to more or less agree. Alongside the optical SED patterns, such frequent IDV establishes AO 0235+164 as a low-energy peaked BL Lac (LBL) object. High energy peaked BL Lacs (HBL) show significantly less optical intraday variability than the LBLs (Heidt & Wagner, 1998; Romero et al., 1999).

The differences in IDV behavior have been attributed to the strength of magnetic fields present in the jet of HBLs. A higher axial magnetic field ($B$) than a critical value ($B_{c}$) may prevent the generation of any bends and Kelvin-Helmhotz instabilities in the jet-base responsible for creating intraday microvariabilities. This indicates the presence of a weaker magnetic field than $B_{c}$ in the jet of AO 0235+164. The critical magnetic field ($B_{c}$) is given in Romero (1995) as

In this work, we conducted a study of long-term and short-term (intraday) variability in the optical multiwaveband observations of the blazar AO 0235+164. Here we summarize our results and the probable physical scenarios.

1. 1.

   We observed a variation of about six magnitudes between the quiescent and flaring episodes, or over two orders of magnitude variation in the SEDs.

2. 2.

   $U\!BV\!I$ lightcurves are highly correlated with the $R$-band lightcurve with zero time lag.

3. 3.

   A significant bluer-when-brighter trend is observed in the $(B-I)$ color variation with $R$-magnitude.

4. 4.

   All the optical $BV\!R$-band SEDs show convexity. These observations indicate that the optical emission is dominated by jet radiation.

5. 5.

   AO  0235+164 frequently shows statistically significant intraday variability in optical wavebands. This implies that AO 0235+164 is an LBL and probably has a weak magnetic field in the jet environment.

6. 6.

   From the analysis of a broad Mg II emission line in a spectrum of AO  0235+164 taken at a low state, we estimate a central black-hole mass of $\sim 7.9\times 10^{7}M_{\odot}$.