Precession-induced Variability in AGN Jets and OJ 287

To verify the applicability of the jet-precession+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), and Britzen 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-squares fitting with the concatenated χ²-statistic, ${\chi }^{2}=({\chi }_{\beta }^{2},{\chi }_{\eta }^{2})$, where ${\chi }_{\beta }^{2}$ is the χ²-statistic that corresponds to the apparent velocities, and ${\chi }_{\eta }^{2}$ is the χ²-statistic that corresponds to the position angles. The fitted parameters converge well, however, with a certain degeneracy close to the minimum χ², i.e., more solutions are possible with different χ² 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 yr.

To better illustrate this, we include two solutions, one with the reduced ${\chi }_{{\rm{r}}}^{2}={\chi }_{\nu ,\beta }^{2}+{\chi }_{\nu ,\eta }^{2}$ of 18.9 and the other with ${\chi }_{{\rm{r}}}^{2}=167.0$. The number of 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 Figure 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 ${\chi }_{{\rm{r}}}^{2}$, have consistent precession parameters within 1σ uncertainties.

Zoom In Zoom Out Reset image size

In Figure 1 (top-left panel), we use Equation (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 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 Figure 3 in more detail. The best fit ${\chi }_{{\rm{r}}}^{2}$ 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 favor the half-opening angle of the nutation cone of at least ∼3°, potentially as much as 5°, see the model with Ω_n = 5° shown in Figure 3 (blue-dotted line). To further constrain the precession-nutation parameters, the longer evolution of the component "a" will be beneficial and/or for detecting the nutation-like motion for more components.

Zoom In Zoom Out Reset image size

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 two-dimensional likelihood contours, which provide information about how tightly each parameter is constrained. To this goal, we maximize the likelihood function ${ \mathcal L }$ using the form

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\left[\displaystyle \frac{{[{q}_{i}^{\mathrm{obs}}-{q}_{i}^{\mathrm{theor}}({\boldsymbol{p}})]}^{2}}{{s}_{{\rm{i}}}^{2}}+\mathrm{ln}(2\pi {s}_{i}^{2})\right],\end{eqnarray} \tag{ 1 }$$

where ${q}_{i}^{\mathrm{obs}}$ are observed quantities (here the position angles and apparent velocities), ${q}_{i}^{\mathrm{theor}}({\boldsymbol{p}})$ is the theoretical model that depends on precession and nutation parameters p , ${s}_{i}^{2}={\sigma }_{i,q}^{2}+\exp {(2\mathrm{ln}f)({q}_{i}^{\mathrm{theor}})}^{2}$ is the variance, which contains an underestimation factor $\mathrm{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 ${{ \mathcal L }}_{\max }$ and infer p , we apply the python-based MCMC solver emcee.

The nonzero 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 $\mathrm{ln}{{ \mathcal L }}_{\max }$, the number of the model parameters k, the data set size N, and Akaike information criterion (AIC) and Bayesian information criterion (BIC) defined as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{AIC} & = & -2\mathrm{ln}{{ \mathcal L }}_{\max }+2k,\\ \mathrm{BIC} & = & -2\mathrm{ln}{{ \mathcal L }}_{\max }+k\mathrm{ln}N.\end{array}\end{eqnarray} \tag{ 2 }$$

Table 2. Summary of the MCMC-inferred Precession and Precession-nutation Parameters Based on the Observed Evolution of the Position Angles η, Apparent Velocities β_app, and Position Angle Evolution of the Quasi-stationary Component "a"

Parameter	η (°)	βapp (c)	ηa (°)
t0 [yr]	${1995.40}_{-0.99}^{+0.81}$	${2000.95}_{-29.68}^{+3.27}$	${2010.91}_{-0.07}^{+1.18}$
Pp [yr]	${23.83}_{-1.81}^{+2.42}$	${34.52}_{-7.58}^{+4.13}$	${31.61}_{-1.25}^{+2.82}$
Pn [yr]	⋯	⋯	${1.58}_{-0.02}^{+5.88}$
γ	⋯	${8.96}_{-0.67}^{+2.89}$	⋯
Ωp [°]	${6.47}_{-1.55}^{+1.77}$	${12.36}_{-4.17}^{+8.37}$	${11.73}_{-6.19}^{+11.06}$
Ωn [°]	⋯	⋯	${2.23}_{-1.04}^{+1.97}$
Φ0 [°]	${15.77}_{-3.70}^{+2.98}$	${11.93}_{-1.53}^{+2.30}$	${156.25}_{-29.38}^{+10.63}$
η0 [°]	$-{120.22}_{-3.97}^{+3.87}$	$-{124.48}_{-17.69}^{+17.41}$	$-{99.45}_{-1.97}^{+3.59}$
$\mathrm{ln}f$	$-{6.74}_{-2.19}^{+2.36}$	$-{1.96}_{-0.27}^{+0.30}$	$-{2.98}_{-0.18}^{+0.21}$
$\mathrm{ln}{{ \mathcal L }}_{\max }$	12.41	−35.85	33.72
Parameter number k (model parameters + $\mathrm{ln}f$)	5+1	6+1	7+1
Data point number N	21	21	92
AIC	−12.81	85.70	−51.44
BIC	−6.55	93.01	−31.26

Note. We also include information on the $\mathrm{ln}{{ \mathcal L }}_{\max }$, AIC, and BIC statistical parameters.

Download table as:  ASCIITypeset image

Overall, the inferred precession-nutation parameters are consistent for all three data sets: the precession period is ∼20–30 yr, 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, Φ₀ ∼ 10°–20°, or a 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 yr and the half-opening angle of the nutation cone of Ω_n ∼ 22. When we inspect the likelihood distributions in corner plots in Appendix A, we see that the viewing angle Φ₀ is not well determined for the position-angle fitting. Both the viewing and 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 Figures 15, 17, and 19 when fitting position angles, apparent velocities, and the quasi-stationary component position angles, respectively.

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

$$\begin{eqnarray}&&\delta (t,\gamma ,{\rm{\Phi }})={\left[\gamma \left(1-\beta \cos \left({\rm{\Phi }}(t)\right)\right)\right]}^{-1}.\end{eqnarray} \tag{ 3 }$$

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, 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 Figure 1, bottom-left panel).

First, we adopt the kinematical parameters corresponding to the best-fit precession-only model (with ${\chi }_{r}^{2}=18.9;$ see Table 1) as inferred from the joint least-squares 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}_{\nu ,0}^{\mathrm{prec}}$, α^prec, and ${t}_{0}^{\mathrm{prec}}$ are listed in Table 3 for the three frequencies, while the precession-model flux densities versus observed flux densities are shown in Figure 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 ∼1980–1990 and ∼2005–2015. However, there are residuals with mean values of 0.64, 0.99, and 1.16 Jy for the 4.5, 8.0, and 14.5+15.0 GHz data, respectively, that can be attributed partly to additional nutation-like jet motion and/or stochastic accretion- or instability-driven red-noise variability.

