© The Authors 2024

1. Introduction

Blazars are among the most energetic objects in the universe. They are generally believed to consist of a super massive black hole (SMBH), referred to as the central engine, surrounded by an accretion disk and, in most cases, a dusty torus, as well as two symmetrical jets of matter emanating from the innermost vicinity of the black hole and the accretion disk. Particles in the jet are accelerated and collimated through a variety of mechanisms (a subject of numerous research studies currently underway), thereby reaching speeds close to the speed of light. This results in the highly energetic emission of radiation across the entire electromagnetic spectrum when these particles interact with the jet itself, the magnetic fields, and the surrounding medium. In the case of blazars, the jet is pointing towards us, thus bringing on some relativistic effects related to light aberrations. Such effects include light travel-time delays that lead to (apparent) superluminal motions or a Doppler boosting of radiation that makes them appear several orders of magnitude brighter than non-blazar jets.

Blazars usually present a spectral energy distribution (SED) with two bumps: the first extending from the radio to optical wavelengths or even X rays in the case of high synchrotron-peaked (HSP) blazars; the second extends from X rays, or γ-rays, to very high energy γ-rays. Synchrotron emission from the interaction of the charged relativistic particles of the jet with the magnetic fields in the medium is meant to account for the first bump. Several scenarios exist to explain the second bump. In the leptonic scenario, the second bump is explained by an inverse Compton effect of relativistic electrons interacting with ambient photons. A distinction is made based on whether these photons originate from the synchrotron emission inside the jet, in the case of which the mechanism is labeled as a “synchrotron self-Compton” (SSC). On the other hand, if the photon field is originated in a region external to the jet (typically, the broad line region or the dusty torus), the mechanism is labeled as “external Compton” (EC). There is an ongoing debate underway about the relevance of other different mechanisms, such as so-called “hadronic scenarios”. Frequently, combination of more than one emission mechanism is necessary to explain the observed SEDs and variability properties of the sources, even if the exact ratio of their contributions, as well as the origin and location of photon fields and particles involved are not sufficiently well-established.

The study of the variability of blazars across the spectrum, combined with the analysis of sequences of ultra-high-resolution very-long-baseline interferometry (VLBI) images, has proven to be an effective way of constraining the different emission models at work in these objects (Blandford et al. 2019). In particular, knowledge about the exact regions around the supermassive black hole and the relativistic jet where the γ-ray emission originates is essential in terms of discarding or supporting different models.

Regarding the location of the γ-ray emission, two main possibilities have been under discussion, differing with respect to the distance to the central black hole (BH). The first one is the so-called “close-zone” scenario, very close to the BH (0.1 − 1 pc), which was frequently used to explain the short time scales of high energy (HE) variability. However, this contradicts the coincidence of γ-ray and mm-wave outbursts that are associated to strong superluminal jet features seen in VLBI image sequences much further ( ≫ 1 pc) from the BH. In the second one, the so-called “far-zone” scenario, the emission region is located farther from the central engine, but multi-zone jet models are needed to explain the short time scales of variability reported at high- and very-high-energy γ-ray emission.

AO 0235+164 is an extragalactic BLLac-type blazar located at redshift z = 0.94 (Cohen et al. 1987). It shows strong variability across all the electromagnetic spectrum and has also displayed an interesting flaring behavior, with the most recent flares occurring in 2008 and 2015 studied with multi-wavelength (MWL) and VLBI. The source typically appears extremely compact at ultra-high-resolution 7 mm VLBI scales (showing the whole of the emission spanning < 0.5 mas) and kinematic and geometrical parameters obtained from VLBI images confirm a highly compact, narrow jet geometry pointing closely towards the observer’s line of sight with a very small opening angle (< 2.4°) at high speed (Doppler factor δ > 24). Altogether, this can explain the violent outbursts reported so far (Jorstad et al. 2001, Weaver et al. 2022). Agudo et al. (2011) reported a detailed analysis of all measurements available up to the 2008 flare, which we extend to 2020 in this paper. Here, we compare the two flaring episodes to shed further light about the origin and mechanisms involved in these extreme flares. The source has also been the subject of several previous observational campaigns that have produced light curves showing flares in previous years, for instance, 1992 and 1998 (see Raiteri et al. 2005). This points to the possibility of certain level of quasi periodicity with a characteristic time scale of about 6 years in the behavior of the source, which can serve as a guidance when developing models of the source, even when the data are not conclusive enough to settle this hypothesis, especially as examples of non-periodic wobbling of blazar jets exist (Agudo et al. 2012).

The 2008 flaring episode has received extensive coverage in the literature. Agudo et al. (2011) analyzed the flare from a multi-wavelength point of view, including the polarimetric data and VLBI imaging of the source. Their results favored a SSC scenario over EC to explain the γ-ray emission and constrained the location of the emitting region at > 12 pc from the central engine. Ackermann (2012) also analyzed the 2008 flare and produced a fit for the SED in the peak of the flare. In their model, EC was the dominating emission mechanism in the γ-rays. However, the EC mechanism fails to explain the observed variability and the correlations between γ-ray and optical emission. Baring et al. (2017) managed to reproduce the SED of AO0235+164 during the peak, including the X-ray excess. Wang & Jiang (2020) concluded in their study that the γ-ray and mm-wave emitting zones coincided within the acceptable errors and were located several parsecs from the central engine. These authors proposed a helical model for the jet to explain the observed polarization, without discarding other possibilities such as the shock-in-jet scenario. For this work, we have used a standard flat ΛCDM cosmological model with Hubble constant H₀ = 67.66 km Mpc⁻¹, as given by Planck Collaboration VI (2020).

2. Observations

We have obtained and compiled time-dependent data in most available ranges of the electromagnetic spectrum, including polarimetry whenever possible and VLBI polarimetric images with submilliarcsecond resolution.

Our observations include 7 mm Very Long Baseline Array (VLBA) images from the blazar monitoring program at Boston University, which were reduced both for total flux and polarization using AIPS (see Weaver et al. 2022 for details on the data reduction and calibration). Single-dish data at 1 mm and 3 mm were obtained from the POLAMI (Polarimetric Monitoring of AGN at Millimeter Wavelengths)¹ program at the IRAM 30 m Telescope (Agudo et al. 2017a,bThum et al. 2017). The optical (R-band) data were taken from Calar Alto (2.2 m Telescope) under the MAPCAT program, Yale University SMARTS blazar program, Maria Mitchell, Abastumani and Campo Imperatore observatories, Steward Observatory (2.3 and 1.54 m Telescopes), Perkins Telescope Observatory (1.8 m Telescope), Crimea Observatory AZT-8 (0.7 m Telescope), and the St. Petersburg State University LX-200 (0.4 m Telescope).

Ultraviolet (UV) measurements were obtained by the Swift-UVOT instrument. The data set also includes X-ray data in the 2.4–10 keV range from the RXTE satellite and in the 0.2–10 keV energy range from Swift-XRT. The light curves and spectral indices were derived from these data using a broken power-law model and the appropriate corrections for extinction. More details on the data reduction procedure from Swift is provided in Appendix B. The γ-ray data in the 0.1–200 GeV range come from the Fermi Large Area Telescope (LAT).

The ±180° polarization angle ambiguity in our R-band measurements was circumvented following the procedure described in Blinov & Pavlidou (2019), which minimizes the difference between successive measurements taking also into account their uncertainty. Clusters of close observations were then shifted by an integer multiple of 180° to match the angle reported at 3 mm. This allows us for a visual comparison of the joint evolution of the optical and millimeter range polarization angles.

Data from the infrared (IR) to the UV bands were corrected following the prescription by Raiteri et al. (2005) and the updated values by Ackermann (2012). This correction accounts for the local galactic extinction at z = 0 and the intervening galaxy ELISA at z = 0.524, as well as for ELISA’s contribution to the observed emission. When these corrections are applied, a UV bump appears in the final spectra for some epochs; this was shown in Raiteri et al. (2005), although the results are in disagreement with the SEDs presented in Ackermann (2012). It must be noted that applying different correction factors available (NED², Junkkarinen et al. 2004, etc.) also produce bumps (albeit of different intensity) but the UV bump is present in every case. Here we have followed Raiteri et al. (2005) when producing the final, extinction-corrected SEDs and used the updated values in Ackermann (2012) for the extinction factors, together with the magnitudes for ELISA reported by Raiteri. A comparison with the older values by Junkkarinen et al. (2004) can be seen in Fig. 15.

