Multiwavelength variability of the blazar AO 0235+164

ABSTRACT

We present a study of the multiwavelength (MW) variability of the blazar AO 0235|$+$|164 based on the radio-to-|$\gamma$|-ray data covering a long time period from 1997 to 2023. The radio data are represented by the 1–22 GHz measurements from the SAO RAS RATAN-600 radio telescope, the 5 and 8 GHz data from the IAA RAS RT-32 telescopes, and the 37 GHz data from the RT-22 telescope of CrAO RAS. The optical measurements in the R-band were collected with the SAO RAS 1-m Zeiss-1000 and 0.5-m AS-500/2 telescopes. Additionally, we used the archive data at 230 GHz from the Submillimetre Array and the |$\gamma$|-ray data in the 0.1–100 GeV band from the Fermi-LAT point source 4FGL-DR2 catalogue. The variability properties during four epochs containing major flares and one epoch of relatively low activity were analysed using the fractional variability indices, discrete correlation functions, Lomb–Scargle periodograms, and structure functions. A significant correlation (⁠|$\ge \!2\sigma$|⁠) between the radio, optical, and |$\gamma$|-ray bands is found for all these periods with time delays from 0 to 1.7 yr. The relation between time delay and frequency is described by a linear law with a negative slope of |$-10$| d GHz|$^{-1}$|⁠. The discovered properties of MW variability for the low-activity period and for flaring states suggest that the mechanisms dominating the radio–|$\gamma$|-ray variations are not substantially different. The detected quasi-periodic oscillations of about 6 and 2 yr are tentative, as the time span of the observations includes fewer than four full cycles for the radio and optical data and only about three cycles for the Fermi-LAT data. These results should be interpreted with caution, given the limited number of observed cycles and the influence of red noise. We used cluster analysis to reliably separate the high and low-activity states and determined statistical differences in the main properties of AO 0235|$+$|164 non-thermal emission. The physical parameters of the radio jet were obtained using the Hedgehog model applied to the average radio spectrum of AO 0235|$+$|164 in the range 0.1–300 GHz. The effectiveness of replacing electrons with protons in the synchrotron radio emission of relativistic jets is shown for describing the nature of blazars and the generation of high-energy neutrinos.

1 INTRODUCTION

Blazars are a specific type of active galactic nuclei (AGNs) known for their relativistic jets aligned with the observer’s line of sight, which causes the radiation to be strongly Doppler boosted (Urry & Padovani 1995). Blazars encompass two categories: BL Lacertae objects (BL Lacs) and flat spectrum radio quasars (FSRQs). BL Lacs display a continuous, featureless emission or weak narrow emission lines in their optical/UV spectra with an equivalent width of |$\leqslant 5$|Å (e.g. Stickel et al. 1991; Marcha et al. 1996), while FSRQs exhibit prominent broad emission lines. The radiation of blazars has a bimodal spectral energy distribution (SED). It has been reliably established that the low-frequency peak is caused by the synchrotron emission of electrons moving in the magnetic field of the jet. The low-frequency peak is most often located in the region from the infrared to the optical part of the spectrum, but for some objects of the BL Lac type it is in the X-ray part. The high-frequency peak extends into the |$\gamma$|-energy region up to the TeV region and is caused by the inverse Compton scattering of thermal or synchrotron photons. Based on the location of the synchrotron peak (⁠|$\nu _{\rm peak}$|⁠) in a SED, blazars are classified as high, intermediate, and low synchrotron-peaked types, i.e. HSP, ISP, and LSP. According to the common criteria, HSP blazars have |$\nu _{\rm peak}\gt 10^{15}$| GHz, the LSP have |$\nu _{\rm peak}\lt 10^{14}$| GHz, and ISP class blazars are with interim values of the peak frequency (Abdo et al. 2010). The flux, polarization, and spectra of blazars are highly variable across the electromagnetic (EM) spectrum from radio to |$\gamma$|-rays (e.g. Gupta et al. 2017 and references therein). The flux variability for blazars is found on different time-scales, ranging from minutes to several decades (Miller, Carini & Goodrich 1989; Gopal-Krishna, Sagar & Wiita 1993; Wagner & Witzel 1995; Gupta et al. 2004). The jet, presumably originating from accretion on to a supermassive rotating black hole (SMBH) surrounded by an accretion disc, contains relativistic electrons, which produce soft photons from radio up to UV (or even soft X-rays) through synchrotron emission and high-energy photons up to GeV and TeV energies via the inverse Compton process that involves scattering of synchrotron photons as well as externally produced soft photons.

One of the effective ways to understand relativistic jet emission is studying the changes in the blazar physical properties that cause the observed multiwavelength (MW) light-curve variations. Investigation of the time-scales of flare events and their correlation across different frequency bands makes it possible to study a connection between spatially different emission regions in the AGN central parts (D’Ammando et al. 2019; Larionov et al. 2020; Abe et al. 2023; Krishna Mohana et al. 2024).

The BL Lac blazar AO 0235|$+$|1644 at z = 0.94 (Cohen et al. 1987) is a good candidate for a MW study of variability processes in AGNs, exhibiting extreme variability of non-thermal radiation across all EM spectrum. It has been extensively studied during the last four decades, across radio to |$\gamma$|-rays, displaying variability at time-scales from less than an hour to several years (Webb et al. 2000; Romero et al. 1997; Romero, Cellone & Combi 2000; Peng & de Bruyn 2004; Raiteri et al. 2001, 2005, 2008; Hagen-Thorn et al. 2008; Ackermann et al. 2012; Wang 2014; Volvach et al. 2015; Fan et al. 2017).

Ground-based VLBI interferometric observations reveal its extreme compactness on submilliarcsecond scales (⁠|$\le \!0.5$| mas) with superluminal apparent speeds and indications of a broad projected jet opening angle (Jorstad et al. 2001; Piner et al. 2006; Jorstad et al. 2017). The RadioAstron space VLBI mission suggests the presence of two scales in the inner compact structure of AO 0235|$+$|164: the first is the ‘core’, which is resolved by ground-based VLBI and equals to |$\sim \!350$||$\mu$|as, the second is ultra-compact, less than about 10 |$\mu$|as or 0.1 pc, and remains unresolved even with the longest ground–space baselines (Kutkin et al. 2018).

The compact nature and extreme flaring along with a small viewing angle (⁠|$\theta =1.7^{\circ }$|⁠) suggest a favourable geometry for observing the source’s jet activity (Hovatta et al. 2009). AO 0235|$+$|16 has undergone significant changes in the direction of its jet, which distinguishes it from most other blazars (Hagen-Thorn et al. 2018).

The first attempt to investigate the time-domain relationship between radio and |$\gamma$|-ray emission in order to localize the |$\gamma$|-ray emission site in blazars was done by Max-Moerbeck et al. (2014). If both types of emission are triggered by shocks propagating along a relativistic jet, the time delay between flares in the two bands depends on their separation. The |$\gamma$|-ray flare becomes observable after crossing the surface of unit |$\gamma$|-ray opacity. In the same way, the radio flare becomes observable upon crossing the surface of unit radio opacity (hereafter ‘radio core’). The time lag between these wavebands provides an estimate of the interval between the emergence of |$\gamma$|-ray and radio radiation. The blazar AO 0235|$+$|164 was the only case with a significance |$\ge 3\sigma$| amongst the 86 objects in the study of Max-Moerbeck et al. (2014); a correlation at a delay of 150 d was demonstrated. Later Wang & Jiang (2020) defined this delay as 45 d at a 3|$\sigma$| significance level. In Cheong et al. (2024), the 2|$\sigma$||$\gamma$|-ray–radio correlation was detected with time lag about 50 d for the period 2013–2019.

The analysis of the correlation between the optical and |$\gamma$|-ray flux variations in AO 0235|$+$|164 highlights their close relationship, emphasizing the interplay between different emission components (Rajput et al. 2021). Wang & Jiang (2020) show that the |$\gamma$|-ray and optical emitting regions are roughly the same and are located |$\sim \!6.6$| pc upstream of the core region of 15 GHz. The estimation of the radio core size predicts that the optical and |$\gamma$|-ray emitting regions are far away from the broad-line region (BLR).

For AO 0235|$+$|164, the quasi-periodic oscillations (QPO) at 3–6 yr time-scale on base of data taken between 1975 and 2000 in the radio and optical light curves were found (Raiteri et al. 2001; Fan et al. 2002), indicating a possible binary black hole system at its centre (Romero, Fan & Nuza 2003; Ostorero, Villata & Raiteri 2004). The analysis of radio data at 37 GHz proposed a model of close binary SMBH with components of similar masses around |$10^{10} \,{\rm M}_\odot$| (Volvach et al. 2015). Roy et al. (2022) present a comprehensive analysis of a long-term optical light curve from 1982 to 2019, revealing a periodicity of |$\sim \!8$| yr. The five major flares were followed by minor flares roughly after |$\sim \!2$| yr. This double-peaked periodicity is consistent with the hypothesis of a binary supermassive black hole system.

