Precession-induced Variability in AGN Jets and OJ 287

The combined study of the flaring of Active Galactic Nuclei (AGN) at radio wavelengths and pc-scale jet kinematics with Very Long Baseline Interferometry (VLBI) has led to the view that i) the observed flares are associated with ejections of synchrotron blobs from the core, and ii) most of the flaring would follow a one-to-one correlation with the component ejection. Recent results have provided mounting evidence that the quasi-regular component injections into the relativistic jet may not be the only cause of the flux variability. We propose that AGN flux variability and jet morphology changes can both be of deterministic nature, i.e. having a geometric/kinetic origin linked to the time-variable Doppler beaming of the jet emission as its direction changes due to precession (and nutation). The physics of the underlying jet leads to shocks, instabilities, or to ejections of plasmoids. The appearance (morphology, flux, etc.) of the jet can, however, be strongly affected and modulated by precession. We demonstrate this modulating power of precession for OJ 287. For the first time, we show that the spectral state of the Spectral Energy Distribution (SED) can be directly related to the jet's precession phase. We model the SED evolution and reproduce the precession parameters. Further, we apply our precession model to eleven prominent AGN. We show that for OJ 287 precession seems to dominate the long-term variability ($\gtrsim 1\,{\rm yr}$) of the AGN flux, SED spectral state, and jet morphology, while stochastic processes affect the variability on short timescales ($\lesssim 0.2\,{\rm yr}$).