The correction of X-ray spectral data was performed using a single absorbed power law with density N_H = 2.8 × 10²¹ cm⁻² (Madejski et al. 1996; Ackermann 2012), which accounts both for galactic extinction and the z = 0.524 absorber. This value agrees with the result obtained when letting N_H vary as a free parameter.

3. Results

3.1. Millimeter, optical, and high energies

The light curves at millimeter wavelengths (VLBA 7 mm, 1 mm, 3 mm), optical bands (R, U, B, V), UV bands (UVW1, UVW2, UVM2), X-rays (0.2–10 keV), and γ-rays (0.1–200 GeV) of AO 0235+164 are presented in Fig. 1. Polarization degree and polarization angles at optical (R-band) and millimeter wavelengths (1 mm, 3 mm, and VLBA 7 mm) are shown in Figs. 2 and 3, respectively. Figure 1 shows that flaring episodes happen almost simultaneously across all the electromagnetic spectrum. Variability is much more pronounced at HE, while it is milder at optical and UV wavelengths and even more so in the millimeter and radio bands.

Fig. 1. Light curves of AO 0235+164 at different wavelengths. The vertical lines mark the epochs whose SEDs were analyzed in Sect. 4.4. The top panel shows both the X-ray and γ-ray fluxes in different axes (left for X-rays and right for γ-rays, respectively).

Fig. 2. Polarization degree evolution of AO 0235+164 at different wavelengths. Together with experimental data, a bayesian block representation is shown superimposed for R (black line) at a 99.9% confidence and for 1 mm (red) and 3 mm (blue) at a 90% confidence level. The vertical lines are the same as in Fig. 1 and correspond to the epochs whose SED was analyzed in Sect. 4.4.

Fig. 3. Polarization angle evolution of AO 0235+164 at different wavelengths. The vertical lines are the same as in Fig. 1 and correspond to the epochs whose SEDs were analyzed in Sect. 4.4. All points in the first box correspond to R band, the colors denote the clusters that were shifted by n × 180°, as noted in Sect. 3.2.

3.2. Polarization

The Bayesian block representation (Scargle et al. 2013) of the polarization degree light curves makes it easier to discern the different behavior between quiescent and flaring states (Fig. 2) because it represents significantly different evolution states of the source. The source exhibits lower polarization degree at both optical and mm wavelengths during the quiescent period in between flares (p_L, R = 9.5 ± 6.0%, p_L, 3mm = 2.5 ± 1.4 % from 2010 to 2014) than during flares (p_L, R = 14.5 ± 8.5%, p_L, 3mm = 3.34 ± 1.3% from 2014 to 2017). The 3 mm polarization angles also varies more slowly during the quiescent period: from 2010 to 2014, the polarization angle at mm wavelengths remains more or less stable, while from 2014 to 2017, it performs three full 180° rotations (as seen in Fig. 3). Rotations in the optical R-band also follow mm rotations, with a stronger variability, sometimes performing several 180° cycles while the 3 mm only varies a full cycle or a partial rotation. There is also an apparent delay of approximately a hundred days between 3 mm and R-band, as can be seen in Fig. 3.

The direction of the electric vector position angle (EVPA) of the VLBA components (indicated in Figs. 4–6 as black lines segments overlaid with the images) coincides with the momentary direction of the jet. This alignment is in agreement with the shock-in-jet model (Marscher et al. 2008), where the compression of the magnetic field in the plane perpendicular to the direction of propagation (slightly askew of the direction of the observer) makes the electric field align with the jet direction. This supports the association of the superluminal components ejected during flares with plane-perpendicular moving shock-waves.

Fig. 4. Selected epochs illustrating the evolution of each identified knot, showing total (contours) and polarized (color scale) intensity. The beam size is indicated as a green ellipse in the first row. Horizontal black lines indicate the position of the core A0, black line segments within the image indicate the direction of polarization (EVPA). The red line in each row is the linear fit to the knot position. It is present for every knot except B7, whose flux was too low to accurately perform a fit. For each row, the spacing between plots is proportional to time, and the total time span is different and indicated in brackets.

Fig. 5. Epoch 102 (2016-09-05) VLBI 7 mm image, showing total flux intensity (contours) and polarized flux intensity (color scale). Black line segments overlaid in the image represent the EVPA. The green ellipse in the lower left corner represents the beam size. The image showcases component B5 close to the peak of the 2015 flaring episode, demonstrating how the polarization angle is aligned with the direction of propagation.

Fig. 6. Epoch 110 (2017-04-16) VLBI 7 mm image. The image shows the trailing component B6 moving in the same direction of B5, maintaining the alignment of the polarization angle (black line segments) with the momentary direction of the jet (northeast).

3.3. VLBI imaging

Our study includes all available 7 mm (43 GHz) VLBA total flux and polarimetric images from the Boston University Blazar Group of the sourced from 2008 to 2020 (from the VLBA-BU-BLAZAR and BEAM-ME programs³). After reducing the data with AIPS (see Weaver et al. 2022), the most prominent jet features were fitted to Gaussian components with Diffmap and then cross-identified along the observing epochs. This was done for a total of 142 observing epochs.

The VLBA images show a compact, stationary component at all epochs, A0, referred to here as the “core”. Other features can be tracked at different epochs and their evolution is traced over time. Figure 4 shows some selected epochs with the identified knot features to give a general idea of the behavior of the source in time. Evolution curves in total and polarized flux intensity and the polarization degree for the total emission and single components were later produced from the images (Figs. 1–3) using the aforementioned identification.

The flux evolution shown in Fig. 1 at all wavelengths, also containing the light curves from the integrated VLBA 7 mm maps, allows us to distinguish two clear “flaring” periods, whose peaks of activity occurred in October 2008 and July 2015, respectively. The 2008 flare is associated with the B2 jet feature that developed southwest of the core (A0). In contrast, the 2015 flare is associated with jet components B5 and B6 that developed northwest. Other weaker components not associated with the main outbursts (e.g., B4), also propagate in different directions. This hints at a possible rotation or wobbling of the jet and supports a helical jet model, and might be associated to a pseudo periodic behavior as proposed by Raiteri et al. (2005). All VLBI jet components have lifetimes lasting several years. During their lifetimes, we observe them propagating quasi-ballistically in the same direction relative to the core, with trailing components maintaining the same direction of propagation as their leading component as well as for its EVPA alignment (parallel to the direction of propagation in the plane of the sky). It is therefore clear that jet nozzle changes direction with time since, in each flaring episode, the direction of propagation of the associated superluminal components is radically different for every on of these episodes.

3.4. Differences and similarities between 2008 and 2015 flares

The flare in 2008, reported by Agudo et al. (2011), featured a superluminal component ejected from the core during June 2008, which kept separating from the core in the southeast direction during the following months, until the last months of 2009, when the component went practically extinct. Agudo et al. also reported correlations between all wavelengths, from radio to γ-rays. The flare in 2008 reached a maximum in October 2008, peaking at 4.7 Jy in 3 mm and a magnitude of 14.2, having increased its brightness by more than a factor of 60 in optical R band with respect to its quiescent state. Although a general correlation was reported between optical and γ-ray bands, the correlation was poorer and less detailed during the main burst.

The second flaring activity begins approximately in autumn 2014, with a brightening of the core visible at 7 mm. The flux densities at all wavelength peak a year later, around autumn 2015, at all wavelengths. Comparatively, this flare is dimmer than the previous one of 2008, reaching 4.2 Jy in 3 mm and lower by 1.3 in magnitude in the optical R band. A plausible explanation for this will be given below through the interpretation of the flare in terms of emission zones and shocks.

The VLBI images reveal that the component responsible for the brightening of the core during the 2015 flare, named B5, originates in the core less than a few weeks before the brightening in total flux density begins. After it separates from the core, the component travels outwards at a position angle of ∼45° while increasing its brightness, and peaks around Fall 2016. Figure 5 shows the total and polarized intensity image of the source during nearest to this event. This is around a year later than the total flux peak. The interpretation for this is that the initial brightening of the total flux is due to the interaction of the component with the core, while the ensuing brightening of the component is due to acceleration or a change in the viewing angle.