Zoom In Zoom Out Reset image size

Table 3. Best-fit Parameters Based on Fitting the Precession-only Model to the Observed Flux Densities of OJ 287

Parameter	4.8 GHz	8.0 GHz	15.0 GHz
${S}_{\nu ,0}^{\mathrm{prec}}$ [Jy]	1.73 ± 0.05	2.48 ± 0.05	2.60 ± 0.04
αprec	−1.73 ± 0.05	−1.76 ± 0.01	−1.69 ± 0.01
${t}_{0}^{\mathrm{prec}}$ [yr]	1993.8 ± 0.2	1992.6 ± 0.1	1993.3 ± 0.1
${\chi }_{{\rm{r}}}^{2}(\mathrm{dof})$	271.3 (975)	683.9 (1299)	680.0 (1722)

Download table as:  ASCIITypeset image

Second, we add the second-order nutation-like motion to the bulk precession motion of the jet, see Figure 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 Figure 4. In addition, the quasi-stationary radio component "a" of OJ 287 exhibits signs of a second-order, short-period motion, see Figure 3. For the precession kinematics, we fix the parameter values close to the best-fit precession model according to Table 1, in particular γ = 9.0, P_p = 22.9 yr, Ω_p = 103, Φ₀ = 15°, and η₀ = −1308. 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 Figure 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}_{\nu ,0}^{\mathrm{pn}}$, spectral index α^pn, nutation period P_n, and the reference epoch ${t}_{0}^{\mathrm{pn}}$. 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 Figure 5, including the residuals, for each frequency band. From the best-fit models, we obtain the inverted spectral index with the mean value of ${\overline{\alpha }}^{\mathrm{pn}}=-2.38\pm 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 ${\overline{P}}_{{\rm{n}}}=1.66\pm 0.14$ yr, which is consistent with the nutation period of 1.6 ± 0.1 yr 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 in 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 precession+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 2. This is especially apparent toward 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.

