Extreme photometric and polarimetric variability of blazar S4 0954+65 at its maximum optical and $\gamma$-ray brightness levels

In 2022 the BL Lac object S4 0954+65 underwent a major variability phase, reaching its historical maximum brightness in the optical and $\gamma$-ray bands. We present optical photometric and polarimetric data acquired by the Whole Earth Blazar Telescope (WEBT) Collaboration from 2022 April 6 to July 6. Many episodes of unprecedented fast variability were detected, implying an upper limit to the size of the emitting region as low as $10^{-4}$ parsec. The WEBT data show rapid variability in both the degree and angle of polarization. We analyse different models to explain the polarization behaviour in the framework of a twisting jet model, which assumes that the long-term trend of the flux is produced by variations in the emitting region viewing angle. All the models can reproduce the average trend of the polarization degree, and can account for its general anticorrelation with the flux, but the dispersion of the data requires the presence of intrinsic mechanisms, such as turbulence, shocks, or magnetic reconnection. The WEBT optical data are compared to $\gamma$-ray data from the Fermi satellite. These are analysed with both fixed and adaptive binning procedures. We show that the strong correlation between optical and $\gamma$-ray data without measurable delay assumes different slopes in faint and high brightness states, and this is compatible with a scenario where in faint states we mainly see the imprint of the geometrical effects, while in bright states the synchrotron self-Compton process dominates.