A few months later another, a weaker component, named B6, originated from the core. This component follows the path of B5, peaking around April 2017 (Fig. 6), while B5 is still visible. This component is the one responsible for the subflare of 2017, visible at all wavelengths in total flux. The behavior of this component is compatible with that of a trailing component of B5 as described by Agudo et al. (2001). By 2019, all activity had ended, with a minimum reached around March 2019.

The two flares present different time profiles at γ-ray energies. A comparison of these profiles can be found in Fig. 7, where also the 7mm flux of each identified VLBI component is shown. We have modeled the γ-ray profiles by fitting to standard exponential shapes given by (Abdo et al. 2010)

$$\begin{array}{c}\hfill y{(t)}_{N_{exp}}={\displaystyle \sum _i^{N_{exp}}}{A}_i·{(e^{(t-{t_c}_i)/{t_r}_i}+e^{(t-{t_c}_i)/{t_d}_i})}^{-1},\end{array}$$(1)

Fig. 7. Fits to the profiles of the 2008 (top) and the 2005 (middle) γ-ray flares by exponential functions as described by Eq. (1). The selected number of terms Nexp were chosen so that the reduced χ2-statistic was minimized (bottom). The best-fit values are given in Tables 1 and 2. An alternative fit is given for the 2008 flare with a similar χ2, that accounts for the double peak at the beginning with a single exponential term. The vertical lines mark the epochs whose SED was analyzed in Sect. 4.4, as in Fig. 1.

which allowed us to derive the rising and decaying times of each subflare, t_ri and t_di. Then, their asymmetry factor could be computed, defined as:

$$\begin{array}{c}\hfill {\xi }_i=\frac{{t_r}_i-{t_d}_i}{{t_r}_i+{t_d}_i}.\end{array}$$(2)

An asymmetry factor close to zero corresponds to the case of a perfectly symmetric flare. There is some uncertainty in the number of exponential terms to use, since the source shows strong variability in timescales shorter than our binning allows us to track. The wide binning applied, (i.e., 7 days) was necessary to accommodate periods of low flux. Still, it can be seen that the source displays significant variations of flux even in these intervals. The value of N_exp chosen was the one that minimized the reduced χ²-statistic. The results of both fits are shown in Fig. 7, and the corresponding parameters in Tables 1 and 2. The distance from the fitted γ-ray subflare maximum to the 7mm maximum in 2015 is 52 days (±8 days), a similar delay to the one found in the DCF analysis in Sect. 4.2 (τ_R, γ = 2 days, τ_R, 7mm = (64  ±  4) days).

For the 2015 flaring episode, the secondary flares in γ-rays are contemporaneous with the appearance of 7mm VLBI components. In particular, the first, second and third maxima happen at approximately the same time A0 rebrightens and the B5 and B6 components appear. This suggests that these emissions are spatially related. The failure in finding components responsible for the subsequent fitted subflares might be due to the difficulty in fitting low-flux VLBI components, the uncertainty in the γ-ray light curve, or a combination of both.

The same can be said about the first and second subflares of the 2008 flaring episode, where it seems that the brightening of the core A0 and the appearance of the B2 component are related to the first and second maximums at γ-rays, taking into account the aforementioned delay.

The γ-ray subflare that can be seen in Fig. 8 occurring one year before the start of the 2014 flaring episode might be related to the B4 component in the same way, but this relation and the associated delay are less clear. This could potentially be the case also with the B1 component and the flare that can be seen in 2006 (before the 2008 episode) in all wavelengths except γ-ray (due to the lack of observational data). Altogether, it seems that there is a direct relationship between the appearance of the 7 mm VLBI components and the successive γ-ray subflares.

Fig. 8. Correlations between fluxes across all wavelengths. Horizontal lines show significance levels for 1σ, 2σ, and 3σ, computed using N = 2000 synthetic light curves, as described in Sect. 4.2. The DCFs here are computed using the whole period of available data, from 2007 to 2020.

For the 2015 flaring episode, the rising and decaying times shown in Table 2, which are as low as ∼18 days (taking into account only the two strongest subflares), are compatible with the sizes found in the SED modeling of Sect. 4.4 (Table 7), which sets the shortest timescale where significant variations of flux can occur at a limit of 15 days (see Eq. (8)). These other values were obtained only from the modeling of the SEDs and the two analyses are completely independent. In the case of the 2008 flaring episode, however, the shorter times (∼4 days) reveal some tension with the region sizes. This might be explained by unaccounted sources of γ-ray variability originating in smaller regions. However, it could also be caused by a wrong estimation of the rising and decaying times. It is possible to produce a fit with a single term accounting for the initial double peak that has a similar χ²-square statistic, but rising and decaying times of 18 and 36 days (dashed line in Fig. 7). In any case, the observed delays between γ and mm might be explained by a combination of adiabatic expansion and cooling times (Tramacere et al. 2022).

4. Analysis

4.1. Kinematic parameters of the VLBI jet components

From the VLBI imaging data, some kinematics parameters associated to the different visible emission zones were computed following the procedure described in Weaver et al. (2022). These include t₀, the ejection time, which is the time where the extrapolated trajectory of the component crosses the core; t_var, the timescale of variability, which is the timescale of the dimming of the component; β_app, the apparent speed in units of c; δ_var, the variability Doppler factor; Γ, the Lorentz factor; and Θ, the viewing angle of the jet component.

The identified knot features in every epoch were traced and their positions adjusted to a linear fit, from which their speeds were obtained and their flux was also fitted to a decaying exponential function of F = F₀exp( − t/t_var), obtaining a timescale of variability (Fig. 9).

Fig. 9. Observed distance from core (top) and flux density (bottom) for every one of the identified component as a function of time, together with a linear weighted fit to knot distance and logarithmic fit to flux. The fit to the flux is done taking into account only the points after the peak of emission of each component. This was done for all components except B7, which due to its low flux, did not have enough points after the peak with low enough uncertainty to perform a fit.

The Doppler factor and apparent speed were then computed as (Jorstad et al. 2005, Casadio et al. 2015):

$$\begin{array}{cc}& {\delta }_{\mathrm{var}}=\frac{1.6a_{S_{max}}{d}_L}{c{t}_{\mathrm{var}}(1+z)},\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\beta }_{\mathrm{app}}=\frac{v_r{d}_L}{c(1+z)},\hfill \end{array}$$(4)

where a_Smax is the FWHM of the component measured at its maximum flux, v_r is the radial velocity of the knot, and d_L is the luminosity distance, but following the more robust approach found in Weaver et al. (2022) and using the value of t_var obtained from the fit. From these, the Lorentz bulk factor,

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{1}{2{\delta }_{\mathrm{var}}}({\beta }_{\mathrm{app}}^2+{\delta }_{\mathrm{var}}^2+1),\end{array}$$(5)

and the viewing angle,

$$\begin{array}{c}\hfill tan\mathrm{\Theta }=\frac{2{\beta }_{\mathrm{app}}}{{\beta }_{\mathrm{app}}^2+{\delta }_{\mathrm{var}}^2-1},\end{array}$$(6)

can be computed.

Our results for these parameters (Table 3) agree with those of Weaver et al. (2022) within the expected margin of error associated with the identification of the components in the VLBA images.

The results agree with the observed behavior of the flares. The estimated viewing angle for the component responsible for the 2008 flare (B2) is 0.2°, between three and four times smaller than the 0.7° of the component responsible for the 2015 flare (i.e., B5). This consistently explains the lower brightness observed in 2015 as being caused by a weaker Doppler boosting of the emission. The viewing angle for the secondary component B6 is similar to the one of B2; however, it is also apparent that its speed is much lower than for any of the others.

4.2. Correlations across the spectrum

Correlations among the different light curves were computed using MUTIS⁴. In particular, since we are dealing with irregularly sampled signals (light curves), we computed the discrete correlation function (DCF) following the method from Welsh (1999), which is a normalized and binned DCF.