Zoom In Zoom Out Reset image size

Table 4. Best-fit Parameters Based on Fitting the Precession+Nutation Model to the Observed Flux Densities of OJ 287

Parameter	4.8 GHz	8.0 GHz	15.0 GHz
${S}_{\nu ,0}^{\mathrm{pn}}$ [Jy]	3.34 ± 0.05	4.46 ± 0.04	5.46 ± 0.04
αpn	−2.36 ± 0.02	−2.30 ± 0.01	−2.39 ± 0.01
Pn [yr]	1.86 ± 0.01	1.56 ± 0.01	1.40 ± 0.01
${t}_{0}^{\mathrm{pn}}$ [yr]	1987.27 ± 0.07	1986.43 ± 0.04	1987.12 ± 0.02
${\chi }_{{\rm{r}}}^{2}(\mathrm{dof})$	278.9 (973)	625.9 (1297)	659.5 (1720)

Download table as:  ASCIITypeset image

Since the orbital period of a putative SMBBH (e.g., 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 nutation—is generally expected, e.g., in analogy to the well-known binary systems SS 433 (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 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 short-term nodding motions in the precessing accretion flows in close binary systems. Here, we apply their formalism to SMBBHs. The nutation frequency may be expressed as (see also Caproni et al. 2013),

$$\begin{eqnarray}&&{\omega }_{{\rm{n}}}=2({\omega }_{\mathrm{orb}}-{\omega }_{{\rm{p}}}),\end{eqnarray} \tag{ 4 }$$

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:

$$\begin{eqnarray}&&{P}_{\mathrm{orb}}^{\mathrm{obs}}={P}_{{\rm{p}}}^{\mathrm{obs}}{\left(\displaystyle \frac{{P}_{{\rm{p}}}^{\mathrm{obs}}}{2{P}_{{\rm{n}}}^{\mathrm{obs}}}-1\right)}^{-1}.\end{eqnarray} \tag{ 5 }$$

For the case of OJ 287, the observed precession period is in the range of ${P}_{{\rm{p}}}^{\mathrm{obs}}\sim 20\mbox{--}30$ yr and the nutation period is in the range of ${P}_{{\rm{n}}}^{\mathrm{obs}}\sim 1\mbox{--}1.6$ yr, see Table 1, which yields the observed SMBHB period in the range ${P}_{\mathrm{orb}}^{\mathrm{obs}}\simeq 2.1\mbox{--}3.8$ yr. This period is about 3–6 times shorter than the orbital period of 12 yr 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 yr and the optical-flare periodicity of 12 yr (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 days 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 polarimetric monitoring. The amplitude of the nutation motion is then expected to be smaller than the precession amplitude by a factor of ${P}_{{\rm{p}}}^{\mathrm{obs}}/{P}_{\mathrm{orb}}^{\mathrm{obs}}\simeq 5.3\mbox{--}14.3$, applying the model of Katz et al. (1982) to OJ 287.

In Figure 6(a) we display the observed period of the SMBBH system OJ 287 that yields the observed precession period of 23 yr, 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 yr in the observed frame is marked as a shaded black rectangle in Figure 6, while the previously discussed orbital period of 11.7 yr related to optical outbursts is marked as a dotted black line. We see that for the assumed total mass of 4 × 10⁸ M_⊙ of the binary the primary mass fraction x_p is much broader for the shorter orbital period of 2–4 yr, x_p ∈ (0.5, 0.80–0.88), while for 11.7 yr, the range shrinks to essentially two equal-mass components, which is unlikely. This situation is the same for the larger total mass of 10¹⁰ M_⊙, see panels (c) and (d) in Figure 6. In addition, for the total mass of m_tot ≃ 10¹⁰ M_⊙ of the SMBH binary, whose components have the observed period of ${P}_{\mathrm{orb}}^{\mathrm{obs}}=10$ yr, the merger timescale decreases just to ≲10⁴ yr according to Peters (1964),

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{merge}} & = & \displaystyle \frac{5{c}^{5}}{256{G}^{3}}\displaystyle \frac{{a}_{0}^{4}}{{x}_{{\rm{p}}}{x}_{{\rm{s}}}{m}_{\mathrm{tot}}^{3}},\\ & = & \displaystyle \frac{5{c}^{5}}{{(4{\pi }^{2})}^{4/3}256{G}^{5/3}}\displaystyle \frac{{P}_{\mathrm{orb}}^{8/3}}{{x}_{{\rm{p}}}{x}_{{\rm{s}}}{m}_{\mathrm{tot}}^{5/3}},\\ & & \sim 12886{\left[\displaystyle \frac{{P}_{\mathrm{orb}}^{\mathrm{obs}}}{10\,\mathrm{yr}(1+z)}\right]}^{8/3}{\left(\displaystyle \frac{{x}_{{\rm{p}}}}{0.5}\right)}^{-1}{\left(\displaystyle \frac{{x}_{{\rm{s}}}}{0.5}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{m}_{\mathrm{tot}}}{{10}^{10}\,{M}_{\odot }}\right)}^{-5/3}\,\mathrm{yr},\end{array}\end{eqnarray} \tag{ 6 }$$

which appears to be statistically less likely, while for the total SMBH binary mass of 4 × 10⁸ M_⊙, τ_merge ≃ 10⁶ yr.

Zoom In Zoom Out Reset image size

For the above reasons, while doing estimates related to OJ 287, we have assumed a total SMBBH mass of 4 × 10⁸ M_⊙, which is preferred, e.g., in Liu & Wu (2002), based on the SMBH-bulge mass correlation. A similar mass, 1.3 × 10⁸ M_⊙, was derived by Neronov & Vovk (2011), based on the timescales of fast γ-ray flares. For completeness, in Figure 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 gray lines as in Figure 6(a).

Kushwaha (2020) discusses new spectral features that appeared during the 2015–2017 multiwavelength high activity state of OJ 287. In particular, a break in the near-IR-optical spectrum and hardening of the megaelectronvolt– gigaelectronvolt emission accompanied by a shift in the location of the peak, were observed. Figure 1(b) 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).

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. 2022), version 0.0.8, ¹⁵ 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²⁷ cm. The radius of the radiating blob is fixed at R_b = 3 × 10¹⁶ 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 Figure 7. We note that these results are possible realizations of a model with a large number of parameters.

