The Repeating Flaring Activity of Blazar AO 0235+164

Context. Blazar AO 0235+164, located at redshift z = 0.94, has undergone several sharp multi-spectral-range flaring episodes during the last decades. In particular, the episodes peaking in 2008 and 2015, that received extensive multi-wavelength coverage, exhibited interesting behavior. Aims. We study the actual origin of these two observed flares by constraining the properties of the observed photo-polarimetric variability, those of the broad-band spectral energy-distribution and the observed time-evolution behavior of the source as seen by ultra-high resolution total-flux and polarimetric Very-long-baseline interferometry (VLBI) imaging. Methods. The analysis of VLBI images allows us to constrain kinematic and geometrical parameters of the 7 mm jet. We use the Discrete Correlation Function to compute the statistical correlation and the delays between emission at different spectral ranges. Multi-epoch modeling of the spectral energy distributions allows us to propose specific models of emission; in particular for the unusual spectral features observed in this source in the X-ray region of the spectrum during strong multi spectral-range flares. Results. We find that these X-ray spectral features can be explained by an emission component originating in a separate particle distribution than the one responsible for the two standard blazar bumps. This is in agreement with the results of our correlation analysis that do not find a strong correlation between the X-rays and the remaining spectral ranges. We find that both external Compton dominated and synchrotron self-Compton dominated models can explain the observed spectral energy distributions. However, synchrotron self-Compton models are strongly favored by the delays and geometrical parameters inferred from the observations.