A uniform bin size of 20 days was used for all correlations. The choice of a uniform bin size was done so that the results of different correlations could be easily compared, the specific value of 20 days was done so that it was large enough to have statistics in any bin, but short enough that the correlations were not overly smoothed and peak positions could still be determined. To confirm the robustness of our choice, we also reproduced our analysis with bin-sizes from 10 to 30 days, obtaining similar results; with the exception of some bins that did not have enough points to compute the correlation, as we discarded bins where the number of pairs was less than 11.

To estimate the significance of the correlations, a Monte Carlo approach was used, generating N = 2000 synthetic light curves for each signal. The randomization of the Fourier transform was used for mm-wavelengths, while for the optical and γ-ray data, we modeled the signals as Orstein-Uhlenbeck stochastic processes (Tavecchio et al. 2020). The uncertainties of the correlations were estimated using the uncertainties of the signals again with a Monte Carlo approach.

4.2.1. Correlation of the whole period (2007–2020)

The results from DCFs of the whole available period of data (2007–2020) show a clear correlation between the flux at almost all wavelengths (> 3σ) with most peak positions close to zero (Fig. 8). The X-ray band is an exception to this, showing no significant (> 3σ) correlation close to zero with some of the other bands. More hints about this will be provided in the following.

The correlation between the polarization degree and total flux is clear for the R band, where it shows a statistically significant maximum near zero, but it is not certain in the other bands. This may possibly be explained by the sparser sampling and larger errors (Fig. 10).

Fig. 10. Correlations between fluxes and polarization degree across all wavelengths. Horizontal lines show significance levels for 1σ, 2σ, and 3σ, computed using N = 2000 synthetic light curves, as described in Sect. 4.2. The DCFs here are computed using the whole period of available data, from 2007 to 2020.

4.2.2. Correlations of flaring episode (2014 and 2017)

The results from DCFs of the flaring episode (2014 to 2020) again show significant (> 3σ) correlation between flux at all wavelengths, with most peaks positions close to zero (Fig. 11). The clear exception to this general correlation is again the X-ray band. The absence of > 3σ correlation close to zero for the X-ray emission with the other bands (Fig. 11) hints at other emission mechanisms located in a different emission zone. This suggests that a different, separated processes could be responsible (at least partially) for the emission in the X-ray. This hypothesis was also favored by the analysis of the spectral energy distribution of the source (Sect. 4.4).

Fig. 11. Correlations between fluxes across all wavelengths. Horizontal lines show significance levels for 1σ, 2σ, and 3σ, computed using N = 2000 synthetic light curves, as described in Sect. 4.2. The DCFs here are computed using only the flaring episode from 2014 to 2017.

The interpretation of the peak positions in the DCFs is not straightforward. A debate on how accurately they represent the timing between different emission episodes is ongoing. This is especially the case when longer periods are taken into account, as a greater number of different processes and regions can be involved in the correlation. For the DCF of the whole period, the correlation only tells us about the probability that the processes causing the emissions are related. However, if we consider only the flaring episode, the relation between the peak position and timing of emissions will be more direct. Even then, the presence of several correlation peaks makes the interpretation of the results difficult. These peaks are the consequence of the low, non-uniform sampling of available data, and the complex structure of the flares. This is a fundamental flaw of any correlation analysis, since this correlation noise might result in peaks that do not correspond to the real delay between the signals. Several ways of dealing with these have been proposed, such as using the centroid instead of the maximum, but they are not free from biases and flaws, such as those discussed in Welsh (1999). Here, we follow Welsh (1999), using just the absolute maximum of the DCF and justifying the decision with the consistency of our results, as shown in the following.

If the position of the peak is to represent the real delay between the signals at different wavelengths, these delays should be more or less compatible between themselves when computed using different sets of correlations, for instance, the delay between A and B plus the delay between B and C should be close to the delay between A and C. In this sense, it is possible to build a compatibility chart, showing the relations between the different positions.

This is shown in Fig. 12 (Table 4) with the DCFs computed for the signals between 2014 and 2017 (the flare period, Fig. 11) and the band R as a reference. Our choice of R as the reference band is motivated by the fact that is the most densely sampled band during the periods of high variability. In this graph, each row corresponds to a band, i. The delays or peak positions between the row band i and any other band j, ${\tau }_{\mathrm{p}}^{i,j}$, are plotted along the x-axis, shifted by the delay between the band, i, and the reference band, ${\tau }_{\mathrm{p}}^{i,R}$, so that they fall aligned on the same positions.

Fig. 12. Correlations compatibility chart for τp. Each row corresponds to a band i. The delays or peak positions between the row band i and any other band j, ${\tau }_{\mathrm{p}}^{i,j}$, are plotted along the x-axis, shifted by the delay between the band i and the reference band, ${\tau }_{\mathrm{p}}^{i,R}$. The peaks selected are the one with highest values of the DCF, marked with red dots in Fig. 11. The DCFs here use the flaring episode from 2014 to 2017.

We see that indeed the positions fall more or less aligned in most cases, justifying our interpretation, and our choice of the maximum. The points for 1 mm and γ are very dispersed, but the correlation between R and γ presents a prominent peak, with high confidence and without spurious peaks close (Fig. 11), so we take ${\tau }_{\mathrm{p}}^{R,\gamma }\sim +2\mathrm{days}$ as the correct delay, τ_R, γ.

Discarding the 1 mm and γ rows, we estimated that the mean delays with respect to R for 7 mm, 3 mm, and X-ray, ${\tau }_{i,R}=⟨{{\tau }_{\mathrm{p}}}^{i,R}⟩$ are 64 ± 4 days, 42 ± 6 days, and 73 ± 4 days, respectively. The delays for 7 mm, 3 mm, and R are consistent with expectations if the mechanism of emission is synchrotron cooling, which should result in delays of the form ${\tau }_s\propto {\nu }_s^{-1/2}$, where ν_s is the synchrotron frequency, as can be seen in Fig. 13.

Fig. 13. Fit to the function y = ax−1/2 + b, which is the expected form if the time delays are due to different synchrotron cooling times. Fit results are a = 1.341 ± 0.009 × 107, b = −6 ± 1 × 10−1, p(χ2) = 0.909, and r2 = 0.9998. Parameter b is close to zero and accounts for an arbitrary reference delay, which is R in our case.

The delay obtained for X-ray with respect to R is much larger than for any other band. This strengthens the hypothesis that emission at X-ray energies might involve a different mechanism, as is already suggested by the lower level of correlation found and as the analysis of the spectral energy distributions revealed (see Sect. 4.4).

4.3. Geometry of the emission regions

We can also compute the corresponding sizes implied by the variability timescales, since they are constrained due to causality and special relativity, according to the following formula:

$$\begin{array}{c}\hfill d=\frac{c\beta \mathrm{\Delta }t}{(1-\beta cos\theta )(1+z)}=\frac{{\beta }_{\mathrm{app}}c\mathrm{\Delta }t}{(1+z)sin\theta }.\end{array}$$(7)

This same formula can also be used to compute the relative distances between emission regions implied by the time delays obtained in the correlation analysis, under the hypothesis that they result from the distances (although this may not necessarily be the case if they arise from synchrotron cooling, as shown in the previous section).

The sizes of the emitting region can be constrained with:

$$\begin{array}{c}\hfill R_b=\frac{c{t}_{\mathrm{var}}{\delta }_D}{1+z}.\end{array}$$(8)

For the moving component B2 corresponding to the 2008 outburst, using the timescale of variability in Table 3 and the δ_D ≃ 67, we obtain sizes of around 2 pc, consistent with the angular measure of VLBI images. Using the γ-ray variability one obtains much lower sizes, of around 0.2 pc, since variability at these energies is observed in timescales as short as 8 days (Ackermann 2012). This smaller high-energy emitting region is in agreement with the expected result of synchrotron cooling in the proposed model which explains the longer duration of flaring activity in mm. The maximum viewing angle of this jet is cited to be ≲2.4° (Agudo et al. 2011). Via the relations between this angle and the true speed and Doppler factor, namely,

$$\begin{array}{cc}& \beta =\sqrt{1-1/{\mathrm{\Gamma }}^2},\hfill \end{array}$$(9)