The investigations of Otero-Santos et al. (2023) suggest a potential QPO in the R-band light curve with a period of |$\sim \!8.2$| yr, in accordance with the previous reports of Raiteri et al. (2006) and Roy et al. (2022). But the authors noted that the temporal coverage of the |$\gamma$|-ray and radio data does not allow one to determine the existence of this period reliably.

This work presents a systematic study of the multiband variability of AO 0235|$+$|164. We used the |$\gamma$|-ray, optical, and radio observations covering the period of 1997–2023. The primary objective of our collaboration is to achieve dense data sampling to identify potential correlations between the blazar emission in different spectral ranges in order to determine the variability time-scale and possible quasi-periodicity. Such a vast data base for this specific source provides an opportunity to investigate variability characteristics on both temporal and spectral domains.

The paper is organised as follows: the details of multiband observations are outlined in Section 2. The MW light curves and some features of the most prominent flares are presented in Section 3; and in Section 4, the results of the used statistical methods, namely the structure function (SF), discrete correlation function, and Lomb–Scargle (L–S) periodogram, are discussed. The intraweek variability of AO 0235|$+$|164, studied from the RATAN-600 radio data obtained on a daily basis, is discussed in Section 5. The average radio spectrum in the range between 0.1 and 300 GHz is constructed in Section 6 with the aim to model it by combining the low-frequency component and high-frequency component (HFC). A set of quasi-simultaneous broad-band radio spectra allows us to discriminate the high- and low-activity states of the blazar in Section 7. Finally, the general results of our work are discussed in Sections 8 (Discussion) and 9 (Summary).

2 MULTIWAVELENGTH DATA

We have collected MW data sets from four radio and two optical instruments and the |$\gamma$|-ray Fermi-LAT telescope (Table 1). Part of the data are from public archives, and some recent observations are published for the first time, see the description of the instruments and data reduction below. At 230 GHz (1.3 mm), we utilized measurements from the Submillimetre Array¹ (SMA; Gurwell et al. 2007), obtained in the period from 2002 October to 2024 January.

2.1 RATAN-600

Flux densities at 1–22 GHz were obtained with the RATAN-600 radio telescope in 1997–2023. Broad-band radio spectra at six frequencies 0.96/1.2, 2.3, 4.7, 7.7/8.2, 11.2, and 21.7/22.3 GHz were measured quasi-simultaneously within 3–5 min (secondary mirrors N1 and N2). The RATAN-600 antenna and radiometer parameters are given in Table 2. We use the word ‘instantaneous’ to such RATAN-600 spectra in order to differ them from the ‘quasi-simultaneous’ spectra, measured often during months at different frequencies. Detailed descriptions of the RATAN-600 antenna, radiometers, and data reduction are given in Parijskij (1993), Verkhodanov (1997), Kovalev et al. (1999), Tsybulev (2011), Udovitskiy et al. (2016), Tsybulev et al. (2018), and Sotnikova (2020). The flux densities at 1−22 GHz |$S_{\rm \nu }$|⁠, their errors |$\sigma$|⁠, and average observing epochs (MJD and yyyy.mm.dd) are presented in Table 3.

Additional measurements of AO 0235|$+$|164 at 2.3 and 4.7 GHz were conducted with daily cadence from 2021May to 2022 June in the radio survey mode on the RATAN-600 Western Sector (secondary mirror N5). The parameters of four four-channel radiometers with a central frequency of 4.7 GHz and a radiometer at 2.3 GHz are given in Table 2. The measurements from the four 4.7 GHz radiometers were averaged. The total data set consists of 370 observing epochs. The averaged flux density errors are equal to 7 per cent and 3 per cent at 2.3 and 4.7 GHz, respectively. The average observing epochs (MJD and yyyy.mm.dd), flux densities and their errors are presented in Table 4. Details of the observations and data processing are described in Majorova, Bursov & Trushkin (2023) and Kudryashova, Bursov & Trushkin (2024).

2.2 RT-32

Flux densities at frequencies of 4.8 and 8.6 GHz were measured at several observing epochs using three RT-32 radio telescopes: Svetloe (Sv), Zelenchukskaya (Zc), and Badary (Bd) (Shuygina et al. 2019). All the antennas and receivers have similar parameters: bandwidth |$\Delta f_{0}=900$| MHz for both receivers with the central frequencies specified above; beam width at half-power level |${\rm HPBW}=7{_{.}^{\prime}} 0$| and |$3{_{.}^{\prime}} 9$|⁠, respectively; flux density limit |$\Delta F$| reaching about 20 mJy per scan with a time constant of 1 s for both frequencies under optimal observing conditions. Observations were performed in the elevation (Bd, Zc) or azimuth (Sv) drift scan mode.

During the first observing epoch in 2015–2017, the intraday variability (IDV) monitoring was carried out at 4.8 and 8.6 GHz with the Zc and Bd telescopes – about two 24-h experiments per month, previously published in Kharinov et al. (2020). In this paper, we use flux densities averaged over a day. The results for AO 0235|$+$|164 show weak flux density variations at levels of about 2–3 Jy on a year-long time-scale. For the calibration source PKS 0316+152, we used accepted flux densities from the ‘CATalogs supporting System’ (CATS; Verhodanov, Andernach & Chernenkov 2005): 3.09 and 1.72 Jy at 4.8 and 8.6 GHz, respectively.

The second observing epoch consists of a single experiment on 2019 August 4 at 8.6 GHz with the Zc telescope: eight measurements during 00:17-07:14 UTC with an average flux density of 0.34 Jy for AO 0235|$+$|164. The source J0521|$+$|1638 with an accepted flux density of 2.49 Jy was used as a reference point.

On the third observing epoch, from 2020 January to 2023 October, the Sv RT-32 had been performing weekly observations at both available frequencies. The source 3C295 was used as a reference with the flux densities from Perley & Butler (2017). Each data point represents an average over about a 30-min exposure (⁠|$\sim \!30$| drift scans).

The flux densities |$S_{\nu }$|⁠, their errors |$\sigma$|⁠, and average observing epochs (MJD and yyyy.mm.dd) are presented in Table 3.

Since the observing frequencies of the RT-32 (5.05, 8.63 GHz) and RATAN-600 (4.7, 8.2) are close, we use their rounded values, 8 and 5 GHz, in the further analysis. The frequencies 21.7/22.3, 11.2, 2.3, and 0.96/1.2 GHz are also used in rounded form: 22, 11, 2, and 1 GHz.

2.3 RT-22

The data on the flux densities at the 36.8 GHz frequency (further 37 GHz) for AO 0235|$+$|164 were obtained with the 22-m RT-22 radio telescope. The beam-modulated receivers were used to acquire the data. The antenna temperature from the source was measured as a difference between the signals from the radiometer output in two antenna positions, when the radio telescope was pointed at the source alternately with one or the other receiving horns (‘on-on observation’ method). Observations of the blazar consisted of 5–20 such measurements to achieve a required signal-to-noise ratio. Details of the observations and data reduction are presented in Volvach, Volvach & Larionov (2023) and Sotnikova et al. (2022a).

The total period of observations with the RT-22 spans the time interval since 2002 February till 2023 October. The flux densities, their errors and average observing epochs are presented in Table 5. The data up to 2012 have previously been used and published in Volvach et al. (2015).

2.4 Zeiss-1000 and AS-500/2 reflectors

The optical study (mainly in the R band) was carried out with the 1-m Zeiss-1000 (2002 August–2023 December) and 0.5-m AS-500/2 (2021 January–2023 December) optical reflectors.

First optical observations of AO 0235|$+$|164 in 2002 were mostly performed using the Zeiss-1000 telescope with a CCD (Charge Coupled Device) photometer based on a small (⁠|$520 \times 580$| pixels) front-illuminated CCD chip with rectangular pixels: |$18 \times 24$| |$\mu$|m. The photometer provided a total field of view of about |$2{_{.}^{\prime}} 4\times 3{_{.}^{\prime}} 7$|⁠. The maximum quantum efficiency of the system was only 55 per cent.

Due to these reasons, further observations with the optical telescope Zeiss-1000 till the end of 2023 were implemented with another CCD photometer in the same Cassegrain focus, equipped with a |$2048\times 2048$| px back-illuminated E2V chip CCD 42-40. Details of the instrumental setup are described in our previous paper (Vlasyuk et al. 2023).