Zoom In Zoom Out Reset image size

Table 5.  a: MJD 57359–57363, b: 57786. Γ = 10 (β = 0.995c) Fixed Based on the Present Work

Spectral State		a (${\chi }_{r}^{2}\approx 30$)	b (${\chi }_{r}^{2}\approx 2$)
Spectral normalization factor	ke (cm−3)	$({4.8}_{-2.2}^{+4.0})\times {10}^{-3}$	$({3.3}_{-1.1}^{+1.6})\times {10}^{-7}$
Low-energy electron spectral index	p1	1.91 ± 0.09	2.67 ± 0.06
High-energy electron spectral index	p2	4.56 ± 0.06	4.05 ± 0.09
Minimum e-Lorentz factor	${\gamma }_{\min }$	102.3	103.0
Lorentz factor at break energy	γb	103.5	104.3
Maximum e-Lorentz factor	${\gamma }_{\max }$	105.0	105.5
Magnetic fields strength	B (G)	${4.56}_{-1.80}^{+2.98}$	${0.16}_{-0.03}^{+0.04}$
Jet power in particles	Pp(erg s−1)	1.8 × 1045	2.5 × 1043
Jet power in the magnetic field	Pm(erg s−1)	1.4 × 1046	1.8 × 1043
Doppler factor	δ	3.7 ± 1.1	45.1 ± 6.0
Viewing angle of the jet	Φ (deg)	∼12.1	<2

Download table as:  ASCIITypeset image

In Table 5 we compare the parameters obtained by SED modeling of the data at the end of 2015 (from 2015 December 13 to 2015 December 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 (∼121 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 :

• 1.  
  Case 0: 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);
• 2.  
  Case 1: 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);
• 3.  
  Case 2: 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 );
• 4.  
  Case 3: 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 χ² between the SED_a model and SED_b data.

The reduced chi-squares are the following: ${\chi }_{r}^{2}=2.06$ (case 0, Table 5), ${\chi }_{r}^{2}=27.7$ (case 1), ${\chi }_{r}^{2}=377.5$ (case 2), and ${\chi }_{r}^{2}=450.39$ (case 3). It seems that AGNPY does not find a good fit when either the electron energy parameters or the Doppler factor is fixed, such that case 2 is worse compared to case 0, much worse than case 1 to case 0. Since the observed synchrotron and SSC flux depend on δ⁴, 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 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 in 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.

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 Figure 8 we show the position angle, viewing angle, and Doppler factor of a plasmoid in the jet stream of OJ 287, with the fitted parameters presented in Table 1.