$$\begin{array}{cc}& {\mu }_s=\frac{1}{\beta }(1-\frac{1}{\mathrm{\Gamma }{\delta }_D})=cos\mathrm{\Theta },\hfill \end{array}$$(10)

and Eq. (8) limits us to a minimum of 1 pc the sizes of the mm-emitting regions.

The relative distances of the core and knots to the base of the jet can be ascertained using a model for the geometry of the jet. Following Wang & Jiang (2020) and using a conical geometry, we have:

$$\begin{array}{c}\hfill r_{\mathrm{core}}=\frac{r_{\perp }}{\phi }=\frac{0.5{\theta }_d{d}_L}{{(1+z)}^2\phi },\end{array}$$(11)

where φ is the half-opening angle of the jet and θ_d is the angular diameter.

With our knot identification, we can estimate the half-opening angle as φ =(Θ_0,max − Θ_0,min)/2 ≃ 0.4°. However, this way of estimating the half-opening angle is very sensitive to the weakest components. A second way to estimate this angle is (Weaver et al. 2022) φ = θ_p sin Θ₀, where θ_p is the projected opening semi-angle of the jet and is taken to be twice the standard deviation of the jet position angle, or of the visible component in the case of a wobbling jet direction. With our parameters, this gives about 0.78°, closer to the more widely cited (Weaver et al. 2022, Wang & Jiang 2020) value of about ≃1° for B2, the brightest component and the responsible for the 2008 flare.

For a core size at 43 GHz of θ_d ≃ 0.059 mas (similar to that obtained by Kutkin et al. 2018), this gives r_{core, 43 GHz} ≃ 17 pc. This would situate the distance from the base of the jet to the 43 GHz core much closer than the r_{core, 15 GHz} ≃ 29 pc obtained by Wang & Jiang (2020) at 15 GHz, consistent with opacity effects. The result is also compatible with the constraint > 12 pc from Agudo et al. (2011).

4.4. Spectral energy distribution

We produced complete SEDs for the two epochs of flaring and quiescent state related to the 2008 outburst, where the MWL coverage was highest: MJD 54761 (2008-10-22), which corresponded to the peak of the flare, and MJD 55098 (2009-09-14). Analogously, we built the SED for two epochs related to the 2015 outburst: MJD 56576 (2013-10-11), which was taken as a quiescent epoch, and MJD 57293 (2015-09-28), as the flaring epoch. The last epoch is the closest one to the peak of the flare with observations in enough bands to perform an accurate modeling. These four epochs were marked with vertical lines in the MWL flux plot (Fig. 1) and their SEDs are represented together in Fig. 14 for comparison.

Fig. 14. Four epochs for which the SEDs were analyzed, represented together for comparison. MJD 54761 and MJD 57293 correspond to flaring epochs of the 2008 and the 2015 outbursts, respectively, while MJD 55098 and MJD 56576 correspond to quiescent epochs. The appearance of an X-ray bump is evident in the 2008 flaring epoch. It is also visible, albeit dimmer, in the 2015 flaring epoch. Both quiescent epochs lack this X-ray feature.

It can be seen that both the 2008 and 2015 flaring epochs, MJD 54761 and MJD 57923, exhibit a softening of the spectrum between the hard UV and soft X-ray ranges (Fig. 14). The feature manifests itself as a increase of the flux from optical to UV wavelengths, with the slope in the SED in the UV becoming positive, and as a increased flux in the X-ray region, with the slope becoming negative. The UV increase is much higher when using the extinction values by Junkkarinen et al. (2004), as seen in Fig. 15, but this probably overestimates the correction in the hard UV. The feature is still present in both epochs when using the values for extinction given by Ackermann (2012), specially for MJD 57923 and considerably dimmer for MJD 54761, and we have used these values to built the final SEDs. The unexplained feature seems to disappear when the source is quiescent, for both the 2008 and the 2015 flares.

Fig. 15. Resulting UV bump applying the extinction correction from Raiteri et al. (2008) and Ackermann (2012). An increase is present in both cases, but it is much dimmer with the values by Ackermann (2012).

The origin of this feature is still under debate, although its presence has been reported before in the literature, also for previous flares of this source. Raiteri et al. (2008) reported the presence of the UV feature in the peak of the 2006-2007 flare, and also in some other earlier epochs where the source was fainter. The change of slope in X-rays was also present in the SED reported in Ackermann (2012) for the MJD 54761-3 epoch, even though the UV increase was not evident in their SED plots. In contrast, our analysis shows that in epoch MJD 54761 that the bump is visible both in the UV and the X-ray. Raiteri et al. (2008) emphasizes that the fact that the bump is visible during flaring states is unusual for quasars. In most cases, similar features are only visible in the faintest states and are attributed to thermal emission from the disk. In contrast, 0235+164 exhibits this feature even during the brightest epochs, thus ruling out such an explanation. The thermal origin of the feature is invalidated further by the high temperatures that would be necessary to result in a bump at these energies and by the fact that the thermal emission from the disk should be approximately stable, while the difference in flux when the feature becomes visible is of more than one order of magnitude. Raiteri et al. (2008) also reported the presence of the feature in the UV for a quiescent epoch related to the 2007 flares, and some intermediate states. This, together with the lower values for the correlation of the X-rays found in our DCF analysis (Sect. 4.2) hints at a different process involved at least partially in the emission at these energies.

In the remainder of this section, we briefly review previous existing models and perform a comparison between them and ours, summarizing the results in Table 5. A schematic representation of the physical setup can be found in Fig. 16.

Fig. 16. Schematic representation of the physical setup. In the proposed scenario, an electron distribution (blob) in the jet, characterized by its geometrical properties, its energy distribution, and its magnetic field (Tables 7 and 8), emits synchrotron radiation. Infrared and optical photons from the disk, the dusty torus and the broad line region reach the blob and are up-scattered to high energies by inverse Compton. The observer, narrowly aligned with the jet, sees the emission boosted by relativistic effects. In the case of bulk-Compton emission, an additional, different distribution would exist, closer to inner region of the blazar (Table 6).

Agudo et al. (2011) postulated that the mechanism of emission was predominantly SSC from the joint analysis of VLBI images, long-term multi-wavelength light curves from mm to γ-ray energies, including polarization, and time delays. They interpreted the outburst as the result of “the propagation of a disturbance, elongated along the line of sight by light-travel time delays” that passes through a standing recollimation shock in the core and propagates down the jet, thereby creating the “superluminal knot”. They also demonstrated the general correlation between the MWL flux at different bands and the appearance of the 43 GHz VLBA superluminal features and obtained the associated time delays. They argued that the variability in γ rays could not be explained within the EC scenario. Instead, they favored a model where the stronger variability in γ-rays is explained by the delayed variability in a multi-zone turbulent cell model (Marscher et al. 2010). This was supported by the general multi-wavelength correlation, as well as the variability of the polarization and the parameters derived from the superluminal components in VLBI images (see Table 5).

Ackermann (2012) produced a model of the SEDs for epochs 54761-3 (2008-10-22 – 24) and 54803-5 (2008-12-03 – 05). The high state epoch 54761 presented a secondary soft X-ray bump, which was modeled as a bulk-Compton feature, although no hint of a bump was present in the hard UV region in the SED. For both of the epochs, ERCIR (Compton emission from infrared radiation from the dusty torus) was the dominating component at higher frequencies. The bulk-Compton feature was not present in the quiescent state. Ackermann argues that EC must dominate SSC for any reasonable covering factor of the broad-line region. The model presents an emission zone located outside the BLR close to the BH (1.7 pc) with a Lorentz factor Γ = 20, opening angle 2.9°, magnetic field B′ = 0.22 G, and viewing angle 2.3°. The electron energy distribution was modeled by a doubly broken power law. The bulk-Compton feature is modeled by a population of cold electrons.