The main characteristics of the instrumental complex of an AS-500/2 reflector are presented in Valyavin et al. (2022a, b). In order to reduce the system readout noise and significant dark current of the FLI (Finger Lake Instrumentation, Inc.) camera used, we installed in July 2023 a smaller but more sensitive detector: the back-illuminated electron-multiplying CCD camera Andor IXon|$^{\rm EM}$| + 897 with high quantum efficiency (⁠|$\sim \!90$| per cent in the 450–700 nm range).

This CCD camera with |$512 \times 512$| pixels, when mounted in the Cassegrain focus of AS-500/2, provides a |$7^{\prime }$| field of view with |$0{_{.}^{\prime\prime}} 82$|/pixel data sampling (1 pixel |$=$| 16 |$\mu$|m). The system readout noise is about 6 |$e^-$|⁠, the CCD operating temperature has been chosen to be |$-70^\circ$|C to minimize the dark current.

Both CCD photometers are equipped with similar sets of filters, which are close to the standard broad-band Johnson–Cousins system given the sensitivity of both CCDs.

The typical integration time for observations of the blazar was 300 s for Zeiss-1000 and 120 s for AS-500/2. During the period of high activity of the object, the integration time was less: up to 30 s for higher time cadence.

In order to obtain blazar flux estimates, we performed the standard reduction steps, which were described earlier (Vlasyuk 1993). We collected the optical data as the instrumental magnitudes of the source and comparison stars in the same field. Magnitude calibration was performed using preferably stars 8–11 from González-Pérez, Kidger & Martín-Luis (2001).

A comparison with the data obtained on the same nights with different telescopes has shown general agreement; in the cases where variations have been detected, they are consistent with the usual variations of the source, as it is shown in Fig. 1. The individual flux estimates are marked by their epochs. The mean deviation of the data from the linear law is |$0\rm{.\!\!^{ {\rm{m}}}}035$|⁠.

In order to provide a joint analysis of the optical and radio data, we averaged the measurements over individual nights and transformed the resulting values into fluxes according to the constant from Mead et al. (1990). We did not subtract from the optical fluxes the addition from ELISA (ELusive Intervening Southern AGN), the faint southern companion of AO 0235|$+$|164 with |$R=20\rm{.\!\!^{ {\rm{m}}}}50$| (Raiteri et al. 2008), supposing that this addition is insufficient for our following analysis.

The optical data are collected from 568 observing nights spanning between 2002 August and 2023 December. The average observing epochs, number of observations for each night, and obtained fluxes with their errors are presented in Table 6. Further in this paper we use the R-band fluxes not corrected for the Galactic extinction.

2.5 Fermi-LAT

We have made use of the |$\gamma$|-ray light curves available in the Fermi Large Area Telescope (Fermi-LAT) public Light Curve Repository (LCR²; Abdollahi et al. 2023). The LCR is a data base of multicadence flux-calibrated light curves for over 1500 sources considered as variable in the 10-yr Fermi-LAT point source (4FGL-DR2) catalogue (Ballet et al. 2020). The analysis has been performed with the standard Fermi-LAT science tools (version v11r5p3), which use maximum likelihood analysis (Abdo et al. 2009), considering a region with a radius of |$12^{\circ }$| centred on the location of the source of interest and the |$\gamma$|-rays energy range between 100 MeV and 100 GeV. A limit of Test Statistics |${\rm TS}\ge 4$| (approximately |$2\sigma$|⁠) was applied to compute the fluxes for each bin of the light curve, while only the measurements with |${\rm TS}\ge 9$| were selected for further MW analysis. We adopted the monthly binned light curves to have more evenly-sampled time series with higher time cadence in the period from 2008 August to 2023 December.

3 MW LIGHT CURVES

Fig. 2 shows the MW light curves of AO 0235|$+$|164 collected from 1997 March until 2023 December. The most complete light curves covering this period are those at radio frequencies of 5–22 GHz. The data for 37 GHz, 230 GHz, and the R-band cover a shorter interval: 2002–2023. The low-frequency radio curve at 2 GHz contains extended gaps due to strong radio interference. There are also seasonal gaps in the optical data, caused by the periods when the object is not visible in the spring–summer seasons. During the years 2009–2014, AO 0235|$+$|164 had been showing quiescent broad-band behaviour with a remarkable feature: very low activity lasting almost 4.5 yr, from 2009 November to 2014 March. We consider this time period as the low-state epoch with the minimum variability level (⁠|$F_{\rm var}$|⁠, see Section 4.1) across the light curves in all the bands. The light curves for this epoch are presented separately in Fig. 3.

We divided the light curves into four epochs (denoted by the blue vertical lines in Fig. 2 and listed in Table 7) in a way that each epoch includes one of the four major flares. As the major flares, we adopt the most prominent flares (⁠|$\gt \!3$| Jy) at 37 GHz, because for this frequency we have the most uniform and abundant measurements. The first epoch covers the period from 2002 August to 2007 October, during which the blazar underwent quasi-periodical flares in the optical range, separated by 480–620 d and having amplitudes between |$\sim \!0.6$| and |$\sim \!1.6$| mJy. The typical duration of these flares lies between 100 and 300 d. This value may be incorrect because the measurements do not cover the total period. In epoch 1, the maximum flux density was observed at 37 GHz, achieving 5.5 Jy in 2007 April, |${\rm MJD}\sim 54200$|⁠.

The second epoch spans the period of 2007 October–2013 April, when the blazar demonstrated only one strong optical flare with a complex shape, located near epoch |${\rm MJD}=54774$|⁠, 2008 November, with an R-band amplitude near 2.6 mJy. This optical outburst corresponds to powerful flares with flux densities between 5 and 6 Jy in the sub-mm and radio waves at a close time interval. The most powerful event in |$\gamma$|-rays over all the analysed period, |$1.2 \times 10^{-6}$| ph cm|$^{-2}$| s|$^{-1}$|⁠, is also located near this date.

Within |$\sim \!1600$| d after this flare, AO 0235|$+$|164 had been in a low state with optical fluxes between 0.1 and 0.3 mJy. The radio flux densities had been varying from |$\sim \!0.5$| Jy to |$\sim \!1.3$|–2 Jy at 1–22 GHz within this epoch. The |$\gamma$|-ray emission was also at minimal levels during this time.

The third epoch spans the period of 2013 April–2019 March, when AO 0235|$+$|164 demonstrated four pronounced flares in the optical range with maxima between |$R \sim 0.52$| mJy (⁠|${\rm MJD}=56952$|⁠, 2014 October) and 1.2 mJy (⁠|${\rm MJD}=57042$|⁠, 2015 January). These outbursts do not show any significant periodicity, and their duration did not exceed 200 d for the longest flare (⁠|${\rm MJD}=57221$|⁠, 2015 July). The flares at the radio frequencies seem to be more smooth and may be a superposition of few events with amplitudes between 1.5 and 3 Jy at |${\rm MJD} \sim 57250$|–57500, 2015 August – 2016 April.

The |$\gamma$|-ray emission is represented here by serial flares. The most powerful one with an amplitude of |$\sim \!7 \times 10^{-7}$| ph cm|$^{-2}$| s|$^{-1}$| is also located near |${\rm MJD}=57230$|⁠, which is close to one of the brightest optical events.

The last epoch covers the observations between 2019 March  and 2023 December. The optical light curve is represented here by a set of |$\sim \!15$| fast flares with amplitudes from 0.29 mJy (⁠|${\rm MJD}=59938$|⁠, 2022 December) to 1.64 mJy (⁠|${\rm MJD}=60179$|⁠, 2023 August). The longest outburst had a duration of about 70 d (the peak near |${\rm MJD}=59197$|⁠, 2020 December), other flares are much shorter: the most intensive flare mentioned above had a duration of about 10 d. The fastest flares span about 6–8 d (the outbursts near |${\rm MJD}=59630$|⁠, 2022 February, and |${\rm MJD}=59938$|⁠). Unfortunately, the real durations of some flares are unknown due to the absence of optical data within the periods when the blazar was unobservable (e.g. the flare with a maximum at |${\rm MJD}=59287$|⁠, 2021 March). The most powerful (and much slower) two or three radio flares with amplitudes between 2 and 2.5 Jy are located near epochs |${\rm MJD}\sim 59200$|⁠, 2020 December, and |${\rm MJD}\sim 59600$|⁠, 2022 January, demonstrating slight delays as the frequency decreases. The |$\gamma$|-ray emission in this period seems to be 2–3 times fainter than in the previous epochs, approaching levels of about |$2 \times 10^{-7}$| ph cm|$^{-2}$| s|$^{-1}$|⁠.

3.1 Parameters of the most prominent flares