Zoom In Zoom Out Reset image size

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 by about 80°. 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 that corresponds to the epoch closest to the lower projected precession reversal point is more than twice the peak Doppler factor derived from the jet precession+nutation model (≈18, see Figure 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 yr) still allow comparing 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.

We test if the three radio light curves could be modeled using the ARIMA model, 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 it 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 quasi-periodic behavior in AGN light curves, especially when only a 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}^{{\prime} }={f}_{t}-{f}_{t-1}$, can generally be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{t}^{{\prime} } & = & c+{\alpha }_{1}{f}_{t-1}^{{\prime} }+{\alpha }_{2}{f}_{t-2}^{{\prime} }+\ldots +{\alpha }_{p}{f}_{t-p}^{{\prime} }\\ & & +{\epsilon }_{t}+{\theta }_{1}{\epsilon }_{t-1}+{\theta }_{2}{\epsilon }_{t-2}\ldots {\theta }_{q}{\epsilon }_{t-q},\end{array}\end{eqnarray} \tag{ 7 }$$

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 yr. 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 B, in particular Table 13 and Figure 20. Using (partial) autocorrelation functions as well as AIC and BIC values obtained from fitting different ARIMA models (p = [0,...,15], d = [0, 1], and 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 differentiation of all 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 Figure 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 yr 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 Figure 21. The small increase in p and q order thus merely reflects the time step decrease by a factor of 2, while the timescale of the stationary stochastic process remains in the range of 0.05–0.2 yr.

Furthermore, we apply fractional differencing to all 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 yr, the term d reaches values intermediate between 0 and 1 and less than 0.5. For a smaller time step of 0.05 yr, 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 Figure 22. For the time step of 0.1 yr and the model with the smallest BIC, p and q are consistently equal to one, while for the time step of 0.05 yr, p is in the range of 1–3 and q = 0. Hence, again the characteristic timescale for the stationary stochastic process is ∼0.05–0.15 yr.

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 yr. 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 Figure 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 yr for stationary stochastic processes. In Figure 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 half of the corresponding radio light curves, tend to converge toward 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 yr; hence, they depend on at most two lagged flux density and white noise values going back in time by ≲0.2 yr. 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 disk-jet system, beyond which they do not affect significantly the emission processes.

Within the 2σ confidence interval, the forecast values are consistent with the observed radio outburst around ∼2010–2011, see Figure 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-curve trends, see Figure 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 yr 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 Figure 9. In contrast, on the timescales ≥1 yr the system is dominated by deterministic periodic processes related to precession and nutation jet motions, which is indicated by precession+nutation model fits and periodogram peaks, see Figures 10 and 11.

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.

Table 7. Previous Table Continued

AGN	TXS 0506+056	3C84	PKS 1502+106
Type	BL Lac	Radio galaxy	Quasar
z	0.3365	0.0176	1.838
Pp,obs [yr]	10.2 ± 1.1	40.5 ± 0.2	20.2 ± 2.7
γ	${2.1}_{-1.1}^{+1.2}$	1.05	1.6 ± 0.4
η0 [deg]	−93.7 ± 94.3	30	113.6 ± 11.4
Φ0 [deg]	10.3 ± 24.3	130	1.90 ± 1.4
Ωp [deg]	19.5 ± 31.9	45	0.7 ± 0.4
Considered period	2009–2018	2000–2018	1997–2020
Model	SMBBH or LT precession	SMBBH	SMBBH or LT precession
References	(10)	(11)	(12)

Note. References are (10) Britzen et al. (2019a), (11) Britzen et al. (2019b), and (12) Britzen et al. (2021).

Download table as:  ASCIITypeset image

As noted earlier, jet precession is intimately connected to the radio flux variability, via the Doppler factor. For favorable configurations, the Doppler beaming can significantly boost the nonthermal jet emission. For comparison, we calculate and plot in Figure 13(a) the time-variable (and periodic) Doppler boosting factor δ for most of the sources in Tables 6 and 7, using the tabulated values of the model parameters. These AGNs, which are classified as BL Lac objects (0735+178, 2200+420, PG1553+113, OJ 287, TXS 0506+056), quasars (3C 273, 3C 279, 3C 345, 1308+326, PKS 1502+106), Seyfert 1 galaxy (3C 120), and radio galaxy (3C 84), have estimated precession periods of the order of 1–40 yr (as deduced from the literature). Recently, jet precession for the nearby low-luminosity radio galaxy M81 was confirmed with a period of ∼7 yr (von Fellenberg et al. 2023). However, another solution is also discussed, specifically a precession+nutation model with the 7 yr period corresponding to nutation, while precession with the period of ∼800 yr 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.

Zoom In Zoom Out Reset image size

Observationally, perhaps a more meaningful parameter than the Doppler factor is the ratio of the maximum to the minimum Doppler factor, $\xi \equiv {\delta }_{\max }/{\delta }_{\min }={({S}_{\max }/{S}_{\min })}^{1/(p+\alpha )}$, which is computed from the temporal profile of the calculated Doppler factor inferred from the observational radio VLBI data, see Equation (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 Doppler factor ratio ${\xi }_{\max }$ (the minimum is unity, which implies zero or 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. (2019b), 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 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 Figure 2, bottom panels in Caproni et al. 2017).

Table 8. Derived Maximum Doppler Factor Ratio ${\xi }_{\max }$ for the AGNs Compared in This Paper

Source	${\xi }_{\max }$
3C 84	1.44
3C 120	1.35
TXS 0506+056	1.80
PKS 0735+178	1.63
OJ 287	13.61
3C 273	3.37
3C 279	18.69
1308+326	1.47
PKS 1502+106	1.003
PG1553+113	19.63
3C 345	1.59
2200+420	1.07

Download table as:  ASCIITypeset image

Homan et al. (2021) derive Doppler factors for the MOJAVE survey sources based on the assumption that all sources in their sample have the same intrinsic median brightness temperature T_B-int. The Doppler factor for each source is then derived by dividing the source's median observed brightness temperature by T_B-int . The majority of the derived Doppler factors for the MOJAVE sources lie in the range between 0 and 40 (Figure 7(a) in their paper). Hovatta et al. (2008) determined the variability brightness temperatures of the fastest flares in 22 and 37 GHz light curves. By assuming the same intrinsic brightness temperature for each source, they calculated the Doppler boosting factors. We compare the values obtained by Homan et al. (2021) and Hovatta et al. (2008) with the range of values for the Doppler factor (and the mean value) obtained by applying our precession model in Table 9.

Table 9. Range of Doppler Factors $({\delta }_{\min },{\delta }_{\max })$ Derived via the Precession Model (This Paper)

Source	${\overline{\delta }}_{\mathrm{prec}}$ $({\delta }_{\min },{\delta }_{\max })$	${\delta }_{{T}_{{\rm{B}}}}$ (Homan et al. 2021)	δvar (Hovatta et al. 2008)
3C 84	1.16 (0.95, 1.37)	4.3	0.3
3C 120	10.24 (8.73, 11.75)	5.7	5.9
TXS 0506+056	2.81 (2.01, 3.61)	1.8	⋯
PKS 0735+178	5.81 (4.42, 7.19)	5.1	3.8
OJ 287	9.40 (1.29, 17.51)	46.3	17.0
3C 273	6.04 (2.77, 9.32)	21.8	17.0
3C 279	5.19 (0.53, 9.86)	140.2	24.0
1308+326	26.12 (21.15, 31.08)	5.1	15.4
PKS 1502+106	2.93 (2.93, 2.94)	5.1	12.0
PG 1553+113	4.33 (0.42, 8.24)	1.4	⋯
3C 345	18.81 (14.51, 23.11)	47.7	7.8
2200+420	9.08 (8.79, 9.38)	18.3	7.3

Note. In addition, we list the mean Doppler factor calculated as $\overline{\delta }=0.5({\delta }_{\min }+{\delta }_{\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.

Download table as:  ASCIITypeset image

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 between the precession-derived Doppler factors and the other estimates are found for 3C 273, 3C 279, 1308+326, and 3C 345. In addition to 3C 273, the values obtained by Homan et al. (2021) also disagree with those obtained by Hovatta et al. (2008) for the three other sources.

Since Doppler factors cannot be measured directly, the derived values rely on model assumptions. For the Doppler factors obtained by 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, instability, or magnetic reconnection and does thus not necessarily serve as an 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 precession 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.

Precession and nutation of a jet basically reflect 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 nonaxisymmetric forces are present. In general terms, this requires a multicomponent system with misaligned rotation axes. These components could be a BH-disk combined or an SMBBH. Precession could then arise from

• (a)  
  Torques on the accretion disk induced by a secondary BH, or
• (b)  
  Torques on the jet launching primary BH induced by a secondary BH,
• (c)  
  The misalignment of the BH spin axis from the angular momentum vector of the accretion inflow, leading to LT precession of the disk that may lead to jet precession.

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 yr, 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 AGNs.

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 an SMBBH to form during the final merger phase (e.g., Mayer 2017, and references therein), or 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 frame dragging 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 disk-like 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 (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 that 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.

Mounting evidence for jet precession has been found for an increasing number of AGNs, as discussed above. Additional support to this scenario comes from observations of microquasars, 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). SS 433 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, 2015, 2017). Another such example is the microquasar 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 & Zimmermann2010).

Since three-dimensional (3D) general relativistic magnetohydrodynamic simulations are very CPU expensive, only recently has the evolving precession (that is essentially 3D) of such systems numerically been 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 an LT precession. These simulations may serve as clear evidence that the jets can be used as probes of disk precession. The precession timescales 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 simulations, 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 timescales (after several precession periods) and the system gradually realigns. Furthermore, Liska et al. (2019) have shown by highly resolved GRMHD 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 nonrelativistic 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 timescales to less than an orbital timescale of the secondary (however thousands of disk orbital periods at the jet footpoint).

In this paper, we focus on the 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 multiwavelength 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 multiwavelength observations of PG1553+113 to investigate if the observed modulation (γ-rays, optical) can be explained by 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 multiwavelength flaring in PG1553+113 has been presented by Caproni et al. (2017). Tavani et al. (2018) also discuss the possible γ-ray (from Fermi-LAT observations) signatures of an SMBBH in PG 1553+113. In the case of 3C 84, Britzen et al. (2019b) find evidence for precession of the central radio structure and a correlation between the radio and γ-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.

Since precession results in a time-variable Doppler factor, this should also influence the temporal evolution 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 multiwavelength 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 in the SED state.

Britzen et al. (2019a) 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, the 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 2021 April 18 by Baikal-GVD (Belolaptikov et al. 2022). On 2022 September 18, 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 an SMBBH model (Britzen et al. 2019a; 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).

Four independent neutrino detectors—including IceCube—reported neutrino emission in 2021 December (see, e.g., Belolaptikov et al. 2022; Sahakyan et al. 2023). PKS 0735+178 is located just outside the 90% point-spread function containment of the IceCube 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).

In the 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 the 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 & Gabányi 2014), LT precession (e.g., Britzen et al. 2019a), etc., the result is the same, namely, a constant pitch.