Baring et al. (2017) modeled the same epoch MJD 54761. They do so with a Lorentz factor Γ = 35. These authors modeled the energy distribution of electrons by simulating their acceleration process through diffusive shock acceleration (DSA). The bulk-Compton feature is also not present in the quiescent state. First- and second-order SSC also contributes to the second bump but is dominated at all energies by the external Compton. However, the authors notice that the Lorentz factor required by this source is significantly higher than the usual value (Γ ∼ 10 − 20) for EC-dominated sources. The Swift-XRT excess is modeled as IC of a seed radiation field of T ∼ 1000 K, postulated to be a dusty torus. Dreyer & Böttcher (2021) also presented a modeling of the SED for the same epoch, based on Baring et al. (2017), where the X-ray bump is explained by bulk-Compton emission. The second bump is also explained by external Compton from the dusty torus. If bulk Comptonization is responsible for the X-ray bump, a prediction is made that it should result in partial polarization in the X-ray bands.

In this work, we have modeled the emission of AO 0235+164 using the JetSeT framework⁵ (Tramacere 2020; Tramacere et al. 2011, 2009), using a SSC + EC scenario. The accretion disk spectrum is modeled as a multi-temperature black body as described in Frank et al. (2002), with a luminosity fixed to L_Disk = 5 × 10⁴⁵ erg s⁻¹, an accretion efficiency (η) fixed to the standard value of 0.08, and a BH mass fixed to 5 × 10⁸ M_⊙, with an external radius of the order of a few hundreds of Schwarzschild radii. The BLR is modeled as a thin spherical shell with an internal radius determined by the phenomenological relation provided by Kaspi et al. (2007), $R_{\mathrm{BLR},\mathrm{in}}=3\times {10}^{17}{L}_{\mathrm{Disk},46}^{1/2}$cm. The external radius of the BLR is assumed to be 0.1R_BLR, in, with a coverage factor of τ_BLR = 0.1. The dusty torus (DT) is assumed to be described by spherical uniform radiative field, with a radius $R_{\mathrm{DT}}=2\times {10}^{19}{L}_{\mathrm{Disk},46}^{1/2}$cm (Cleary et al. 2007), and a reprocessing factor τ_DT = 0.1. The emitting region is modeled as a single spherical zone with a radius, R, located at a distance, R_H, from the central black hole. The jet has a conical geometry, with an half opening angle of ϕ ≈ 3 deg, with the emitting region size determined by R = R_Htanϕ. The emitting region moves along the jet axis with a bulk Lorentz factor, Γ, oriented at a viewing angle, θ, and a consequent beaming factor of $\delta =1/(\mathrm{\Gamma }\sqrt{1-{\beta }_{\mathrm{\Gamma }}cos(\theta ))}$. For the relativistic emitting electron distribution (EEE), we tested a broken power law (BKN) distribution:

$$\begin{array}{c}\hfill n(\gamma )=N\{\begin{array}{cc}{\gamma }^{-p}\hfill & {\gamma }_{min}\le \gamma \le {\gamma }_b,\hfill \\ {\gamma }^{-p_1}{\gamma }_b^{p-p_1}\hfill & {\gamma }_b<\gamma <{\gamma }_{max},\hfill \end{array}\end{array}$$(12)

with an index of p and p₁ below and above the break energy, γ_b, respectively, and a power law distribution with a cut-off (PLC) distribution of

$$\begin{array}{c}\hfill n(\gamma )=N{\gamma }^{-p}exp\frac{\gamma }{{\gamma }_{\mathrm{cutoff}}},{\gamma }_{min}\le \gamma \le {\gamma }_{max}.\end{array}$$(13)

The initial values of L_Disk and T_Disk were determined by JetSeT during the pre-fit stage and L_Disk is frozen to the value of L_Disk = 5 × 10⁴⁵ erg s⁻¹. The model minimization was performed using the JetSeT ModelMinimizer module plugged into the iminuit Python interface (Dembinski et al. 2020). The errors were estimated from the matrix of second derivatives, using the HESSE method. We fit the data above 30 GHz, excluding data below the synchrotron self-absorption frequency. To avoid the small errors in the UV-to-radio frequencies biasing the fit toward the lower frequencies, we added a 20% systematic error to data below 10¹⁶ Hz. We find that the PLC model provides a slightly better fit to the data, thus, in the following, we present only the results for this model. All the states presented in this analysis can be modeled by a single-zone EC-dominated (see Figs. A.4, A.6, and A.8 as well as Table 7) or an SSC-dominated scenario (see Figs. A.5 and A.7 and Table 8), with the SSC-dominated scenario resulting in systematically lower values of B, needed to accommodate for the proper U_e/U_B ratio able to match the peak flux and frequency of the IC emission. On the contrary, for the flaring state on MJD 54761, the presence of a strong and soft bump in the X-ray makes both the SSC and EC unable to model the data. As suggested by Celotti et al. (2007), Ackermann (2012), this spectral feature can be explained by the Comptonization of the external radiative fields by a population of cold electrons. We have introduced such bulk-Compton (BC) component, modeled as a spherical region with a radius, R_BC, moving with corresponding bulk factor of Γ = 10, at a distance, r, from the BH and with a total number of particles, N_BC.

We notice that for a purely cold population, namely, for electrons with γ_min = γ_max = γ = 1, the resulting shape of the BC radiation was always too steep to reproduce the observed data (see e.g., Celotti et al. 2007); on the contrary, we found that a reasonable fit to the data was provided by increasing the fit range of γ_max to 5, and setting r = 1.5 × 10¹⁶ cm. With this model configuration, the fit converged with a resulting value of γ_max ≈ 4 and a resulting total number of cold electrons of N_BC ≈ 1.8  ×  10⁵⁴ (see Fig. A.2 and left column in Table 6). These values are compatible with those reported in Ackermann (2012; N_BC = 2.4  ×  10⁵⁴ and r = 5 × 10¹⁵ cm), anyhow we stress that in Ackermann (2012) the BC spectral shapes is assumed to be a PL, whilst in the present analysis, it is obtained by the actual Comptonization of the cold electrons. We also applied the BC model to an SSC-dominated scenario (see Fig. A.3 and Tables 6 and 8), we notice that even though the overall agreement of the model with the data is still reasonable, the model shows some tension in the optical-IR and X-ray data. At the high-energy branch of the X-ray data, the excess of flux in the model is due to the broader spectrum of IC emission compared to the EC case, originating from the broader spectrum of the synchrotron seed photons compared to the narrower seed photon spectrum of the external fields.

Another possible option that would be able to produce a PL shape for the BC could be obtained by assuming a purely cold electron population with a truncated conical geometry. The higher energy part of the BC would be produced by the low-number electrons closer to BH and the higher energy would be produced by the larger number of electrons in the upper part of the truncated cone. To mimic such a geometry we implemented a BC model with two spherical regions. The radius of the two regions is obtained so that the two spheres match the volume of the upper and lower parts of the truncated cone. We find that a reasonable modeling of the BC emission is obtained assuming a truncated cone, with an opening angle of 45° and an height of ≈9 × 10¹⁵ cm. The smaller spherical region corresponds to the segment of the truncated cone with an height of ≈5 × 10¹⁵ cm and the larger spherical region corresponds to the segment with an height of ≈8 × 10¹⁵ cm The total number of cold electrons is of N_BC ≈ 1.4 × 10⁵⁴ (see Fig. A.2 and left column in Table 6). Since the introduction of this extra component introduces new parameters, first, we used the ModelMinimizer to fit the model to the data without the BC component and excluding the X-ray data (statistics are reported in Tables 7 and 8) and, in a second step, we added the X-ray data and we proceeded to a qualitative fitting of the BC conical component (the values of the BC component are reported in Table 6).

The flaring epoch MJD 57293 could also be modeled using a single-zone model, although the observed softening of the X-ray spectrum could be explained by bulk-Compton emission in two-zone model. This could be done in a similar manner to the case of MJD 54761, as the DCF analysis for 2014-2017 indicates.

5. Discussion and conclusions

We have presented new and updated multi-wavelength photometric and polarimetric data of AO 0235+164, all across the spectrum from radio cm and mm wavelengths up to γ-ray energies. The analysis of the correlations have shown that the emission at different wavelengths is statistically correlated, linking their emission mechanisms, with the notable exception of the X-ray band.

We have analyzed and shown the compatibility between the positions of the peaks of the different correlations, strengthening their interpretation as the delay between emissions. In this context, we have also shown that the obtained delays are compatible with the proposed emission mechanisms: from mm to optical wavelengths, the delays agree with what it is to be expected for synchrotron emission.

In addition, we have also seen that the γ-ray light curve is indeed correlated with the mm and R-band emissions, which is to be expected if the dominating emission mechanism is SSC or EC. Furthermore, the γ-ray subflares seem to be related to the appearance of identifiable VLBI components.

On the other hand, we have not found a significant correlation between the X-ray light curve and the rest of the bands. This is explained by the presence of the X-ray bump in the SED. This bump can not be accounted for by a closely correlated emission (SSC or EC) with the rest of the bands. Instead, it is proposed that it corresponds to bulk-Compton emission from a different population of particles. The large obtained delays imply that this emitting zone is separated by a large distance from the main emission component and this is further confirmed by the results from the SED modeling.

Understanding how our observational data and results fit in the current landscape of existing blazar models is a difficult task. There is the rebrightening of knot features, which could be explained by successive recollimation shocks with the jet, and the difference in Doppler factor and speed between different components, which could be explained by different energies of a shock wave, points toward a shock-in-jet model. The observed post-maximum subflares in 3 mm and γ-ray can be explained by less energetic recollimation of the same dulled shockwave, analogously to the rebrightning of knot features farther from the jet as seen in the VLBA images, they even appear to be more or less simultaneous. The observed longer duration of the flare in mm wavelengths is explained in this model by the longer cooling of synchrotron electrons. This smears out the peak in the correlation and shifts the correlation shape to show a delay of mm emission.

The question about whether SSC or EC dominates the high energy bump does not have a clear, definite answer. EC-dominated SED models seem to be favored by literature (Ackermann 2012, Dreyer & Böttcher 2021). However, as we present in this paper, SSC-dominated models are also possible (as shown in Sect. 4.4). It is generally easier and more common to produce a fit with dominant EC, however, the model is harder to explain physically and the obtained delays in the correlation analysis and the results from VLBI observations favor SSC-dominated models.

The delays between signals are not directly interpretable as the relative time at which emissions at different wavelengths start; this interpretation would be valid only if the signals had the same shape but were shifted with respect to each other, which is not the case here. However, the correlation between R and γ shows a clear peak whose position is ${{\tau }_{\mathrm{p}}}^{R,\gamma }$ of 2 days, which corresponds to a distance of less than 1 pc after accounting for relativistic effects. Meanwhile, the large delay obtained between R and X-ray place the emission regions at tens of parsecs away, which aptly fits the obtained distances in the SSC scenario where the X-ray is produced by bulk-Compton emission.

The results from the kinematic analysis of VLBI components show that the 43 GHz core is located at distances from 12 pc to 17 pc downstream from the central BH, assuming a conical jet geometry. The best-fit distances obtained in SSC-models (Table 8) are in better agreement with the ones obtained from the VLBI kinematic analysis and, in any case, since the SSC emission is less dependent on the distance to the BH, other distances are easier to accommodate. This is not the case in the EC-scenario.

Scenarios where the γ-emitting zone is close to the central BH are ruled-out by the long-term and highly significant correlation (Fig. 10) between γ, R, and mm light curves, since the emissions must be close enough and from analysis of VLBI images we know this is more than ten parsecs away from the central engine. SED models also help us discard these scenarios.

The presence of IC flares after the synchrotron flares has already ceased, as in the case in some instances between the 2008 and 2015 flares, is also an indicator of SSC (Sokolov et al. 2004). They can be explained by the time-delays and crossing times, specially for small viewing angles such as AO 0235+164. However, this is not valid in a EC scenario. Also, the observed stronger variability in γ-rays with respect to low energies is harder to explain in the EC scenario, where there is no reasonable source of increased variability.

A good test to determine whether the emission is SSC or EC might be the polarization of the γ-rays. Generally, EC is not expected to have significant polarization, while SSC is expected to have a polarization degree about half of the corresponding synchrotron emission. While the X-ray polarization is already being measured by some instruments (IXPE), γ-ray polarization is still not possible, although recent technological developments have opened the possibility up in the next decade.

Acknowledgments

The IAA-CSIC team acknowledges financial support from the Spanish “Ministerio de Ciencia e Innovación” (MCIN/AEI/10.13039/501100011033) through the Center of Excellence Severo Ochoa award for the Instituto de Astrofísica de Andalucía-CSIC (CEX2021-001131-S), and through grants PID2019-107847RB-C44 and PID2022-139117NB-C44. 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. J.Y.K. was supported for this research by the National Research Foundation of Korea (NRF) grant funded by the Korean government (Ministry of Science and ICT; grant no. 2022R1C1C1005255). IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). The VLBA is an instrument of the National Radio Astronomy Observatory, USA. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This study was based (in part) on observations conducted using the 1.8 m Perkins Telescope Observatory (PTO) in Arizona (USA), which is owned and operated by Boston University. The BU group was supported in part by US National Science Foundation grant AST-2108622, and NASA Fermi GI grants 80NSSC20K1567 and 80NSSC22K1571.