At high radio frequencies above 10 GHz, flux density variations in AGNs can be well modelled by a superposition of flare components (Valtaoja et al. 1999; Hovatta et al. 2009). According to these papers, the ascending and descending wing of each flare may be modelled by an exponential law with the same amplitude and epoch of the flare and different time-scales. We used independent time-scales for the rise and decay of a flare in the approximation procedure, whilst the authors of the papers cited above supposed a close connection between them.

To pre-process the light curves, we used cubic Hermite interpolation. The flares were then detected using the Bayesian block method described in Scargle et al. (2013) and implemented in python by Wagner et al. (2022). As a result, for all the light curves we obtained a picture of decomposition into exponential components (Fig. 4).

The behaviour of the flares along the spectral ranges shows systematic variation: they become slower and more smooth as the frequency decreases, depicting less energetic processes. The most abundant ranges are the optical band and |$\gamma$|-rays: they demonstrate the presence of tens of strong and faint flares, both fast and slow. With emission energy falling down, these flares become wider, exhibiting a time delay with respect to the optical/|$\gamma$|-ray flares, the greater as the wavelength grows.

As a result, for the light curves at frequencies of 5 and 8 GHz, we can see only few strongest of them: the fainter events that leave only slight smooth waves in the total curve. For example, the delay of the most powerful outburst (⁠|${\rm MJD} \sim 54730$|⁠, 2008 September for |$\gamma$|-ray and optical emission) increases from 30 d for radio frequencies of 37 and 230 GHz to 70 d at 5 and 8 GHz.

The amplitudes of the most powerful flares vary through radio bands within 20–30 per cent of the mean values with only one exception: the leftmost outburst in Fig. 4 (MJD near 54200, 2007 April). Its amplitude, being near 5–6 Jy at high frequencies (230 and 37 GHz), becomes 3 times lower at 8 GHz and is only slightly noticeable at 5 GHz. A similar behaviour of flare components may be noticed for the group of outbursts at epoch |${\rm MJD} \sim 57000$|–57300, 2015 January–September. A complete analysis of event delays over different wavebands is presented in Section 4.3.

The behaviour of the flare in epoch 1 differs from the flares in epochs 2–4 by a steeper decrease of its amplitude with frequency decrease. As a result, the flare practically does not appear at 2 and 5 GHz in epoch 1. This behaviour corresponds to an optically thick emission area, where the spectrum is inverted at 2–22 GHz (⁠|$\alpha _{2-5}=+0.59$| and |$\alpha _{5-22}=+0.21$| in the epoch around |${\rm MJD}= 54267$|⁠).

Within the framework of the ‘generalized shock model’ of Valtaoja et al. (1992), which supposes the presence of growing and decaying shocks in relativistic jets, we can attribute this outburst as ‘high-peaked,’ and the others should be attributed as ‘low-peaked’ flares. Among possible explanations, there are higher magnetic fields or a smaller size of the emitting region at the moment of ‘high-peaked’ flare initiation.

Also in the case of this flare in epoch 1, the relevant explanation is the possibility of higher absorption of low-frequency radio emission because of a larger region where it originates compared to the flares at other epochs, which also appears as a more stretched variability time-scale. The analysis of AO 0235|$+$|164 long-term variability in Kutkin et al. (2018) found a similar behavior of some flares observed in 1980–2012, which corresponds to extreme compactness of the blazar core, although the flare in epoch 1 stands out between them as possessing the steepest spectral index.

4 ANALYSIS OF LONG-TERM VARIABILITY

4.1 Fractional variability

Quantification of the degree of variability in each band can provide information about the particle population that dominates in the emission in a certain energy band. For this purpose, we used fractional variability |$F_{\rm var}$| (Vaughan et al. 2003):

where |$V^2$| is variance, |$\bar{S}$| is the mean flux density, and |$\sigma _{\rm err}$| is the root mean square error. The uncertainty of |$F_{\rm var}$| is determined as

Fractional variability was computed for the light curves shown in Fig. 2 for the whole period as well as for each of four flaring states and the low-state period. The |$F_{\rm var}$|–|$\log \nu$| plot is given in Fig. 5. The |$F_{\rm var}$| errors are from 0.2 to 6 per cent in different wavelength bands. Fractional variability increases with increasing frequency for the whole time period and for each selected epoch.

|$F_{\rm var}$| reaches maximum levels in all bands in epoch 2 (above 150 per cent for the optical and |$\gamma$|-ray bands). In the low-state epoch fractional variability is from 7 to 20 per cent at frequencies of 2–230 GHz and reaches 35 and 52 per cent in the optical and |$\gamma$|-ray bands.

4.2 Structure function analysis

The SF is a method of searching for typical time-scales and periodicities in non-stationary processes (Simonetti, Cordes & Heeschen 1985; Hughes, Aller & Aller 1992). SF analysis provides time variability quantification and gives information on the nature of the process that caused the variations. A characteristic time-scale in a light curve, defined as the time interval between a maximum and an adjacent minimum or vice versa, is indicated by a maximum of the SF, whereas periodicity in the light curve causes a minimum (Heidt & Wagner 1996).

The SF of the first order normalized to the variance of the signal |$\sigma ^{2}$| is usually determined as

where |$f(t)$| is the signal at a time t, and |${\tau }$| is the lag. The slope of the power part of the curve is estimated as

An ideal SF consists of two plateaus and an upward-sloping curve with a slope b between them. One of the important characteristics of the SF is |$\tau$|⁠, the point where the SF reaches its higher plateau with an amplitude equal to 2|$\sigma ^{2}$|⁠. This time-scale gives the maximum time-scale of correlated signals or, equivalently, the minimum time-scale of uncorrelated behaviour. The value of the slope b determines the nature of the variable process. If the light curve can be modelled as a combination of white (flicker) and red (shot) noises, then the slope value is between 0 and 1 (Hughes et al. 1992). For a single dominating outburst in the light curve, the slope is usually steeper than 1. If the slope |$b \sim 2$|⁠, then there is a strong linear trend or a strong periodic oscillation (Hufnagel & Bregman 1992).

In this work, for the non-uniform and finite data series |$f(i)$| with |$i=1,2,...,N$|⁠, obtained from the original temporally non-uniform observed data series by partitioning into intervals with the values averaged over the interval, the first-order SF is calculated as

where |$N_{1}(k)=\sum {w(i)w(i+k)}$| and the weight factor |$w(i)$| is equal to 1 if the observations exist in the interval i and is 0 if there are no observations.

For |$k=1,2,...,L$|⁠, we construct our own set of intervals. The initial time interval |$k=1$| is chosen in such a way that it is equal to or slightly larger than the average time interval between observations (ignoring very long gaps in them). For radio frequencies of 2–22 GHz, intervals between 10 (for the low-state epoch) and 50 (epochs 1–4) d were taken as the initial ones. The corresponding intervals for 37 GHz, 230 GHz, and optical data lie between 4 and 20 d, whereas |$\gamma$|-ray data are described by an SF with an initial time interval of 30 d. The final interval |$k=L$| is calculated based on the length of the time-scale of the original sequence so that all the values of the original series are divided into two intervals.

We estimated the uncertainties on the SF through an independent bootstrap technique. A model light curve was generated for each source and frequency by applying a 2-d boxcar smoothing and a subsequent averaging process to the original light curve. The resulting 2-d model was then subtracted from each light curve, yielding a set of residuals. To preserve the time sampling of the original data, a new simulated light curve was created by adding a randomly selected residual to each point in the model. This process ensured that each residual was utilized once. Repetition of this procedure 1000 times allowed us to generate multiple simulated light curves, allowing the SF to be recalculated each time. The confidence limits were then determined, such as the 99 per cent confidence limit on a given time-scale, which was identified as the point where only 5 data points deviated from the value. This approach does not mandate symmetric confidence limits.

The computed time-scale values for the four epochs and for the low-state epoch are given in Table 8 and in Figs A1–A5. The dashes in Table 8 for some frequencies denotes the absence of a clear structure in the obtained SFs. The absence of the variability scale at 37 GHz for the low state epoch is due to the impossibility to determine its value reliably. We have indicated its lower limit (⁠|$\ge 0.5$| yr) as well as for epoch 3 at 37 GHz (⁠|$\ge 0.9$| yr).

As it was mentioned by Emmanoulopoulos, McHardy & Uttley (2010), results of SF analysis must be interpreted with some caution. Strong SF breaks frequently occur in data sets without any sort of characteristic time-scales. The position of these artificial breaks depends on the length of the observations and the nature of the variation process. On the other side, intensive modelling shows that SFs derived from typical blazar observations resemble those coming from the shot-noise model. However, in the frequency domain the same model does not describe correctly the observations and thus may not be physically realistic.