INTRODUCTION
With the term "blazar" we indicate a jetted active galactic nucleus (AGN) with one jet directed towards us. Leptons moving at relativistic speeds along the magnetic field lines inside the jet produce low-energy synchrotron radiation and high-energy radiation through inverse-Compton scattering of soft photons. Processes involving hadrons may also be responsible for the high-energy emission (e.g. Böttcher et al. 2013). Because of the jet orientation, this radiation is relativistically Doppler beamed (e.g. Urry & Padovani ★ E-mail:claudia.raiteri@inaf.it 1995). Consequences of the Doppler beaming are that the flux that we observe is enhanced in comparison to what is emitted by the source, and the variability time scales are shortened. This is why blazars often show extreme variability at all wavelengths, from the radio to the rays, on time scales ranging from years to minutes (e.g. Wagner & Witzel 1995;Aharonian et al. 2007;Albert et al. 2007;Shukla & Mannheim 2020;Weaver et al. 2020). The origin of such multiscale flux changes is still debated, but it is clear that different processes must intervene to account for the variety of observed variability events. Flares suggest that particles get accelerated in the jet. The two main acceleration mechanisms that are commonly invoked are shock waves propagating in the jet (e.g. Hughes et al. 1985;Marscher & Gear 1985), and magnetic reconnection, possibly triggered by kink instabilities (e.g. Sironi et al. 2015;Zhang et al. 2018;Bodo et al. 2021;Zhang et al. 2022). Moreover, turbulence is likely to play an important role (Marscher 2014). But since Doppler beaming depends on the viewing angle, strong flux variations are expected if the jet emitting regions change their orientation with respect to the line of sight (e.g. Raiteri et al. 2017). This can happen because of jet precession, or rotation induced by orbital motion in a supermassive black hole (BH) binary system, or jet twisting due to kink instabilities developing inside the jet.
The jet physics is determined by the magnetic field, and this can be studied by means of polarimetric observations. In blazars, both the polarization degree and the polarization angle are very variable (e.g. Raiteri & Villata 2021, and references therein). Their behaviour is often uncorrelated with the total flux density, which makes the interpretation of the polarization observations very difficult.
The above variability features are common to both the blazar subclasses, i.e. flat-spectrum radio quasars (FSRQs) and BL Lacertaetype objects (BL Lacs). The two types have originally been distinguished on the basis of the strength of the emission lines in their spectra (Stickel et al. 1991;Stocke et al. 1991), the former showing broad emission lines with equivalent width greater than 5Å in the rest frame, the latter exhibiting (nearly) featureless spectra.
Polarimetric observations show that the polarization degree undergoes rapid changes at both radio (Gabuzda et al. 2000) and optical wavelengths (Morozova et al. 2014;Raiteri et al. 2021). Wide rotations of the polarization angle have been observed (Hagen-Thorn et al. 2015).
The source was detected at GeV energies by the Compton Gamma Ray Observatory (CGRO; Mukherjee et al. 1995) and at TeV energies by the Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescopes during an exceptionally bright optical state in 2015 (MAGIC Collaboration et al. 2018). Raiteri et al. (2021) analysed well-sampled optical light curves obtained in 2019-2020 by the Whole Earth Blazar Telescope 1 (WEBT; e.g. Villata et al. 2002Villata et al. , 2006Raiteri et al. 2017;Larionov et al. 2020;Jorstad et al. 2022;Raiteri et al. 2023) Collaboration, with the addition of 2-min cadence data acquired by the Transiting Exoplanet Survey Satellite (TESS) during three observing sectors of about 1 month each. They detected several characteristic times scales of variability, ranging from 6 hours in the TESS light curves, to several weeks in the whole data set. Moreover, they identified quasi-periodic oscillations (QPOs) with period of about 1 month, which were interpreted as produced by a rotating inhomogeneous helical jet, whose pitch angle changes in time.
A QPO of 1.52 d was detected in the TESS light curves by Kishore et al. (2023). They suggested that the QPO may originate from the orbital motion of some blob or flare in the innermost part of the accretion disc, and estimated a BH mass of ∼ 2×10 8 ⊙ or ∼ 10 9 ⊙ for a Schwarzschild or Kerr BH, respectively.
QPOs with periods of 66 d and 210 d were recognised by Gong  al. (2023) in the -ray light curves from the Fermi satellite. The most plausible scenario was found to be a plasma blob following a helical path inside the jet. In April 2022, S4 0954+65 was observed in very high optical (Bachev & Strigachev 2022;Vlasyuk et al. 2022a,b) and -ray (Rani et al. 2022) states. This triggered intensified observations by the WEBT (Marchini et al. 2022), whose members have been regularly monitoring the source since 2007 through the GLAST-AGILE Support Program (GASP, e.g. Villata et al. 2008Villata et al. , 2009.
In this paper we present the results of the optical monitoring by the WEBT in the period 2022 April 6 to July 6, together with theray observations by the LAT instrument onboard the Fermi satellite. We explore the details of the flux variability at both low and high energies and analyse their correlation. We also investigate the optical polarization behaviour.

OPTICAL PHOTOMETRY
Optical observations were carried out in the framework of the WEBT Collaboration at the observatories listed in Table 1.
The source magnitude was derived using the photometric sequence published by Raiteri et al. (1999). Data belonging to 31 datasets from 25 observatories in 16 countries around the northern hemisphere were carefully assembled and processed to get a homogeneous and precise optical light curve in band. No significant systematic offset of the data points of individual datasets was found with respect to the others. Particularly noisy datasets from the same telescope were binned over time intervals of a few minutes and clear outliers were eliminated. The final light curve is shown in Fig. 1; it includes 3628 data points acquired during 91 days, from April 6 to July 6 2022 (JD = 2459676.0-2459767.0). In this period, the source showed wide brightness variations, with a maximum amplitude of 3.28 mag, and reached its maximum brightness level = 12.83 ± 0.01, exceeding the levels observed during the 2015 outburst (Bachev 2015;MAGIC Collaboration et al. 2018).

INTRA-DAY VARIABILITY
As mentioned in the Introduction, S4 0954+65 has often shown intense intra-day variability. One of the most extreme episodes was reported by Bachev (2015), who observed an almost 0.7 mag brightness decrease in 5 h. In the period analysed here, we found even more dramatic fast variability.
The brightness rise that brought the source to its historical maximum = 12.83 ± 0.01 on May 9 (JD ≈ 2459708.61) involved a variation of roughly 2.8 mag in less than 2 d (47 h). Three days later we observed a 0.73 mag brightness increase in 1.2 h.
Another extreme IDV episode was detected on May 19-20 (JD 2459719), with a ∼ 0.92 mag fading in 4.8 h. Fig. 2 displays the source light curve in this period.
On June 6-7 (JD 2459737) we observed a ∼ 1.4 mag brightness increase in about 5.7 h, followed by a ∼ 0.67 mag brightness decrease in about 1.4 h, which was part of a longer dimming phase of ∼ 1.2 mag in less than 10 h. This fading trend continued in the following two days, with a total variability amplitude of about 2.5 mag in less than 2 d (47 h). An enlargement of the optical light curve in this period is shown in Fig. 3.
Overall, on 17 occasions we detected mag changes of more than 0.5 mag in less than 12 h.

SIZE OF THE OPTICAL EMITTING REGION
As detailed in the previous section, unprecedented intraday variability was observed in the period analysed. From causality arguments based on light travel time, the observed minimum variability time scale Δ min can put an upper limit to the size of the emitting region: where is the speed of light, is the Doppler factor and the source redshift. The value of Δ min can be obtained from the wellsampled, extreme IDV episode on JD 2459737, which is shown in Fig. 4. Flux densities have been obtained from magnitudes using the calibrations by Bessell et al. (1998) and correcting for Galactic absorption ( = 0.259 mag from the NASA/IPAC Extragalactic Database 2 ).
We modelled the flare according to Valtaoja et al. (1999): where 0 is the base level, the flare amplitude, peak the time of 2 https://ned.ipac.caltech.edu/ the flare peak, and Δ 1 and Δ 2 the time scales before and after the peak, respectively.
The least-squares best-fit model of the flare is obtained by the following parameters: 0 = (9.5 ± 0.5) mJy, = (8.6 ± 1.1) mJy, peak = (2459737.459 ± 0.003), Δ 1 = (0.028 ± 0.007) d, and Δ 2 = (0.012 ± 0.005) d. Setting Δ min ≈ 17 min, and = 13.6 (see Section 6 and Fig. 8) we found < 3 × 10 14 cm, i.e. about 10 −4 parsec. This upper limit to the emitting region size responsible for the flare is in general smaller than the typical size assumed for the blazar jets, and in particular it is more than 66 times smaller than that assumed by Raiteri et al. (1999) when applying the homogeneous model by Ghisellini et al. (1998) to the SEDs of S4 0954+65 during a faint state observed in 1994-1998. This suggests that we may be seeing flux fluctuations in a jet subregion.
We note that blazar microvariability with as short as a few minutes time scale was also detected at rays, in both GeV and TeV energy domains. In the case of 3C 279 observed by the Fermi satellite, a very fast flare in 2018 was ascribed to magnetic reconnection in a region of about 8 × 10 14 cm (Shukla & Mannheim 2020).   Table 1. Uncertainties are plotted in grey and are typically smaller than the symbol size.

OPTICAL POLARIMETRY
Optical polarimetric data for this work were acquired at the Belogradchik, Calar Alto, Crimean, Skinakas, and St. Petersburg observatories. The degree of polarization and the electric vector position angle (EVPA) are shown in Fig. 5 together with the optical flux densities.
Both and EVPA display a flickering behaviour. The values of range from 0.03% to 39.6%, those of EVPA span a ∼ 87 • interval. A general anticorrelation between and the flux density in band, , is recognisable, which is confirmed by the plot of versus shown in Fig. 6. Values of higher than 23% are found only when < 8 mJy, and the two values greater than 39% correspond to flux densities as low as ∼ 4 mJy. We note that values of of ∼ 40% are close to the maximum values observed in blazars. Here they are found just before and just after the very narrow flare on JD 2459737 that seems to conclude the preceding period of strong activity.
Another way of assessing the general anticorrelation between and is through the discrete correlation function (DCF; Edelson & Krolik 1988;Hufnagel & Bregman 1992). A strong correlation results in a positive DCF peak with value close to one, while an anticorrelation gives negative DCF values. As visible in Fig. 7, the DCF between and assumes negative values for time lags around zero and in general does not show positive peaks greater than 0.25.
It is interesting to compare the source polarization behaviour in 2022 with that observed in 2019-2020 and described by Raiteri et al. (2021). At that time, the source flux density was lower, oscillating between ∼ 1 mJy and ∼ 10 mJy, while was ranging from about   Table 1.  2% to 29%. The authors commented that did not seem to correlate with flux, thus the general anticorrelation seen in the current dataset was previously not detected. We plotted these earlier data in Fig. 6. We first note that in the common range of flux densities, the values of in 2019-2020 were in average lower than in 2022, which could be due to a more turbulent magnetic field in that period. Then we infer that the anticorrelation was not detected at that time because it shows up clearly only when the source exceeds the brightness values that were observed in 2019-2020.

INTERPRETATION OF THE FLUX AND POLARIZATION BEHAVIOUR
In several papers by the WEBT Collaboration we have proposed that the long-term flux variability at a given frequency is due to orientation changes of the corresponding jet emitting region, which produce a variation of the Doppler beaming. In contrast, the fast variability superimposed on the long-term trend is likely produced by intrinsic processes, such as turbulence, shocks or magnetic reconnection. We can obtain the optical long-term trend by performing a cubic spline interpolation through the light curve shown in Fig. 5, after binning it over variable time intervals that decrease with increasing brightness in order to take into account the effect of Doppler beaming on time scales (see discussion in Raiteri et al. 2017). The long-term trend derived in this way is shown in Fig. 8, and in our view it represents the amount of variability that can be ascribed to changes of the viewing angle. Because the flux density depends on the Doppler factor as ∝ + , where = 2 for a continuous jet (Urry & Padovani 1995) and = 1.8 is the source mean optical spectral index (Raiteri et al. 2021), from the long-term trend of the flux (i.e. the spline) we can derive the behaviour in time of , which is shown in Fig. 8. Moreover, we know that depends on the viewing angle as: is the bulk Lorentz factor and is the plasma bulk velocity in units of the speed of light. Therefore, we can also infer the trend of in time, once reasonable assumptions are made on the values of the other parameters. The trends of and shown in Fig. 8 were obtained by fixing Γ = 10 and min = 1.5 • . The first value is the same adopted by Raiteri et al. (2021) and is comparable with that obtained by Jorstad et al. (2017); in contrast, min was halved with respect to that in Raiteri et al. (2021) because the brightness levels in the current dataset are much higher. Indeed, at this minimum viewing angle, the Doppler factor reaches its maximum value ≈ 18.7, and consequently we observe the maximum brightness level of the long-term trend. In contrast, the lowest value of ≈ 10.6 corresponds to the largest value of ≈ 5.4 • and to the faintest long-term flux level.
What is the prediction of this geometrical scenario for the behaviour of the polarization? We investigated different possibilities. Lyutikov et al. (2005) analysed different models of relativistic jets characterised by helical magnetic fields. The authors assume a cylindrical shape for the jet, with the emitting plasma moving parallel to the jet axis. They do not consider bulk rotation of the jet. In cases where the number density of relativistic particles scales with the square of the intrinsic magnetic field, it is possible to infer (see also Raiteri et al. 2013 where sin ′ = sin . The values of max can be derived from the comparison with the observational data. In Fig. 8 we report the results obtained with max = 7, 25, and 43%. The case with max = 25% fairly reproduces the average observed polarization level. However, the data points show stronger variability. There are phases where the observed is lower than predicted, suggesting that the magnetic field was less ordered. This occurs in particular during the first stage of activity, before the flare at JD ≈ 2459690, and also before the flare at JD ≈ 2459730. In contrast, there are other phases where reaches quite high values, as just after these two flares. In this case, something must have happened in the jet to order the magnetic field. We note that the observed values of are roughly included between the models with max = 7% and 0 = 43%. We next considered the model of a rotating relativistic jet with helical magnetic field by Pariev et al. (2003). The jet is cylindrical and the magnetic field has a uniform poloidal component, while the toroidal component decreases from the centre outwards, until it disappears at the jet boundary. When assuming a value = 3 for the index of the power-law energy distribution of the emitting particles, Figure 8. Top: -band flux densities corrected for Galactic extinction (grey dots); the red line represents the long-term trend obtained as a cubic spline interpolation through the light curve binned over time intervals depending on brightness. Middle: Behaviour of the Doppler factor (red line) and viewing angle (blue line) in time derived from the long-term trend. Bottom: Observed degree of polarization (grey diamonds); the blue, red, and green thick lines represent predictions by the models in Eq. 3 with max = 25%, Eq. 4 with Ω = 7, and Eq. 5 with = 1.4, respectively. The yellow area is delimited by the models in Eq. 3 with max between 7% and 43%, which are very similar to those by Eq. 4 with Ω = 5.4-9.7, and to those by Eq. 5 with = 1.10-1.85.
the authors could derive an analytical solution: where Ω is a dimensionless parameter defining the angular rotational velocity of the magnetic field lines. Again, we derive the values of Ω from the data. As shown in Fig. 8, values of Ω = 5.4, 7, and 9.7 produce a trend of that is almost equal to that obtained with Eq. 3 and max = 7, 25, and 43%. The variation of Ω can come either from an effective change in the angular rotational velocity of the jet or possibly in its radial dependence inside the jet. Finally, we tried the shock-in-jet model (Hughes et al. 1985, see also Larionov et al. 2013 andRaiteri et al. 2013), where a random magnetic field is compressed by the passage of a shock wave. In this case: where 0 = ( + 1)/( + 5/3) is the degree of linear polarization for particles with a power-law energy distribution with index = 2 +1, and the parameter represents the degree of compression of the magnetic field by the shock wave. As in the previous cases, is set by comparison with the data. If we set = 1.4, we find the same average behaviour as in the two previous models (see Fig. 8), while by setting = 1.10 and 1.85 we can reproduce fairly well the lower and upper bounds, respectively. In summary, the three models considered above lead to the same results for reasonable choices of their parameters. They all imply a dependence of the polarization degree on the viewing angle, which anticorrelates with the Doppler factor and thus with the flux, and therefore can explain the general anticorrelation between and seen in Section 5.

-RAY OBSERVATIONS
The Fermi satellite has been monitoring the -ray sky since 2008. Light curves produced with data from its Large Area Telescope (LAT) instrument (Atwood et al. 2009) with 3-day, 1-week, and 1-month binning can be retrieved from the Fermi LAT Light Curve Repository 3 (Abdollahi et al. 2023). We show the complete 1-week binned light curve of S4 0954+65 in Fig. 9 to put our analysis into context. With this time resolution, the source -ray flux exceeded the 10 −6 ph cm −2 s −1 level only twice, at the beginning of 2015 and in 2022, and in the latter period it reached the historical maximum. This is the period we are considering in the present paper. Given the high number of -ray photons detected, we expect to be able to obtain a -ray light curve with a time resolution better than those available in the Repository, for an optimal comparison with the densely sampled optical light curve. Therefore, we first built -ray light curves with fixed integration time intervals, starting from 5 d and then decreasing the time bin down to 1 h. By combining the data from these light curves, we could build a composite light curve, with better sampling during the source high states. We checked for spectral variations, which were found to be negligible. Then we adopted an adaptive binning method with constant spectral shape, which showed to be superior than the fixed binning method from the sampling point of 3 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ LightCurveRepository/about.html view. Details on the two methods are given below. We believe that the comparison between their results can be instructive.

Fixed binning method
The -ray data in the time period considered in this paper were analysed in the 0.1-300 GeV energy range using the FermiTools package version 2.2 installed with Conda 4 , with instrument response function P8R3_V3, Galactic diffuse emission model gll_iem_v07, and isotropic background model iso_P8R3_ SOURCE_V3_v1. We performed a binned likelihood analysis, adopting a region of interest of radius 30 • , a maximum zenith angle of 90 • and only 'Source' class events (evclass=128, evtype=3). Fluxes of the sources within a 10 • radius were set as free parameters of the model, whereas fluxes of more distant sources were fixed to their mean values according to the 4FGL catalogue. As in the 4FGL catalogue, for S4 0954+65 (4FGL J0958.7+6534) we used a log-parabola model of the type: with a "break" energy b = 674.472 MeV and normalisation 0 . When integrating over the whole period, the maximum likelihood gives a test statistic TS = 4426.11, with an integrated average flux of (2.78 ± 0.025) × 10 −7 ph cm −2 s −1 and spectral parameters wp = 2.033 and wp = 0.060, very close to their catalogue values. We calculated light curves with fixed integration time bins of 5, 4, 3, 2, 1 days and 12, 6, 3, 1 hours and spectral parameters fixed to wp and wp . The source was assumed to be detected if the TS exceeded 25. Then we merged the data to get a densely-sampled composite light curve, shown in Fig. 10, starting from the 1-h binned light curve and filling the time gaps by adding data from light curves with progressively larger bins, however respecting a time distance between the data points that depends on their light curve bin. In order to investigate possible spectral changes that are known to occur in blazars, we repeated the whole procedure letting the spectral parameters free to vary. The corresponding composite light curve is shown in Fig. 10 and is very close to the previous one. The behaviour in time of the spectral parameters is displayed in Fig. 11. They show a large dispersion, sometimes reaching the boundaries set by the procedure (0 and 5 for , −5 and 10 for ). We conservatively set TS = 80 as the limit above which the values of the parameters can be trusted. By plotting these versus flux, we see that they are mostly in agreement with wp and wp , and we cannot recognise any spectral trend.

Adaptive binning method
The adaptive binning -ray light curve was computed using the standard Fermitools software 5 version 2.2.0 packaged within a FermiBottle container 6 . During the analysis we used the same instrument response function, Galactic diffuse emission model, and isotropic background model as in the fixed binning method described above. The computations were performed in the unbinned likelihood regime in the 0.1-200 GeV energy range.
The background model includes all the sources from the 4FGL catalogue that fall within 15 • radius around the target location. The fluxes of the background objects within 10 • radius were set as free parameters of the model if their significance is higher than 5.0 according to the 4FGL catalogue. The fluxes of more distant or less significant sources were fixed to their mean values according to the 4FGL catalogue. The flux of the target itself was modelled using a log-parabolic spectral energy distribution (Eq. 6), with the spectral parameters fixed to their catalogue values ( cat = 2.125, cat = 0.052, b = 674.472 MeV), and only the normalisation factor 0 was set as a free parameter to compute the flux. In order to obtain the highest possible temporal resolution of the light curve, we used an adaptive temporal binning strategy such that the periods of active state are covered with shorter time bins to obtain a fine structure of the time variability, while quieter states with low flux are covered with wider bins to obtain a signal-to-noise ratio not lower than some predefined value. In order to do so, we start the integration with a time bin as short as one hour and increase it gradually by 15 min until we reach the test statistic value TS = 25 (which roughly corresponds to = 5). When the desired TS value is reached, we stop the integration, save the current flux parameters and start a new bin until the whole time range is covered.
The adaptive binning -ray light curve is shown in Fig. 10; the time bins range from 1 h to about 6 d. The number of epochs is 163, while the composite light curve includes 102 epochs in the case of fixed parameters and 108 in the case of the variable parameters. Therefore, for the following analysis we will use the adaptive binning light curve.

COMPARISON BETWEEN THE OPTICAL AND -RAY FLUXES
The correlation between the blazar flux variations at -ray and optical frequencies has been extensively investigated, especially after the launch of Fermi, which is scanning the sky every ∼ 3 hours. The topic is of great interest, since it can provide clues on the origin of the -ray radiation, i.e. whether a leptonic or hadronic mechanism is more plausible (de Jaeger et al. 2023) and, in the first case, what is the nature of the soft seed photons (Cohen et al. 2014). However, the results are not always in compliance with the theoretical predictions, and can vary from source to source, and even for the same source observed in different periods (Raiteri et al. 2012). In most cases, a strong correlation is found, with (nearly) simultaneous variations in the two bands, as predicted by leptonic models (e.g. Raiteri et al. 2011Raiteri et al. , 2013Hovatta et al. 2014;Larionov et al. 2016;Carnerero et al. 2017;D'Ammando et al. 2019). However, in several cases the -ray variations were observed to lead those in the optical band, especially in FSRQs (Hayashida et al. 2012;Cohen et al. 2014;Carnerero et al. 2015;Larionov et al. 2020), but there were also cases where the -ray changes appeared delayed (Jorstad et al. 2011;Cohen et al. 2014).
A visual comparison between the optical and -ray fluxes of S4 0954+65 in Fig. 10 shows a general good agreement, with the exception of the first optical flaring phase at JD ∼ 2459680-2459699, which has only a minor -ray counterpart.
We investigate the -optical correlation with the DCF. The result is shown in Fig. 12, where the optical data have previously been averaged over 1 h and the DCF is calculated over 1 d bins. The main peak at r ∼ 0.7 indicates fair correlation with no time delay (time lag = 0 d). This suggests that -ray photons can indeed be produced by inverse-Compton scattering of soft photons off the same relativistic electrons that are responsible for the optical synchrotron photons, as predicted by leptonic models. A deeper analysis can reveal further details on the source variability and on the nature of the soft photons seeds. We transformed the -ray fluxes from counts to physical units taking into account that where d /d is given in Eq. 6. The -ray fluxes at 1 GeV (log = 23.383) were then compared to the optical fluxes in the band (log = 14.670). We paired each -ray data point of the adaptive binning light curve with the mean of the optical fluxes included in the same bin interval. In this way we obtained 97 -optical pairs, where for each -ray data point we averaged from 1 to 470 optical points. The result is shown in Fig. 13. The linear regression line on all pairs has a slope of 1.48 ± 0.05, but the scatter is rather large. This is actually expected, as two different mechanisms are likely acting. One is the geometrical effect due to the variation of the viewing angle discussed in Section 6, which would give a slope equal to 1. The other is the synchrotron self-Compton (SSC) process, where the soft photons that are inverse-Compton scattered up toray energies are the same optical photons (Maraschi et al. 1992;Bloom & Marscher 1996). This would yield a slope of 2. Because of the presence of these two mechanisms, the distribution of data points would approximately lie within a parallelogram, as shown by Larionov et al. (2016). Therefore, a cubic regression (see Fig. 13) is a more suitable fit, since it can follow the different slopes at the various brightness states.
To check for the double nature of the -optical correlation, we considered high and low brightness states separately. The slope of the linear regression on the data points corresponding to the brightest optical states only (log > −10.5, see Fig. 10) is 1.98±0.12, while in the case of the faintest states the slope decreases to 1.05 ± 0.19. This matches well the trend of the cubic fit and can be understood as follows. In faint states, longer integration intervals are required to get significant -ray signals, and also the optical data to be paired with them are consequently averaged on long time bins. This is highlighted by the large horizontal bars in Fig. 13, which represent the standard deviations of the optical data averaged in the time bin of the corresponding -ray data point. Therefore, the fast, intrinsic flux changes are smoothed out, and the slope becomes close to 1, in agreement with a geometrical origin of the long-term flux variations Figure 13. Gamma-ray fluxes at 1 GeV versus the -band fluxes. The red/blue data points correspond to optical data lower/higher than log = −10.5. Cyan-filled symbols represent data before JD = 2459699. The grey, blue, red, and cyan lines are linear regressions to: all the data points, the brightest data points, the faintest data points, and the data points up to JD = 2459699, respectively. The green line shows a cubic fit to all the data points. Vertical grey bars represent errors on the -ray fluxes; horizontal bars indicate standard deviations on the mean of the optical fluxes paired with the -ray ones. (Raiteri et al. 2017). We note that a similar slope is found during the first optical active phase, which has only a minor -ray counterpart, which smooths the variability. Indeed, the optical-correlation before JD = 2459699 has a power-law index 0.97 ± 0.18 (see Fig. 13). In contrast, during the brightest phases, the time bins become shorter and it is possible to appreciate the squared dependence of the -ray fluxes on the optical ones due to the SSC process, whose signature prevails over the geometric effect.

DISCUSSION AND CONCLUSIONS
We have presented the photometric and polarimetric optical data obtained by the WEBT Collaboration during a very active phase of the BL Lac-type object S4 0954+65 in 2022, together with the -ray data from the Fermi satellite. In this period the source reached its historical brightness maxima in both optical and -ray band.
Many unprecedented, extreme episodes of optical intraday variability were detected. A model fit to one of the fastest flares implies a variability time scale of about 17 min and thus a size of the optical emitting region less than 3 × 10 14 cm (about 10 −4 parsec). This means that we are likely observing emission coming from jet subregions, confirming earlier suggestions arising from the detection of minute-scale variability also at GeV (Shukla & Mannheim 2020) and TeV energies (Aharonian et al. 2007;Albert et al. 2007).
We have analysed the behaviour of the optical polarization, which shows large fluctuations in both the polarization degree and angle, and a general anticorrelation between and the optical flux density . The presence of this defined trend together with strong dispersion around it may indicate that we are seeing the combination of different processes.
We have shown that if the long-term trend of the optical flux density is due to a variation of the Doppler factor caused by orientation changes of the jet emitting region, then a simple helical magnetic field in a possibly rotating jet can well describe the observed average behaviour of . The same result can also be obtained in the frame-work of a shock-in-jet model. All the polarization models that we considered predict that the trend of follows that of the viewing angle . Therefore, since anticorrelates with , which in our view determines the long-term behaviour of the source brightness, these models naturally explain the general anticorrelation between and .
However, there are periods in which the observed is lower than the average trend predicted by the models, and periods where it is higher. One way to explain the whole range of values is to assume a variation of the model parameters, such as the jet angular rotational velocity or the strength of the shock waves. More realistically, when becomes much lower than predicted, we can imagine that some process is leading to a less ordered magnetic field, most likely increasing its turbulent component. This would also be consistent with our finding that the emission comes from jet sub-regions, because turbulence can be seen as the overlapping of multiple stochastic contributions (e.g. Marscher 2014). In contrast, the observation of values of higher than predicted suggests that something has produced a more ordered field. This can occur when shock waves propagate in the jet, compressing the magnetic field lines (Hughes et al. 1985;Marscher & Gear 1985), or in any other case when the field lines appear more parallel along the line of sight through the jet.
An alternative mechanism that can produce flux and polarization variability is magnetic reconnection. According to the simulations by Zhang et al. (2020), who explored the magnetic reconnection flux and polarization signatures in relativistic jets, the outcome of magnetic reconnection strongly depends on the model parameters, i.e. on the physical conditions. In particular, the value of the polarization degree is very sensitive to the strength of the guide field, which represents the component of the magnetic field which is perpendicular to the reconnecting magnetic field. This mechanism could also explain fast flux variability from small size emitting regions, since reconnection can lead to the formation of a large number of plasmoids (Giannios 2013).
The comparison between optical and -ray data leads to results that are consistent with our geometrical interpretation of the long-term trend. The optical and -ray flux variations are well correlated with no appreciable time delay, which means that the observed optical and -ray photons come from the same jet region. Moreover, the power-law dependence of the -ray fluxes on the optical ones has an index ≈ 1 during faint states, where the variability is dominated by the long-term trend because of the low -ray statistics. An index one is what is expected if the long-term trend is due to orientation changes. In contrast, the index becomes ≈ 2 in bright states, when the higher statistics allow us to detect what is most likely the signature of the intrinsic SSC process. We note that an SSC nature of the -ray emission in S4 0954+65 is in agreement with previous results that favour an SSC mechanism for BL Lacs in general, and explains the -optical correlation with no appreciable time delay (e.g., Cohen et al. 2014;Hovatta et al. 2014, and references therein).
In conclusion, we interpret the long-term optical flux and polarization behaviour as the result of variations in the viewing angle of the optical emitting region. The short-term variability would instead be produced by energetic processes occurring inside the jet. Polarimetry suggests that there are periods where turbulence dominates the optical emitting region, while in other periods the magnetic field becomes more ordered maybe because of the passage of a shock wave. Magnetic reconnection could also be a viable explanation for the short-term photometric and polarimetric variability. The -ray emission, which comes from the same emitting region, seems to confirm the geometric nature of the long-term trend, and during the brightest states reveals its SSC nature. vices. HG and KW acknowledge financial support from the Department of Science and Technology, India, through an INSPIRE Faculty Award IFA17-PH197 at ARIES, Nainital, India. The -band photometric data from the University of Athens Observatory (UOAO) were obtained after utilizing the robotic and remotely controlled instruments at the facilities (Gazeas 2016 -2010-5, No. 265772), the observing and financial grant support from the Institute of Astronomy and Rozhen NAO BAS through the bilateral SANU-BAN joint research project "GAIA astrometry and fast variable astronomical objects", and support by the SANU project F-187. MDJ thanks the Brigham Young University Department of Physics and Astronomy for continued support of the ongoing AGN monitoring program at the West Mountain Observatory. JAP acknowledge financial support from the Spanish Ministry of Science and Innovation (MICINN) through the Spanish State Research Agency, under Severo Ochoa Program 2020-2023 (CEX2019-000920-S). This paper is partly based on observations made with the IAC-80 telescope operated on the island of Tenerife by the Instituto de Astrof'isica de Canarias in the Spanish Observatorio del Teide and on observations made with the LCOGT 0.4 m telescope network, one of whose nodes is located in the Spanish Observatorio del Teide. EB acknowledges support from DGAPA-PAPIIT GRANT IN119123. This work is partly based upon observations carried out at the Observatorio Astronómico Nacional on the Sierra San Pedro Mártir (OAN-SPM), Baja California, Mexico. VJ acknowledges the support provided by the Department of Science and Technology (DST) under the "Fund for Improvement of S & T Infrastructure (FIST)" program (SR/FST/PS-I/2022/208).

DATA AVAILABILITY
Data acquired by the WEBT collaboration are stored in the WEBT archive and are available upon request to the WEBT President Massimo Villata (massimo.villata@inaf.it). Fermi data can be downloaded from the National Aeronautics and Space Administration (NASA) site (https://fermi.gsfc.nasa.gov/ ssc/data/access/) and analysed with the Fermitools software ( https://fermi.gsfc.nasa.gov/ssc/data/analysis/ documentation/); the -ray light curves published here can be obtained from the authors under request.