Introduction
Blazars are among the most energetic objects in the universe.They are widely accepted to consist of a super massive black hole, referred to as the central engine, surrounded by an accretion disk and usually a dusty torus, and 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, still under research, to speeds close to the speed of light.This results in highly energetic emission of radiation across the entire electromagnetic sion from the interaction of the -charged-relativistic particles of the jet with the magnetic fields in the medium is accepted to account for the first bump.Several scenarios exist to explain the second bump.In the leptonic scenario, the second bump is explained by inverse Compton effect of relativistic electrons interacting with ambient photons, and distinction is made whether these photons originate from the synchrotron emission inside the jet, in which case the mechanism is labeled as 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 about the relevance of other different mechanisms, such as the 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, and the origin and location of photon fields and particles involved is not sufficiently well established.
The study of the variability of blazars across the spectrum, combined with the analysis of sequences of ultra-high-resolution 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 to discard or support different models.
Regarding the location of the γ-ray emission, two main possibilities have been under discussion, differing in the distance to the central black hole (BH).The first one is the so called "closezone" scenario, very close to the BH (0.1 − 1pc), that 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 shown interesting flaring behavior with the most recent flares occurring in 2008 and 2015 that have been studied with multi-wavelength (MWL) and Very Large Baseline Interferometry (VLBI).The source typically appears extremely compact at ultra-high resolution 7mm 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) which can together 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 in this paper to 2020, where 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, which have produced light curves showing flares in previous years, e.g.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 ∼ 6 years in the behavior of the source, which can serve as a guidance while developing models of the source, even if data is not conclusive enough to settle this hypothesis and 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 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 et al. (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 at γ-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 errors and were located several parsecs from the central engine, and 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 0 = 67.66km/Mpc as given by Planck Collaboration et al. 2020.

Observations
We have obtained and compiled time dependent data in most available ranges of the electromagnetic spectrum, including polarimetry whenever it was possible, and VLBI polarimetric images with submilliarcsecond resolution.
Our observations include 7 mm VLBA images from the blazar monitoring program at Boston University; reduced both for total flux and polarization using AIPS (see Weaver et al. 2022 for details about data reduction and calibration).Single dish data at 1mm and 3mm were obtained from the POLAMI (Polarimetric Monitoring of AGN at Millimeter Wavelengths)1 program at the IRAM 30m Telescope (Agudo et al. 2017a, Thum et al. 2017, Agudo et al. 2017b); optical (R-band) data from Calar Alto (2.2 m Telescope) under the MAPCAT program, the Yale University SMARTS blazar program, Maria Mitchell, Abastumani and Campo Imperatore observatories, Steward Observatory (2.3 and 1.54 m Telescopes), the Perkins Telescope Observatory (1.8 m Telescope), the Crimea Observatory AZT-8 (0.7 m Telescope) and St. Petersburg State University .
Ultraviolet measurements were obtained by the Swift-UVOT instrument.The dataset 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; from where light curves and spectral indices were derived using a broken-power law model and the appropriate corrections for extinction.More details about the data reduction procedure from Swift is provided in Appendix A. Gamma-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 3mm.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 ultraviolet (UV) bands were corrected following the prescription by Raiteri et al. (2005) and the updated values by Ackermann et al. (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 ultraviolet bump appears in the final spectra for some epochs; as shown in Raiteri et al. (2005) although in disagreement with the SEDs presented in Ackermann et al. (2012).It must be noted that applying different correction factors available (NED2 , 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 et al. (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 21 cm −2 (Madejski et al. 1996, Ackermann et al. 2012), which accounts both for galactic extinction and the z = 0.524 absorber.This value agrees with the value obtained by letting N H vary as a free parameter.
Figure 1 show that flaring episodes happen almost simultaneously across all the electromagnetic spectrum.Variability is much more pronounced at HE, milder at optical and UV wavelengths, and softer in the millimeter and radio bands.

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.30 % from 2014 to 2017).The 3mm 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 can be 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 3mm only varies a full cycle or a partial rotation.There is also an apparent delay of approximately a hundred days between 3mm and R-band as can be seen in Fig. 3.
The direction of the EVPA (Electric Vector Position Angle) of the VLBA components (indicated in Figs. 4, 5 and 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 off the direction of the observer, makes the electric field to align with the jet direction.This supports the association of the superluminal components ejected during flares with plane-perpendicular moving shock-waves.

VLBI imaging
Our study includes all available 7mm (43GHz) 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 programs3 ).After reducing the data with AIPS (see Weaver et al. 2022), most prominent jet features were fitted to Gaussian components with Diffmap and then cross-identified along 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 here as the core.Other features can be tracked at different epochs, and their evolution traced in 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 polarization degree for the total emission and single components were later produced from the images (Figs. 1, 2 and 3) using the aforementioned identification.
The flux evolution shown in Fig. 1 at all wavelengths, also containing the light curves from the integrated VLBA 7mm 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, as in each flaring episode the direction of propagation of the associated superluminal components is radically different for every on of these episodes.

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, ergs/cm 2 /s) 2 0 0 8 -1 0 -2 2 2 0 0 9 -0 9 -2 4 2 0 1 3 -1 0 -1 1 2 0 1 5 -0 9 -2 8     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 99.9 % confidence and for 1mm (red) and 3mm (blue) at 90 % confidence level.The vertical lines are the same as in Fig. 1 and correspond to the epochs whose SED was analyzed in sec.4.4.
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 3mm 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 around Fall 2014, with a brightening of the core visible at 7mm.The flux densities at all wavelength peak a year later, around Fall 2015, at all wavelengths.Comparatively, this flare is dimmer than the previous one of 2008, reaching 4.2 Jy in 3mm and 1.3 magnitudes less in 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.
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 of viewing angle.Some months later another, a weaker component, named B6, originates 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) which allowed us to derive the rising and decaying times of each subflare, t ri and t di .Then, their asymmetry factor, defined as could be computed.An asymmetry factor close to zero corresponds to the case of a perfectly symmetric flare.There exists some uncertainty in the number of exponential terms to use, since the source shows strong variability in timescales shorter than our binning allows to track.The large binning used, of 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 χ 2 -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 sec.4.2 (τ R,γ = 2 days, τ R,7mm = 64 ± 4 days).
During 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, to the uncertainty in the γ-ray lightcurve, 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 Table 1.Parameters for the fit to the γ-ray lightcurve of the 2008 flaring episode to functions of shape given by eq. ( 1).The resulting reduced χ 2 -statistic for the fit is shown, and also the computed symmetry factor ξ i for each subflare.Some values could not be computed.Upper limits are indicated with '<'.The result can be seen in Fig. 7.  7 (0.4)  1.5 (0.4)  0.5 (0.1)  0.2 (4 × 10 3 )  0.09 (4 × 10 1 )   25.7 (2.1)  5.3 (5.0)  4.6 (2.9)  < 17.2  < 290   < 4.4  27.9 (4.7)  37.7 (21 −0.7  −0.8  −1.0  −1.0   Table 2. Parameters for the fit to the γ-ray lightcurve of the 2015 flaring episode.The resulting reduced χ 2 -statistic for the fit is shown, and also the computed symmetry factor ξ i for each subflare.The result can be seen in Fig. 7.
2) 0.7 (0.1) 0.15 (0.07) 58 ( 12) 49 ( 35 to the first and second maximums at γ-rays, taking into account the aforementioned delay. The γ-ray subflare that can be seen in Fig. 10 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 is 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 7mm VLBI components and the successive γ-ray subflares.For the 2015 flaring episode, the rising and decaying times shown in Table 2, as low as ∼ 18 days (taking into account only the two strongest subflares) are compatible with the sizes found in the SED modeling of sec.4.4 (Table 7), which limits to 15 days the shortest timescale where significant variations of flux can occur (eq.8).These other values were obtained only from the modeling of the SEDs and the two analysis are completely independent.In the case of the 2008 flaring episode, however, the shorter times (∼ 4 days) are in 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 χ 2 -square statistic, but rising and decaying times of 18 and 36 days (dashed line in Fig.The selected number of terms N exp were the ones such 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 sec.4.4, as in Fig. 1.
7).In any of the ways, the observed delays between γ and mm might be explained by a combination of adiabatic expansion and cooling time (Tramacere et al. 2022).

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 0 , 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 speed were obtained, and also their flux was fitted to a decaying exponential function F = F 0 exp (−t/t var ), obtaining their timescale of variability (Fig. 8).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 for 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) where a S max 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, and the viewing angle, could 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; but it's apparent speed is much less than any of the others.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 to have statistics in any bin but short enough that the the correlations were not smoothed too much and peak positions could still be determined.To confirm the robustness of our choice, we have also reproduced our analysis with bin-sizes from 10 to 30 days, obtaining similar results (except for some bins disappearing due to not having 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.Randomization of the Fourier transform was used for mm-wavelengths, while for 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.

Correlation of the whole period (2007 to 2020)
The results from DCFs of the whole available period of data (2007 to 2020) show a clear correlation between flux at almost all wavelengths (> 3σ), with most peak positions close to zero (Fig. 10).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 is not certain in the other bands, possibly explained by the sparser sampling and larger errors (Fig. 12).

Correlations of flaring episode (2014 and 2017)
The results from DCFs of the flaring episode (2014 to 2020) show again 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 at a different emission zone.This suggests that a different, separated processes could be responsible at least partially for the emission in X-rays, hypothesis also favored by the analysis of the spectral energy distribution of the source (section 4.4).
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 on-going.Specially as longer periods are taken into account, since more and 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, and use simply the absolute maximum of the DCF, justifying the decision by 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, e.g. the delay between A and B plus the delay 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 was done in Fig. 13 (Table 4) using 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, τ i, j p , are plotted along the x-axis, shifted by the delay between the band i and the reference band, τ i,R p , so that they fall aligned on the same positions.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 1mm 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 τ R,γ p ∼ +2 days as the correct delay τ R,γ .Discarding the 1mm and Gamma rows, we estimate that the mean delays with respect to R for 7mm, 3mm, and X-ray, τ i,R = τ p i,R are 64 ± 4 days, 42 ± 6 days and 73 ± 4 days respectively.The delays for 7mm, 3mm, and R are consistent with expectations if the mechanism of emission is synchrotron cool-ing, which should result in delays of the form τ s ∝ ν −1/2 s , where ν s is the synchrotron frequency, as can be seen in Fig. 9.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 already suggested by the lower level of correlation found, and as analysis of the spectral energy distributions in section 4.4 reveal.
Table 4.Estimated delays (in days) obtained from peaks of the DCFs (Fig. 11) and represented in the compatibility chart (Fig. 13), with their average and dispersion.The delay between R and γ can be extracted directly from the DCF in Fig. 11

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 formulas 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 it needs not be the case as seen in the previous section if they arise from synchrotron cooling).
The sizes of the emitting region can be constrained with For the moving component B2 corresponding to the 2008 outburst, using the timescale of variability in Table 3 and the δ D ≃ 67, one obtains sizes of around 2 pc, consistent with the angular measure of VLBI images.Using the gamma-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 et al. 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) , which, through the relations between this angle and the true speed and Doppler factor and relation ( 8) limits to a minimum of 1 pc the sizes of the mmemitting 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, 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 Θ 0 , 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. For a core size at 43GHz 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 15GHz, consistent with opacity effects.The result is also compatible with the constraint > 12 pc from Agudo et al. 2011.