INTRODUCTION
The term blazar encompasses BL Lac Objects and Flat-Spectrum Radio Quasars into a subgroup of active galactic nuclei (AGN). The identification of blazars as jets of relativistic plasma oriented close to our line-ofsight has helped to establish the current paradigm of AGN variability (Rees 1978). Many fundamental ques-tions regarding the jet phenomena, however, remain unanswered. In particular, key questions remain concerning the physical connections between the innermost constituents of an AGN -the compact jet, the accretion disk, and the central SMBH (see Karas et al. 2021 for a recent review and references therein).
The central engine is directly linked to the issue of the origin of the temporal flux variability and the emission observed in different wavebands across the electromagnetic spectrum. In order to probe this, radio monitoring campaigns of blazars and other AGN with singledish telescopes have been increasingly combined with detailed radio interferometric monitoring of the morphological changes witnessed in parsec-scale jets. While the former yields extensive information on radio flux variability on different time scales and frequencies, the latter reveals the kinematics of the radio-emitting components in the parsec-scale jet. It is indeed necessary to combine the temporal and the structural information to make significant claims on the deterministic nature and/or the quasi-periodicity of the radio flares since light curves alone often exhibit false periodicities that are purely of stochastic nature rather than deterministic (see e.g., Vaughan et al. 2016;Ait Benkhali et al. 2020;Tarnopolski et al. 2020).
OJ 287 is one of the most promising candidates for harbouring a supermassive binary black hole (SMBBH) (see e.g. Sillanpaa et al. 1988;Valtonen et al. 2009;Britzen et al. 2018;Gómez et al. 2021), although there are other plausible scenarios to account for its characteristic periodicities in the optical domain ( 11 years). Britzen et al. (2018) find that its radio jet is precessing on a time scale of 22 − 23yr. The data is also consistent with an additional jet-axis rotation on an approximately yearly time scale. Besides jet rotation, nutation of the jet could also explain the observed second-order motion of the jet axis. The bulk jet precession explains the variability of the total continuum radio flux density via viewing angle changes, which leads to variable Doppler beaming. An SMBBH or Lense-Thirring precession (misaligned disc around a single BH) seem to be required to explain the time scale of the precessing motion. Such a scenario intrinsically connects the radio and the optical periods, with the radio period being twice the optical period when viewed close to the rotation axis. Zhao et al. (2022) recently published the first GMVA and ALMA observations of OJ 287 observed at 3.5 mm. The radio structure resembles a precessing jet in projection which seems to confirm the precessing jet proposed by Britzen et al. (2018).
In this paper, we focus on investigating the radio flux variability in conjunction with the observed structural changes in the parsec-scale radio jets of OJ 287 and other prominent AGN. This has led us to propose a possible alternative to the currently popular paradigm that the dominant cause of blazar flux variability is the occurrence of internal shocks within the relativistic nuclear jets (the so-called "shock-in-jet" model, see e.g. Marscher & Gear 1985;Hughes et al. 1985;Qian et al. 1991).
Here we argue that a major, if not dominant contribution to flux variations of most blazars on longer timescales, may come from the precession effects (which are, in principle, of deterministic type). A number of observational findings as well as recent results of general relativistic magneto-hydrodynamic simulations reinforce this proposed explanation for AGN radio variability.
The present work argues for jet precession, which, however, does not explain by itself the physical nature of jet components. Radio knots could manifest themselves as shocks in jets (see e.g. Marscher & Gear 1985;Hughes et al. 1985;Qian et al. 1991, instabilities (see e.g. Ferrari et al. 1978Hardee 1987), or plasmoid injections due to magnetic reconnection (see e.g. Comisso et al. 2017;Vourellis et al. 2019;Nathanail et al. 2020a;Ripperda et al. 2022). Both mechanisms -shock-in-jet and precession -can operate at the same time. However, shocks are more likely if injections are always in the same direction. In case of precession, components are moving in different directions and do not collide and thus are less likely to produce shocks. While the jet might be made of shocks or plasmoids, its wiggly structure and fluxdensity variability of the various jet features are very likely modulated (boosted) by precession. Without precession, jets would appear more straight (ignoring for the moment jet interactions or other phenomena that could make a jet turn), AGN variability aperiodic, and less pronounced. Stochastic processes will add to the variability, which predominantly appears organized by a deterministic process.
The approach of a precessing jet model was already applied to several sources to explain the radio variability as well as smoothed jet morphological changes (Abraham 2000;Abraham & Carrara 1998;Abraham & Romero 1999;Caproni & Abraham 2004a,b;Caproni et al. 2007Caproni et al. , 2013Caproni et al. , 2017. We have successfully applied this model to the "Rosetta stone" blazar OJ 287 (Britzen et al. 2018), the central galaxy of the Perseus cluster 3C84 (Britzen et al. 2019a), as well as to blazars exhibiting high-energy neutrinos (Britzen et al. 2019b(Britzen et al. , 2021. In addition, jet precession has been invoked for jetted tidal disruption events, where the infalling star is generally misaligned with respect to the equatorial plane of a rotating supermassive black hole (Stone & Loeb 2012;Franchini et al. 2016). This shows that jet precession is likely an integral part of AGN engines, which is not surprising since the main trigger for disk-jet precession is a secondary massive black hole (BH) or a misaligned accretion disk (Abraham 2018), both of which are integral and recurring features of galaxy evolution.
We will first briefly describe our precessing "nozzle" (jet base) model originally developed for the BL Lac object OJ 287. The model directly relates the evolving VLBI structure (radio components in the parsec-scale jet) to the flux variability seen in the integrated lightcurve. Thus, in this framework, both phenomena -flux variability and jet evolution -are coupled and represent two manifestations of the same physical process. For OJ 287, we show that the crossing of the projected reversal point is the most important precessional phase of the VLBI jet and it is related to a new phase in the Spectral Energy Distribution (SED).
In the following, we shall compare and discuss the most prominent reported cases of AGN where a coupling between flux variability and jet evolution may be at work (our own work, supplemented with data taken from the literature). We shall also examine the precession scenario from the perspective of the multi-wavelength spectral energy distribution (SED).
We will then summarise our results for the other blazars in our sample and recapitulate some mechanisms proposed for the jet precession -the most prominent of them being the presence of a secondary BH in the proximity of the central engine. Finally, we discuss some possible implications of our findings to the AGN physics and potential follow-up studies. While our focus here remains on cm-wavelength observations (VLBI and single-dish), we shall also mention some recent multiwavelength observations of OJ 287, which offer support to our model.
The geometric modelling of the jet evolution has parallels in other contexts related to blazars. For instance, we recall the population of Extreme High-frequency Peaked BL Lac objects (EHBLs) which emit a substantial fraction of their power in the GeV-TeV range (e.g., Costamante et al. (2001)). The synchrotron component in the SEDs of these relatively low-power blazars peaks in the EUV, X-ray, or even the MeV band (Ghisellini 1999). The origin of this high-energy emission is yet to be well understood before the basic difference of these extreme blazars from the more typical ones can be grasped. For the jets of these blazars, VLBI measurements commonly reveal low (i.e., at most mildly superluminal) apparent speeds (e.g., Piner & Edwards 2004, 2018Giroletti et al. 2004;Saugé & Henri 2004). Such speeds are in stark contrast to the large bulk Lorenz factors (≥50) commonly deduced for their jets, e.g., from gamma-ray transparency arguments (Tavecchio et al. 1998;Henri & Saugé 2006). In order to explain this huge discrepancy, a geometrical argument, namely stratified (spinesheath) jet structures (e.g., Swain et al. 1998;Chiaberge et al. 2000;Laing et al. 1999;Giroletti et al. 2004;Britzen et al. 2021) has often been invoked. It has also been shown that the inferences about the jet kinematics and viewing angle, made from observations of apparent speed and brightness of the VLBI knots and flux variability, can depend quite sensitively on another geometrical factor, namely a conical shape of the jet (in particular uniform or stratified, see e.g., Gopal-Krishna et al. 2004Boutelier et al. 2011).
Thus, it seems quite plausible that geometric effects, such as those explored here, can play a vital role in defining the general phenomenology of blazar jets. The manuscript is structured as follows. In Section 2, we introduce the precession-nutation model and apply it to radio jet components of OJ 287 as well as to its radio light curves. In Section 3, we show that the detected change of the SED shape of OJ 287 can be addressed by the Doppler factor variation as the jet precesses. In Section 2, we analyze power density spectra and periodograms of OJ 287. We can detect potentially significant periods in periodograms that correspond approximately to precession-nutation periods as well as to the putative orbital period of the SMBH binary or its multiples. In addition, we list precession-model parameters for 12 sources in Subsection 4.2 and we compare Doppler factors inferred from precession model with those based on the "shock-in-jet" model in Subsection 4.3. We discuss our results in the broader context in Section 5 and conclude with Section 6.

MODELLING PRECESSION AND NUTATION
Precession of AGN jets can generally be caused by a misaligned binary supermassive black-hole system or the Lense-Thirring (hereafter LT) effect (Thirring 1918;Lense & Thirring 1918;Abraham 2018). The first effect is by nature classical -the accretion disk orbiting the primary SMBH is perturbed by a secondary black hole whose orbital plane is misaligned with respect to the accretion disk plane by angle Ω. Due to induced torques, the accretion disk precesses which is followed by the jet precession along the precession cone with a half-opening angle of Ω. The Lense-Thirring effect is by nature relativistic and is induced by the frame-dragging effect acting on the accretion disk misaligned with respect to the equatorial plane of the central Kerr black hole. Both effects can cause precession as well as nutation of the disk-jet system. By jet precession, we understand the  periodic bulk motion of the jet with respect to the ambient medium, while nutation is understood as (i) the "nodding motion" of the jet components with respect to the jet frame or (ii) as a secondary bulk motion of the jet with respect to the primary precession motion. The "nodding motion" of the nutation and rotation may be hard to disentangle (in time and space). Kinematically, these two cases of secondary, nutation-like motion are difficult to distinguish as they are manifested by the same modulation of the time-dependent Dopplerbeaming factor with respect to the observer. The physical distinction depends on the type of intrinsic kinematics of jet components -in case-(i) nutation components rotate in the jet frame, while in case-(ii) nutation components move rather linearly in the jet frame but they are rotated around the mean jet motion by an external per-turbation, which leads to the same "deferent-epicycle" motion as in case (i). Accreting binary systems composed of a compact object and a normal star can also be sources of jets (Mirabel & Rodríguez 1999). The microquasar SS433, with an orbital period of 13.1 days (Crampton et al. 1980;Cherepashchuk 1981), is of particular interest in the context of our work, because the accretion disk and the jet in this source are known to exhibit both precession and nutation (see, e.g., Fabrika 2004 for a review). Precession was first detected in the periodically Doppler shifted Balmer and He I emission lines (Abell & Margon 1979).
The precession-nutation model we will apply in this manuscript is of general kinematic nature and can be applied independently of the origin of the precession-nutation motion. Here we only discuss the kinematic model and do not investigate the reason of the precession nor do we consider different effects due to a BZ- (Blandford & Znajek 1977) or BP- (Blandford & Payne 1982)jet. We assume that the disk-jet system is intrinsically connected, and in particular that the disk precession leads to the jet precession. This assumption of the precessing disk-jet connection has recently been justified by the 3D general relativistic magnetohydrodynamic simulations (Liska et al. 2018).
A concise mathematical description of our jet precession and nutation modelling was first applied in Britzen et al. (2018). We here summarize the basic concept and refer the reader to Britzen et al. (2018) for details of the modeling. For a schematic view, we refer the reader to Fig. 1, which illustrates the basic parameters of the model. In order to model the jet precession and nutation in the observer's frame we apply a set of eight physical quantities (see Table 1), specifically the precession period P p , the nutation period P n , the Lorentz factor γ, a half-opening angle Ω p of the precession cone, a halfopening angle Ω n of the nutation cone for case-(ii) nutation or it coincides with the jet half-opening angle for the jet rotation, the angle Φ 0 between the line of sight and the precession-cone axis (viewing angle), and the position angle η 0 of the precession-cone axis projected on the sky plane. Finally, the jet position is expressed with respect to the reference epoch t 0 that is inferred from fitting the model to the data.
The parameters defining the geometry of the system -in particular the angles Ω p , Ω n , Φ 0 , η 0 -are depicted in Fig. 1 (top right panel). The bottom-right panel is a zoom-in into the central engine that depicts the two most common drivers of the jet precession -(a) a secondary massive black hole whose orbital plane is misaligned with respect to the accretion-disk plane around the primary black hole and (b) the Lense-Thirring precession due to a frame dragging acting on a misaligned accretion disk with respect to the spin of a (single) SMBH (Britzen et al. 2018).
The nutation of the jet as a second-order precession effect can be observed as an extra rotation-like motion of jet components in the jet frame. Therefore, the possible existence of nutation in precessing jets was motivated by a rotation-like motion of the so-called stationary component a (see Fig. 5a in Britzen et al., 2018, for the plot on the motion perpendicular to the jet ridge line). The term stationary is usually attributed to jet components which do not show any motion over longer periods of time, as in Mrk 501 (see, e.g., Edwards & Piner (2002), Britzen et al. (2023). For clarity, we call component a quasi-stationary component in the following text.

Modelling of the jet components in OJ 287
To verify the applicability of the jetprecession/nutation model, we refit the OJ 287 jet component viewing and position angles Φ and η inferred for the component ejection times by Gabuzda et al. (1989), Tateyama et al. (1999), andBritzen et al. (2018). Specifically, we fit β app (t) and η(t) simultaneously to the temporal evolution of the jet apparent velocities and the position angles at the time of the component ejection. The best-fit parameters are inferred using the global least-square fitting with the concatenated χ 2 -statistic, χ 2 = (χ 2 β , χ 2 η ), where χ 2 β is the χ 2 -statistic corresponding to the apparent velocities and χ 2 η is the χ 2 -statistic corresponding to the position angles. The fitted parameters converge well, however, with a certain degeneracy close to the minimum χ 2 , i.e. more solutions are possible with different χ 2 values. The parameter values are, however, comparable within the uncertainties, which justifies the application of the precession model for describing the temporal evolution of the OJ 287 VLBI data for a time period spanning 40 years.
To illustrate this better, we include two solutions, one with the reduced χ 2 r = χ 2 ν,β + χ 2 ν,η of 18.9 and the other with χ 2 r = 167.0. The number of the degrees of freedom for the apparent velocities is 15, while for the position angles, it is 16. The application of the precession model to the data is graphically depicted in Fig. 2, where we also highlight the two exemplary fits according to the legend. The best-fit parameter values for both solutions are listed in Table 1. These two fits, although differing considerably in terms of χ 2 r , have consistent precession parameters within 1σ uncertainties.
In Fig. 1 (top left panel), we use Eq. (8) from Britzen et al. (2018) to fit the observed position-angle evolution of the quasi-stationary component a and we obtain P n and Ω n (see Table 1) in addition to the precession parameters that are consistent within uncertainties with the precession-only model (see also Britzen et al. 2018), except for the precession period and the viewing angle of the precession axis that differ by more than 1σ. These differences can be interpreted by the lack of the velocity information that is not included in the fit since this component appears to be stationary with regard to the radial distance from the core. The positional-angle fit is shown along with the residuals in Fig. 3 in more detail. The best-fit χ 2 r is 2.3 for 86 degrees of freedom. When the evolution of the quasi-stationary component is fitted with the precession-only model, the residuals are clearly larger. The data appear to favour the half-opening angle of the nutation cone of at least ∼ 3 • , potentially as much as 5 • , see the model with Ω n = 5 • in Fig. 3 Figure 2. The best-fit solutions of the precession model applied to the apparent velocities (left panel; expressed in terms of the speed of light) and the position angles (right panel; expressed in degrees) of the jet components of OJ 287. One model solution with the reduced χ 2 r = 18.9 (χ 2 r = χ 2 r,β + χ 2 r,η = 12.7 + 6.2 with 15 and 16 degrees of freedom for β and η fitting, respectively) is depicted by a solid black line, while the other solution with χ 2 r = 167.0 (χ 2 r = χ 2 r,β + χ 2 r,η = 25.4 + 141.6 with 15 and 16 degrees of freedom for β and η fitting, respectively) is represented by a solid blue line. The bottom panels show the corresponding residuals (data-model). The best-fit parameters are listed in Table 1. While the black solution provides a formally better fit, the residuals are comparable as well as the best-fit parameters within the uncertainties, see Table 1 . Table 1. List of parameters for the precession-only and the combined precession-nutation jet model of OJ 287 based on the refitted model in this work and the original analysis presented in Britzen et al. (2018). For the precession-nutation model, we include the specific values of the parameters inferred for the quasi-stationary component a in parantheses.
Complementary to the combined least-square fitting, we apply the Bayesian analysis using Markov Chain Monte Carlo (MCMC) inference of precession and precession-nutation parameters. The advantage of the method is that based on the prior flat distributions of relevant parameters, we obtain marginalized posterior one-dimensional likelihood distributions and twodimensional likelihood contours, which provide information about how tighly each parameter is constrained. To this goal, we maximize the likelihood function L using the form, where q obs i are observed quantities (here the position angles and apparent velocities), q theor i (p) is the theoretical model that depends on precession and nutation parameters p, s 2 i = σ 2 i,q + exp (2 ln f )(q theor i ) 2 is the variance, which contains an underestimation factor ln f apart from the measurement errors σ i,q . The quantity q theor (p) depends on the precession-model parameters p, which we already introduced. To find L max and infer p, we apply the python-based MCMC solver emcee.
The non-zero flat priors and posterior likelihood distributions for each parameter are shown in Appendix A separately for the case of the position angle and apparent velocity precession models. We also apply the MCMC method to infer precession-nutation parameters based on the temporal evolution of the quasi-stationary component a. The inferred posterior parameters for each case with 1σ uncertainties are summarized in Table 2. We also list the obtained values of ln L max , the number of the model parameters k, the dataset size N , and Akaike and Bayesian information criteria defines as, Overall, the inferred precession-nutation parameters are consistent for all three datasets: the precession period is ∼ 20 − 30 years, the half-opening angle of the precession cone is ∼ 10 • , the viewing angle of the precession cone axis is close to the line of sight, Φ 0 ∼ 10 • −20 • or the similar value when subtracted from 180 • , and the position angle of the precession cone axis less than −100 • . From the apparent velocities, we obtain a tight constraint on the Lorentz factor of the moving jet components, γ 9. From fitting the quasi-stationary component a, we obtain the nutation period of P n ∼ 1.6 years and the half-opening angle of the nutation cone of Ω n ∼ 2.2 • . When we inspect the likelihood distributions in corner plots in Appendix A, we see that the viewing angle Φ 0 is not well determined for the position angle fitting. Both the viewing and the position angles are degenerate while fitting apparent velocities. On the other hand, when fitting the quasi-stationary component a, all the parameters are reasonably well determined. The prior values for the parameters were iteratively narrowed down during the comparison of median maximum-likelihood models with the actual data in Figs. 15, 17, and 19 when fitting position angles, apparent velocities, and the quasi-stationary component position angles, respectively.

Modelling radio light curves of OJ 287
The time-variable jet kinematics that is predicted by the precession and nutation translates into the observed flux density variations via the Doppler-boosting factor Assuming that the underlying intrinsic jet emission can be described by synchrotron emission, S ν,0 ∝ ν −α with the spectral index α, the time-variable Doppler-boosted flux density can be expressed as S ν (t) = S ν,0 δ p+α , where p is the geometry factor for boosting, which adopts the value of p 2 for a uniform cylindrical jet emission and p 3 for discrete, isotropically radiating jet components (Abraham 2000(Abraham , 2001; for the original derivation of the Doppler beaming, see McCrea (1972) and Blandford & Königl (1979). The modulation of S ν (t) by the Doppler-boosting factor δ(t, γ, Φ) leads to the continuum outbursts with the characteristic periodicity of P p and P n that can potentially be traced in the radio light curves of OJ 287 (see Fig. 1, bottom left panel).
First, we adopt the kinematical parameters corresponding to the best-fit precession-only model (with χ 2 r = 18.9; see Table 1) as inferred from the joint leastsquare fitting of the apparent velocities and position angles. Then we use the time-variable flux-density function S ν = S ν,0 δ p+α to fit the temporal evolution of the flux density of OJ 287 at 4.5, 8.0, and 14.5+15.0 GHz (combined UMRAO+OVRO radio data), where we fix p = 2 since p = 3 leads to unrealistically large negative spectral indices (steep inverted spectrum). The best-fit parameters S prec ν,0 , α prec , and t prec 0 are listed in Table 3 for the three frequencies, while the precession-model flux densities vs. observed flux densities are shown in Fig. 4. The precession-only model reproduces well the main observed radio flares, i.e. the intrinsic jet flux density is clearly modulated by the bulk jet precession or at least it can address the increased radio flux density at   Table 2. Summary of the MCMC-inferred precession and precession-nutation parameters based on the observed evolution of the position angles η, the apparent velocities βapp, and the position angle evolution of the quasi-stationary component a. We also include the information about ln Lmax, AIC, and BIC statistical parameters.   Second, we add the second-order nutation-like motion to the bulk precession motion of the jet, see Fig. 1, top right, for an illustration. This is motivated by significant residuals of the order of 1 Jy between the observed flux densities and the precession-only model, see Fig. 4. In addition, the quasi-stationary radio component a of OJ 287 exhibits signs of a second-order, short-period motion, see Fig. 3. For the precession kinematics, we fix the parameter values close to the bestfit precession model according to Table 1, in particular γ = 9.0, P p = 22.9 yr, Ω p = 10.3 • , Φ 0 = 15 • , and η 0 = −130.8 • . The nutation-cone half-opening angle is set to Ω n = 3 • , which is within the uncertainties consistent with the best-fit value obtained from fitting the position-angle temporal variation of the quasi-stationary component a, see Fig. 3. The geometric boosting parameter is set to p = 2.0 as in the precession-only model, which corresponds to the continuous (cylindrical) jet emission model. We fit the precession+nutation flux density modulation, S ν = S ν,0 δ p+α , to the observed continuum flux densities to infer the intrinsic jet base flux density S pn ν,0 , spectral index α pn , nutation period P n , and the reference epoch t pn 0 . The best-fit parameters for 4.5, 8.0, and 14.5+15.0 GHz data are listed in Table 4 and the best-fit models are graphically depicted in Fig. 5, including the residuals, for each frequency band. From the best-fit models, we obtain the inverted spectral index with the mean value of α pn = −2.38 ± 0.07, which is consistent with the self-absorbed, inverted radio spectrum for the precession-only model. The mean period of the short-term nutation motion is P n = 1.66 ± 0.14 years, which is consistent with the nutation period of 1.6 ± 0.1 years inferred from fitting the position angle of component a, see Table 1. Adding the nutation to the bulk jet precession kinematics can temporally reproduce the short-term flux-density variability; however, the decrease of the mean values of residuals -0.61 (4.5 GHz), 0.93 (8.0 GHz), and 1.18 Jy (14.5 GHz) -is only mild or essentially negligible, which is due to the fact that the nutation does not reproduce well the observed short-term flare amplitudes. Therefore, the short-term radio variability is not yet fully captured by the pre-cession+nutation model. This may be due to the fact that the current model does not include stochastic radio variability, nor do we perform radiative transfer that can also modify the radio variability as such. Most of the residuals are due to a few large radio flares that exceed the Doppler-boosted flux density level at the peaks by about a factor of two. This is especially apparent towards higher frequencies.
In Section 4, we analyze in more detail the coexistence of purely stochastic radio variability with the deterministic, quasi-periodic kinematical effects, such as the precession and the nutation.

Precession and nutation in blazar jets: a case for binary systems
Since the orbital period of a putative supermassive binary black hole (SMBBH) (for example OJ 287) is not negligible as compared to the precession period of the accretion disk and the jet, the presence of a smaller nodding motion of the jet -here referred to as nutationis generally expected, e.g., in analogy to the well known binary systems SS433 (Margon 1984) and the accretion disk of the X-ray pulsar Hercules X-1, as proposed by Shakura et al. (1999).
Essentially, the gravitational torque of the secondary BH acts on the accretion disk around the primary BH (with a jet). When averaged over the binary orbital period, this yields a steady torque resulting in a mean counter-precession of the accretion disk (opposite sense to the disk rotation). This steady torque is accompanied by a nutation torque that is periodic at the synodic period derived from the precession frequency and the binary's orbital frequency. This nutation torque has the same magnitude as the precession torque, but the  , which is depicted from the left to the right panels, respectively. The legends list the reduced χ 2 values with the corresponding number of degrees of freedom. The best-fit parameters for each frequency band are listed in Table 3. The best-fit spectral indices are negative for all three frequencies, i.e. the spectrum is inverted, which is consistent with the increasing flux density for higher frequencies as is visible from the left to the right panels in the sequence of the increasing frequency.  5 GHz data -they are marked as orange points), which are shown in the panels from the left to the right, respectively. The legends list the reduced χ 2 values with the corresponding number of degrees of freedom. The best-fit parameters for each frequency band are listed in Table 4. The best-fit spectral indices are negative for all three frequencies, i.e. the spectrum is inverted, which is consistent with the increasing flux density for higher frequencies as is see from the three panels above.
amplitude of its effects is smaller by the ratio of the precession frequency and the binary's orbital frequency. Katz et al. (1982) developed a formalism to calculate the frequency as well as the amplitude of shortterm nodding motions in the precessing accretion flows in close binary systems. Here we apply their formalism for SMBBHs. The nutation frequency may be expressed as (see also Caproni et al. 2013), where the angular frequency of the precession ω p should be negative with respect to the orbital angular frequency ω orb . Equation (4) can then be employed to derive the observed binary orbital period: For the case of OJ 287, the observed precession period is in the range of P obs p ∼ 20 − 30 years and the nutation period is in the range of P obs n ∼ 1−1.6 years, see Table 1, which yields the observed SMBHB period in the range P obs orb 2.1 − 3.8 years. This period is about three to six times shorter than the orbital period of 12 years that is adopted in most studies of OJ 287, directly from the observations of (periodic) optical flares (see Valtonen et al. 2016, and references therein).
We note that the inferred orbital period of ∼ 2 − 4 years and the optical-flare periodicity of 12 years (e.g., Sillanpaa et al. 1988) are not mutually incompatible. It is quite plausible that the optical flares of OJ 287 occur at certain specific phases of the SMBBH orbital motion. Such a possibility has indeed been suggested to explain the 35-day periodicity of X-ray flares of the X-ray pulsar Her X-1, even though its binary has an orbital period of just 1.7 days. The flares would occur when the rapid nodding motion combines with the slower precession, effectively removing the absorbing matter from the line of sight (Katz et al. 1982). More detailed investigation of this effect is needed for OJ 287, in combination with the polarimetric monitoring. The amplitude of the nutation motion is then expected to be smaller than the precession amplitude by a factor of P obs p /P obs orb 5.3 − 14.3, applying the model of Katz et al. (1982) to OJ 287.
In Fig. 6 [a] we display the observed period of the SMBBH system OJ 287 that yields the observed precession period of 23 years, as a function of the component distance (in milliparsecs) and the mass ratio x p of the primary to the total SMBH binary mass, i.e. x p = m p /(m p + m s ). The SMBHB orbital period range of 2.1 − 3.8 years in the observed frame is marked as a shaded black rectangle in Fig. 6, while the previously discussed orbital period of 11.7 years related to optical outbursts is marked as a dotted black line. We see that for the assumed total mass of 4 × 10 8 M of the binary the primary mass fraction x p is much broader for the shorter orbital period of 2-4 years, x p ∈ (0.5, 0.80−0.88), while for 11.7 years, the range shrinks to essentially two equal-mass components, which is unlikely. This situation is the same for the larger total mass of 10 10 M , see Panels [c] and [d] in Fig. 6. In addition, for the total mass of m tot 10 10 M of the SMBH binary, whose components have the observed period of P obs orb = 10 years, the merger timescale decreases just to 10 4 years according to Peters (1964), which appears to be statistically less likely, while for the total SMBH binary mass of 4 × 10 8 M , τ merge 10 6 yr. For the above reasons, while doing estimates related to OJ 287, we have assumed a total SMBBH mass of 4×10 8 M , which is preferred, e.g., in Liu & Wu (2002), based on the SMBH-bulge mass correlation. A similar mass, 1.3 × 10 8 M , was derived by Neronov & Vovk (2011), based on the time-scales of fast γ-ray flares. For completeness, in Fig. 6[b] we also show the merger timescale as a function of a component separation and the SMBBH mass ratio. The same two orbital periods are depicted by thick grey lines as in Fig. 6[a].

PRECESSION-INDUCED (TEMPORAL)
VARIABILITY OF THE OJ 287 SED As discussed in this paper, radio variability of blazars is substantially, even predominantly, influenced by deterministic processes (i.e., precession). It seems plausible that the precession would then also influence the overall Spectral Energy Distribution (SED) of a blazar in a deterministic fashion. Turning this argument around, a predictable shift in the SED corresponding to a certain phase of the precession should strengthen the case for precession playing a prominent role in blazar variability. We examine this in this section. In essence, precession means a change of the Doppler factor due to a changing viewing angle, see Eq. (3). A timevarying Doppler factor should also be reflected in temporal evolution of the SED. Both synchrotron emission [a] [b] [c] [d] Figure 6. Panel [a] shows the observed orbital period P obs orb of the binary system of OJ 287 as a function of the component distance in milliparsecs rps and the mass ratio xp of the primary BH to the total SMBH binary mass (expressed as logarithms), i.e. xp = mp/(mp + ms). The black shaded rectangle marks the range of the putative orbital period of 2.1 − 3.8 years as derived based on the nutation period and the dotted black line represents the orbital period of 11.7 years based on the optical flares. The total BH mass was assumed to be 4 × 10 8 M . The upper white region represents the parameter space where the outer disc radius is larger than the primary-secondary BH distance, r out d > rps, and hence rigid precession cannot proceed. The lower white region represents the parameter space where the merger timescale is shorter than the orbital period of the SMBHB. Panel and Inverse-Compton emission should get boosted (deboosted) in tandem.

OJ 287 reveals an atypical SED-state
Kushwaha (2020) discusses new spectral features which appeared during the 2015-2017 multi-wavelength high-activity state of OJ 287. In particular, a break in the NIR-optical spectrum and hardening of the MeV-GeV emission accompanied by a shift in the location of the peak, were observed. Fig. 1b in Kushwaha (2020) shows a flaring SED (MJD: 57359-57363), a typical SED during the VHE phase (MJD: 57786) and the quiescent state SED after the VHE activity (MJD: 57871). According to Kushwaha (2020), the new spectral features support the disk-impact binary SMBH model over the geometrical class of models. The 2015-2017 multiwavelength flare coincides with the time of the projected lower reversal point as deduced from VLBA studies of the jet kinematics and identified as the sign of precession by Britzen et al. (2018). In the following, we explore whether the (geometrical) approach of changing just the viewing angle can reproduce the new SED-state discussed by Kushwaha (2020).

The SED variation of OJ 287 during the passage of the lower projected precession reversal point
In the following, we assume a single emission zone generating the SED of OJ 287. We fit the SED observed at two epochs close in time to see if the change is consistent with a jet component following a typical trajectory reflecting the precession + nutation model.
The fitting to the SED of OJ 287 with self-absorbed synchrotron and synchrotron self-Compton (SSC) contributions is done by employing the AGNPY package (Nigro et al. 2021), version 0.0.8 1 , written in python for numerically computing the photon spectra produced by leptonic radiative processes in radio-loud AGN. The electron energy spectrum (described by the spectral normalization factor), the low and high-energy spectral index, the minimum and maximum electron Lorentz factors, the Lorentz factor at the break energy, as well as the jet Doppler factor, inclination and the magnetic field strength were free in the fitting. The luminosity distance of OJ 287 (z = 0.306) is d L = 4.9 × 10 27 cm. The radius of the radiating blob is fixed at R b = 3 × 10 16 cm. We present the parameters of the fitted leptonic SEDs for MJD 57359-57363 and at MJD 57786 in Table 5 and plot the computed SEDs in Fig. 7. We note that these results are possible realizations of a model with a large number of parameters.
1 https://github.com/cosimoNigro/agnpy In Table 5 we compare the parameters obtained by SED modeling of the data at the end of 2015 (from 2015-12-13 to 2015-12-17, 2015.92, SED a ), in the so-called flaring state, with the values obtained for the VHE phase (SED b ). The main difference lies in the viewing angle of the jet (∼ 12.1 • at the flaring phase, ≤ 2 • at the VHE phase) and the corresponding Doppler factors (3.7 and 45.1 respectively).
It seems that the fit by AGNPY produces different values of the Doppler factor for SED a and SED b , but in addition, other important quantities describing the synchrotron and the Compton emission might have changed according to the SED fits. We explore whether variations in: (i) both the component parameters and the global value of the jet's Doppler factor, (ii) just the global Doppler factor, or (iii) just the component parameters, can best model the observed SED. We made several SED fitting-runs for SED b : • case0: the synchrotron and SSC parameters, as well as the Doppler factor are allowed to change during the fit (see the results in Table 5, column 'b'); • case1: the synchrotron and SSC parameters are allowed to change during the fit (initial parameters from SED a ), the Doppler factor is fixed to the value of SED a (see Table 5, column 'a'); • case2: the synchrotron and SSC parameters of SED b are fixed to the values of SED a , the Doppler factor is allowed to change (initial value from SED b ); • case3: the synchrotron and SSC parameters of SED b as well as the Doppler factor are fixed to the values of SED a . In this case there is no extra fitting, we will only calculate the χ 2 between SEDa model and SEDb data.
The reduced chi-squares are the following: χ 2 r = 2.06 (case0 , Table 5.), χ 2 r = 27.7 (case1), χ 2 r = 377.5 (case2), χ 2 r = 450.39 (case3). It seems that AGNPY does not find a good fit when either the electron energy parameters or the Doppler factor is fixed, such that case2 is worse compared to case0, much worse than case1 to case0. Since the observed synchrotron and SSC flux depend on δ 4 , we find that the changes in the Doppler factor have played a more decisive role between SED a and SED b , compared to the parameters describing the electron energy spectrum. We further note that the increase of the Doppler factor between SED a and SED b cannot be evaded by increasing the spectral normalization factor (k e ) of SED b , since the shape of the SED changes both along the frequency and the energy axes. Moreover, changing k e affects differently the synchrotron and SSC contributions, and the Doppler factor can only raise the SED amplitude by a factor that is independent of the frequency.
By fitting the SED of OJ 287 for two epochs close in time and under the single-zone (leptonic) approximation, we find that the change of the shape of the SED of OJ 287 is compatible with a variation of the Doppler factor of a blob following a typical trajectory in the jet reflecting the precession and nutation motions.
3.3. New OJ 287 SED state correlates with special precessional phase By fitting the jet kinematics of OJ 287, we derived a precession-only model with the precession period of P p ≈ 23 yr. In addition, we derived a precession + nutation model with a similar precession period (P p ≈ 32 yr) together with a much shorter nutation period (P n ≈ 1.6 yr). Changes in the observed apparent velocities of the identified jet components are attributed to changes in the viewing angle, hence changes in the Doppler factor. In Fig. 8 we show the position angle, viewing angle, and Doppler factor of a plasmoid in the jet stream of OJ 287, with fitted parameters presented in Table 1.
SED a (2015.92) was observed close to a local minimum of the Doppler factor, while SED b (2017.02) was observed shortly after a global maximum of the Doppler factor. Meanwhile the position angle changed about 80 degrees. This suggests that during the time elapsed between the observation of SED a and SED b , the VLBI structure of the inner jet changed drastically within a short time.
During the observation time of SED a , the viewing angle derived from the SED fitting as well as the viewing angle derived from the jet kinematics were larger than the corresponding values for SED b . Thus, the two independent methods give consistent timing.
The peak value of the Doppler factor from modeling SED b (≈ 45), i.e. the SED which corresponds to the epoch closest to the lower projected precession reversal point, is more than twice of the peak Doppler factor derived from the jet precession+nutation model (≈ 18, see Fig. 8). It is important to note that we derived an upper limit for the viewing angle of the jet for SED b (Φ < 2 • ) based on the lower limit on the Doppler factor. The viewing angle is probably closer to 0, which means extreme beaming of the light. Though there is a numerical tension between the Doppler factors from the kinematical and the SED modeling, we can infer from the analyses that the jet was pointing very close to the observer's line of sight at that time. We also note that the errors on the nutation angle (i.e. half-angle of the jet) and precession angle (i.e. half-angle of the precession cone), combined with the short nutation period (1.6 years) still allow to compare the Doppler factors. We thus infer that the viewing angle changes are able to explain the new SED state in the framework of the precession+nutation model presented in this paper.

GENERALIZED APPROACH TO MODELLING RADIO LIGHT CURVES
Given our analysis that can explain the (quasi)periodic radio outbursts by jet precession and nutation, radio light curves could in general be modelled as a superposition of periodic and non-periodic, stochastic processes. In the power-density spectrum S(f ) or periodograms, periodic processes show up as peaks at specific frequencies f per , while the non-periodic, stochastic processes are represented by a power-law function of the frequency S(f ) ∝ f γ (Timmer & Koenig 1995), with γ = −1 corresponding to the flickering or "red noise" process, γ = 0 stands for the white noise, and γ = −2 stands for the random-walk noise. In our model of the radio flux variability, based on the bulk jet precession and nutation, the power-law noise is complemented by the periodic signal that consists of a precession with a frequency ω p and a nutation with a generally higher frequency ω n = 2(ω orb − ω p ), where ω orb is the orbital frequency of the BH binary.
To demonstrate this specifically, we take the radio continuum light curves at 4.8, 8.0, and the combined 14.5+15 GHz light curve and calculate their power spectral densities (PSDs) using the fast Fourier transform algorithm. The PSDs are shown in Fig. 9. It is seen that the PSDs decrease approximately as simple powerlaw functions. To analyze this, we binned the Fourier powers into 10 equal-sized bins between the minimum and the maximum frequencies. The Fourier power per bin was calculated as an average value of the powers that fall into the bin and the uncertainty was determined as the standard deviation, see ten points in Fig. 9. The binned power density spectra were then fitted with simple power-law functions using the form log PSD = γ log f + β. The best fit power-law functions are shown as lines in Fig. 9. For all the three radio frequencies, the slope is the same within uncertainties, γ ∼ −1.7, while the intercept increases from β = 2.67 to β = 3.61 from 4.8 to 14.5 + 15 GHz, respectively. This demonstrates the presence of the stochastic process with the red-noise and the damped random-walk nature (Timmer & Koenig 1995;Kelly et al. 2009). The value of γ falls into the interval (−2, −1) inferred by Tarnopolski et al. (2020) for a limited blazar sample, Table 5. a: MJD 57359-57363, b: 57786. Γ = 10 (β = 0.995c) fixed based on the present work. . Pure leptonic SEDs for OJ 287 in different spectral states, where the SED points are adopted from Kushwaha (2020).
Left: The dashed blue curve shows the (self-absorbed) synchrotron, the dotted red curve shows the SSC contribution to the total leptonic SED between MJD 57359-57363, that is shown by purple continuous line. Right. The dashed blue curve shows the (self-absorbed) synchrotron, the dotted red curve shows the SSC contribution to the total leptonic SED at MJD 57786, that is shown by black continuous line. The main difference between the models applied in the figures is the viewing angle of the jet (left plot Φ ≈ 12.1 • , right plot Φ ≤ 2 • ) and the corresponding Doppler factor (left plot δ ≈ 3.7, right plot δ ≈ 45.1).  though in the γ-ray domain. The presence of periodic processes is not apparent in PSDs. The potential periodicities appear when one expresses PSDs considering the red-noise power-law process, such as in the periodogram of the orthogonal multi-harmonic analysis of variance (MHAOV; Schwarzenberg-Czerny 1996), which is suitable for unevenly sampled light curves. In Fig. 10, we show calculated MHAOV periodograms for the three frequencies -4.8, 8.0, and the combined 14.5+15 GHz -from the top to the bottom panel, respectively. In the periodograms, we show 1σ, 2σ, and 3σ significance levels that are inferred from fitting the generalized extreme value distribution functions to the histograms of periodogram peak distributions. The periodogram peak distributions are constructed from a few hundred bootstrap realizations. In the left panel of Fig. 10, we infer confidence levels for the whole frequency range, while in the right panel, we calculate 1σ, 2σ, and 3σ levels for different frequency bins with the logarithmic increment of log 2.5 0.4 whose central frequencies are represented by points. For 4.8 GHz light curve, we find the periodogram peaks at 0.97, 10.62, and 32.47 years that are at least at 1σ level. For 8.0 GHz light curve, we find similar periodicity peaks at 0.97, 9.56, and 30.45 years with the significance of 1σ. For the combined 14.5+15.0 GHz light curve, the periodicity peak distribution is more complex with the peaks at 1.03, 1.78, 2.59, 3.84, 6.82, and 31.43 years that are close to at least 1σ significance level, see Fig. 10. For all the three light curves, the most significant peak is at ∼ 30 − 32 years. The periodicity analysis demonstrates the stochastic nature of the radio variability for OJ 287 in conjunction with periodic, deterministic processes at different periods. The periodic process close to one year (though at only 1σ significance) may be associated with the jet nutation, while the longer periods (∼10 and ∼30 years) that have ∼ 2−3σ significance could be associated with the precession in the framework of our geometric precession-nutation model. The intermediate candidate periods at 14.5 + 15.0 GHz, i.e. 1.78, 2.59, 3.84, and 6.82 years, could be related to the orbital motion of the putative SMBH binary.
Interestingly, the same periodicities are also reproduced using the Lomb-Scargle (LS) periodogram, see Fig. 11. In the left panels of Fig 11, we show the histograms of LS periodogram maximum constructed from 2000 bootstrap realizations. From the fitted Gaussian function, we inferred 1, 2, and 3σ confidence levels for the whole frequency range. In the right panel, we show LS periodograms for 4.8, 8.0, and 14.5+15.0 GHz light curves, where we also depict 1, 2, and 3σ confidence levels calculated for different frequency bins in a similar way as for the MHAOV periodograms (using 1000 bootstrap realizations for each frequency bin). In this case, the periodocity peak significance is generally lower, with the peaks at ∼ 30 years and ∼ 1 year being above the 1σ level (see the 8.0 GHz light curve), while the other peaks that are close to and below ∼ 10 years are at or mostly below 1σ level.
The combined model light curves versus the observed light-curves at 4.8, 8.0, and 14.5+15.0 GHz are shown in Fig. 12. The model light curves are constructed as the sum of periodic components and can reproduce well the main features of the observed radio light curve, with the potential to predict future radio outbursts of OJ 287, which can be utilized for testing the model. For 4.8 and 8.0 GHz radio light curves, the main features can be reproduced by summing three signals with periods of ∼ 30, ∼ 10, and ∼ 1 year, which is roughly consistent with the longer precession period, short nutation timescale, and an additional periodic signal that could be associated with the orbital period of the binary black hole. For the higher frequency at 14.5 + 15 GHz, the periodogram is more complex, with the ∼ 30-and ∼ 1year signals present, but additional four periodic signals (6.82, 3.84, 2.59, 1.78 years) are needed to reproduce the full variability. We note that at 14.5+15 GHz, the model has also predicted correctly the radio flare around 2017.4. A smaller radio flare is expected around 2024.8-2024.9. In the following subsection, we approach the light curve modelling from a different perspective, making use of the autoregressive intergrated moving average (ARIMA) models.

Light curve modelling using ARIMA model: determining timescale of stochastic processes
We test if the three radio light curves could be modelled using the autoregressive integrated moving average model (ARIMA), which can be applied in astronomical time series analysis (Feigelson et al. 2018;Caceres et al. 2019) to describe time series using stationary stochastic processes (see e.g. Tarnopolski et al. 2020, and references therein for the application on blazar light curves). Training the ARIMA model on the observed light curve or its part, one can make a prediction for the future evolution of the system and compare with the actual observations to see whether the source behaves in a stationary stochastic way as described by ARIMA on shorter or longer timescales. Further motivation for such an analysis is given by Vaughan et al. (2016) who show that stochastically driven flux variability may mimic quasiperiodic behaviour in AGN light curves, especially when only few cycles are captured (≤ 2; see also Ait Benkhali et al., 2020 for the periodicity detection in γ-ray detected Fermi -LAT blazars).
ARIMA(p, d, q) process contains p autoregressive terms, d denotes the order of time series differencing, and q stands for the number of moving-average terms. The autoregressive part AR(p) calculates the flux density at time t as the linear combination of p weighted past values of the flux density. The moving average MA(q) part calculates the current value as the sum of the current white noise and q weighted past white-noise values . The integrated order of differencing d ensures that the studied light curve is stationary. The final ARIMA(p, d, q) model equation for the once-differenced time series, i.e. light curve transformed by f t = f t −f t−1 , can generally be expressed as, where the parameters c, α i , and θ j can be determined using the Box-Jenkins method as implemented e.g. in python statsmodels. For each radio light curve (4.8, 8.0, 14.5 + 15.0)GHz, we first performed a linear interpolation to regular time intervals separated by 0.1 years. Subsequently, with the help of the python module statsmodels, we performed the analysis to determine the order of p, d, and q parameters. We provide details about the ARIMA model fitting in Appendix Subsection B, in particular Table 13 and Fig. 20. Using (partial) autocorrelation functions as well as AIC and BIC values obtained from fitting different ARIMA models (p = [0, . . . , 15], d = [0, 1], q = [0, . . . , 15]), we determined the preferred ARIMA model with p = 2, d = 1, and q = 1 for the lower frequency light curves at 4.8 and 8.0 GHz. This class of ARIMA models has the smallest BIC, while the AIC indicates the higher p order of ≥ 5. For the 14.5+15.0 GHz light curve, the preferred model is ARIMA(1,1,2) corresponding to the smallest BIC. The Augmented Dickey-Fuller (ADF) test indicates that its ADF p-value drops below the threshold value of 0.05 after the first differenciation of all the three light curves, which supports setting d = 1 for reaching stationarity. The best-fit coefficients for the corresponding AR and MA terms are listed in Table 13 for each light curve. In Fig. 23, we show the best-fit to the light curves for the corresponding ARIMA(p, d, q) model.
We also perform additional tests using the light curves that were regularly sampled to 0.05-year time steps. For cases d = 0 and d = 1, the number of AR terms p is in the range between 1 and 3, while the number of MA terms q is in the range between 0 and 4, see Table 14 and Fig. 21. The small increase in p and q order thus merely reflects the time-step decrease by a factor of two, while the timescale of the stationary stochastic process remains in the range of 0.05 − 0.2 years.
Furthermore, we apply fractional differencing to all the three light curves using the python library fracdiff (de Prado 2018). This ensures that light curves can be described by a stationary stochastic process, while still keeping the largest possible signal. For a larger time step of 0.1 years, the term d reaches values intermediate between 0 and 1 and less than 0.5. For a smaller timestep of 0.05 years, d decreases to close to zero for 4.8 and 14.5+15 GHz light curves, while it is d = 0.41 for the 8.0 GHz light curve. The p and q values for fractionally differenced light curves are displayed in Fig. 22. For the time step of 0.1 years and the model with the smallest BIC, p and q are consistently equal to one, while for the time step of 0.05 years, p is in the range if 1 − 3 and q = 0. Hence, again the characteristic timescale for the stationary stochastic process is ∼ 0.05 − 0.15 years.
We also fit the ARIMA model to the temporal evolution of the quasi-stationary component a, which was interpolated with the time-step of 0.15 years. Using fractional differencing and the ADF test, we determined the difference degree of d = 0.555 to reach the stationarity of the temporal evolution, see Fig. 25. Based on the minimum AIC and BIC values, we found that the optimal ARIMA models have a low order of autoregressive and moving average terms with p ≤ 3 and q ≤ 2, which indicates the timescale of ∼ 0.15 − 0.45 years for stationary stochastic processes. In Fig. 26, we show the best-fit ARIMA(1,0.555,0) model (left panel) and the forecast testing (right panel).
The fitted ARIMA models and the ADF test suggest that without differencing, radio light curves are nonstationary, i.e. they exhibit a long-term trend, which is apparent by eye. The long-term jet precession and nutation are viable candidate mechanisms to address this trend as we showed in Section 2. Forecast mean ARIMA models, which were trained on approximately the half of the corresponding radio light curves, tend to converge towards a constant value, which reflects the short-term memory of the preferred ARIMA models, i.e. p and q are at most 2 for the time step of 0.1 years, hence they depend on at most two lagged flux density and white noise values going back in time by 0.2 years. The p order of ≤ 2 was also found by Tarnopolski et al. (2020) while fitting blazar γ-ray light curves using the ARIMA model. The low order implies the presence of a damping timescale of the propagating disturbances in the compact disc-jet system, beyond which they do not affect significantly the emission processes. In the left panel, we show the histogram of MHAOV periodogram peaks constructed from 50 generated boostrapped light curves with 10 best identified peaks. From the fitted Gumbel probability density function, we inferred 1σ (0.3173, red), 2σ (0.0455, green), and 3σ (0.0027, orange) confidence levels for the whole frequency interval where the smallest frequency is given by the light curve duration, while the highest frequency is set to (2 × mean sampling separation) −1 , i.e. the Nyquist frequency. In the right panel, the confidence levels depicted in the MHAOV periodogram are calculated for specific frequency bins whose central frequencies are marked by points and they exhibit a decreasing tendency towards higher frequencies.
Within 2σ confidence interval, the forecast values are consistent with the observed radio outburst around ∼ 2010-2011, see Fig. 24. However, the sum of three main deterministic periodic processes indicated by the periodograms, in contrast to ARIMA stationary stochastic processes, can reproduce the main long-term light curves trends, see Fig. 12, with a clear trend forecasting in comparison to the ARIMA model.
Overall, the preferred ARIMA models with the weighted terms lagged by < 0.2 year imply that the OJ 287 system is dominated by stochastic processes on short timescales, which is also consistent with the simple power-law shape of PSDs, see Fig. 9. In contrast, on the timescales ≥ 1 year the system is dominated by deterministic periodic processes related to precession and nutation jet motions, which is indicated by precessionnutation model fits and periodogram peaks, see Figs. 10 and 11.

Comparing precession parameters for a sample of Blazars
In this section, we extend our analysis to a larger number of systems in order to identify sources where jet precession/nutation seems to be present.
In order to explain the structural changes on parsec scales and the concurrent flux-density changes observed for individual blazars, several authors have employed bulk precession of the jet. In Tables 6 and 7 we list sources for which the entire set of precession parameters has been successfully fitted or estimated. The tabulated parameter values have been obtained from the literature and go back to studies obtained in various frequency bands. The information is based on data sets widely differing in quality and quantity, hence the uncertainties in the precession parameters also vary substantially.
As noted earlier, jet precession is intimately connected to the radio flux variability, via the Doppler factor.
Recently, jet precession for the nearby low-luminosity radio galaxy M81 was confirmed with the period of ∼ 7 years (von Fellenberg et al. 2023). However, another so-lution is also discussed, specifically a precession-nutation model with the 7-year period corresponding to nutation, while precession with the period ∼ 800 years can account for the observed long-term trend of the position angle. Hence, because of the ambiguity, this source is not included in our sample of precessing jets.
Observationally, perhaps a more meaningful parameter than the Doppler factor is the ratio of the maximum to the minimum Doppler factor, ξ ≡ δ max /δ min = (S max /S min ) 1/(p+α) , which are computed from the temporal profile of the calculated Doppler factor inferred from the observational radio VLBI data, see Eq.
(3). This ratio ξ is a more realistic measure of variability, being proportional to the ratio of the maximum to the minimum flux density. However, since the values of the Doppler factor ratio depend on the ratio between the maximum and minimum flux densities in the light curves, the observationally inferred values of ξ would strongly depend on the considered epochs of monitoring. If an outburst is included, the ratio would increase, while epochs of nearly quiescence would lead to a low ratio. We therefore included information on the considered epochs in Tables 6 and 7. Figure 13 [b] displays the values of ξ for the sources listed in Table 6.
We list the computed values of the maximum Dopplerfactor ratio ξ max (the minimum is unity that implies zero or a negligible variability due to the Doppler boosting of a precessing jet) in Table 8. For several sources listed in this table, the precession model fits the light-curves well. In Britzen et al. (2019a), the correlation between the radio flux-density evolution and jet precession is shown. For 3C 84, the radio flux-density is modeled and nicely fits the OVRO and UMRAO data (15 GHz). A superposition of the observed optical magnitude in the B band and the Doppler factor of the precession model is shown for 3C 279 (Abraham & Carrara 1998). For 3C 345 and 3C 120 the authors show the contribution of the precessing jet to the flux-density in the B-band (3C 345: Caproni & Abraham, 2004a;3C 120: Caproni & Abraham, 2004b). For PG 1553+113, the Doppler boosting factor is superposed to the Fermi γ-ray light curve (see In the left panels, we show the distributions of LS periodogram peaks inferred from 2000 bootstrap realizations of the radio light curves. The vertical lines mark 1, 2, and 3σ confidence intervals calculated for the whole frequency range. In the right panels, we show LS periodograms with marked 1, 2, and 3σ confidence levels that were inferred for frequency bins whose central points are represented by points. The dotted horizontal gray line marks 5% false alarm probability level. [a] [b] Figure 13. [a] shows the time-dependent Doppler-boosting factor δ(t, γ, Φ) for the radio galaxies listed in Table 6. In [b] the ratio of maximum to minimum Doppler factor ξ = δmax/δmin derived for one precession phase is displayed. The dot-dashed line denotes a ratio ξ = 2; most of the sources lie below this limit (8/12), while four lie above.  Table 9. The values for 3C84, 3C 120, TXS 0506+056, PKS 0735+178, OJ 287, PKS 1502+106, PG 1553+113 and 2200+420 calculated based on the precession model fall in a similar range as those obtained by Hovatta et al. (2008) and Homan et al. (2021). The largest discrepancy Table 6. A selection of AGN with variability and/or jet morphological changes that can be modelled applying a jet-precession model. References are: 1) Abraham & Carrara (1998); 2) Abraham & Romero (1999) 1964-1991 1960-1990 1995-2009 1980.93-2007.682 2009-2016 1977-1995 1976-1999 1995-2017 1995- Homan et al. (2021), the apparent speeds of components and the brightness temperature based on these speeds are used to derive a Doppler factor. To our knowledge, the MOJAVE team assumes the shock-in-jet scenario and does not account for any geometric effect (e.g. precession) which can lead to different apparent speeds. In the case of the values obtained by Hovatta et al. (2008), the brightness temperatures of the fastest flares are used to determine the Doppler factor. The fastest flare might be due to a shock, an instability or magnetic reconnection and does thus not necessarily serve as a indicator of the boosting. Both approaches (apparent speed, fastest flare) thus rely on model assumptions.
Both teams provide one single value for the Doppler factor. The Doppler factor based on the observed preces-sion takes decades of measurements into account and can predict the evolution of the Doppler factor. This predicted Doppler factor evolution can be confronted with future measurements.

DISCUSSION
Many observational studies of AGN radio variability have been reported in the literature (e.g., Hughes et al. 1992;Aller et al. 1999;Fuhrmann et al. 2016).
The kinematics of AGN jets on the parsec scale of AGN is usually probed using VLBI observations. Brighter parts in the jet -so-called components -can be traced and are seen to move at apparent speeds typically of the order of 5-10 times the speed of light (see e.g. CJF -Caltech-Jodrell Bank Flat-Spectrum sample: Britzen et al., 2008;for MOJAVE: Lister et al., 2009). The physical nature of these apparently moving jet components is still unclear, in particular the question whether the apparent motion merely reflects a pattern speed, or a real material object like a blob of synchrotron plasma is still open. The most often commonly invoked scenario to explain the component motion in jets is that of a shock travelling within the jet stream -the so-called shock-in-jet model. This directly relates to the most common explanation for variability as outlined in the Introduction.
Despite the prominence of the shock-in-jet scenario in the literature, several studies have indicated that the jet kinematics may in fact be largely dictated by geometrical processes. In particular, jet precession has been invoked to explain the observed phenomena (e.g., Abraham & Romero 1999; Kudryavtseva et al. 2011;Abraham 2018).
Precession as a potential cause of AGN flux variability in AGN has been discussed by a number of authors in the past (see e.g., Caproni et al. 2013;Abraham 2018).
Precession is not a property exclusively associated with blazars only; Seyfert 1 galaxies as well as radio galaxies exhibit signs for precession -3C 84 being a prominent example (Britzen et al. 2019a).

Origin of jet precession and nutation
Precession and nutation of a jet basically reflects the motion of its base. Depending on the dominant mechanism for jet launching, the base can be close to the central BH (Blandford-Znajek jet, Blandford & Znajek 1977) or to the surrounding accretion disk (Blandford-Payne jet, Blandford & Payne 1982). Essentially, both types of outflow can be set into precession, so long as requisite non-axisymmetric forces are present. In general terms, this requires a multi-component system with misaligned rotation axes. These components could be a Table 9. The range of Doppler factors (δmin, δmax) derived via the precession model (this paper). In addition, we list the mean Doppler factor calculated as δ = 0.5(δmin + δmax). We list Doppler factors derived from the literature (determined without assuming precession) for the AGN compared in this paper. No Doppler factor has been determined for "-" entries. For the case of OJ 287 (for details see Britzen et al. 2018), all three scenarios have been considered for the periodic motion inferred for the disk-jet system; however, considering timescales of the order of 10 years, option (b) could be excluded, leaving for the probable cause of precession an efficient disk-jet coupling. We further refer to, e.g. Caproni & Abraham (2004b) for the application of the SMBBH model to the jet kinematics, and to Caproni et al. (2004) for the observational analysis of the presence of LT precession in several AGN. Both processes, (a) and (b), can be expected to occur in the aftermath of galaxy mergers, and are thus a natural outcome of hierarchical galaxy formation scenario. The relevance of SMBBH formation has been highlighted in numerous studies (e.g., Begelman et al. 1980;Wilson & Colbert 1995;Mayer 2017). One expects a SMBBH to form during the final merger phase (e.g. Mayer 2017, and references therein), or a misalignment of disk and BH spin axes as a direct result of two supermassive BHs merging. Process (a) and (b) requires two massive BHs with sub-parsec separation, thus in the final merger stage. Process (c) arises due to the framedragging in the vicinity of a spinning BH and requires the offset of the accretion flow angular momentum vector from the BH spin. Such an offset can arise in the case of either: (i) post-merger systems due to their random orientations of the original orbital plane of the merging BHs with respect to the accretion flow around a (primary) BH, (ii) accretion flows that are fed by winds from a disklike stellar system as expected in common models for the hot accretion flow for Sgr A* (Cuadra et al. 2006;Shcherbakov & Baganoff 2010;Yalinewich et al. 2018). For this example the misalignment between the accretion disk axis and the BH spin is specifically discussed in Dexter & Fragile (2013).
Note that in scenario (c) the BH-driven jet is launched without precession as the BH axis remains stable in space, while it is the accretion disk which is precessing. However, disk winds and jets launched from the disk will follow the disk precession and thus set their jet axis in precession. As a result, the environment of the central spine jet will change in time, and also the spine jet will structurally follow the precession of the disk jet and disk (see e.g. McKinney et al. 2013). Thus, there is an overall consistency between the precessing jet model for AGN, presented in our work, and the large-scale picture of galaxy formation.

Numerical simulations of jet precession
Mounting evidence for jet precession has been found for an increasing number of AGN, as discussed above. Additional support to this scenario comes from observations of micro-quasars, i.e. stellar binary systems in which a compact object accretes matter from a donor star and a jet is launched (Mirabel & Rodríguez 1999).
As a typical example, the GRAVITY collaboration recently confirmed evidence for precession in the microquasar SS433 (Gravity Collaboration et al. 2017), as already reported in several studies (see e.g. Blundell et al. 2011). SS433 is known to eject a two-sided radio jet precessing with a period of 162.5 days (Margon 1984). Recent hydrodynamical simulations are capable of reproducing the observed precession signatures of SS433 very well (see e.g. Monceau-Baroux et al. 2014, 2017. Another such example is the micro-quasar GX 339-4, an X-ray detected BH binary. It also shows variability in the IR and X bands, which could be modeled via jet precession driven by the LT mechanism (Malzac et al. 2018). Yet another case of jet precession is the microquasar V404 Cygni (Miller-Jones et al. 2019). Also the jet precession in LS I +61 • 303 (Wu et al. 2018, and references therein) is likely the result of LT precession (Massi & Zimmermann 2010).
Since 3-dimensional (3D) general relativistic magnetohydrodynamic simulations are very CPU expensive, only recently the evolving precession (that is essentially 3D) of such systems could be numerically investigated (see reviews by Hawley et al. 2015;Martí 2019). Recent progress has been reported by Liska et al. (2018) who investigated the dynamical evolution of a tilted thick accretion disk around a rapidly spinning BH. In fact, the authors find that the disk-jet system as a whole undergoes a LT precession. These simulations may serve as clear evidence that the jets can be used as probes of disk precession. The precession time scales obtained within their simulations are consistent with the timescales that were derived for precession and nutation in OJ 287 (Britzen et al. 2018). However, we caution that despite the clear precession seen in these simula-tions, the application of these simulations to the observationally inferred precession of the AGN central engine is somewhat limited. The reason is that the precession amplitude found in these simulations decreases on longer time scales (after several precession periods) and the system gradually realigns. Furthermore, Liska et al. (2019) have shown by highly resolved GR-MHD simulations of the BH-disk-jet system that any misalignments may disappear due to the Bardeen-Petterson effect. Nonetheless, the onset of disk warping and subsequent wind launching directed along the disk rotation axis seems confirmed (White et al. 2019).
The evolution of jet launching in binary systems has also been investigated recently (Sheikhnezami & Fendt 2015, 2018, albeit not in the relativistic regime. The authors investigate using non-relativistic 3D MHD simulations the evolution of jets launched from the disk around a primary compact object. The disk starts warping and with that the direction of the disk jet changes. The jet axis starts tracing a circle, indicating strongly the formation of a precession cone. In particular, it could be shown how the spiral structure of the accretion disk subsequently affects the mass loading and the magnetic field of the disk outflow, namely leading to a spiral magnetohydrodynamic structure of the jet (Sheikhnezami & Fendt 2022). Also here, CPU constraints as well as the physical constraints like the mass loading of the disk limit simulations time scales to less than an orbital time scale of the secondary (however 1000s of disk orbital periods at the jet foot-point).

Implications of jet precession for multi-wavelength flaring
In this paper we focus on radio variability of AGN which seems to arise from jet precession. However, several authors have presented evidence that jet precession may also be responsible for the variability observed in other bands of the electromagnetic spectrum. A deeper examination of the multi-wavelength aspect of jet precession is planned for a forthcoming paper, but some examples will be mentioned here. Since precession is a geometric phenomenon, it would affect the synchrotron flux density at a wide range of frequencies. For OJ 287 as a representative of BL Lac sources, the optical emission is dominated by the synchrotron radiation of the jet and its observed optical periodicity is expected to be at least partly related to the jet precession (e.g., Abraham 2000;Britzen et al. 2018). In the case of OJ 287 we found a correspondence between the radio periodicity and the periodicity seen in the optical (Britzen et al. 2018). The time interval of the optical flaring is roughly half of the radio flaring. The authors suggest that most likely both periodic phenomena originate from jet precession.
The precession phase in OJ 287 relates to the SED variability, as discussed in the following subsection. Stamerra et al. (2016) use multi-wavelength observations of PG1553+113 to investigate if the observed modulation (gamma-rays, optical) can be explained with geometrical variations in the jet, possibly pointing to jet precession. They discuss PG1553+113 as a probe for a geometrical periodic modulation. Independent evidence for jet precession being responsible for the periodic multi-wavelength flaring in PG1553+113 has been presented by Caproni et al. (2017). Tavani et al. (2018) also discuss the possible gamma-ray (from Fermi -LAT observations) signatures of a SMBBH in PG 1553+113. In the case of 3C 84, Britzen et al. (2019a) find evidence for precession of the central radio structure and a correlation between the radio and gamma-ray light-curves with the higher energy emission lagging the low-energy emission. This delayed correlation is difficult to explain in terms of the shock-in-jet scenario. A plausible scenario, however, is that the correlation arises from jet precession and the time delay due to the emissions originating in different parts of the radio structure. Bhatta et al. (2020) studied the deterministic aspect of the γray variability of 20 γ-ray bright blazars dependent on the source class. They find that the dominant physical processes in FSRQs are more of deterministic nature. The following four sources from our paper are in common with Bhatta et al. (2020): 3C 273, 3C 279, PKS 1502+106, and 2200+420. In PKS 1502+106, they find the strongest deterministic case.

Precession phase related to time variable SED
Since precession results in a time-variable Doppler factor, this should also influence the temporal evoluton of the SED, affecting the synchrotron as well as the Inverse-Compton emission in tandem. Indeed, we show in this paper that the new spectral features seen in the 2015-2017 multi-wavelength high activity state of OJ 287 (Kushwaha 2020) could correspond to a special precession phase. In this paper, we demonstrate that the Doppler factors derived from SED-fitting and precession-fitting peak at the same time. We conclude that the viewing angle changes caused by precession also induce variations of the SED state. Britzen et al. (2019b) have reported evidence that the neutrino emission observed from TXS 0506+056 might be related to a collision of jetted material. In their model, a special viewing angle of the precessing inner jet plays a major role in obtaining the right circumstances for neutrino generation. The timescale inferred from the observed precession signatures is of the order of 10 yr. The observational clue presented in that work appears promising for addressing the issue of neutrino production, intricacy of which is highlighted, e.g., by Rodrigues et al. (2019). Most recently, further neutrino emission has been observed from the direction of TXS 0506+056 on April 18, 2021 by Baikal-GVD (Belolaptikov & Dzhilkibaev 2021). On September 18, 2022, the IceCube Collaboration indicated that another neutrino arrived from the direction of TXS 0506+056 (Becker Tjus et al. 2022). Both events support the prediction for repeated neutrino-emission based on the precession model for this source and a supermassive binary black hole model (Britzen et al. 2019b;de Bruijn et al. 2020). For another IceCube candidate emitter, the AGN PKS 1502+106, Britzen et al. (2021) find evidence for precession (and nutation).

Precessing jets as potential neutrino emitters
Four independent neutrino detectors -including Ice-Cube -reported neutrino emission in December 2021 (see e.g., Belolaptikov & Dzhilkibaev 2021;Sahakyan et al. 2022). PKS 0735+178 is located just outside the 90% point-spread function containment of the Ice-Cube event, within the expected systematic uncertainty of the determination of the arrival direction. For PKS 0735+178, evidence for precession has been presented in Britzen et al. (2010a).

Lighthouse effect versus Precession
In case of well observed and sufficiently long monitored VLBI jets there is a simple phenomenological way to discriminate between geometric precession and internal helicity. In case of geometric precession the pitch of the jet helix (which is the height of one complete helix turn, measured parallel to the axis of the helix) is nearly constant. The helix we see is only a pattern, the components move essentially along straight or ballistic trajectories. Whether the cause of the periodic jet launching direction is the spin-orbit precession, Newtonian binary motion (e.g. Kun et al. 2014), LT precession (e.g. Britzen et al. 2019b), etc., the result is the same, a namely constant pitch.
In case of internal helicity the pitch is not constant, and the jet components move along the helix (e.g. Steffen et al. 1995;Hardee 1987). This is the lighthouse model (see below), and the Kelvin-Helmholtz (or any magnetic field/instability related) model. While Kelvin-Helmholtz instabilities can be expected from almost any hydrodynamic shear flow, it remains unclear whether the flow pattern can actually develop into a large-scale (meaning pc-long) structure. Numerical simulations show that these instabilities develop rather quickly, resulting more in a turbulent cocoon between the jet beam and the ambient gas than in a wave-like, bent jet morphology (see e.g. Martí et al. 1997).
If one is able to measure and de-project the pitch, one can discriminate between external (changing jet launching direction) or internal helical motion (the jetcomponents follow MHD-determined paths).
As an alternative explanation for the observed (intraday) variability of blazars, Camenzind & Krockenberger (1992) proposed the so-called lighthouse effect. In this model, radiating plasma blobs are injected along the magnetized relativistic jet and then follow the overall jet kinematics. The radiation from the relativistically moving blobs is beamed in the direction of motion and thus when the blob follows a helical trajectory, its light beam sweeps across the observer repeatedly, giving rise to a quasi-periodic emission flares. Due to the overall acceleration and collimation of the bulk outflow, the geometry of the helical trajectory probably changes along the jet. We note that the lighthouse model is a kinematic model which neglects the dynamics of the jet flow. In fact, the nature and the dynamics of the postulated jet blobs is not explained. Naturally, any MHD jet is rotating, however, the poloidal motion dominates rotation for super-Alfvénic flows. Also, due to the rotation, MHD jets naturally imply the existence of a jet helical magnetic field. This is different from jet precession where the motion of the ejected blobs can be purely poloidal -being ejected ballistically into different directions over time. The nozzle itself is rotating (i.e. precessing), but the jet velocity is expected to exceed the precession speed substantially, so in principle even a jet with purely poloidal speed and magnetic field will change its viewing angle. The latter does not happen in the light-house model.
From an observational point of view, the periodicity time scale and the number of events detected in a given time may vary from source to source and also from time to time for a given source. Britzen et al. (2000) showed that the simultaneous detection of radio/optical/gamma-ray flaring in the blazar PKS 0420-014 is in agreement with the predictions of the lighthouse model. Similarly for 3C 454.3, Qian et al. (2007) found periodic variations in the radio flux that are consistent with the lighthouse scenario. 5.7. OJ 287: Unsolved questions regarding the "Impact-Model" There seems to be a need to explore an alternative to the currently popular model of the double-peaked optical outbursts appearing every 12 years (Sillanpaa et al. 1988;Valtonen et al. 2009) which invokes periodic piercing of the orbiting secondary BH through the accretion disk of the primary BH. The disk-piercing model has scored well in predicting the times of the giant optical flares (with ∼ 12 year periodicity), by suitably updating of the input dynamical parameters. However, Komossa et al. (2023) did not detect the outburst predicted for October 2022 by the "Impact-model". Neither was the predicted thermal bremsstrahlung spectrum observed (at any epoch). The Impact-Model also predicts a very large mass of the primary black hole (∼ 10 10 M ). Komossa et al. (2023) estimate the mass of the primary black hole in OJ 287 to be ∼ 10 8 M and thus confirm the mass value proposed by Britzen et al. (2018), based on the precession model.
The Impact-Model scenario is a hybrid in the sense that it identifies the periodic giant optical flares to thermal radiation arising from the disturbed/ejected gas in the disk as the latter is impacted by the accretion disk of the primary BH by the orbiting secondary BH. An expected corollary of this is that the optical flaring event would be followed by an enhancement of jet activity (i.e., synchrotron radiation) on account of gradual accretion of the disturbed disk gas on to the primary BH. Thus, conceivably, jet precession and resulting variability of its synchrotron and IC flux could take place even in this scenario as well. However, as the immediate major outcome of the impact (identified with the giant optical flare), the model invokes free-free radiation from the impact-heated gas belonging to the accretion disk. One would then also expect to see bright optical/UV emission lines from the same cloud of thermal plasma, although the required cooling of the gas cloud to an appropriate temperature might delay the onset of bright recombination line emission. Such emission features, to our knowledge, have never been reported for this source and therefore it has steadfastly retained its BL Lac classification. Oddly, this has happened despite its high brightness (m V ∼ 12 mag) as well as the fact that the free-free cooling rate is very rapid in the temperature regime of 10 5 -10 6 K. An intensive optical spectral monitoring would play a critical role in placing this model on a firm footing (or, otherwise). A few observational pointers relevant to this issue are: • We are NOT dealing with a thin accretion disk (BL Lac). Hence the putative "impact" is a nebulous concept (see Villforth et al. 2010).
• For all the flares, the primary peak lacks a radio counterpart, while some secondary peaks have shown radio counterparts.
• Optical polarisation can be very high even for a primary peak (Villforth et al. 2010). This is a problem for the thermal interpretation (Valtaoja et al. 2000;Villforth et al. 2010).
• The amplitude of the optical flares has been steadily declining from flare to flare. This holds for both primary and secondary peaks in a flare.
At present, the key arguments in support of the thermal optical interpretation for the giant optical flares comes from the observed drop in the fractional linear polarisation during the giant optical flares, a possible indicator of enhanced thermal contamination of the optical synchrotron radiation (Valtonen et al. 2017;Kushwaha et al. 2018). While such a contamination may well be occurring, we would like to recall that the drop in a fractional polarisation with a rise in brightness may also occur in case of a mildly swinging synchrotron jet (Gopal-Krishna & Wiita 1992). Thus, we believe that a clinching evidence in favour of the thermal interpretation of the giant optical flare awaits an unambiguous detection of optical/UV emission lines following the flares. Given the high brightness of this blazar, even a 1-metre class optical telescope should suffice.
Recently, Suková et al. (2021) performed GRMHD simulations of an orbiting perturbed (secondary black hole) punching through the hot flow surrounding the primary SMBH. Thus, they correctly accounted for the ADAF typical for BL Lac sources. The model can address the 12-year optical periodicity by the enhanced accretion that occurs every time the secondary component passes through the accretion flow (twice per orbit), while the 24-year periodicity in the radio light curves could be attributed to the ejected plasmoids due to the viewing angle close to the jet axis (the counter-jet and its components are deboosted). However, the model does not provide a prediction for the light curves in different bands yet. In addition, the duty cycle of the modelled flares appears rather short in their model, while the large-amplitude radio flares seen in OJ287 have a long duty cycle of ∼ 10 years more consistent with the smooth viewing angle changes due to jet precession.
Finally, we may recall that in the literature, there are claims that the giant optical flares are accompanied by a Seyfert-like behaviour, e.g., a soft X-ray excess (Pal et al. 2020). However, the central engine in these objects is known to accrete at a high Eddington rate (e.g. Tortosa et al. (2022), unlike BL Lacs. Moreover, in order to justify the analogy, one expects to observe optical emission lines which are ubiquitous in case of Seyferts. A more detailed analysis of these issues is planned for a future publication.

Some observational challenges to the shock-in-jet scenario
While a coincidence is often seen between flaring at millimetre wavelengths and the injection of a radio knot into the base of the parsec-scale VLBI jet (see e.g. Savolainen et al. 2002), this is not always the case. As an example, Bach et al. (2006) studied the connection using the flares in the radio and optical light-curves of BL Lac, and found that only some of the new VLBI components had an associated event in the light curves. Several AGN, where no component ejections could be detected in VLBI observations coinciding or following radio outbursts have been reported (e.g., PKS 0735+178, Britzen et al., 2010a;3C 454.3 Britzen et al., 2013).
Stationary radio components, which maintain a constant separation from the core, have been observed in the VLBI images of 1803+784 (Britzen et al. (2005)) and several other AGN. For Gamma-ray bright blazars, an exhaustive list of stationary components (which do not move along the jet) can be found in Jorstad et al. (2017).
Such stationary radio knots have been identified with "recollimation shocks" envisioned to form within the relativistic jet flows (e.g., Stawarz et al. 2006). However, (Britzen et al. 2010b) shed new light on them, by demonstrating that while the "stationary" knots do not appear to move along the jet, they do move (rather rapidly) orthogonal to the main jet ridge line in 1803+784. For OJ 287 this kind of motion is shown to be consistent with jet nutation (Britzen et al. 2018). Throughout this paper we call these components quasi-stationary components. They need to be distinguished from those stationary components, where no motion at all is observed over a longer time span. This is the case for, e.g., Mrk 501 (Edwards & Piner 2002;Britzen et al. 2023).
AGN feedback models struggle to steadily quench star formation. Kinetic jets are less efficient in quenching star formation in massive galaxies unless they have wide opening or precession angles (Su et al. 2021). According to this study, AGN feedback models require a jet with an optimal energy flux and a sufficiently wide jet cocoon with long enough cooling time at the cooling radius.
In contrast to straight jets, precessing jets sweep through a larger volume of the interstellar medium (ISM). Therefore they have the potential for enhanced feedback into the ISM as well as the intergalactic medium, i.e. quenching or sometimes triggering star for-mation (e.g., Aalto et al. 2016, for a molecular jet) or providing radio-mechanical heating feedback that prevents hot X-ray-emitting atmospheres in galaxy clusters from runaway cooling and thus stabilizes them (see e.g. Zajaček et al. 2022).
On a different scale, in Henize 2-10, a dwarf starburst galaxy with a potential central massive black hole, a precessing bipolar outflow connecting the region of the black hole with a star-forming region has been reported (Schutte & Reines 2022). 5.9. Why the shock-in-jet model is not sufficient to explain variability The past decades have seen a predominant interpretation of variability and jet morphology in terms of the shock-in-jet scenario. While the shock-in-jet scenario might explain radio knots, this scenario is insufficient to explain the observed flux-density and structural changes of many AGN jets. The dominance of the shock-in-jet scenario over the past decades is no argument against exploring mechanisms that might unravel more basic details about the way the jet and the central engine of AGN work. As noted above and in other studies, the precession model can explain certain observational aspects of the radio jets and of AGN variability that are difficult to comprehend relying solely within the framework of the shock-in-jet scenario. We think that precession therefore modulates the appearance of the jet. Identifying deterministic patterns within variability data will enable us to discriminate between geometric (precession) phenomena and stochastic (e.g., flaring) events. This will provide a deeper and more fundamental understanding of AGN physics. It will also allow us to predict times of distinct states, e.g., flaring states due to predictable precession events (reversal points, see Fig. 4[b] in Britzen et al., 2018). This is progress compared to multiple multi-wavelength campaigns which, in the past, had to be planned and conducted "blindly" due to the lack of knowledge of predictable, correlated and inter-connected processes. Applying this previous knowledge to predict the deterministic behaviour will allow us to much more carefully explore also those processes which are not deterministic but of stochastic nature. Magnetic reconnection might play an important role in the launching of jets as well as the injection of components (Vourellis et al. 2019;Ripperda et al. 2020;Nathanail et al. 2020b).

CONCLUSIONS
In this paper we have revisited the key problem of flux variability/flaring of AGN at radio frequencies. Conventionally, it is interpreted in terms of the "shock-in-jet" scenario and its variants. Recently this paradigm has encountered problems in explaining certain instances of radio knot ejection unaccompanied by a flare and vice versa. A similar difficulty is posed by the correlated flux variability of 3C 84 at gamma-ray and radio frequencies, albeit the gamma-ray profile exhibiting a very large time lag of 300-400 days.
We have therefore proposed an alternative (deterministic/geometric) interpretation according to which the observed structure and flux density of the parsec scale radio jet predominantly reflects the precession (and nutation) of the relativistic jet, leading to a time variable Doppler beaming of its emission. We note that the detection of the signatures of nutation (a second-order precession effect) in the light curves of certain blazars provides further support to this interpretation. Hence, the monitoring of radio flux and parsec-scale structure of blazars is likely to facilitate extraction of more information about the jet ejection process as compared to its subsequent evolution via shocks, etc. The recent confirmation of repeated neutrino emission from TXS 0506+056 fosters the jet precession model. We also find that the jet precession scenario is in accord with the concept of hierarchical evolution of the universe wherein mergers of galaxies (along with their central BH) are supposed to be a key ingredient. As a caveat, the applicability of our model in connecting the jet kinematics to features in the light-curve is unclear for M87 (e.g. Britzen et al. 2017a) and Sgr A (Gravity Collaboration et al. 2018). For M87, while indication for jet precession has been found (Britzen et al. 2017a) and a tilted disk/jet system has been considered (Chatterjee et al. 2020), its link to radio flux variability needs to be explored further. In the case of Sgr A , the presence of a jet is still to be demonstrated (e.g., Falcke 1999;Li et al. 2013;Yusef-Zadeh et al. 2020;Fragione & Loeb 2020;Cecil et al. 2021). For Sgr A*, asymmetric X-ray flares seem to be associated with nearly edge-on orbits (Mossoux et al. 2015;Eckart et al. 2017;Karssen et al. 2017), which could be the sign of the accretion-flow precession or the disk warping due to the Bardeen-Petterson effect. Dexter & Fragile (2013) modelled tilted thick accretion flows around Sgr A* and obtained a good agreement with the near-infrared and the mm-domain observations in terms of the variability and the spectral evolution.
In the following we recapulate our main arguments and the results obtained.
(1) Precession of AGN jets can generally be triggered in a binary black hole system or by the Lense-Thirring effect on the accretion disk. Simultaneously, both ef-fects can also lead to nutation of the disk-jet system. In this paper, we have applied a kinematic precession model that is of general nature, and can, thus, be applied independently of the origin of the precession. We have provided a concise mathematical description of our modeling of jet precession and nutation, extending the approach of Britzen et al. (2018).
(2) We apply the kinematic precession model to twelve AGN. We model their radio light curves and determine the long-term evolution of the Doppler factor. From our comparison of the derived Doppler factors with those derived by the well-known shock-in-jet scenario, we conclude that precession is likely to play a dominant role in causing the (radio) variability, and the apparent speeds of the jet components.
(3) A corrollary of the above is that a time-varying Doppler factor should also be reflected in the temporal evolution of the SED. Indeed, we show that the Doppler factor derived from SED-fitting and the one derived independently from precession-fitting (based on the derived VLBI jet kinematics) peak at the same time taking OJ 287 as an example. We thus propose that the change of viewing angle over time, as envisioned in our precession model, can explain the variation in the SED state.
(4) The point that needs to be emphasized is that precession and nutation of a jet basically reflect the varying orientation of its base and at the same time, they can be detected by studying the pc-scale motion of the VLBI knots. Signatures of precession are also frequently observed on large-scales in Karl G. Jansky Very Large Array (VLA) and MERLIN radio maps of powerful jet sources (e.g., Krause et al. 2018). For a complete sample of 33 3CR radio sources the authors find strong evidence for jet precession in 24 cases (73%). We argue in case of a collimated disk wind or jet, the torques needed for causing the precession and nutation can be induced on the accretion disk by a secondary BH, or by Lense-Thirring precession by an inclined BH axis. In case of a precessing spine jet launched by the Blandford-Znajek process, torques on the jet launching primary BH may be induced by a secondary BH. State-of-the art numerical simulations confirm our scenario.
(5) Fitting radio light curves using the ARIMA model implies that the OJ 287 system is dominated by stochastic processes on rather short timescales ( 0.2 years), while the clearly present long-term trends are dominated by deterministic periodic processes related to precession and nutation jet motions on timescales ≥ 1 year.
(6) While we mainly focus on radio variability in this paper, we briefly discuss other bands of the electromagnetic spectrum as well, finding evidence for variability caused by jet precession at these higher frequencies as well. We further argue that jet precession and the subsequent strong dynamical changes in the jet kinematics may eventually have also led to neutrino emission in three AGN. The three known radio sources (TXS 0506+056, PKS 1502+106, PKS 0735+178) have been identified by IceCube as neutrino emitters.
(7) Earlier, we had pointed out problems with the popular model of the quasi-periodically flaring blazar OJ 287, which invoked a secondary black hole repeatedly impacting the accretion disk surrounding the primary black hole (see Britzen et al. 2018). Here we have noted additional difficulties with that scenario, accentuated due to the lack of evidence for optical/UV recombination line emission associated with the giant optical flares which are purported to arise from free-free emission (e.g., Valtonen et al. 2009).
(8) We hope that the various observational evidences for the AGN jet precession, highlighted and modelled in this work, would encourage long-term VLBI and SED monitoring of bright blazars.
The authors appreciate the significant and important help provided by the anonymous referee. The insightful suggestions enhanced the presentation of the manuscript's results considerably. The authors thank N. MacDonald for carefully reviewing the manuscript and for providing many valuable comments and thoughts that significantly improved the paper. We thank A. Witzel, S.-J. Qian, J.-P. Breuer, C. Raiteri, M. Villata, Z. Abraham, and M. Krause for very helpful discussions. Special thanks go to P. Kushwaha for providing data and a plot modified for this manuscript. M.Z. acknowledges the financial support by the GAČR EX-PRO grant No. 21-13491X "Exploring the Hot Universe and Understanding Cosmic Feedback". G-K. acknowledges a Senior Scientist fellowship from the Indian National Science Academy. E.K. thanks the Hungarian Academy of Sciences for its Premium Postdoctoral Scholarship. This research has made use of data from the OVRO 40-m monitoring program (Richards et al. 2011)  . This research has made use of data from the University of Michigan Radio Astronomy Observatory which has been supported by the University of Michigan and by a series of grants from the National Science Foundation, most recently AST-0607523. This research has made use of data from the MOJAVE database, which is maintained by the MOJAVE team (Lister et al. 2018). In Fig. 14, we show the corner plot with marginalized one-dimensional likelihood distributions and two-dimensional likelihood contours for the precession parameters t 0 , P p , Ω p , φ 0 , and η 0 inferred from the temporal evolution of the ejected component position angles η. We also display the posterior distribution for the variance underestimation factor log f , which is included in the definition of the likelihood function. Hence we deal with 6 free parameters. In Fig. 15, we compare the posterior median model with the VLBI data, including 1000 random samples of parameter chains. This way the precession model uncertainties around the model median are apparent. We summarize the non-zero flat prior parameter ranges in Table A.1.

A.2. Precession parameters based on the apparent-velocity evolution
In Fig. 16, we show the corner plot with marginalized one-dimensional likelihood distributions and two-dimensional likelihood contours for the precession parameters t 0 , P p , γ, Ω p , φ 0 , and η 0 inferred from the temporal evolution of the ejected component apparent velocities β app . We also display the posterior distribution for the variance underestimation factor log f , which is included in the definition of the likelihood function. Hence we deal with 7 free parameters. In Fig. 17, we compare the posterior median model with the VLBI data, including 1000 random samples of parameter chains. This way the precession model uncertainties around the model median are apparent. We summarize the non-zero flat prior parameter ranges in Table A Figure 14. Marginalized one-dimensional likelihood distributions and two-dimensional likelihood contours at 1σ, 2σ, and 3σ confidence levels for the parameters related to the position-angle temporal evolution due to the bulk precession of the jet. The dashed vertical lines in diagonal panels stand for 16%, 50%, and 84% percentiles of the inferred posterior parameter likelihood distributions.

A.3. Precession-nutation parameters based on the stationary-component evolution
In Fig. 18, we show the corner plot with marginalized one-dimensional likelihood distributions and two-dimensional likelihood contours for the precession-nutation parameters t 0 , P p , P n , Ω p , Ω n , φ 0 , and η 0 inferred from the temporal evolution of the position angle of the stationary component a. We also display the posterior distribution for the variance underestimation factor log f , which is included in the definition of the likelihood function. Hence we deal with 8 free parameters. In Fig. 19  random samples of parameter chains. This way the precession-nutation model uncertainties around the model median are apparent. We summarize the non-zero flat prior parameter ranges in Table A.3. Here we provide additional plots and supporting statistical verification related to the radio light curve modelling using stationary stochastic processes represented by ARIMA(p, d, q) models described in Subsection 4.1.
We apply the python function ARIMA in statsmodels that requires a regular binning of the light curves. We interpolate the data into two regular time grids -one with the time step of 0.1 years and the other with the time step of 0.05 years. The mean time sampling of the light curves is ∼ 0.03 years, hence we do not apply a smaller time step than 0.05 years. Figure 16. Marginalized one-dimensional likelihood distributions and two-dimensional likelihood contours at 1σ, 2σ, and 3σ confidence levels for the parameters related to the apparent-velocity temporal evolution due to the bulk precession of the jet. The dashed vertical lines in diagonal panels stand for 16%, 50%, and 84% percentiles of the inferred posterior parameter likelihood distributions.
In Table 13, we describe basic results of fitting different ARIMA(p, d, q) models to light curves at 4.8, 8.0, and 14.5+15 GHz that were interpolated to regular 0.1-year time steps. Here f t denotes the light curve that is differenced once, i.e. f t = f t − f t−1 . The table lists Augmented Dickey-Fuller (ADF) test statistic values and the corresponding p-values for no differencing and the first-order differencing. Then we include the preferred ARIMA model for the first-order differencing including its Akaike and Bayesian information criterium (AIC and BIC) values as well as its best-fit parameters. For completeness, we also list minimum AIC and BIC values for d = 0 and d = 1 cases as well as corresponding ARIMA models, which we inferred from AIC and BIC values as functions of different p and q parameters, see Fig. 20. For completeness, we evaluate the degree of fractional differencing using the python library fracdiff (de Prado 2018), for which the light curves become stationary, i.e. they pass the Augmented Dickey-Fuller test at the 95% confidence level, and at the same time they keep the maximum possible signal for the stationary case.
In Table 14, we list the same quantities and parameters for the case when the radio light curves are interpolated to regular 0.05-year time steps. We may notice small differences of the p and q order for the ARIMA models corresponding to the smallest BIC values. However, considering a smaller time step by a factor of two, p and q values increased at most by a factor of two, which implies that stochastic stationary processes are relevant on the timescales of ∼ 0.05−0.2 years, i.e. significantly smaller timescales in comparison with the precession and nutation timescales. The AIC and BIC values for different p and q parameters are graphically shown in Fig. 21 for the case of 0.05-year time step.
In Fig. 22, we fit the fractionally differenced radio light curves with ARIMA models, which again indicates that p = 1 and q = 1 for the minimum BIC cases for the 0.1-year time step. For the 0.5-year time step, p is equal to 3, 1, and 2 for 4.8, 8.0, and 14.5+15.0 GHz light curves, respectively, while q = 0 for the minimum BIC cases.
In Fig. 23, we present ARIMA(p, d, q) fits of preferred models to 4.8, 8.0, and 14.5+15.0 GHz light curves (from the left to the right panels; the regular time step is 0.1 years). For 4.8 and 8.0 GHz datasets, the preferred model based on the minimum BIC value is ARIMA(2,1,1), while for 14.5+15.0 GHz dataset, the preferred model is ARIMA (1,1,2).
In Fig. 24, we present forecasts of the trained ARIMA(p, d, q) models. The models were first fit to an approximately first half of light curves and then the forecasts of the ARIMA means, 1σ and 2σ confidence intervals were performed for the second half of the light curves. The test light curve data are denoted by dashed lines.
In addition to the ARIMA modelling of the radio light curves, we apply the model to the position-angle evolution of the stationary component a, see Fig. 25. The original evolution with 92 measurements is interpolated to the regular grid with the time step of 0.15 years. This evolution is not stationary and we have to apply the fractional differencing with d = 0.555 to remove the long-term trend, see the left panel of Fig. 25. Subsequently, we fit a series of ARIMA(p,q) models to the differenced evolution, see the middle and the right panels of Fig. 25 for the AIC and BIC distribution. The minimum AIC value is reached for p = 3 and q = 2 and the minimum BIC value is reached for p = 1 and q = 0, i.e. just the autoregressive polynomial seems to be necessary. Overall, the p and q order is rather small, implying the timescale for the stationary stochastic process in the range of ∼ 0.15 − 0.45 years.
We fit the ARIMA(1,0.555,0) model, which corresponds to the minimum BIC value, to the differenced temporal evolution of the position angle, see Fig. 26. The best-fit model is η(t) = (0.963±0.020)η(t−1)+(−27.851±16.033)+ (t), i.e. having just the autoregressive term. Next, we first fit ARIMA(1,0.555,0) to the first half of the evolution and   Table 14. Summary of the radio light curve analysis using ARIMA models. The light curves were interpolated to a regular time-step of 0.05 years.  Figure 20. Distribution of Akaike and Bayesian Information Criteria (AICs and BICs) for ARIMA(p,0,q) (upper two rows) and ARIMA(p,1,q) (lower two rows) models fitted to 4.8, 8.0, and 14.5+15.0 GHz light curves (from the left to the right columns). AIC and BIC minima are denoted by red-framed squares. The radio light curves were interpolated to the regular time-step of 0.1 years. The white squares in the middle upper two panels represent p and q parameters, for which the ARIMA model fit did not converge.
perform the forecast for the rest, see the right panel in Fig. 26. The forecast mean clearly deviates from the measured position angles, while within 1σ confidence region, the measured detrended evolution is nearly consistent with the forecast model.  Figure 21. Distribution of Akaike and Bayesian Information Criteria (AICs and BICs) for ARIMA(p,0,q) (upper two rows) and ARIMA(p,1,q) (lower two rows) models fitted to 4.8, 8.0, and 14.5+15.0 GHz light curves (from the left to the right columns). AIC and BIC minima are denoted by red-framed squares. The radio light curves were interpolated to the regular time-step of 0.05 years.