References

Appendix A: Spectral Energy Distribution models

This is a collection of the spectral energy distribution (SED) models presented in Sect. 4.4.

Fig. A.1. SED model for epoch MJD 54761. The model includes the usual synchrotron plus SSC components, but the high-energy bump is dominated by EC emission from a dusty torus (Table 7). The X-ray bump is modeled by bulk-Compton emission from the disk by a secondary particle distribution much closer to the central engine (Table 6), consistent with the much lower correlation and higher delays in the DCFs between X-ray and the other bands (Fig. 12).

Fig. A.2. SED model for epoch MJD 54761. The EC-dominated model includes a bulk-Compton component with a conical shape.

Fig. A.3. SED model for epoch MJD 54761 in the SSC-dominated scenario. The X-ray bump is modeled as bulk-Compton emission from the disk in a similar way to the EC-dominated model (Table 6).

Fig. A.4. SED model for epoch MJD 55098 in the EC-dominated scenario (Table 7).

Fig. A.5. SED model for epoch MJD 55098 in the SSC-dominated scenario (Table 8).

Fig. A.6. SED model for epoch MJD 56576 in the EC-dominated scenario (Table 7).

Fig. A.7. SED model for epoch MJD 56576 in the SSC-dominated scenario (Table 8).

Fig. A.8. SED model for epoch MJD 57293 in the EC-dominated scenario (Table 8). Unlike the models for the older flaring epoch 54761 (Figs. A.1, A.2), this model does not include a bulk-Compton component and can be explained with only the usual SSC+EC components. The source exhibits an important softening in the X-ray spectrum that could also be explained by bulk-Compton emission.

Appendix B: Swift observations

The Neil Gehrels Swift observatory satellite (Gehrels et al. 2004) carried out 195 observations of AO 0235+164 between 2005 June 28 (MJD 53549) and 2016 February 11 (MJD 57429). The observations were performed with all three on-board instruments: the X-ray Telescope (XRT; Burrows et al. 2005, 0.2–10.0 keV), Ultraviolet/Optical Telescope (UVOT; Roming et al. 2005, 170–600 nm), and the Burst Alert Telescope (BAT; Barthelmy et al. 2005, 15–150 keV).

All XRT observations were performed in photon counting mode (for a description of XRT read-out modes, see Hill et al. 2004). The XRT spectra were generated with the Swift-XRT data product generator tool at the UK Swift Science Data Centre⁶ (for details, see Evans et al. 2009). Spectra having count rates higher than 0.5 counts s⁻¹ may be affected by pile-up. To correct for this effect, the central region of the image was excluded and the source image has been extracted with an annular extraction region with an inner radius that depends on the level of pile-up (see e.g., Moretti et al. 2005). We used the spectral redistribution matrices in the Calibration database maintained by HEASARC. The X-ray spectral analysis was performed using the XSPEC 12.13.0c software package (Arnaud 1996). Data were grouped for having at least 20 counts per bins with grppha and the chi square statistics is used. All XRT spectra are fitted with an absorbed log-parabola model, except for cases with low number of counts, and a HI column density fixed to 2.8×10²¹ cm⁻² for taking into account the absorption effects of both our own Galaxy and an intervening z = 0.524 system (see e.g. Madejski et al. 1996).