In the case of internal helicity, the pitch is not constant, and the jet components move along the helix (e.g., Hardee 1987; Steffen et al. 1995). 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 parsec-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 jet components 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 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 that neglects the dynamics of the jet flow. In fact, the nature and dynamics of the postulated jet blobs are 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 lighthouse model.

From an observational point of view, the periodicity timescale 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/γ-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.

There seems to be a need to explore an alternative to the currently popular model of the double-peaked optical outbursts appearing every 12 yr (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 yr periodicity), by suitably updating the input dynamical parameters. However, Komossa et al. (2023) did not detect the outburst predicted for 2022 October 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 BH (∼10¹⁰ M_⊙). Komossa et al. (2023) estimate the mass of the primary BH in OJ 287 to be ∼10⁸ 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 the gradual accretion of the disturbed disk gas onto 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⁵–10⁶ 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:

• 1.  
  We are NOT dealing with a thin accretion disk (BL Lac). Hence, the putative impact is a nebulous concept (see Villforth et al. 2010).
• 2.  
  For all the flares, the primary peak lacks a radio counterpart, while some secondary peaks have shown radio counterparts.
• 3.  
  Optical polarization 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).
• 4.  
  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 come from the observed drop in the fractional linear polarization 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 contamination may well be occurring, we would like to recall that the drop in a fractional polarization with a rise in brightness may also occur in the case of a mildly swinging synchrotron jet (Gopal-Krishna & Wiita 1992). Thus, we believe that definitive evidence in favor 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 m class optical telescope should suffice.