Also it was emphasized in the paper cited above that the data gaps affect severely the SF estimates in an unpredictable way, introducing systematic deviations. Even the bootstrap method cannot yield statistically meaningful errors depicting the deviations between the gappy and continuous SFs. This effect as well as the presence of some powerful outbursts may be responsible for difficulties with SF estimates in the optical data, where SF was successfully built only at the quiet epoch between 2009 and 2014.

The SF analysis reveals different variability time-scales: for epochs 1 and 2 the |$\tau$| value is between 0.3 and 1 yr; for epochs 3 and 4 |$\tau$| is almost two times higher, between 0.9 and 1.7 yr, respectively. The variability time-scale grows steadily from |$\gamma$|-rays (about 0.4–0.5 yr) to longer radio waves (between 0.9 and 1.7 yr for different epochs), but significantly differs from one epoch to another. The variability time-scales for the quiet epoch are similar at all wavelengths: |$\tau =0.4$|–0.5 yr.

The given trends of |$\tau$| are consistent with the findings of (Raiteri et al. 2001; Fan et al. 2002), which also indicate time-scales of up to a year that are related to accretion processes and jet activity.

4.3 Multiband correlations

Cross-correlation functions (CCFs) are commonly used in the multifrequency analysis of light curves observed over approximately the same period. The CCF analysis described in many digital signal processing texts is usually applied to the regular time series, which is not always possible in the case of astronomical observations. To solve this problem, discrete correlation function (DCF) method for irregularly sampled time series began to be used several decades ago to find the time lags between the light curves at different frequencies (Edelson & Krolik 1988).

In this paper the DCF method was applied in the same way as in our earlier work (Vlasyuk et al. 2023). The python-based software described in Robertson et al. (2015) was used to calculate cross-correlation functions between all the frequency bands mentioned above. Bin widths |$DT=10$|–40 d were taken.

To estimate the confidence level of a resulted DCF, we utilized Monte-Carlo simulations. The method is outlined in Emmanoulopoulos, McHardy & Papadakis (2013), where related software is also described. The point is to simulate a large number of artificial light curves that have the same statistical properties (probability density function, power spectrum density) as the real light curve has. To find the probability of getting a given DCF value purely by chance, one can calculate cross-correlation functions for an ensemble of paired artificial light curves and analyse the DCF distribution for each lag.

Since the quality of the data was quite different at various frequencies and in different epochs, sometimes a DCF had very large errors that made impossible to assert confidently its values, such cases were skipped. The DCF examples are presented in Figs A6–A9.

The |$2\sigma$| confidence level DCFs between the pairs of frequencies and the corresponding time lags over the whole period and for epochs 1, 2, 3, and 4 are displayed in Table 9. For this analysis we changed the definition of the epochs due to the poor number of data points for the first and second flares and complicated DCF calculation. Therefore, for epoch 1, we combined the first two flares, the second epoch is the low state, the third and fourth epochs are the same.

The analysis of the DCF data shows different time delays and behaviour of multiband variability for epochs 1–4. Evidently, all the analysed light curves exhibit flux variations, which are typical for blazars, where the flux changes at higher frequencies lead those at lower frequencies, as it was noted in Raiteri et al. (2001), Gaur et al. (2015), and references therein. For AO 0235 |$+$| 164, the value of such lags differs for epochs under analysis. A minimum delay, 0–40 d, is observed in epoch 1 with the two brightest flares for the pairs ‘37–2 GHz’ and ‘5–2 GHz.’ The time lags for epochs 3 and 4 are higher, hundreds of days, and approach 450 d for the pairs ‘230–5 GHz’ and ‘230–2 GHz.’ This observational result confirms the scenarios where time lags at lower frequencies are due to different opacities of the matter, as greater synchrotron self-absorption could occur at the lower frequencies compared to the higher ones.

4.4 Search for long-term periodicity

The Lomb–Scargle (L–S) periodogram method (Lomb 1976; Scargle 1982) is widely used in astronomy to search for periodicities in unevenly sampled time series, like the light curves. The result of this method is a periodogram similar to the power spectrum estimate in the standard Fourier analysis, where the most prominent peaks show the periods of oscillations in given data.

To find the long-term periodicity in our data, we used the astropy L–S implementation³ based on the algorithms described in detail in Ivezić et al. (2014). For a cross-check, we also used the generalized L–S periodogram module from pyastronomy,⁴ implemented as described by Zechmeister & Kürster (2009), and found that these two implementations give very similar results.

The interpretation of the L–S method results should be made with care. First of all, special attention should be paid to the window function. In most cases it has a shape of a ‘comb’ with a minimum interval of 1 d between the successive observations and with large gaps when the source is not observable. The convolution of the window function with the true light curve spectral image in the spectral domain can lead to false peak detection.

It was found that for our data the most prominent window function feature was the peak around a period of 365 d in the R-band light curve, caused by the yearly nature of time intervals when the source is unobservable during daytime.

At other frequencies the window function produces mild comb-like features at the scales of interest (e.g. at long-period scales).

The L–S periodograms show peaks with significance above |$2\sigma$| (with the exception of |$\gamma$|-rays and the 37 and 230 GHz radio bands; Table 10, Columns 2–3, Fig. A10) related to periods of 5.0–6.9 yr for all the bands, with a mean value of about 6 yr. With the exception of three frequencies, 2, 22, and 37 GHz, other light curves give very close periods between 5.5 and 6.4 yr.

Although the QPO of about 6 yr appears in the periodogram, the time span covers fewer than four full cycles for the optical and radio bands and fewer than three cycles for the Fermi-LAT data. This presents a challenge for making a robust claim of a periodic signal. Additionally, the influence of red noise has been considered, as it significantly affects the periodogram’s power at low frequencies. To address this, we conducted a simulation of red noise (following the procedure described in Roy et al. 2022) and corrected the periodogram by subtracting the mean red noise power. The results indicate that the detected peaks, while present, should be interpreted with caution.

This is consistent with the periods found by Raiteri et al. (2001, 2006) (quasi-periodicity with periods of about 5–6 yr in the optical and radio bands) and Fan et al. (2002, 2007) (about 5.8 yr in the optical and radio bands), but differs from the newest results with periods between 8 and 9 yr (Roy et al. 2022; Otero-Santos et al. 2023).

In order to detect the hidden periodicity for AO 0235|$+$|164, we performed an analysis of its light curves within the ‘quiet’ epoch that was not affected by the influence of powerful outbursts. The L–S periodograms reveal the presence of peaks with significance above |$2.5\sigma$| for almost all spectral bands except |$\gamma$|-rays and the 22 and 230 GHz radio bands(Fig. A11). All of the periods lie in the range from 1.3 to 2.3 yr (Table 10; Columns 4–5), which may indicate the presence of a single quasi-periodic process responsible for the manifestation of brightness variability throughout the studied epoch. The variation of the periods over the EM spectrum may be explained by a more significant effect of the |$1/f$| red noise at lower periodogram frequencies, where the peaks concentrate.

This result is in good agreement with the conclusions of Tripathi et al. (2021), who analysed the AO 0235|$+$|164 radio emission variability over three decades. The authors did not find significant variability at periods of 6–8 yr, but identified periods with a time-scale from 1.8 to 2.6 yr.

5 VARIABILITY AT SHORTER TIME-SCALES

Using the 2 and 5 GHz data obtained in daily observations from 2021 May to 2022 June (Fig. 6), we searched for the short-term temporal and spectral variation. The observations coincide with the period of the prolonged activity in epoch 4 (Fig. 2) and cover 376 d. The spectral index |$\alpha _{2-5}$| correlates with the 5 GHz flux density and varies from |$+0.20$| to |$+0.75$|⁠, which means that the corresponding emission area is optically thick over the whole period of observations in 2021–2022. We define the spectral index |$\alpha$| from the power-law |$S_{\nu }\sim \nu ^{\alpha }$|⁠, where |$S_{\nu }$| is the flux density at a frequency |$\nu$|⁠. The SF reveals a variability time-scale of the order of 100–120 d (Fig. 7) at 5 GHz, and similar value at 2 GHz (the low level). This probably indicates the presence of a shorter variability scale than it was determined for the long-term light curves in epochs 2 and 3 at these frequencies.

The DCF analysis indicates significant correlation between the 2 and 5 GHz flux densities with a peak |${\rm DCF}=0.65\pm 0.05$| (at a confidence level |$\gt 2\sigma$|⁠) and with a time lag closer to zero, |${\rm lag}\approx 0\pm 1.5$| d (Fig. 6). Using light curve interpolation, we obtain a similar result: |${\rm lag}\approx 7\pm 1$| d (⁠|${\rm DCF}=0.86\pm 0.06$|⁠, |$\gt 3\sigma$|⁠). The L–S periodogram at 5 GHz shows a few peaks at significance levels between 1.5 and |$2\sigma$|⁠, which may hint at QPOs with characteristic times of about 0.25 and 0.4 yr (Fig. 8).