Spectral energy distribution
We have 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 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.
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 et al. (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.
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 et al. (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, hence ruling out such an explanation.The thermal origin of the feature is further discarded 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 (section 4.2), hints at a different process involved at least partially in the emission at these energies.
In the remaining of this section we will briefly review previous existing models and perform a comparison between them and ours, the result of which is summarized in Table 5.A schematic representation of the physical setup can be found in Fig. 16.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 "a consequence 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 to create 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 multiwavelength correlation, the variability of the polarization and the parameters derived from the superluminal components in VLBI images (see Table 5).Ackermann et al. (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 Xray 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 et al. 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) models the same epoch MJD 54761.They do so with a Lorentz factor Γ = 35.They model 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 (Γ ∼ 10 − 20) for ECdominated 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 framework6 framework (Tramacere 2020;Tramacere et al. 2011Tramacere et al. , 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 45 , erg/s, an accretion efficiency (η) fixed to the standard value of 0.08, and a BH mass fixed to 5×10 8 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 Disk,46 cm.The external radius of the BLR is assumed to be 0.1R BLR,in , with a coverage factor τ BLR = 0.1.The dusty torus (DT) is assumed to be described by spherical uniform radiative field, with a radius Disk,46 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 ϕ ≈ 3 deg, with the emitting region size determined by R = R H tan ϕ.The emitting region moves along the jet axis with a bulk Lorentz factor Γ, oriented at a viewing angle θ, and a consequent beaming factor δ = 1/(Γ 1 − β Γ cos(θ)).For relativistic emitting electron distribution (EEE) we tested a broken power law (BKN) distribution with an index of p and p 1 below and above the break energy γ b , respectively, and a powerlaw distribution with a cut-off (PLC) Disk temperature Notes. (a) The components considered are Synchrotron (Synch), Synchrotron Self-Compton (SSC), External Compton (EC) from the Dusty Torus (DT), and Bulk Compton (BC). (b) The values for Agudo et al. (2011) do not come from a SED model, but from the DCF analysis and kinematic parameters from VLBI images, assuming a SSC scenario.The value for Γ is not cited in the paper, it is the one obtained by Weaver et al. 2022 for the same component in VLBI. (c) Parameters for the blazar zone (their model includes a second population of relatvistic cold electrons to account for the secondary soft X-ray bump, whose parameters are indicated between parenthesis). (d) The secondary soft x-ray bump is modeled by bulk Comptonization of a background seed field from a dusty torus); the H.E. bump by External Compton of the electron population. (e) The energy distribution is simulated from Diffusive Shock Acceleration (DSA), resulting parameters are not explicitly indicated. distribution The initial values of L Disk and T Disk , are determined by JetSeT during the pre-fit stage, and L Disk is frozen to the value of L Disk = 5 × 10 45 erg s −1 .The model minimization is performed using the JetSeT ModelMinimizer module plugged to iminuit python interface (Dembinski & et al. 2020).The errors are estimated from the matrix of second derivatives, using the HESSE method.We fit data above 30 GHz, excluding data below the synchrotron self-absorption frequency.To avoid that the small errors in the UV-to-radio frequencies biasing the fit toward the lower frequencies, we add a 20% systematic error to data below 10 16 Hz.We find that the PLC model provides a slightly better fit to the data, hence 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 22, 24 and Table 7) or SSC-dominated scenario (see Figures: 21,23 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 et al. 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 Γ = 10, at a distance of r from the BH, and with a total number of particle N BC .
We noticed that for a purely cold population, i.e. for electron with γ min = γ max = γ = 1, the resulting shape of the BC radiation was always to 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 16 cm.With this model configuration, the fit converged with a resulting value of γ max ≈ 4 and a resulting total number of cold electrons N BC ≈ 1.8 × 10 54 (see Figure 18 and left column in Table 6).These values are compatible with those reported in Ackermann et al. 2012 (N BC = 2.4 × 10 54 and r = 5 × 10 15 cm), anyhow we stress that in Ackermann et al. 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 SSCdominated scenario, (see Fig. 19 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 a 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 to the broader spectrum of the synchrotron seed photons compared to the narrower seed photon spectrum of the external fields.
Another possible choice, able to produce a PL shape of the BC, can be obtained assuming a purely cold electron population with a truncated conical geometry, with the higher energy part of the BC being produced by the low-number electrons closer to BH, and the higher energy being 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 in order that the two spheres match the volume of the upper and lower part 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 15 cm, with the smaller spherical region corresponding to the segment of the truncated cone with an height of ≈ 5 × 10 15 cm, and the larger spherical region corresponding to the segment with an height of ≈ 8 × 10 15 cm The total number of cold electrons is of N BC ≈ 1.4 × 10 54 (see Figure 18 and left column in Table 6).Since the 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 (the statistics are reported in Table 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 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 Xray spectrum could be explained by bulk Compton emission in two-zone model, in a similar manner to MJD 54761, as the DCF analysis for the 2014-2017 points towards to.

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 gamma-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 seen that indeed also the γ-ray light curve is 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 SED modeling.
Understanding how our observational data and results fit in the current landscape of existing blazar models is a difficult task.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 postmaximum subflares in 3mm 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.ECdominated SED models seem to be favored by literature (Ackermann et al. 2012, Dreyer & Böttcher 2021).However, as we present in this paper, SSC-dominated models are also possible, as shown in section 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 correlation analysis and the results from VLBI observations favor SSCdominated 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.But the correlation between R and γ show a clear peak whose position is τ p R,γ 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 nicely 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 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; which 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. 12) 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 ended, such as some of those 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, but not 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 not a reasonable source of increased variability.
A good test to determine whether the emission is SSC or EC might be polarization of the gamma-rays.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 X-ray polarization is already being measured Notes.Uncertainties for the best-fit values were automatically obtained using the HESSE method of second derivatives and are indicated between parenthesis, parameters without them were frozen during the fit.For uncertainties smaller than the third significant digit, an upper limit is given.The rest of the parameters for the models can be found in Tables 7 and 8, together with their uncertainties and fit statistic. (a) Only the parameters of the Bulk Compton emission are shown here.See Tables 7 and 8 for the rest of the parameters.
by some instruments (IXPE), gamma-ray polarization is still not possible, although recent technological development open the possibility in the next decade.
Table 7. Parameters in the External Compton (EC) scenario for epochs MJD 54761 (Figs. 17 and 18), MJD 55098 (Fig. 20), MJD 56576 (Fig. 22) and MJD 57293 (Fig. 24).Notes.Uncertainties for the best-fit values were automatically obtained using the HESSE method of second derivatives and are indicated between parenthesis, parameters without them were frozen during the fit.For uncertainties smaller than the third significant digit, an upper limit is given.The degrees of freedom and the χ 2 statistic for the model fit are indicated in the last rows, the residuals are shown in the figures. (a) The model for this epoch includes an additional component which is independently modeled as Bulk Compton emission from the disk, the reported values refer to the SSC/EC components alone, with the exclusion of the X-ray data.Two possible geometries where considered for the EC-dominated scenario, and their parameters can be found Table 6.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 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. 13).8).Unlike the models for the older flaring epoch 54761 (Figs. 17,18), this model does not include a bulk Compton component and can be explained with only the usual SSC+EC components.The source exhibits however an important softening the X-ray spectrum that could be explained also by bulk Compton emission.

Fig. 1 .
Fig. 1.Light curves of AO 0235+164 at different wavelengths.The vertical lines mark the epochs whose SED was analyzed in sec.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.
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 99.9 % confidence and for 1mm (red) and 3mm (blue) at 90 % confidence level.The vertical lines are the same as in Fig.1and correspond to the epochs whose SED was analyzed in sec.4.4.

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

Fig. 4 .
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 .
Fig. 5. Epoch 102 (2016-09-05) VLBI 7mm image, showing total flux intensity (contours) and polarized flux intensity (color scale).Black line segments overlaid in the image represent the Electric Vector Position Angle (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, and demonstrates how the polarization angle is aligned with the direction of propagation.

Fig. 6 .
Fig. 6.Epoch 110 (2017-04-16) VLBI 7mm 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).

Fig. 7 .
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 N exp were the ones such that the reduced χ 2 -statistic was minimized (bottom).The best-fit values are given in Tables1 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 sec.4.4, as in Fig.1.

Fig. 8 .
Fig.8.Observed distance from core (up) and flux density (down) 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 for B7, which due to its low flux, did not have enough points after the peak with low enough uncertainty to perform a fit.
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.

F
Fig. 14.The 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, and also visible, although dimmer, in the 2015 flaring epoch.Both quiescent epochs lack this X-ray feature.

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

Fig. 17 .
Fig. 17.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 (Table7).The X-ray bump is modeled by bulk Compton emission from the Disk by a secondary particle distribution much closer to the central engine (Table6), consistent with the much lower correlation and higher delays in the DCFs between X-ray and the other bands (Fig.13).

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

Fig. 19 .
Fig.19.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 (Table6).

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

Table 5 .
Summary of the comparison between different models for the flaring epoch MJD 54761.

Table 6 .
Parameters for the models of the Bulk Compton (BC) emission both with simple and conical geometries, shown inFigs.19, 17, 18 for  epoch MJD 54761. ).