Recently, Suková et al. (2021) performed GRMHD simulations of an orbiting perturbed (secondary BH) 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 yr optical periodicity by the enhanced accretion that occurs every time the secondary component passes through the accretion flow (twice per orbit), while the 24 yr 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 de-boosted). However, the model does not provide a prediction for the light curves in different bands yet. In addition, the duty cycle of the modeled flares appears rather short in their model, while the large-amplitude radio flares seen in OJ 287 have a long duty cycle of ∼10 yr 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 behavior, 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 the case of Seyferts. A more detailed analysis of these issues is planned for a future publication.

While a coincidence is often seen between flaring at millimeter 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 AGNs. For γ-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 with Mrk 501, for example (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 a 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 formation (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-0, a dwarf starburst galaxy with a potential central massive BH, a precessing bipolar outflow connecting the region of the BH with a star-forming region has been reported (Schutte & Reines 2022).

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 Figure 4(b) in Britzen et al. 2018). This is progress compared to multiple multiwavelength campaigns, which, in the past, had to be planned and conducted blindly due to the lack of knowledge of predictable, correlated and interconnected processes. Applying this previous knowledge to predict the deterministic behavior will allow us to much more carefully explore also those processes that 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; Nathanail et al. 2020; Ripperda et al. 2020).

In Figure 14, we show the corner plot with marginalized one-dimensional likelihood distributions and two-dimensional likelihood contours for the precession parameters t₀, P_p, Ω_p, ϕ₀, and η₀ inferred from the temporal evolution of the ejected component position angles η. We also display the posterior distribution for the variance underestimation factor $\mathrm{log}f$, which is included in the definition of the likelihood function. Hence, we deal with six free parameters. In Figure 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 nonzero flat prior parameter ranges in Table 10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 10. Summary of Nonzero Flat Prior Parameter Ranges for the Precession Parameters Inferred from the Temporal Evolution of the Jet Position Angle