6 BROAD-BAND RADIO SPECTRUM

In this section, we consider the properties of the quasi-simultaneous and average broad-band radio spectra of AO 0235|$+$|164.

6.1 Quasi-simultaneous radio spectra

To derive quasi-simultaneous broad-band radio spectra at 2–230 GHz (Fig. 9), we divided the observing period into equal time intervals of 30 d and averaged flux densities within these intervals. A total of 270 spectra have been obtained in such a way, although flux densities for some frequencies may be missing due to the differing gaps in the multiband light curve (see Fig. 2).

This sample of spectra was investigated using the method of cluster analysis (clustering) with an implementation similar to that described in Kudryavtsev et al. (2024). The analysis aims at identifying possible different states of non-thermal radiation properties during selected activity epochs. Cluster analysis allows one to combine investigated objects into groups based on their similarity, which is evaluated based on a number of natural or specially designed characteristics of the objects, named features. The features constitute axes of a multidimensional feature space, and the distance in this space, after proper scaling, is the metric of similarity. The obtained groups can be further investigated as different types/classes of the considered objects. Cluster analysis is a multiparametric technique that makes it possible to study a sample using the entire set of its available properties and their relationships.

To select the necessary features to form the clustering feature space, we first had to trim the quasi-simultaneous radio spectra sample to drop the dates earlier than MJD 54700 (roughly the year 2008.6) because a number of the observed frequencies had no measurements for these earlier epochs, particularly the |$\gamma$|-ray measurements that are pretty abundant for later dates. Having then the flux densities at different radio frequencies and in the optical R and |$\gamma$|-ray bands as well as radio indices |$\alpha$| calculated for different radio frequency pairs, we evaluated the multicollinearity and the percentage of missing values for particular characteristics. The rule of thumb is to drop the characteristics with more than 40 per cent of missing values and to choose one of the multicollinear characteristics (correlation greater than 0.7) with the least number of missing values.

As of the multicollinearity, we made an exception for the flux densities at radio frequencies of 5, 8, 11, 22, and 37 GHz. All these characteristics have the percentage of missing values lower than 40 per cent and are strongly correlated, but the gaps in the obtained measurements may be located at different epochs. The emission at radio frequencies is also roughly of the same physical nature. We took all this into consideration and, instead of dropping all but one radio flux densities, combined them into one metafeature, thus preserving a greater number of data points for the clustering. The metafeature, let us call it radio flux, was calculated as the first principal component in the standard principal component analysis (PCA) for the 5, 8, 11, 22, and 37 GHz flux densities. In order to perform the standard PCA, we, nevertheless, had first to impute missing values in the combined features. This was made using probabilistic PCA (PPCA; Tipping & Bishop 1999), particularly its implementation⁵ based on Porta, Verbeek & Krose (2005).

After performing the above mentioned steps of feature selection, we ended up with a total of five features in the clustering data set: combined radio flux |$S_{\rm radio}$|⁠, |$\gamma$|-ray flux |$S_{\gamma }$|⁠, and radio indices |$\alpha _{{11}-{37}}$|⁠, |$\alpha _{{11}-{22}}$|⁠, |$\alpha _{{5}-{11}}$|⁠. The number of the quasi-simultaneous spectra in the dataset after dropping those with missing feature values amounted to 187.

The following k-means clustering can be performed both in the original five-dimensional feature space or in a space of primary components after PCA dimensionality reduction. We performed both and obtained the same results. The scikit-learn library (Pedregosa et al. 2011) was used. The optimal number of PCA primary components describing the ‘true’ dimensionality of the data according to Minka’s maximum-likelihood estimation (Minka 2000) is two, so the clustering after dimensionality reduction was performed onto a simple two-dimensional plane.

As is seen in Fig. 10, the obtained clusters can be attributed to the low- and high-activity states of the source. Note that this result is obtained based not only on the observed fluxes, the three radio indices mentioned above influenced the clustering in equal measure. The figure shows the light curves in the |$\gamma$|-ray range and radio frequencies of 37 and 5 GHz. In the bottom panel, the variation the |$\alpha _{2-5}$| and |$\alpha _{5-37}$| radio indices are presented. During the high states, the spectrum becomes optically thick (⁠|$\alpha \gt 0$|⁠) and remains in this condition during several years.

The shapes of the quasi-simultaneous spectra of the high- and low-activity states are clearly different (Fig. 9). The distributions of flux densities and radio indices are shown in Fig. 11 as box plots. A box is the interquartile range (IQR, 25th to 75th percentile, or Q1 to Q3), the median of a distribution is shown as a vertical line inside the box, and the ‘whiskers’ extend to show the rest of the distribution, except for points that are determined as ‘outliers,’ locating beyond the median |$\pm 1.5$| IQR range, these outliers are shown by dots.

6.2 Average spectrum

The AO 0235|$+$|164 average spectrum is constructed based on the experimental data described in Section 2 and on the literature data taken from the CATS data base (Verkhodanov, Trushkin & Chernenkov 1997). The spectrum has a peak at millimetre waves in the observer’s frame of reference. AO 0235|$+$|164 is an extremely variable AGN and its broad-band radio spectrum (Fig. 12) show a strong variety of shapes: from steep (⁠|$\alpha \lt -0.5$|⁠) to rising (⁠|$\alpha \,\gt\, 0$|⁠). During the 27 yr of observations, the spectral index |$\alpha _{11-22}$| had been less than |$-0.5$| only during three epochs: February 2010, October 2011, and January 2023.

According to the Hedgehog jet model (Kovalev, Kovalev & Nizhelsky 2000), suggested by Kardashev in 1969, the strong longitudinal magnetic field B of the radio jet depends on the distance r as |$B= B_0\, (r/r_0)^{-2}$|⁠. In Fig. 12, the result of fitting the sum of two components in this model to the averaged observed data is given. Such average model spectra can be observed if the continuous flow |${\rm d}N(t)/{\rm d}t$| of emitting particles across the base of the jet is constant for a long time. The variability of |${\rm d}N(t)/{\rm d}t$| is converted in the model to the variability of the HFC and the total spectrum. A possible excessive emission at frequencies above 80 GHz is still hardly explained by this model and needs to be analysed in future studies.

Similar to Vlasyuk et al. (2023), we obtained the following fitted physical parameters of the radio-emitting jet (the HFC): the angle to the observer’s line of sight |$\vartheta \sim 1.0^\circ$|⁠, the flux |$S_m$| and frequency |$\nu _m$| of the jet spectrum maximum |$S_m \sim 0.53$| Jy at |$\nu _m \sim 60$| GHz, |$\gamma \sim 1.5$|⁠, and supposed |$\gamma _E \sim 300$|⁠. Neglecting the redshift, the next parameters can also be estimated (Sotnikova et al. 2022b):

Here, |$\nu _m$| (Hz), |$\lambda _m$| (m), and |$S_m$| (W m|$^{-2}$| Hz|$^{-1}$|⁠) are the frequency, wavelength, and flux density of the HFC maximum; the gamma-factor |$\gamma _E = E/(Mc^2)$| of emitting particles is supposed to have the same value for electrons and protons; |$B_\perp$|  = |$B \sin \vartheta$| (G), |$M_{2e} = 1$| for electrons, and |$M_{2e} = 1836$| for protons. |$\Theta$| (rad), |$k_B$|⁠, and |$T_b$| (K) are the angular diameter of the emitting jet in the picture plane, the Boltzman constant, and the brightness temperature. The estimated parameters are shown in Table 11.

The results in Table 11 have supported a new idea (Plavin et al. 2020; Kovalev et al. 2020; Sotnikova et al. 2022b; Kovalev et al. 2022) that replacing electrons with protons in the synchrotron radio emission of relativistic jets in the strong longitudinal magnetic field of an active nuclei may be an important assumption to describe both the nature of blazars and the generation of high-energy neutrinos via the high-energy part of the distribution of relativistic protons in blazars.

7 DISCUSSION

In this study we have enhanced the time coverage for the optical and radio data for the blazar AO 0235|$+$|164 up to 2023, achieving a 27-yr-long time period for the MW analysis. This period includes two multiband major flares: in 2007 (marked as epoch 1 in our study) and in 2015–2016 (epoch 3), which matches the |$\sim \!8.2$|-yr period described by the previous authors mentioned above. The flare in 2007 was followed by another massive outburst in 2009; the flare in 2015 was also followed by an outburst in 2016, but it was less prominent. In the latter case, that second outburst even merges with the major flare at 2.3 and 5 GHz. In addition to that quasi-periodic major flares, we spot other minor flares and flux density variations. At higher frequencies the structure of the flares is more complex and the amplitude variations are more notable compared to the low-frequency radio bands.