The hard X-ray flux of this source is usually below the sensitivity of the BAT instrument for daily short exposures. Moreover, the source is not included in the Swift-BAT 157-month catalogue⁷. During the Swift pointings, the UVOT instrument observed the sources in its optical (v, b, and u) and UV (w1, m2, and w2) photometric bands (Poole et al. 2008; Breeveld et al. 2010). The UVOT data in all filters were analysed with the uvotimsum and uvotmaghist tasks and the 20201215 CALDB-UVOTA release. Source counts were extracted from a circular region of 5 arcsec radius centered on the source, while background counts were derived from a circular region with a 20 arcsec radius in a nearby source-free region.

All Tables

All Figures

	Fig. 1. Light curves of AO 0235+164 at different wavelengths. The vertical lines mark the epochs whose SEDs were analyzed in Sect. 4.4. The top panel shows both the X-ray and γ-ray fluxes in different axes (left for X-rays and right for γ-rays, respectively).
In the text

	Fig. 2. Polarization degree evolution of AO 0235+164 at different wavelengths. Together with experimental data, a bayesian block representation is shown superimposed for R (black line) at a 99.9% confidence and for 1 mm (red) and 3 mm (blue) at a 90% confidence level. The vertical lines are the same as in Fig. 1 and correspond to the epochs whose SED was analyzed in Sect. 4.4.
In the text

	Fig. 3. Polarization angle evolution of AO 0235+164 at different wavelengths. The vertical lines are the same as in Fig. 1 and correspond to the epochs whose SEDs were analyzed in Sect. 4.4. All points in the first box correspond to R band, the colors denote the clusters that were shifted by n × 180°, as noted in Sect. 3.2.
In the text

Fig. 4. Selected epochs illustrating the evolution of each identified knot, showing total (contours) and polarized (color scale) intensity. The beam size is indicated as a green ellipse in the first row. Horizontal black lines indicate the position of the core A0, black line segments within the image indicate the direction of polarization (EVPA). The red line in each row is the linear fit to the knot position. It is present for every knot except B7, whose flux was too low to accurately perform a fit. For each row, the spacing between plots is proportional to time, and the total time span is different and indicated in brackets.

In the text

Fig. 5. Epoch 102 (2016-09-05) VLBI 7 mm image, showing total flux intensity (contours) and polarized flux intensity (color scale). Black line segments overlaid in the image represent the EVPA. The green ellipse in the lower left corner represents the beam size. The image showcases component B5 close to the peak of the 2015 flaring episode, demonstrating how the polarization angle is aligned with the direction of propagation.

In the text

	Fig. 6. Epoch 110 (2017-04-16) VLBI 7 mm image. The image shows the trailing component B6 moving in the same direction of B5, maintaining the alignment of the polarization angle (black line segments) with the momentary direction of the jet (northeast).
In the text

Fig. 7. Fits to the profiles of the 2008 (top) and the 2005 (middle) γ-ray flares by exponential functions as described by Eq. (1). The selected number of terms Nexp were chosen so that the reduced χ2-statistic was minimized (bottom). The best-fit values are given in Tables 1 and 2. An alternative fit is given for the 2008 flare with a similar χ2, that accounts for the double peak at the beginning with a single exponential term. The vertical lines mark the epochs whose SED was analyzed in Sect. 4.4, as in Fig. 1.

In the text

	Fig. 8. Correlations between fluxes across all wavelengths. Horizontal lines show significance levels for 1σ, 2σ, and 3σ, computed using N = 2000 synthetic light curves, as described in Sect. 4.2. The DCFs here are computed using the whole period of available data, from 2007 to 2020.
In the text

Fig. 9. Observed distance from core (top) and flux density (bottom) for every one of the identified component as a function of time, together with a linear weighted fit to knot distance and logarithmic fit to flux. The fit to the flux is done taking into account only the points after the peak of emission of each component. This was done for all components except B7, which due to its low flux, did not have enough points after the peak with low enough uncertainty to perform a fit.

In the text

	Fig. 10. Correlations between fluxes and polarization degree across all wavelengths. Horizontal lines show significance levels for 1σ, 2σ, and 3σ, computed using N = 2000 synthetic light curves, as described in Sect. 4.2. The DCFs here are computed using the whole period of available data, from 2007 to 2020.
In the text

	Fig. 11. Correlations between fluxes across all wavelengths. Horizontal lines show significance levels for 1σ, 2σ, and 3σ, computed using N = 2000 synthetic light curves, as described in Sect. 4.2. The DCFs here are computed using only the flaring episode from 2014 to 2017.
In the text

Fig. 12. Correlations compatibility chart for τp. Each row corresponds to a band i. The delays or peak positions between the row band i and any other band j, ${\tau }_{\mathrm{p}}^{i,j}$, are plotted along the x-axis, shifted by the delay between the band i and the reference band, ${\tau }_{\mathrm{p}}^{i,R}$. The peaks selected are the one with highest values of the DCF, marked with red dots in Fig. 11. The DCFs here use the flaring episode from 2014 to 2017.

In the text

	Fig. 13. Fit to the function y = ax−1/2 + b, which is the expected form if the time delays are due to different synchrotron cooling times. Fit results are a = 1.341 ± 0.009 × 107, b = −6 ± 1 × 10−1, p(χ2) = 0.909, and r2 = 0.9998. Parameter b is close to zero and accounts for an arbitrary reference delay, which is R in our case.
In the text

Fig. 14. Four epochs for which the SEDs were analyzed, represented together for comparison. MJD 54761 and MJD 57293 correspond to flaring epochs of the 2008 and the 2015 outbursts, respectively, while MJD 55098 and MJD 56576 correspond to quiescent epochs. The appearance of an X-ray bump is evident in the 2008 flaring epoch. It is also visible, albeit dimmer, in the 2015 flaring epoch. Both quiescent epochs lack this X-ray feature.

In the text

	Fig. 15. Resulting UV bump applying the extinction correction from Raiteri et al. (2008) and Ackermann (2012). An increase is present in both cases, but it is much dimmer with the values by Ackermann (2012).
In the text

Fig. 16. Schematic representation of the physical setup. In the proposed scenario, an electron distribution (blob) in the jet, characterized by its geometrical properties, its energy distribution, and its magnetic field (Tables 7 and 8), emits synchrotron radiation. Infrared and optical photons from the disk, the dusty torus and the broad line region reach the blob and are up-scattered to high energies by inverse Compton. The observer, narrowly aligned with the jet, sees the emission boosted by relativistic effects. In the case of bulk-Compton emission, an additional, different distribution would exist, closer to inner region of the blazar (Table 6).

In the text

Fig. A.1. SED model for epoch MJD 54761. The model includes the usual synchrotron plus SSC components, but the high-energy bump is dominated by EC emission from a dusty torus (Table 7). The X-ray bump is modeled by bulk-Compton emission from the disk by a secondary particle distribution much closer to the central engine (Table 6), consistent with the much lower correlation and higher delays in the DCFs between X-ray and the other bands (Fig. 12).

In the text

	Fig. A.2. SED model for epoch MJD 54761. The EC-dominated model includes a bulk-Compton component with a conical shape.
In the text

	Fig. A.3. SED model for epoch MJD 54761 in the SSC-dominated scenario. The X-ray bump is modeled as bulk-Compton emission from the disk in a similar way to the EC-dominated model (Table 6).
In the text

	Fig. A.4. SED model for epoch MJD 55098 in the EC-dominated scenario (Table 7).
In the text

	Fig. A.5. SED model for epoch MJD 55098 in the SSC-dominated scenario (Table 8).
In the text

	Fig. A.6. SED model for epoch MJD 56576 in the EC-dominated scenario (Table 7).
In the text

	Fig. A.7. SED model for epoch MJD 56576 in the SSC-dominated scenario (Table 8).
In the text

	Fig. A.8. SED model for epoch MJD 57293 in the EC-dominated scenario (Table 8). Unlike the models for the older flaring epoch 54761 (Figs. A.1, A.2), this model does not include a bulk-Compton component and can be explained with only the usual SSC+EC components. The source exhibits an important softening in the X-ray spectrum that could also be explained by bulk-Compton emission.
In the text