Parameter	Prior Range
t0 [yr]	[1970, 2010]
Pp [yr]	[10, 40]
Ωp [°]	[0, 90]
ϕ0 [°]	[10, 20]
η0 [°]	[−180, 0]
$\mathrm{log}f$	[−10, 1]

Download table as:  ASCIITypeset image

In Figure 16, we show the corner plot with marginalized one-dimensional likelihood distributions and two-dimensional likelihood contours for the precession parameters t₀, P_p, γ, Ω_p, ϕ₀, and η₀ inferred from the temporal evolution of the ejected component apparent velocities β_app. We also display the posterior distribution for the variance underestimation factor $\mathrm{log}f$, which is included in the definition of the likelihood function. Hence, we deal with seven free parameters. In Figure 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 nonzero flat prior parameter ranges in Table 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 11. Summary of Nonzero Flat Prior Parameter Ranges for the Precession Parameters Inferred from the Temporal Evolution of the Jet Apparent Velocity

Parameter	Prior Range
t0 [yr]	[1970, 2007]
Pp [yr]	[10, 40]
γ	[1, 50]
Ωp [°]	[0, 50]
ϕ0 [°]	[10, 15]
η0 [°]	[−150, −100]
$\mathrm{log}f$	[−10, 1]

Download table as:  ASCIITypeset image

In Figure 18, we show the corner plot with marginalized one-dimensional likelihood distributions and two-dimensional likelihood contours for the precession-nutation parameters t₀, P_p, P_n, Ω_p, Ω_n, ϕ₀, and η₀ 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 $\mathrm{log}f$, which is included in the definition of the likelihood function. Hence, we deal with eight free parameters. In Figure 19, we compare the posterior median model with the VLBI data, including 1000 random samples of parameter chains. This way the precession+nutation model uncertainties around the model median are apparent. We summarize the nonzero flat prior parameter ranges in Table 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 12. Summary of Nonzero Flat Prior Parameter Ranges for the Precession-nutation Parameters Inferred from the Temporal Evolution of the Stationary component "a"

Parameter	Prior Range
t0 [yr]	[1995, 2016]
Pp [yr]	[10, 40]
Pn [yr]	[0, 10]
Ωp [°]	[0, 90]
Ωn [°]	[0, 10]
ϕ0 [°]	[0,180]
η0 [°]	[−180, −50]
$\mathrm{log}f$	[−10, 1]

Download table as:  ASCIITypeset image