The variability time-scales |$\tau$| range from 110 to 620 d for individual epochs 1–4. In general time-scales decrease smoothly with an increase of frequency, and at each frequency |$\tau$| increases from epoch 1 to epoch 4. This reflects the differences in the physical properties during each flare, such as the size of the new compact components originating from the central region and observed in 2007, 2008, 2014, and 2020. We estimate the upper limit on the physical size of the emitting region (R) from |$\le \!0.9$| pc (for the 100-d variability time-scale) to |$\le$|5.5 pc (for 600 d), assuming

where c is the speed of light, |$t_{\mathrm{var}}$| is the variability time-scale, |$\delta =21$| (Kutkin et al. 2018) is the Doppler factor, and |$z=0.94$| is the redshift.

For comparison, we calculated the brightness temperature and associated Doppler factor of each major flare at 37 GHz according to equation (1) from Lähteenmäki & Valtaoja (1999) using the method of flare decomposition described in Section 3.1. We adopted the Hubble constant |$H_0 = 68$| km s|$^{-1}$| Mpc|$^{-1}$| and the intrinsic brightness temperature |$T_{\rm int} = 7\times 10^{10}$| K (Kutkin et al. 2018), which is comparable to the equipartition brightness temperature |$T_{\rm eq}~\backsimeq 5 \times 10^{10}$| K suggested by Readhead (1994). As we note from the results (Table 12), our estimated Doppler factors are quite conservative: for example, Kutkin et al. (2018) determined |$\delta \sim 21$| for the second flare, while we estimated |$\delta \sim 6$| for it. The difference arises from the adopted variability time-scale: we estimated the flare rise time-scale (⁠|$\tau =189$| d) from the exponential fits (Valtaoja et al. 1999), while Kutkin et al. (2018) obtained a much shorter time-scale of core variability of |$\sim \!33$| d. If we consider our Doppler factors, the physical size of the emitting region at 37 GHz is about |$R \sim 0.4$|–0.5 pc.

A cross-correlation analysis between MW light curves reveals that flares at longer wavelengths follow the short-wavelength flares, with different time delays during individual episodes. In the period of two major flares in 2006–2007, a 40-d lag is found for many pairs of the emission bands considered. The maximum lags up to 450 d have been measured for the two epochs in 2013–2019 and 2019–2023, and they correspond to the delay between the extreme radio frequencies of 230 and 2 GHz.

Cross-correlation between the optical and |$\gamma$|-ray flux variations is not so evident: DCF maxima for epochs 2–4 are rather high, lying between 0.66 and 0.86, but their significance does not exceed 2|$\sigma$|⁠. The time delay at this level is close to 0 d with an uncertainty of about 5 d, but with two major exceptions: the presence of a higher secondary peak (peak value |$= 0.86$|⁠) at 50 d for epoch 3 and a DCF peak equal to 0.74 at |$-70$| d for epoch 4. This result is caused by a complex structure of both light curves, with ‘sterile’ (optical without a |$\gamma$|-ray counterpart) flares at MJD |$\sim \!57050$|⁠, 59150, and 60200 and the lack of data in the optical light curve on yearly basis.

The same pairs of light curves (radio–|$\gamma$|⁠, radio–optical, radio–radio) show different time lags for various epochs, indicating a complex MW connection between them and, likely, the presence of several sources of EM radiation in AO 0235|$+$|164. Based on our calculated time lags, we can evaluate the distance between the emitting regions using the following equation (Pushkarev, Kovalev & Lister 2010):

where c is the speed of light, |$t_{\mathrm{lag,rest}}$| is the time lag in the source’s frame of reference, |$\beta \sim 10$| is the apparent angular jet speed, and |$\theta =1.7^{\circ }$| is the viewing angle (Kutkin et al. 2018). We obtain for a time lag of 10 d a distance |$D\sim 1.4$| pc, and for the greatest time lag of 450 d a distance |$D\sim 65$| pc.

The assumption for these estimates is that the synchrotron opacity in the nuclear region causes the lags, and different production zones in a jet for the radio, optical, and |$\gamma$|-ray emissions become observable at different distances. This theory is consistent with the findings that the flares in higher energy bands leads the lower energy flares (e.g. Pushkarev et al. 2010; Max-Moerbeck et al. 2014; Kramarenko et al. 2022), hence suggesting that the former emission originates upstream with respect to the latter emission. The cases with the shortest time lags (10 d) indicate almost co-spatial (within several pc) origin for some |$\gamma$|-ray and high-frequency radio flares (León-Tavares et al. 2011).

Epochs 3 and 4 show the typical behaviour of time lag for the |$\gamma$|-ray and radio band emission in quasars: the high-energy flare followed by the 230 GHz flare in about 10 d, which in turn is followed by the 37 and 22 GHz flares after 1 month, then the flare subsequently appearing at 11, 8, and 5 GHz in 6–7 months, and, finally, emerging after 10 months at 2 GHz. For example, Kramarenko et al. (2022) found the delays for the |$\gamma$|-ray and 15 GHz core emission of 331 AGNs to lie between 1 and 8 months with the most significant delays lying in the range of 3–5 months.

The relationships between time lag and frequency (Fig. 13) for 37, 230 GHz, R-band, and |$\gamma$|-rays are well fit for epoch 3 by the following straight lines with a negative slope:

For epoch 4,

The negative slope reveals that the lower-frequency emission lags that at higher frequencies. The quantitative estimate of these slopes is about |$-10$| d GHz|$^{-1}$| (five of seven better-defined slopes) and corresponds to the rate at which the lag decreases with increasing frequency. A possible reason of this trend is that synchrotron self-absorption is greater at low frequencies than at the higher ones due to different optical opacities. For the blazar 3C 279, Krishna Mohana et al. (2024) found a slope of about |$-30$| d GHz|$^{-1}$| at the radio frequencies. The differences can be explained by many reasons, mainly the physical conditions in a compact relativistic jet. The redshift of the object should also be taken into account for calculation of the ‘lag versus frequency’ relation in the observer’s frame of reference. Thus, the measurements of time lags in MW light curves give important information about the properties of AGN relativistic jets.

We have carried out an analysis of the multifrequency 1–230 GHz spectra obtained with the radio telescopes RATAN-600 at 1–22 GHz in 1997–2023, RT-32 at 5 and 8 GHz, and RT-22 at 37 GHz and the data taken from the CATS data base. As a result, we have estimated the physical parameters for AO 0235|$+$|164 based on the Hedgehog model of variable synchrotron radio emission in the strong magnetic field of the jet from an active nucleus of the blazar. An HFC with a flux about 0.3–0.5 Jy may emit at frequencies above 100 GHz near the jet base in addition to the jet model used. The results have supported a new idea that replacement of electrons with protons in the synchrotron radio emission model may be important both for the investigation of the nature of AGNs and for the generation of high energy neutrinos by relativistic protons in blazars.

7.1 Low-activity state in 2009–2014

During the extended period of low activity in 2009 November–2014 March, a clear variability in all wavelength bands was observed with the highest flux variations of 35 and 45 per cent occurring in the optical range and |$\gamma$|-rays, respectively. The fractional variability for all spectral ranges is factor 2–3 less compared to epochs 1–4 with the exception of |$\gamma$|-rays at epoch 4. We estimate almost same variability time-scales, |$\tau =150$|–180 d, at all wavelengths (Table 8). For several wavelength pairs, strong positive correlations are found (⁠|$\ge \!2\sigma$|⁠, Table 9, Fig. A7) with time lags of 100–200 d for |$\gamma$|–22 GHz and R–230 GHz, and with zero lags for the 230–37, 22–11, 11–8, and 8–5 GHz light curves. This can mean that the mechanisms that dominate the radio, optical, and |$\gamma$|-ray variations during the low state are not substantially different from those that are responsible for the emission during flaring states. Some indications of flux variation quasi-periodicity is detected in most of the wavelength bands in 2009–2014 (Fig. A11; Table 10) with main periods at 1.4–2.1 yr, which is much shorter than for the whole period of 1997–2023.

A similar result was obtained for the blazar Mrk 501 during a period of its extremely low broad-band activity (Abe et al. 2023). For this epoch the authors found significant flux variations in all wavebands with the highest values in X-rays and very high energy (⁠|$E\gt 0.1$| TeV) |$\gamma$|-rays. A strong correlation between |$\gamma$|-rays and radio with a time delay of 100 d was discovered. These properties are considered as the Mrk 501 basic emission properties.

7.2 Periodicity search: two periods may exist

AO 0235|$+$|164 is known for its MW emission variations on various time-scales, from hours to years. Continuous systematic observations of the source since 1975 in the radio (Aller et al. 1985) and optical (Spinrad & Smith 1975) bands allowed us to reveal long-term QPO patterns in its optical and radio light curves.

Based on the estimates of the |$\sim \!6$| yr period of multiwave variability for AO 0235|$+$|164, Raiteri et al. (2001) predicted subsequent outburst for it, but it did not occur (e.g. Raiteri et al. 2006, 2008). The persistent attempts to detect the quasi-periodicity of its brightness are explained by the idea that this object is a good candidate for the binary black hole scenario (Romero et al. 2003; Ostorero et al. 2004). The binary system would allow one to explain the |$\sim \!6$| yr variability time-scale by the precessing jet model, as it is shown in Escudero Pedrosa et al. (2024). The absence of the predicted flare as well as the fact that the quasi-period for AO 0235|$+$|164 varies significantly for different spectral ranges or analysed data temporal coverage (see Raiteri et al. 2001, 2008; Tripathi et al. 2021; Roy et al. 2022) need to be explained.

The L–S method used in the search for periodicity in light curves primarily takes the fluxes at the moments of increased brightness due to their greater statistical significance. Moreover, the first |$\sim \!6$|-yr period estimate was done by averaging the mutual distances between the peaks of five large-amplitude outbursts in the 8 GHz radio band throughout the analysed time period (Raiteri et al. 2001).

Thus, powerful flares in all ranges of the EM spectrum, which are largely stochastic in nature (Kushwaha et al. 2017), can mask the manifestations of the real periodicity of the brightness due to the regular properties of the active nucleus, like its binary nature or any features of nuclear structures (e.g. a disc or a torus). Obviously, the stochastic behaviour of the active nucleus itself may provoke regular processes in surrounding media, thereby regularizing the observed phenomenon. Moreover, the results of the DCF analysis may be affected by the lack of data for long time intervals, mainly in the optical bands because of the solar conjunctions or inappropriate weather conditions. Another feature of optical data is the presence of fast flux variations, which may affect the results of frequency analysis (VanderPlas 2018).

The reliability of our result is given by the fact that the analysis of the total LCs gives similar periods of variability from the data obtained in different ranges of the spectrum and with various instruments: 6.1 and 6.4 yr for |$\gamma$|-rays and the optical R-band and from 5.0 to 6.9 yr for all radio ranges, respectively. Thus, we have been able to get rid of systematic errors that affect the final result. As it is shown in Fig. A10, the significance of the L–S periodograms has decreased after applying the red-noise model. As a result, we can assign |$6.0\pm 0.5$| yr as the most robust but debatable quasi-period estimate. This value may be related to the quasi-periodicity of the outburst events in AO 0235|$+$|164.

At the same time, the 2-yr period we have suspected throughout the ‘low-state’ epoch can more reliably describe the intrinsic properties of the active nucleus and its environment. The detected periods, |$1.9\pm 0.3$| yr, are also mutually consistent, which strengthens our argument in favour of this result. As it is noted above, taking into consideration the red-noise influence has decreased the significance of the L–S peaks (Fig. A11).

Further note that a systematic analysis of the |$\gamma$|-ray AGN data measured by Fermi-LAT resulted in hints of periodic oscillations of |$\sim \!2$| yr in 7 BL Lacs and 4 FSRQs (Peñil et al. 2020). Supermassive binary black hole systems have been suggested as a plausible origin of periodicity (Prokhorov & Moraghan 2017). Indeed, Rieger (2007) suggested AO 0235|$+$|164 as a supermassive binary black hole candidate.

In conclusion, we express the idea that the general picture of the brightness variability in AO 0235|$+$|164 can be represented by a composition of the well-pronounced variations with an average repetition time of about 6 yr, probably due to the stochastic nature of the active nucleus radiation, and the latent variability with a 2-yr periodicity, which may reflect the real features of the nucleus internal structure and can be revealed reliably within its low-state periods.

8 SUMMARY

We report a comprehensive analysis of the MW behaviour of AO 0235|$+$|164 during its several activity epochs in 1997–2023, including four major flares and relatively low state in 2009–2014. We implemented the long-term broad-band radio measurements with RATAN-600, RT-32, RT-22, and optical observations with Zeiss-1000 and AS-500/2. Additionally we used the literature data at 230 GHz (SMA) and in the |$\gamma$|-rays (Fermi-LAT). From our study, we sum up the following conclusions:

• The inspection of the light curves by cross-correlation analysis has shown that time lag decreases with frequency in the range from 0 to 450 d. The relation between lag and frequency is described by a linear fit with a negative slope of |$-10$| d GHz|$^{-1}$| for the pairs of most radio frequencies versus 230 GHz, 37 GHz, R band, and |$\gamma$|-rays. This behaviour is common for many blazars, and we can assume that the higher energy emission originates closer to the central object upstream of the jet, whereas lower energy synchrotron emission is self-absorbed and becomes visible at farther distances downstream the jet. Additionally, the time lag differs from epoch to epoch for the same pairs of bands. This indicates a complicated relationship for multiband flares and a possible overlap of several emission sources for different flaring processes.

• A smooth decrease of variability time-scale from radio waves to |$\gamma$|-rays is observed for all the epochs. Typical values of |$\tau$| range from 0.3 to 1.7 yr. During the low state, the variability time-scale |$\tau =0.4$|–0.5 yr practically does not vary with frequency, suggesting a similar size of the non-thermal emission region for all the bands. The variability time-scale in 2021–2022, obtained using the data at 2 and 5 GHz with daily time resolution, is about 100–120 d, which is also seen for the epochs with longer time periods and caused by the appearance of a new compact component responsible for the flare in epoch 4 (2019–2024).

• For 4 yr of low activity in 2009–2014, we have measured clear variability throughout the EM spectrum. The highest flux density variations occurred in the R-band and |$\gamma$|-rays and reached 30 and 50 per cent, respectively. This variability might be caused by interaction between the remnants of the shocks followed from major outbursts. A correlation is revealed between the radio, optical, and |$\gamma$|-ray light curves for this period, which means that the mechanisms dominating the radio–|$\gamma$|-ray variations in the low state are not different from those during the active phase.

• The L–S periodogram reveals disputable peaks with close values, near 6 yr, for almost all of the wavelengths. However, because our data span less than 4 full estimated periods, we cannot assert confidently the presence of the 6 yr quasi-period. For the low state in 2009–2014, a shorter period of |$\sim \!2$| yr has also been suspected. The disputable quasi-period of 6 yr may reflects the time between the most prominent outbursts in the light curves. These flares have probably a stochastic nature, and the detected quasi-periodicity does not have tight connection with the characteristics of the active nucleus and relativistic jet. On the other side, we suppose that the periodicity found for the low state reflects general AGN properties.

ACKNOWLEDGEMENTS

The reported study was funded by the Ministry of Science and Higher Education of the Russian Federation under contract 075-15-2022-1227. The observations were carried out with the RATAN-600 scientific facility, Zeiss-1000 and AS-500/2 optical reflectors of SAO RAS, and RT-22 of CrAO RAS. The observations at 5.05 and 8.63 GHz were performed with the Badary, Svetloe, and Zelenchukskaya RT-32 radiotelescopes operated by the Shared Research Facility Center for the Quasar VLBI Network of IAA RAS (https://iaaras.ru/cu-center/).

YYK was supported by the M2FINDERS project which has received funding from the European Research Council (ERC) under the European Union’s Horizon2020 Research and Innovation Programme (grant no. 101018682). VAE and VLN are grateful to the staff of the Radio Astronomy Department of CrAO RAS for their participation in the observations.

SR acknowledges the support from the National Research Foundation, South Africa, through the BRICS Joint Science and Research Collaboration programme with grant no. 150504.

This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration; the CATS data base, available on the Special Astrophysical Observatory website; the SIMBAD data base, operated at CDS, Strasbourg, France. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. The development of the Fermi-LAT Light Curve Repository has been funded in part through the Fermi Guest Investigator Program (NASA Research Announcements NNH19ZDA001N and NNH20ZDA001N). The SMA is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. We recognize that Maunakea is a culturally important site for the indigenous Hawaiian people; we are privileged to study the cosmos from its summit.

DATA AVAILABILITY

The data underlying this article are available in the article and in its online supplementary material. The measured flux densities are distributed in the VizieR Information System.

The SMA data are available at http:// sma1.sma.hawaii.edu/callist/callist.html, the Fermi-LAT data are presented in public Light Curve Repository at https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository/about.html, the RATAN-600 data partly available in the BLcat on-line catalogue at https://www.sao.ru/blcat/.

Facilities: RATAN-600, Zeiss-1000, AS-500/2, RT-22, RT-32, Fermi LAT, SMA.

Footnotes

REFERENCES

APPENDIX A: PLOTS

See Figs A1–A11.