© The Authors 2026

1. Introduction

Blazars are a subclass of active galactic nuclei (AGNs) whose jet axes are closely aligned with the Earth’s line of sight. The small viewing angle leads to a number of unique observational properties due to relativistic Doppler boosting. These include superluminal motion (e.g., Rees 1966; Lister et al. 2019) and violent flux variability across the electromagnetic spectrum (e.g., Raiteri et al. 2017; Sobacchi et al. 2017), as well as visible changes in the jet morphology and orientation (e.g., Britzen et al. 2017; Kun et al. 2018).

The periodic and erratic changes in the position angle (PA) of the inner jet have been detected in many blazars (Lister et al. 2013, 2021; Kostrichkin et al. 2025). Such PA variations could be associated with the jet launching mechanism because they occur near the jet base, yet there is no unique paradigm pertaining to their origin. Some proposed mechanisms are instabilities in the jet flow, disk precession, or jet precession. The latter two can arise for a number of reasons, such as the Lense-Thirring (LT) precession of the disk or the precession of the black hole (BH) spin axis (Lense & Thirring 1918; Lu 1992), as well as the orbital motion in a supermassive black hole binary (SMBHB) in some cases (Begelman et al. 1980). One aspect of special interest is that inner jet PA changes often arise together with flux variability and with the flux periodicity on occasion (Britzen et al. 2009; Liu et al. 2012; Lobanov & Roland 2005; Kun et al. 2014; Britzen et al. 2018; Jiang et al. 2023).

3C 66A is a well-known BL Lac object. Very Large Array (VLA) observations indicate that the source has an extended jet with a PA of 170° and weak radio lobes toward both the north and south (Price et al. 1993; Taylor et al. 1996). Only the southern jet has been detected in very long baseline interferometry (VLBI) images. Cai et al. (2007) find two bendings at ∼1.2 and ∼4 mas from the core in 2.3, 8.4, and 22 GHz VLBI images obtained between 2001 and 2002. The bent structure is also confirmed by Zhao et al. (2015) using data from 2004–2005. 15 GHz VLBI observations from 2003–2013 moreover show complex jet kinematics (Lister et al. 2019). The bright features identified present a wide range of apparent speeds, from apparently stationary to superluminal motions. In addition, some features are characterized by transverse motion, while some are in the direction of the core, as also reported in Zhao et al. (2015). Similar complex kinematics were reported at the mm-regime (Jorstad et al. 2001, 2005, 2017). Using 10 years of 43 GHz data from the VLBA (Very Long Baseline Array) BU-BLAZAR monitoring program, Weaver et al. (2022) also reported a bent jet and found four stationary components, some of which present PA variability of more than 20°. Moreover, they determined six moving knots in the jet that reach apparent superluminal speeds. Recently, Kostrichkin et al. (2025) used an automated algorithm to study jet orientation changes in hundreds of AGNs. From 15 GHz and 43 GHz data, these authors also found a variable PA in the jet of 3C 66A.

In addition, 3C 66A is known for its optical flux variability, found periodic on timescales between approximately 2 and 5 years in a significant number of independent studies (Fan et al. 2002; Belokon & Babadzhanyants 2003; Kaur et al. 2017; Fan et al. 2018; Otero-Santos et al. 2020; Cheng et al. 2022). Cheng et al. (2022) argued that the two-year period is persistent in the historical light curve (LC) of 3C 66A and is not a transient feature. Proposed explanations for the periodicity include instabilities in the accretion disk or orbital motion in a SMBHB. However, there is no consensus yet regarding the origin of this optical periodicity. On the other hand, Otero-Santos et al. (2020) did not find any evidence of periodic behavior in the γ-ray LC of 3C 66A.

In this paper, we confirm that the 3C 66A parsec-scale jet is twisted (i.e., characterized by multiple bendings) using high-resolution KaVA (KVN and VERA array) images. Moreover, we report the detection of oscillation and continuous clockwise shift in the inner jet PA, based on long-term VLBA observations. The same PA periodicity is found in the downstream jet with a phase-shift and amplitude reduction. In addition, we report high radio flux variability, with three possible quasi-periods. These features have motivated us to investigate plausible physical scenarios leading to the observed jet orientation and flux periodicities. In particular, we investigated jet precession that could be triggered by a SMBHB system or the LT effect. We also explore further instabilities and disk-related mechanisms that might be relevant to 3C 66A.

The paper is organized as follows. In Section 2, we describe the datasets and the data reduction for the KaVa and VLBA data. In Section 3, we present our analysis and results concerning the jet structure and the core flux. In Section 4, we model and constrain precession parameters for the jet of 3C 66A. In Section 5, we investigate the SMBHB hypothesis and show that it can explain all the observables reported in this paper. In Section 6, we consider additional scenarios such as LT precession and instabilities to account for some of the properties of 3C 66A. Our findings are summarized in Section 7.

We adopted a standard Λ cold dark matter (ΛCDM) cosmology with parameters Ω_M = 0.3, Ω_Λ = 0.7, and H₀ = 68.0 km s⁻¹ Mpc⁻¹, rounded from the latest Planck Collaboration results (Planck Collaboration VI 2020; Tristram et al. 2024). 3C 66A has an uncertain spectroscopic redshift (0.335 < z < 0.444; Furniss et al. 2013) due to the lack of spectral features. Here, we use the redshift of z = 0.340 that was more recently determined from a cluster the host galaxy of 3C 66A belongs to (Torres-Zafra et al. 2018). At this distance, 1.0 mas corresponds to 5.0 pc.

2. Observations and data reduction

2.1. KaVA observations

At the time of these observations, KaVA consisted of seven antennas at three stations (Yonsei, Ulsan, and Tamna) of the Korean VLBI network (KVN) and four (Mizusawa, Iriki, Ogasawara, and Ishigakijima) of the VLBI Exploration of Radio Astrometry (VERA) array (Niinuma et al. 2014; Hada et al. 2017). KaVA is the precursor of the East Asia VLBI Network (EAVN) (An et al. 2018b; Akiyama et al. 2022). In 2014, 3C 66A was observed with KaVA during three epochs: April 16-17 (DOY = 106–107), September 30-October 1 (273–274), and December 15-16 (349–350). The observations for each epoch were conducted over two successive days, with the first day at 22 GHz and the second at 43 GHz. On April 16-17, the on-source time for 3C 66A was ∼6 hours at both frequencies, while it was ∼2 and ∼3 hours at 22 and 43 GHz, respectively, for the latter two epochs. The KVN-Ulsan station failed to participate in the second epoch observations due to instrumental malfunctions, and during the observation on December 15, several stations, especially those providing long baselines, suffered from severe weather conditions. Single polarization (left-hand circular polarization) was recorded with a total bandwidth of 256 MHz at each frequency: 16 IFs (intermediate frequencies or sub-bands) and 16 MHz per IF. The data were correlated at the Daejeon hardware correlator (Lee et al. 2015b).

A priori amplitude and phase calibration were conducted using the National Radio Astronomy Observatory (NRAO) AIPS software (Schwab & Cotton 1983). We first used opacity-corrected system temperatures and antenna gain curves to perform amplitude calibration. Additionally, we multiplied a factor of 1.35 to all visibility amplitudes in order to compensate for the known amplitude loss due to digital sampling and quantization (Lee et al. 2015a). The source 3C 84 or NRAO 150 was selected as the bandpass and manual fringe fitting calibrator, and KVN-Yonsei was used as the reference antenna. We corrected the parallactic angles for KVN antennas only. The global fringe fitting was done with a solution interval of 1 and 0.5 minutes for the 22 and 43 GHz data, respectively. The output visibilities from AIPS were imported to the DIFMAP package (Shepherd 1997) for imaging and model fitting. We performed an hybrid imaging procedure (Cornwell & Wilkinson 1981), with CLEAN (Högbom 1974) and phase-only self-calibrations iteratively in the early stage and amplitude self-calibration at the end of the process. We modeled the UV visibilities with several circular Gaussian components by modelfit. The Gaussian model fitting determined the flux, separation from the core, PA, and full width at half maximum (FWHM) of each component. The uncertainties of each parameter were calculated following the equations described in Lee et al. (2008).

2.2. VLBA archival data

We used all available 3C 66A VLBA data from observations performed at 43 GHz between 1996 and 2025. Out of the total 109 epochs, 68 were obtained within the VLBA-BU-BLAZAR program observed between October 2008 and June 2020. The detailed description of this program and the data reduction are presented in Jorstad et al. (2017) and Weaver et al. (2022). The five epochs obtained in 2024-2025 are from the BEAM-ME project¹ by the same group. The observations were carried out with 9 to 10 antennas (36 or 45 out of 45 possible baselines) and the typical integration time is of 40 minutes.

With the CLEAN algorithm, by carefully setting CLEAN windows, we imaged the self-calibrated UV data without further self-calibration. The 36 other epochs of the 43 GHz VLBA archival data, spanning the period August 1996–November 2005, were taken from the NRAO Data Archive System. These observations were also carried out with 9 to 10 antennas and have a similar integration time. We performed the data reduction using AIPS and DIFMAP in a similar way to that described in Section 2.1. For the archival VLBA data, an opacity fitting was applied during the amplitude calibration. One of the south-west antennas (“FD”, “PT”, “LA”, or “KP”) was selected as the reference antenna and one of the bright, well-sampled sources (3C 279, 3C 345, 3C 454.3, or 0420−014) was selected as the calibrator source for bandpass calibration and manual pulse-calibration. The solution interval was set to 30 seconds. After AIPS calibration, imaging was performed in DIFMAP and we fitted circular Gaussian components to each self-calibrated VLBA dataset.

3. Results

3.1. Parsec-scale jet morphology

Figure 1 shows the total intensity contour maps obtained from the KaVA observations at 22 and 43 GHz with the Gaussian jet components superimposed. The mapping parameters are summarized in Table E.1. A core-jet structure is found in the north-south direction down to ∼16 mas at 22 GHz (upper panel of Figure 1) and down to ∼3 mas at 43 GHz (lower panel of Figure 1). As the observing circumstances for each epoch are not the same (see Section 2.1), we can see various dynamical ranges and restoring beam shapes (Table E.1). The parameters and errors of the Gaussian components for each epoch are listed in Table E.2. The brightness temperature of each component at the observing frequency ν in the source frame is calculated with the equation (e.g., Condon et al. 1982): T_b = 1.22 × 10¹²Sν⁻²d⁻²(1 + z) K, where S is the flux density of the component in Jy, z is the redshift of the source, and d is the FWHM of the circular Gaussian component in mas. Component K is identified as the core, since it is located at the most upstream of the jet with the largest flux density, the highest brightness temperature, and the flattest spectrum (Zhao et al. 2015). The core separations and PAs of the other components are relative to component K.

Fig. 1. KaVA intensity maps of the 3C 66A radio jet at 22 GHz (top) and 43 GHz (bottom) in 2014. The black contour levels start at three times the rms value, scaling twice. The grey shaded ellipses above each contour show the restoring beam. The detailed imaging parameters are summarized in Table E.1. The horizontal spacing is proportional to the observation time gap. The exact observing time is in the unit of years. The red circles are the model components obtained from modelfit. The continuous curved jet extends downwards to ∼3 mas on both panels.

We found evidence of multiple bendings in the VLBI jet of 3C 66A. At ∼1.5 mas downwards from the core in the 43 GHz maps (bottom panel of Figure 1), the jet curves from the southwest to the southeast. In the 22 GHz maps (top panel of Figure 1), the jet axis further bends to the south direction between 4 and 8 mas. We found additional bending between 8 and 15 mas based on the bright downstream knot visible at ∼15 mas in the southeasterly direction. We define this structure as twisted, in the sense that there are several bendings with changing direction along the jet. This apparent morphology is confirmed by the PA distribution of the Gaussian components as a function of the core separation (Figure 2). The bendings found are consistent with the results of previous studies (Cai et al. 2007; Zhao et al. 2015) and are visible up to ∼20 mas southwards in the lower frequency (2.3 and 8.4 GHz) maps of Zhao et al. (2015). This twisted structure is also visible on smaller scales with the 43 GHz VLBA data, as illustrated by the ridgelines of six different epochs in Figure A.3.

Fig. 2. PA of the Gaussian components of the 3C 66A jet as a function of the core separation obtained with the KaVA observations at 22 (left) and 43 GHz (right) in 2014. Each epoch is marked by green circles, blue triangles, and red squares.

3.2. Evolution of the jet orientation

We studied the changes in the jet orientation for the 43 GHz VLBA archival data. We used two methods to determine the PA of a region in the jet: estimating the PA of the circular Gaussian components obtained from modelfit and estimating the PA from the clean components in annular bins. Since the clean components are the deconvolved solutions of the source structure, we were able to derive the PA of the jet by using the distribution of clean components (e.g., Lister et al. 2013; Rani et al. 2014). We began by analyzing the inner jet PA evolution. Using Gaussian fitting, the inner jet region corresponds to the component J1 (see Figure A.3). This component remains at a core separation of 0.2–0.6 mas throughout all epochs. For the clean components method, we set an annular bin (called annulus A) in the radius range from 0.15 to 0.45 mas away from the (0, 0) central point. We adopted the lower limit of the annulus to prevent distortion of the determined PA by the core component. We defined the flux-weighted average of the PA from the clean components in the bin as the inner jet PA and estimate the statistical error of the PA using the bootstrap method. We randomly resampled the clean components from the bin and calculate the flux-weighted mean PA from the bootstrapped sample. After repeating the calculation 10 000 times, we adopted the standard deviation as the statistical error. Moreover, we considered systematic errors to the PA to be a fraction, f_err, of the beam size. This fraction depends on the bandwidth Δν (that ranges between 20 MHz and 486 MHz) of each observation: f_err = 0.2(Δν_ref/Δν)^1/2, with Δν_ref = 32 MHz. We took the PA uncertainty to include both the statistical and systematic error contributions. The inner jet PA evolution obtained from the Gaussian component J1 and and from the clean components in annulus A are almost identical (see Figure A.2). The upper panel of Figure 3(a) shows the result for the annular bin method, presenting drastic changes in the PA over time. From 1998 to 2004, the PA decreases gradually, and it increases afterwards until 2013. It decreases once again until 2017 and increases until 2020, which is the last epoch observed prior to a four-year gap.

Fig. 3. (a) Inner jet PA versus time at 43 GHz. The dashed curve in the upper panel is the analytic sinusoidal+linear fit. The lower panel shows the residuals. The fitting parameters are summarized in Table 1. (b) Time evolution of the total (blue) and core (black) flux density.

The variations observed in the inner jet PA yield evidence of a periodicity, which we confirmed using three periodicity search methods (details in Appendix B). We applied these methods to the annulus A PA data. For all the methods, the standard deviation on the periodicity was obtained from the FWHM of the peak found in the period domain. Firstly, we applied the well-established Jurkevich method (Jurkevich 1971). We detected a clear minimum for a period of approximately 12 years, which reaches well below the threshold typically used to determine significance, hinting at a reliable periodicity signal. The drop in the total variance is very wide and its gaussian fit indicates a period of P = 14.44 ± 5.07 years. Secondly, we used the Lomb-Scargle (LS) periodogram (Lomb 1976; Scargle 1982) and obtained a periodicity at P = 10.78 ± 2.34 years, largely above both the 5σ significance level assuming red noise and the 4e−5 false alarm level (FAL); namely, the level needed to attain 0.004% false alarm probability (FAP). Finally, we apply the weighted wavelet Z-transform (WWZ) method (Foster 1996; Grossmann & Morlet 1984). We find a peak at P = 16.76 ± 4.25 years in the time-shift averaged WWZ power, once again with more than 5σ significance assuming red noise. We note that, while the WWZ method is characterized by time-frequency resolution and the ability to detect transient quasi-periods, it does not include data error when searching for periodicity. This could explain, in particular, the longer period found with respect to the LS periodogram. The large standard deviations in all the periodicity search techniques are due to the “short” observational time span with respect to the detected period and to the wide data gaps in the historical PA curve. While the possibility of irregular variations cannot be fully ruled out, these independent search methods, which include red noise estimation, all confirm a 5σ periodicity signal. Continuous monitoring of 3C 66A will allow us to better constrain the exact value of the periodicity.

We fit the data to a harmonic function with a global linear trend to characterize the PA evolution (e.g., Martí-Vidal et al. 2011; Kun et al. 2014), expressed as

$$\begin{array}{c}\hfill \theta =Asin\left(\frac{2\pi }{P}(t-t_0)\right)+S(t-t_0)+{\theta }_0,\end{array}$$(1)

where A is the magnitude, P the period, t₀ a reference time, S the slope of the linear trend, and θ₀ the PA at t = t₀. We performed a nonlinear least square fit using the Levenberg-Marquart method that searches the parameter set where χ² is minimized. Table 1 gives the fitting parameters. The result indicates that the PA changes sinusoidally with a magnitude of 7.77 ± 0.79 degrees and a period of 10.94 ± 0.22 years. Additionally, the mean PA drifts linearly by −0.83 ± 0.07 degrees per year. As expected, the reduced chi-square χ_red² is ∼2.25. The residual of the fit is shown in the lower panel of Figure 3(a). Table E.1 also shows the fitting results for the inner jet PA obtained with the Gaussian component J1, which are in very good agreement with the values from the annular bin approach. We therefore extract two evolution timescales from the inner jet PA evolution: the ∼11 years sinusoidal motion and a clockwise shift of ∼ − 0.8 degrees/year.

The inner jet PA oscillation has not been noticed in previous studies (likely due to relatively short time spans) to the spatial resolution and to variability criteria. Jorstad et al. (2001, 2005, 2017) respectively covered 2 (1995–1997), 3 (1998–2001), and 5 years (2008–2013), making it difficult to detect the oscillation, although they used 43 GHz VLBA data as in our study. 3C 66A was also monitored with the VLBA at 15 GHz for 10 years (2003–2013) under the monitoring of jets in AGNs with VLBA experiments (MOJAVE) program (Lister et al. 2013). However, it remains difficult within this duration to detect the periodicity. While Lister et al. (2019) covered the period 1994–2016, observations at 15 GHz have three times lower resolution and probe a relatively outer region compared to the 43 GHz observations. Weaver et al. (2022) presented the results from the 43 GHz VLBA data between 2007–2018 and estimated the PA evolution of eight knots in 3C 66A, classifying the behavior as constant overall (Weaver et al. 2022, Figure 6). This is related to the criterion chosen for the definition of constant evolution and to the time span of ∼11 years. We note that the evolution of their knot A2, which corresponds to our inner jet component (J1 or annulus A), presents a PA increase peaking around 2013 and then a decrease until 2017, in agreement with part of the second period that we observe. The difference between component A2’s minimum and maximum PA is slightly more than 20 degrees, which matches the amplitude of the PA sinusoidal fit that we obtain. Recently, Kostrichkin et al. (2025) analyzed the inner jet direction of 3C 66A as one of the AGNs for which the PA change was studied with an automated algorithm. From 43 GHz observations spanning slightly more than a period, the authors find a similar oscillation as the one reported here, although the data is fitted with a linear regression. The same is done for the 15 GHz PA data, that is also analyzed over approximately 10 years and which reports less variability. In comparison, our analysis covers approximately 30 years with high resolution (typical beam size of 0.4 × 0.2 mas), so that we can detect nearly three periods of the PA swing in 3C 66A.

On top of the evolution of the PA for the inner jet, we study the PA evolution of two components further from the core (details in Appendix A). As described above, the PA of the additional components is also estimated using two methods and their resulting evolution is equivalent (the results are shown in Figures A.1 and A.2). The second component is part of the core-jet structure that extends southwards over ∼1 mas: Gaussian component J2 or annulus B. We call it downstream component, since it is found south of the inner jet one. The third component is the separated region found at ∼2.5 mas from the core: Gaussian component J3 or annulus C. The position of these components can be seen clearly in Figure A.3. The downstream component presents a PA swing apparently similar to that of the inner jet one. We fit the PA evolution with Eq. 1 and find a similar period of 10.50 ± 0.29 yr topped with a slightly slower clockwise shift of −0.68 ± 0.06 degrees/year from the annular bin estimate. Moreover, the amplitude of the oscillations is significantly reduced to 4.35 ± 0.56 degrees and the sine is phase-shifted by ∼1.5 years compared to the inner jet oscillation. The fit to the PA evolution of Gaussian component J2 leads to a similar result as obtained for annulus B. In the case of the third component, there is no periodicity in the PA evolution, but it follows a decreasing linear trend. The annulus C PA yields a continuous shift of −0.31 ± 0.04 degrees/year, identical to the Gaussian J3 PA with a value of −0.31 ± 0.03 degrees/year. The details of the fits to the PA can be found in Table 1. It appears that the components studied all undergo a long-term clockwise PA shift, which is strongest for the inner jet, weaker for the downstream jet, and less than half the rate for the furthest component.

3.3. Flux density variability

To compare the PA evolution with the radio flux density, we consider the sum of all clean components’ flux density as the total flux density, and the sum of clean components within a radius of 0.20 mas as the core flux density. The latter is in excellent agreement with the flux density obtained from the fitted circular Gaussian core, as well as with the peak flux density of the convolved map. We adopt one-tenth of the flux density as its error.

Figure 3(b) shows the total and core flux density over time. The total flux density varies between 0.2 and 1.1 Jy. Except for the overall higher flux, its evolution appears identical to the core flux density’s, as it is core-dominated. The major features in the flux density variability therefore originate at the base of the jet. The outbursts in the flux can be associated with either the emergence of VLBI components or with geometric processes as a result of Doppler boosting.

If the latter is the case, it should be manifested as the correlation between the core flux and the jet axis direction. To check the correlation between the core flux density and the inner jet PA, we apply the discrete correlation function (Edelson & Krolik 1988, DCF). We probe the time lags between −7 and +7 years. The time-lag bins are set as 62 days which is the median cadence of the observations. To test the significance of the correlation in the DCF analysis, we simulate 10 000 LCs having the same power spectral density (PSD) and probability distribution function (PDF) as those of the observations (Emmanoulopoulos et al. 2013).

Figure 4(a) shows the DCF results. The grey dotted and dash-dotted lines mark the 2σ and 3σ significance levels, respectively. The correlation coefficient between the inner jet PA and core flux density slightly exceeds 2σ at a positive time lag of ∼2 years, indicating a weak but possible correlation with PA changes leading the core flux density changes. The correlation is also hinted at by the behavior of the core flux density versus PA (Figure 4(b)). The left and right panels of the plot show the results derived from the observations in 1996-2010, and 2010-2020, respectively (we leave the five 2024-2025 epochs in gray). Each panel shows crescent-like patterns for the two subdivided datasets that match the inner jet PA decrease-and-increase repeated twice, corresponding approximately to the two periods observed in the jet orientation evolution until 2020. This pattern could be the indication of a potential correlation, as argued by Rani et al. (2014) for the BL Lac object 0716+714.

Fig. 4. (a) DCF analysis between the PA and the core flux density. The grey dotted and dash-dotted lines show the 2σ and 3σ significance levels, respectively. (b) Core flux density vs. inner jet PA plot. Each panel displays the data between 1996 and 2010 (left) and between 2010 and 2025 (right). We use changing colors to show the time sequence: grey denotes the values from all observing epochs, purple to dark blue to clear blue denotes the data observed between 1996, 2002, and 2010, while green to orange to red between 2010, 2017, and 2020.

Motivated by the optical flux periodicity detected in a number of studies on 3C 66A (Fan et al. 2002; Belokon & Babadzhanyants 2003; Kaur et al. 2017; Fan et al. 2018; Otero-Santos et al. 2020; Cheng et al. 2022) and by the apparent high variability of the radio LC, we also perform a periodicity search for the core flux density using the three methods mentioned in Section 3.2 (details in Appendix A). The analysis with the Jurkevic method yields a peak at 3.53 ± 0.25 years, at 6.92 ± 0.21 years and at 10.51 ± 0.56 years that reach the F = 0.5 threshold, hinting at a strong periodic signals. The LS periodogram results indicate a peak at 6.22 ± 0.85 years which attains the 3σ significance level considering red noise. Two other peaks, at 3.43 ± 0.30 years and 2.38 ± 0.15 years, reach the 2σ significance level. Using the WWZ method, we find two peaks in the time-shift averaged power that reach the 2σ significance level. The longest period is ∼10 years and merges with the ∼7 years peak. We also find a stronger signal reaching the 3σ significance level and which is a better-defined Gaussian peaking at 4.98 ± 0.49 years. Taken together, these results indicate three possible periodicities in the radio LC of 3C 66A: at ∼10 years, ∼7 years, and ∼3–4 years. Each periodicity may not be exactly sinusoidal, despite the detection with methods such as the LS periodogram. For instance, the radio flux may be characterized by repetitive flares as has been observed in several sources (Jiang et al. 2022; Britzen et al. 2023).

In addition to the high variability of the core and total radio flux densities, we draw attention to their overall decrease, in particular between 1996 and 2020. The baseline is clearly decreasing over time, with a rate of −9.40 ± 1.43 mJy/year, obtained from a linear fit of the core flux density.

4. Jet precession model

4.1. Case for two precession timescales

Radio-interferometric monitoring studies of AGNs find an increasing number of sources showing jet PA variations (e.g., Lister et al. 2013, 2021; Cui et al. 2021; Kostrichkin et al. 2025, and reference therein), often accompanied by flux periodicity at one or more wavelengths (e.g., Kun et al. 2014; Britzen et al. 2018; Dey et al. 2021; Jiang et al. 2023). While several possible scenarios to trigger this phenomenon have been suggested, jet precession is the most frequently considered model to explain periodic PA behavior (e.g., Martí-Vidal et al. 2011; Roland et al. 2015; Caproni et al. 2017; Qian et al. 2019; Britzen et al. 2018; Cui et al. 2023; Jiang et al. 2023). The model predicts that newly-born components close to the core orient at a different PA from older components farther away from the core, resulting in a bent jet as witnessed in 3C 66A. Fendt & Yardimci (2022) show with 2D (M)HD (magnetohydrodynamic) simulations that precessing nozzles lead to S-shaped jets, similar to the 3C 66A twisted jet structure revealed by our KaVa observations. We believe that jet motion in 3C 66A is governed by two precession timescales: the first associated to the ∼11 year PA oscillation and the second to its slow re-orientation.

First, regarding the decade-long precession, the wide range of apparent speeds and the variability of the core flux density in 3C 66A could be explained as a result of the variation in the Doppler factor depending on the viewing angle. The possible correlation between the PA change and the core flux variability discussed in Section 3.3 is consistent with this argument. In addition, Mohana et al. (2021) find a clear baseline flux decrease in the Fermi γ-ray LC of 3C 66A occurring in mid-2011. Their results show a persistent high-baseline flux for 3 years beforehand and a persistent low-baseline one during the subsequent 8 years. They find evidence of a similar baseline flux drop in the optical band. To explain these baselines’ variations the authors argue in favor of Doppler factor or jet inclination change. The flux drop date coincides with the second peak that we observe in the inner jet PA. Given the observational timescale analyzed by Mohana et al. (2021), the 3 years of high baseline match the period in which the PA changes westward to an average of 200°, while the low baseline matches the PA southward shift reaching 180°. We consider this coincident multi wavelength flux drop an additional argument in favor of jet precession for the decade-long PA oscillation. A similar relation between flux density and jet orientation change was reported, for example, by Martí-Vidal et al. (2011) in the case of M81, as evidence of a precessing jet. More recent data continue to support the precession model for this source (Jiang et al. 2023).

Second, the PA evolution in 3C 66A presents a superposed long-term shift of approximately −0.8° per year. This slow re-orientation could constitute a fraction of a longer periodic variation in precession models where two periods coexist. Indeed, this slow PA re-orientation, whether positive or negative (anti-clockwise or clockwise), has been observed in other sources and has been modeled as a precession occurring over periods of the order of magnitude of 10³ years (e.g., Lobanov & Roland 2005; Jiang et al. 2023). In our case, the continuous decrease of the radio core flux density also supports the precession origin for the long-term PA shift. Comparing the core baseline flux decrease to the inner jet PA evolution, we observe approximately a 10 mJy flux decrease for a re-orientation of −0.8°. We chose not to model this long-term precession, given the much longer timescale on which it occurs compared to our observation time span. Indeed, we expect the period to be of approximately 1000 years based on comparison with Lobanov & Roland (2005). We also performed a sinusoidal fitting of the inner jet PA with a long-term sine. While a wide range of parameters allow to approximately recover the −0.83°/year continuous shift over the 29 years of our time span, the period must fall approximately between 500 and 1000 years.

4.2. Model parameters and uncertainties

We propose a simple precessing jet model to reproduce the approximate 11 year period oscillation in the PA of 3C 66A and use the value of −0.83°/year obtained from the sinusoidal plus linear fit of the PA (from the annular bin method) to account for its long-term decrease. With the reduced number of parameters, we obtain in this way a better constraint from the data. The precession model consists of a solid cone of half-opening angle Ω, precessing with period P_p (angular velocity ω = 2π/P_p) about the precession cone axis directed along z_p. We refer to Figure 11 in Britzen et al. (2018) for an equivalent geometrical scheme. In order to express the PA θ in the observer reference frame (Cartesian coordinates x, y, z) as the projected jet axis from the jet reference frame, we must perform a series of rotations, as detailed for example in Cui et al. (2023). The angle between the precession cone axis and the line of sight z₀ is denoted ϕ₀ and the projected precession axis onto the sky plane makes an angle θ₀. At any time t, the observed PA is:

$$\begin{array}{cc}\hfill \theta (t)& =arctan\left(\frac{y(t)}{x(t)}\right),\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill x(t)& =A(t)cos{\theta }_0-B(t)sin{\theta }_0\text{,}\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill y(t)& =A(t)sin{\theta }_0+B(t)cos{\theta }_0,\hfill \end{array}$$(4)

where the quantities A(t) and B(t) are given by

$$\begin{array}{cc}\hfill A(t)& =cos\mathrm{\Omega }sin{\varphi }_0+sin\mathrm{\Omega }cos{\varphi }_{0}sin[\omega (t-t_0)],\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill B(t)& =sin\mathrm{\Omega }cos[\omega (t-t_0)],\hfill \end{array}$$(6)

and t₀ is the reference time. There are five unknown parameters in the model: Ω, P_p, t₀, ϕ₀ and θ₀. We perform a Bayesian analysis and use the Markov chain Monte Carlo (MCMC) Python module emcee² (Foreman-Mackey et al. 2013). The details of the Bayesian analysis are given in Appendix C. The half-opening angle Ω of the precessing cone and the angle between the cone axis and the line of sight ϕ₀ are degenerate, as visible in the purple corner plots shown in Figure 5. These correspond to the case ϕ₀ ∈ [0, 90]°, which is the expected parameter range when no additional information on the viewing angle Θ is known. The viewing angle relates to the precession model quantities as

Fig. 5. Corner plots of the likelihood analysis for the precession model given three ranges of ϕ0. The parameter values shown are the median and 1σ credible intervals of the posterior samples.

$$\begin{array}{c}\hfill \mathrm{\Theta }(t)=arcsin\sqrt{x{(t)}^2+y{(t)}^2}.\end{array}$$(7)

A constraint on the viewing angle reduces the parameter range for ϕ₀ and consequently for Ω. Since 3C 66A is a BL Lac object, we expect the viewing angle to be below 20°. More restrictively, Weaver et al. (2022) find that the viewing angle of BL Lac objects peaks between 2 − 4°. In the case of 3C 66A, Jorstad et al. (2017) find an average viewing angle ⟨Θ⟩ = 1.7 ± 0.5° based on two moving knots. This value has been used for example in spectral energy distribution (SED) modeling by Mohana et al. (2021). Weaver et al. (2022) estimate, for 5 different jet knots in 3C 66A, viewing angles ranging from 2.6 ± 0.7° to 5.0 ± 0.5°. Given these constraints, we perform the likelihood analysis for two additional cases: the relaxed constraint case ϕ₀ < 15° which is equivalent to Θ < 18° approximately (blue plots in Figure 5), and the tight constraint case ϕ₀ < 5° corresponding to Θ < 6° (green plots in Figure 5). Such a joint constraint approach based on the viewing angle or inclination of the source has been employed by a number of authors (Kun et al. 2014; Cui et al. 2023; Jiang et al. 2023). With the additional viewing angle prior, we obtain a good fit to the data in both cases for an equivalent chi-square of 3.43. This is the same value as obtained when fitting the observed PA to a simple sinusoidal plus linear shift, showing that the precessing jet model fits the data as successfully as possible assuming the PA behavior is sinusoidal. We expect the deviation from the ideal value of 1 to be caused by fluctuations of the PA that occur on smaller scales than the precession period, typically between days and months, for example because of instabilities in the jet or component ejections (Zhao et al. 2015; Weaver et al. 2022).

We consider the tight constraint case ϕ₀ < 5° to be the most physically relevant one, due to the small reported viewing angles for 3C 66A. We therefore detail the jet precession parameters for the prior ϕ₀ < 5° in Table 2. They are given as the mean and standard deviation of the MCMC samples. The precessing cone has a very small half-opening angle and the precession axis is closely aligned to the line of sight, by a few degrees. The period obtained from the MCMC mean is of 10.88 years, which is in particularly good agreement with the periodicity detected with the LS periodogram, peaking at 10.78 years, and also agrees with the Jurkevic search method within the uncertainties. We note that this period obtained from the MCMC samples does not depend on the prior chosen for the viewing angle. In the following sections, we take the precession period, P_p = 10.88 ± 0.24 year, to describe the inner jet PA short-term oscillation.

4.3. Possible precession-driven Alfvén waves

Before exploring the origins of the observed precession, we include some considerations on the short-term PA oscillation of the downstream jet. Indeed, a ∼10-year period was also measured for the annulus of 0.45 < r < 1.00 mas using the annular binning method, as well as for the second component using the Gaussian fitting method (see Appendix A). In both cases, the periodic PA pattern shows a reduced amplitude and a phase shift. At the large redshift of 3C 66A, 1.0 mas corresponds to approximately 5.0 pc. If we consider the deprojected distance, taking the approximate viewing angle Θ = 5° (Weaver et al. 2022), 1.0 mas corresponds to ∼60 pc, or more in the case of a smaller viewing angle. Studying the PA change over small annuli along the jet, we determine that the ∼10 year PA oscillation is observed up to 1.3 mas from the core, equivalent to ∼80 pc in deprojected distance. However, reducing the annular size implies that we have less data to measure the PA at each epoch and we cannot fit the periodic behavior reliably. We therefore have two confident measures of the oscillation pattern only: at 0.3 mas and 0.7 mas from the core on average. This makes it challenging to interpret the change in this pattern with distance from the core.

Nonetheless, the reduced PA oscillation amplitude with r implies that the pattern is not ballistic. Regarding the phase-shift, we consider two cases. First, it may be caused by a complex geometrical structure of the jet, in which case the PA oscillation in the downstream jet directly translates the precession motion of that region. Second, we can interpret the phase-shift as a signature of propagation of the periodic PA pattern. If we assume the latter, we can take the shift to be Δτ_shift ≈ 1.5 years between the inner and downstream jet regions (0.3 and 0.7 mas from the core, respectively). The pattern appears to be a downstream-propagating wave-like disturbance with a rough speed of 0.3 mas/year, which corresponds to the apparent speed β_wave^app ≈ 5. Here, we note the similarity with the transverse patterns propagating superluminally on the jet ridge line of BL Lacertae, observed by Cohen et al. (2015) from 15 GHz MOJAVE data. The authors interpret this motion as Alfvén waves that propagate downstream, assuming a well-ordered helical magnetic field. They moreover find a correlation between the waves’ timing and direction and the swinging of a recollimation shock, located at approximately 0.34 pc (projected). This leads the authors to argue that the swinging of the shock component at the jet base excites the MHD waves, which propagate in a way resembling “waves on a whip” along the jet. The approximate propagation speed estimated for the periodic PA pattern in 3C 66A is similar to that of the transverse waves in BL Lacertae. In addition, Cohen et al. (2015) considered the transverse jet motion over a distance of approximately 13 pc deprojected, which is of the same order of magnitude as the one probed by our inner and downstream jet measurements. The periodic PA pattern observed in the downstream jet of 3C 66A may therefore be the signature of an Alfvén wave excited by the precession of the jet nozzle and propagating on the longitudinal component of the magnetic field.

5. SMBHB system

Jet precession can originate from mechanisms such as the misalignment of the central engine spin axis with surrounding accretion disks by spin-induced disk precession (Bardeen & Petterson 1975; Caproni et al. 2004), from disk instabilities (Pringle 1997; Lai 2003), or from a SMBH companion (Begelman et al. 1980). The SMBHB scenario, through the phenomenon of jet precession, has been proposed as an interpretation of the PA periodicity for a large number of sources (Caproni & Abraham 2004; Britzen et al. 2018; Abraham 2018; Dey et al. 2021; Kun et al. 2023; Otero-Santos et al. 2023). The presence of a SMBHB has also widely been brought forward to explain flux periodicity found from radio to high energies, whether accompanied by PA changes or not (Li et al. 2022; Kun et al. 2022; Mao & Zhang 2024). Regarding 3C 66A, the two distinct precession timescales found, as well as the known optical flux periodicity, lead us to argue in favor of the SMBHB model. While it is clear that a SMBHB system can lead to periodic behavior, several mechanisms can explain the observed PA, optical and possible radio periodicities of 3C 66A when assuming a SMBHB system. We show in Section 5.1 that the short-term PA period can be associated with orbital motion in a SMBHB system. We then consider two scenarios to explain the long-term PA period: spin-orbit precession in Section 5.2; and disk-orbit precession in Section 5.3. We discuss the origin of 3C 66A’s flux variability in Section 5.4. We show that orbital motion and disk-orbit precession in a SMBHB can naturally account for all the observed properties of 3C 66A. We consider the implications regarding other detection methods for this SMBHB system in Section 5.5.

5.1. Orbital motion

In the case of the PA, the periodicity is associated with a precessing jet both for short-term (∼11 years) and long-term (∼10³ years) periods. The ∼11 years PA oscillation is only detectable up to approximately 80 pc from the core. On the other hand, a long-term PA shift is also measured for the furthest component (J3 or annulus C) at approximately 2.5 mas from the core, so at the deprojected distance of ∼150 pc. It is therefore clear the the decade-scale jet precession affects only a part of the jet close to the core, while the long-term precession appears to act on the full jet. These findings strongly support the fact that the short-term and long-term periodicities have different physical origins.

We argue that the short-term precession can be associated with the orbital motion of the primary BH in the binary system, while the long-term precession cannot. We can estimate the expected central SMBHB mass M₈, in units of 10⁸ M_⊙, given an orbital period P (Begelman et al. 1980): M₈ ≈ 18^−8/5P^8/5q^−3/5, and compare it to the expected mass. Here q = M₂/M₁ is the mass ratio of the system, with M₁ and M₂ the masses of the primary and the secondary BH, respectively. For 3C 66A, the primary BH mass M₁ can be computed from optical variability timescale measurements (e.g., Kaur et al. 2017; Fan et al. 2018) that give M₁ = (1.20 ± 0.05)×10⁸ M_⊙ (see Appendix D). We consider the typical range of mass ratios 30:1−3:1 derived by Gergely & Biermann (2009) for SMBH mergers. Taking q = 0.18 ± 0.15, we obtain the central BH mass M = (1.42 ± 0.19)×10⁸ M_⊙. If the long-term period of at least 500 years, corrected for redshift, corresponds to the binary’s orbital period, then M ≈ 2 × 10¹⁰ − 4 × 10¹¹ M_⊙ for q = 0.99 − 0.01. These values are two to three orders of magnitude above the expected mass of the central BH in 3C 66A. They also reach the maximum observable SMBHB mass estimated to be ∼5 × 10¹⁰ M_⊙ (King 2016). It therefore appears highly improbable that the observed long-term precession in 3C 66A translates orbital motion in the SMBHB system. On the other hand, we can assume that the decade-scale precession reveals the orbital motion of the system. The inherent short-term precession periodicity is P^PA = P_obs^PA/(1 + z) = 8.12 ± 0.18 years. The estimated SMBHB mass, following Begelman et al. (1980), then lies in the range M ≈ 3 × 10⁷ − 4 × 10⁸ M_⊙ for q = 0.99 − 0.01, which is in agreement with expectations for 3C 66A.

5.2. Spin–orbit precession

In this section, we interpret the short- and long-term periods as precession due to orbital motion versus precession induced by a misalignment between the primary BH spin and the orbital plane axis. This is done in analogy to the quasar S5 1928+738, which presents a similar periodic and monotonic trend in its inclination and PA evolution (Kun et al. 2014, 2023). We take the period of the precession model fit to the inner jet PA to be the observed orbital period, P_obs. The intrinsic binary period is therefore assumed to be P^PA = 8.12 ± 0.18 yr. From the equation of motion of a system with total mass M = M₁ + M₂, the binary separation can be estimated to Newtonian order (Kidder 1995):

$$\begin{array}{c}\hfill r_{\mathrm{sep}}={\left(\frac{P}{2\pi }\right)}^{2/3}{(GM)}^{1/3}\end{array}$$(8)

with G the gravitational constant. Given the central BH mass M = (1.42 ± 0.19)×10⁸ M_⊙, the binary separation between the two putative BHs is r_sep = (1.02 ± 0.05)×10⁻² pc =(2.10 ± 0.10)×10³ AU. In gravitational radii R_g = GM₁/c² = 5.7 × 10⁻⁶ pc, relative to the primary BH, the separation is r_sep = (1.78 ± 0.08)×10³R_g.

The post-Newtonian (PN) parameter, quantifying the validity of the Newtonian approximation, is given by

$$\begin{array}{c}\hfill ϵ=\frac{\mathit{GM}}{r_{\mathrm{sep}}{c}^2}.\end{array}$$(9)

We have ϵ = (6.64 ± 0.94)×10⁻⁴ which lies just outside the expected inspiral phase range 0.001 < ϵ < 0.1 (Gergely & Biermann 2009). The coalescence of such a SMBHB is predicted to occur within the gravitational (merger) timescale following

$$\begin{array}{c}\hfill T_{\mathrm{merge}}=\frac{5GM}{32c^3{ϵ}^4}\frac{{(1+q)}^2}{q}.\end{array}$$(10)

For the assumed mass ratio of 3C 66A, T_merge ≈ 136 Myr.

From the computed PN parameter, we can estimate the spin-orbit precession timescale of the system and establish the validity of this interpretation for the long-term re-orientation of the PA in 3C 66A. We consider that the precession is induced by the misalignment between the primary’s spin S₁ and the orbital motion angular momentum L. In the expected mass ratio range for merging SMBHs, the spin S₂ of the secondary BH can be neglected (Gergely & Biermann 2009). Ignoring 2PN effects, of second order 𝒪(ϵ²), the primary BH precesses around the total angular momentum vector J = S₁ + L with period (Kun et al. 2014):

$$\begin{array}{c}\hfill P_{\mathrm{SO}}=2\pi {\mathrm{\Omega }}_{\mathrm{SO}}^{-1}=2\pi \frac{\mathit{GM}}{2c^3}{ϵ}^{-5/2}\frac{{(1+q)}^2}{q}.\end{array}$$(11)

We obtained a spin-orbit precession period P_SO ≈ 47 000 years, which is largely above the expected period obtained from the fit of the slow PA re-orientation. Even considering the extreme case q = 0.99, the spin-orbit precession period remains an order of magnitude above our expectation. We therefore disfavor spin-orbit precession to account for the long-term jet precession in 3C 66A.

5.3. Disk-orbit precession

We show in Sections 5.1 and 5.2 that orbital motion can account for the decade-scale jet precession, but that spin-orbit precession is unlikely to drive the long-term PA change in 3C 66A. Ignoring the LT effect, discussed in Section 6.1, the only satisfying explanation for the long-term jet precession is that the SMBHB orbital plane is misaligned with the primary BH’s accretion disk. In this case, the ∼11 year period in the PA of 3C 66A is expected to result from a combination of orbital motion and disk-orbit precession (Katz et al. 1982; Bate et al. 2000). Indeed, in a binary system the gravitational torque from the secondary BH affects the disk of the primary and induces an additional oscillation in the accretion disk’s orientation, that depends both on the precession frequency from the disk-orbit misalignment and on the orbital frequency of the system. We therefore interpret our short-term precession as an effect of orbital motion and disk-orbit misalignment (also called “wobbling” or “nutation” in the literature), while the long-term precession results from disk-orbit misalignment alone. A similar case is argued for the double periodicities found in the PA evolution of M81 (Jiang et al. 2023) or of the famous SMBHB candidate OJ 287 (Katz 1997; Britzen et al. 2023).

Bate et al. (2000) found that the short-term precession, or “wobbling”, corresponds to approximately half of the binary period, while the long-term precession should be on the scale of 20 periods. This implies that our observed binary period would be of approximately 20 years and the long-term period of approximately 400 years, which is in agreement with the fit to the PA re-orientation of 3C 66A. We first check the possibility of disk-orbit precession in the case of 3C 66A more rigorously, by determining the expected rate of disk precession from a misalignment with the orbital plane (Katz 1997):

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{DO}}=-\frac{3}{4}\frac{G{M}_2}{r_{\mathrm{sep}}}{\left(\frac{r_{\mathrm{d}}}{r_{\mathrm{sep}}}\right)}^2\frac{cos({\theta }_{\mathrm{d}})}{{(G{M}_1{r}_{\mathrm{d}})}^{1/2}},\end{array}$$(12)

where r_d is the radius of the accretion disk and θ_d the angle between the disk and the orbital plane. Katz (1997) stresses that this corresponds to evaluated torques acting onto the disk and is not expected to match exactly the precession period observed. Nonetheless, the order of magnitude between both should be in agreement. For instance, in the case of the X-ray binary system SS 433, for which this model was applied before SMBHB systems, the derived timescale P_DO is approximately half of the observed precession timescale (Katz 1980). We assume that the accretion disk around the primary BH precesses as a rigid body, implying that r_d < r_sep, and take r_d = 0.5 × 10⁻² pc as roughly half the binary separation. For the disk inclination to the orbital plane, we choose θ_d = 20° in analogy to SS 433. From the separation r_sep = (1.02 ± 0.05)×10⁻² pc, obtained by assuming that the observed short-term period translates directly the orbital period, we get P_DO ≈ 200 years, which is of the right order of magnitude. We therefore consider disk-orbit precession to be a valid mechanism to explain the long-term period observed in 3C 66A.

In a system undergoing disk-orbit precession, we expected the orbital period observed to be modulated by the disk-orbit precession period, so that we must refine our SMBHB parameter estimation with respect to Section 5.2. We applied the results of Katz et al. (1982) to obtain a more accurate value of the observed orbital period of the system (Caproni et al. 2013; Britzen et al. 2023):

$$\begin{array}{c}\hfill P_{\mathrm{obs}}^{\mathrm{orb}}=P_{\mathrm{obs}}^{\mathrm{prec}}{(\frac{P_{\mathrm{obs}}^{\mathrm{prec}}}{2P_{\mathrm{obs}}}-1)}^{-1},\end{array}$$(13)

where P_obs^prec is the long-term precession period that we obtain from the slow PA re-orientation. We take P_obs^prec = 750 ± 250 years, with 750 the central estimate and 250 the uncertainty reflecting the plausible range [500 − 1000] years, obtained from the long-term period fit of the PA. Then, Eq. 13 yields P_obs^orb = 22.41 ± 0.56 years, equivalent to P^orb = 16.72 ± 0.42 years in the source rest frame. Using Eqs. 8 and 9, we determined, in this case, a larger binary separation of r_sep = (1.65 ± 0.08)×10⁻² pc =(3.4 ± 0.2)×10³ AU=(2.89 ± 0.14)×10³R_g and a smaller PN parameter of ϵ = (4.10 ± 0.58)×10⁻⁴. The latter lies again outside the expected inspiral phase range; at the current stage, the bound system is dominated by dynamical friction. From Eq. 10, we also estimate the expected merger timescale, T_merge ≈ 934 Myr. Given the new binary separation from the updated orbital period, Eq. 12 finally gives the disk-orbit precession timescale for the system: P_DO ≈ 846 years, which falls perfectly in the range obtained from the fitting of the long-term period.

Since disk-orbit precession implies the misalignment between the disk and orbital motion, we expect that either the primary BH spin is also misaligned with the orbital motion, or that it is misaligned with the accretion disk. The first would result in spin-orbit precession, which in this case would be negligible because occurring on a timescale of 10⁵ years. The second could lead to LT precession, discussed in Section 6.1.

The above estimates strongly support orbital motion and disk-orbit precession as the driving mechanisms for the observed short- and long-term precession timescales in 3C 66A. In the rest of Section 5, we therefore refer to the orbital parameters computed for this model as those of the putative SMBHB system. We summarize these parameters in Table 3.

5.4. Flux variability

The variability in the flux of 3C 66A, whether optical or radio, can be caused by geometrical effects in the precessing jet or by flaring from the accretion disk because of perturbed accretion rates due to the secondary BH’s influence on the primary’s disk. This is, for example, the most supported case to explain the flares observed in OJ 287 (Lehto & Valtonen 1996; Dey et al. 2018).

For 3C 66A, our favored scenario regarding radio and optical flux variability is Doppler boosted emission from the precessing jet. This is indeed one of our motivations to fit the PA change with a precession model in Section 4. The possible correlation between the radio core flux density and the inner jet PA supports this theory, as well as the periodicity timescales themselves (Rieger 2004). The optical period P_obs^optical = 2.38 ± 0.19 years (Cheng et al. 2022) of 3C 66A also motivates the jet-boosted emission scenario. Indeed, if we assume that the optical flux variability is dominated by the jet rather then the accretion disk and is the result of Doppler boosting, the intrinsic period is P^optical = δP_obs^optical/(1 + z) = 7.92 ± 3.31 years, given the Doppler factor δ = 4.46 ± 0.19 obtained from observations (see Appendix D). We find it interesting that the value of P^optical matches the intrinsic precession period P^PA resulting from the short-term PA oscillation.

The flux variability may also occur due to modulation of the disk’s accretion rate, leading to flaring events (Gold et al. 2014). In this case, the periodicity may or may not need to be corrected for Doppler boosting: if the flaring comes from an inner disk region subjected to relativistic effects, some boosting may occur, but if it comes from a standard Shakura-Sunyaev accretion disk region it is negligible (Shakura & Sunyaev 1973). In either case, the mass accretion rate can increase due to a SMBHB once, twice or more times per orbit (Hayasaki et al. 2013). Given our knowledge of the source, it is difficult to link the mass accretion variability to orbital parameters if we assume this origin for the optical flux periodicity. In the context of disk-orbit precession, considering that the optical flux originates from a standard, isotropic disk, the ∼2 year period would have to correspond to enhanced mass accretion some ten times per binary orbit, which appears unlikely.

We have shown that Doppler boosted emission from the jet, which precesses due to orbital motion, can explain the observed optical periodicity of 3C 66A. The correlation between the radio core flux density and inner jet PA change supports a similar origin for the radio variability. While accretion disk modulation may contribute to the flux variability in 3C 66A, we do not believe it is a dominant effect.

5.5. Implications regarding other detection methods

Inspiraling SMBHBs lose energy through GW emission that may be detected through techniques such as pulsar timing array (PTA) or the Laser Interferometer Space Antenna (LISA) in the upcoming decade (e.g., Bi et al. 2023; EPTA Collaboration 2024; Agazie et al. 2024; Sah et al. 2025; Amaro-Seoane et al. 2023; de Bruijn et al. 2020). Based on the SMBHB parameters obtained for 3C 66A, we can estimate the strain amplitude for this system. The computed central mass, M, and assumed ratio, q, yield the chirp mass, $M_{\mathrm{chirp}}=M^{9/5}{\left[\frac{q}{{(1+q)}^2}\right]}^{3/5}=(4.18\pm 1.52)\times {10}^7{M}_{\odot }$. The strain amplitude is then (Taylor et al. 2016; Feng et al. 2019):

$$\begin{array}{c}\hfill h=\frac{2(G{M}_{\mathrm{chirp}}^{5/3}{(2\pi f)}^{2/3}}{c^4{D}_{\mathrm{L}}},\end{array}$$(14)

where f = 1/P is the orbital frequency in the source rest frame and D_L = 1847.9 Mpc is the luminosity distance. This distance is computed with the CosmoCalc³ tool (Wright 2006). The resulting strain amplitude is of the order of 10⁻¹⁹. Given the frequency and strain estimated for the SMBHB in 3C 66A, no current PTA collaboration has the required sensibility to detect GW from this system (see Fig. 11 in D’Orazio & Charisi 2026). Nonetheless, we encourage PTA investigation of 3C 66A in order to place stringent upper limits on the chirp mass of the putative SMBHB system.

In addition to GWs, detecting the orbital motion of SMBHBs is a very promising method to directly test the presence of these systems at the center of galaxies. The continuous increase in sensitivity and resolution of VLBI arrays and subsequent data analysis prone us to believe that firm detections of SMBHB motions will soon be possible (An et al. 2018a; Zhao et al. 2024; Gurvits et al. 2025). 3C 66A is a bright source and the orbital period of its putative SMBHB is of the order of decades, making it a serious candidate for orbital motion detection with VLBI. Interestingly, 3C 66A is very close in the sky to the radio galaxy 3C 66B (z = 0.0213), which is a famous candidate SMBHB due to a 1-year core orbital motion found by Sudou et al. (2003). This motion was obtained using 3C 66A as the stationary position reference after observing both sources simultaneously in each beam at 2.3 and 8.4 GHz. We note that there is a non-negligible possibility that the motion attributed to the core of 3C 66B is at least partially the result of core motion in 3C 66A. In particular, at lower frequencies the observed core is most likely blended with the inner jet, which was found at higher frequencies in this work. Although Iguchi et al. (2010) further support the SMBHB hypothesis for 3C 66B with a periodicity found at mm wavelengths, no other evidence of a SMBHB has been found. Moreover, PTA measurements refuted the presence of a SMBHB in 3C 66B for the parameter range derived by Sudou et al. (2003) (Jenet et al. 2004). It should be noted, however, that later PTA studies put an upper limit on the chirp mass of the system that is consistent with the one refined by Iguchi et al. (2010) (Feng et al. 2019; Arzoumanian et al. 2020, 2023; Agazie et al. 2024; Tian et al. 2025). If the detected orbital motion corresponds to the core of 3C 66B, this near-merger SMBHB would have an extremely short merger timescale (of the order of years to thousands of years) compared to the typical life-time of these systems, meaning that we caught it at a very rare moment. We believe it possible that 3C 66B hosts a SMBHB, but the strong evidence found for the candidate 3C 66A motivates further core motion observations of each source using different references.

6. Discussion

We have interpreted the twisted jet, PA change, and flux periodicity in 3C 66A with jet precession caused by a SMBHB, since this scenario can easily explain all of the observed features. However, other scenarios could explain some of the properties listed in Section 3. Here, we review additional mechanisms that could account for at least part of the features observed in 3C 66A.

6.1. Lense-Thirring disk precession

The LT effect corresponds to frame dragging inducing the precession of particles orbiting in misalignment with respect to the equatorial plane of a Kerr BH (Lense & Thirring 1918). Equivalently, in a single Kerr BH system where the viscous accretion disk is misaligned with the spin axis of the BH, frame-dragging effects combined with the internal viscosity lead to precession of the accretion disk (also known as Bardeen-Petterson effect) that can be modeled as a rigid body (Bardeen & Petterson 1975; Nelson & Papaloizou 2000). This precession can be extended to the entire disk-jet system, as shown by 3D GRMHD (general relativistic MHD) simulations (Liska et al. 2018). Through the Bardeen-Petterson effect, the accretion disk axis and the BH spin are progressively forced to align. The LT effect is claimed to be responsible for the jet precession observed in a number sources, such as 3C 454.3 (Qian et al. 2014) or M87 (Cui et al. 2023). We investigate whether the LT effect leading to the precession of the disk may cause the observed short and or long-term jet precession in 3C 66A.

We consider a BH with dimensionless parameter −1 ≤ a ≤ 1 such that its angular momentum is |J|=GM²|a|/c. In this part of our study, we take the single BH mass M to be the primary BH mass of M₁ = (1.20 ± 0.05)×10⁸ M_⊙ derived in Section 5 from variability timescales. Following the approach by Caproni et al. (2004), the inner radius of the accretion disk is assumed to be the radius of the marginally stable orbit, expressed as R_ms = ξ_msR_g, where R_g = GM/c² is the gravitational radius and the dimensionless variable, ξ_ms, is

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{ms}}& =3+A_2\mp {[(3-A_1)(3+A_1+2A_2)]}^{1/2},\text{where}\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill A_1& =1+{(1-a^2)}^{1/3}[{(1+a)}^{1/3}+{(1-a)}^{1/3}]\text{and}\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill A_2& ={(3a^2+A_1^2)}^{1/2}.\hfill \end{array}$$(17)

The ∓ sign in the ξ_ms definition is − for prograde motion and + for retrograde motion. The outer radius of the precessing disk is dimensionlessly given by ξ_out = R_out/R_g, in which case some choice on the value of this ratio must be made. The viscosity in the disk is modeled through the disk surface density, that we take to be either a power-law Σ_d(ξ) = Σ₀ξ^S (Nelson & Papaloizou 2000) or an exponential Σ_d(ξ) = Σ₀e^σξ (D’Angelo et al. 2003), where Σ₀, s and σ are constants.

The precession period of the rigid disk is then

$$\begin{array}{c}\hfill P_{\mathrm{LT}}^{\mathrm{disk}}=\frac{2\pi GM}{c^3}\frac{{\int }_{{\xi }_{\mathrm{ms}}}^{{\xi }_{\mathrm{out}}}{\mathrm{\Sigma }}_{\mathrm{d}}(\xi ){[\mathrm{{\rm Y}}(\xi )]}^{-1}{\xi }^{3}d\xi }{{\int }_{{\xi }_{\mathrm{ms}}}^{{\xi }_{\mathrm{out}}}{\mathrm{\Sigma }}_{\mathrm{d}}(\xi )\mathrm{\Psi }(\xi ){[\mathrm{{\rm Y}}(\xi )]}^{-2}{\xi }^{3}d\xi },\end{array}$$(18)

with Υ(ξ) = ξ^3/2 + a and Ψ(ξ) = 1 − (1 − 4qξ^−3/2 + 3a²ξ⁻²)^1/2. In the case of a power-law surface density, we vary the power-law slope from s = −2 to s = 1, fixing Σ₀ = 1. The resulting period as a function of a is given in Figure 6 (upper panels). For the exponential surface density case, we present the results for σ = −2 to σ = 0.5 (lower panels). We have explored the range from R_out = 10R_ms to R_out = 1000R_ms, but only show the most relevant curves corresponding to the cases R_out = 10R_ms (left panels) and R_out = 100R_ms (right panels).

Fig. 6. Period of LT disk precession as a function of the spin parameter a, outer radius Rout and density power law slope s (upper panel) or density exponential factor σ (lower panel). The black dashed line corresponds to PLTdisk = 10 years and the gray one to PLTdisk = 1000 years.

First, we consider the short-term PA oscillation of ∼11 years. If we assume a power-law surface density, the right timescale for the LT precession can be achieved in a disk where the precessing ring has a small outer radius (R_out = 10R_ms) if the BH spin is −0.2 < a < 0.15 approximately, for an expected disk surface density power law slope lower than 1. If the surface density must decrease with radius, then for a disk with such a small precessing ring the central BH must have a spin |a|< 0.05 approximately. This allowed spin range decreases quickly with s, so that in the case of a steep power-law s = −2, the BH spin must be almost null. Conversely, if we allow for a larger part of the disk to precess (R_out = 100R_ms or more) then an ∼11 year LT precession period can be obtained for a wide range of negative values of s and non-null BH spin, which better agrees with our expectations for an accretion disk. If we assume an exponential function for the surface density, it seems less likely that the LT disk precession can account for the short-term PA periodicity. Indeed, considering σ < 0 implies that the BH spin is practically null in any disk configuration, whether R_out is small or large. However, if the disk surface density satisfies an exponential increase with radius, then the right timescale can be achieved for R_out ≈ 10R_ms and a > 0.

Second, we investigate the ∼1000 year period that we expect from the long-term PA re-orientation. For a power-law surface density, whether the slope s is negative or positive, the precessing part of the disk must have an outer radius larger than approximately 30R_ms in order to have a spinning BH with |a|≥0.1. To have a decreasing surface density in the disk, for instance s = −1, R_out must reach the order of 100R_ms. The value of R_out must be even larger for a steeper power-law profile. In the case where the surface density follows an exponential function, even with R_out > 100R_ms the period of ∼1000 years can only be obtained for positive values of σ.

We have shown that LT precession of the disk can account for the observed PA oscillation on short or on long timescales. In addition, there is a wider acceptable parameter space to match the observed periods when assuming that the disk surface density is governed by a power-law function. We note, however, that the ∼11 and ∼1000 year oscillations can only be reproduced for distinct assumptions regarding the size of the precessing ring in the disk and the BH spin. Therefore, if 3C 66A hosts a single SMBH, then LT precession of the disk can either be the cause of the fast PA oscillation or of its long-term re-orientation, but cannot account for both.

6.2. Disk-induced BH precession

Next, we considered how the disk axis to spin misalignment may induce precession of the BH spin axis. If a massive, misaligned outer accretion disk exerts sufficient torque, it can drive both the inner disk region and the BH spin axis to precess. Similar to the case of LT disk precession, this may lead the entire disk-jet system to precess and explain the observed variations in the jet orientation. The timescale for such a configuration is of the order of thousands of years or more. Given the dimensionless viscosity parameter, α_vis (Shakura & Sunyaev 1973), and the mass accretion rate, Ṁ, we can compute the expected precession period for 3C 66A (Lu 1992):

$$\begin{array}{c}\hfill P^{\mathrm{BH}}={10}^{9.25}{a}^{5/7}{\alpha }_{\mathrm{vis}}^{48/35}{\left(\frac{M}{{10}^8{M}_{\odot }}\right)}^{1/7}{\left(\frac{\dot{M}}{{10}^{-2}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{-6/5}\mathrm{yr}.\end{array}$$(19)

Assuming α_vis = 0.1 (Armitage 1998; King et al. 2007), we show in Figure 7 how the period changes with the accretion rate and spin. We first consider the Eddington accretion rate Ṁ_Edd from the estimated Eddington luminosity L_Edd of 3C 66A:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{Edd}}=\frac{L_{\mathrm{Edd}}}{\eta c^2}=\frac{4\pi GM{m}_{\mathrm{P}}}{\eta c{\sigma }_{\mathrm{T}}},\end{array}$$(20)

Fig. 7. 2D plot of the BH spin axis precession period as a function of the spin parameter a and the accretion rate Ṁ for αvis = 0.1. The black dashed line corresponds to PBH = 10 years and the gray one to PBH = 1000 years. The area in between the vertical dotted white lines corresponds to the typical range of accretion rates for BL Lac objects.

where m_P is the proton mass, σ_T is the Thomson cross-section and η translates the efficiency in converting gravitational into radiative energy. We ignore general relativistic effects and take the typical value η = 0.1 expected for a moderately spinning BH (Shakura & Sunyaev 1973). We obtain Ṁ_Edd ≈ 3 M_⊙/yr. BL Lac objects such as 3C 66A are expected to be powered by radiatively inefficient, hot, geometrically thin and optically thick disks, known as ADAF (advection dominated accretion flow; Narayan & Yi 1994). Population studies have shown that BL Lac objects are highly sub-Eddington systems typically accreting with a rate Ṁ < 0.01 Ṁ_Edd (Wang et al. 2002; Ghisellini et al. 2010). We therefore expect 3C 66A to accrete with a rate of, at the very most, Ṁ ≈ 0.1 M_⊙/yr. In Figure 7, the expected range of Ṁ = 10⁻³−10⁻¹ M_⊙/yr for BL Lac objects corresponds to the area in between the vertical dotted white lines.

For the maximal value Ṁ = 0.1 M_⊙/yr, we obtain a precession period of approximately 10⁶ years when a = 0.1. This timescale is orders of magnitude above the ∼11 year period of the inner jet PA. The BH precession could only account for the short-term PA oscillation if the accretion rate were much larger than expected for 3C 66A, or if the BH spin were of the order of 10⁻⁸. Even when decreasing the disk viscosity to a reasonable value of 10⁻², the spin parameter would have to be smaller than 10⁻⁶. It is therefore highly unlikely that precession of the BH can explain the short PA periodicity. A similar conclusion is drawn by Caproni & Abraham (2004) or Britzen et al. (2018) regarding the observed PA oscillations of 3C 345 and OJ 287, respectively, that are also of the order of 10 years. Regarding the long-term PA change occurring over approximately 1000 years, the acceptable parameter range for α_vis, a and Ṁ is more reasonable. In the best scenario, the spin parameter would be approximately 10⁻⁵ which remains very small. While this interpretation for the long-term PA periodicity is possible, we do not favor it due to the very low probability of a single SMBH at the center of 3C 66A having such a negligible spin.

Therefore, the misalignment between the accretion disk and the BH spin axis can explain one of the observed PA periods in 3C 66A only through LT precession of the disk, not through BH spin axis precession. In addition, LT effects alone cannot account for both the short-term jet precession and the long-term one. Additional mechanisms are needed to explain all of the observed features in 3C 66A.

6.3. Jet instabilities

A broad array of instability types are known to potentially occur in relativistic jets, thanks both to early analytical investigation and to the fast development over recent decades of (R)MHD and PIC (particle-in-cell) simulations (Nishikawa et al. 2019; Ortuño-Macías et al. 2022). One of the most discussed instabilities to explain changes in AGN jets is the current-driven kink instability (CDI) (McKinney 2006; Mizuno et al. 2009, 2014; Singh et al. 2016). In the presence of strong toroidal magnetic fields it leads to helical jet distortion. Recently, Jorstad et al. (2022) interpret a ∼13 hours quasi-periodic oscillation in the optical and γ-ray flux of BL Lacertae as CDI propagating near a recollimation shock, at approximately 5 pc from the BH core. CDI is excited by strong magnetic fields (Hu et al. 2025), which is often an argument in favor of this instability to explain jet variability near the base. We verify whether CDI could explain the short-term PA periodic pattern in the jet of 3C 66A. The growth timescale of CDI (m = 1) kink-modes is governed by the crossing time τ_cross = R_j/v_A where R_j is the jet radius and v_A the Alfvén speed. We take R_j ≈ 0.1 mas from the average FWHM of the fitted inner jet circular components and consider v_A ≈ 0.8c for the relativistic jet. The resulting crossing time is τ_cross ≈ 2 years. Mizuno et al. (2009) show that CDI affects the jet structure over a characteristic timescale of ∼100R_j/v_A and during this linear phase the instability grows. To understand whether CDI affects a relativistic jet, we can check whether it satisfies the condition

$$\begin{array}{c}\hfill \frac{\mathcal{D}}{\mathrm{\Gamma }c}>100\frac{R_{\mathrm{j}}}{v_{\mathrm{A}}},\end{array}$$(21)

with Γ the Lorentz factor and 𝒟 the jet distance over which we expect the instability to grow. We take Γ ≈ 20 from the study of 5 jet knots in 3C 66A by Weaver et al. (2022). For the length-scale 𝒟, we choose the distance of 18 pc which corresponds to the average at which the PA of the inner jet component is measured. If CDI is responsible for the periodic pattern observed in the inner jet, then it should have been able to grow significantly by this distance. The left-hand side of eq. 21 yields 𝒟/Γc ≈ 1 × 10³ years, while 100τ_cross ≈ 2 × 10² years. The proper propagation time is approximately one order of magnitude larger than 100 crossing times. It is therefore possible for the inner jet region to be significantly affected by CDI. We can moreover compare the CDI growing timescale to the travel time of the presumed wave in 3C 66A, which is manifested by the PA oscillation shift with distance from the core.

As done in Section 4.3, taking Δτ_shift ≈ 1.5 years between the inner and downstream jet regions, the wave appears to propagate with the apparent speed β_wave^app ≈ 5. Between 0.1 and 1 mas from the core (limits of the inner and downstream jet regions), the travel time of the wave is therefore approximately 3 years, which is larger than the crossing time. Ours is an approximate estimate due to the fact that we can only identify a periodic pattern for two regions in the jet, and their distance from the core falls in a wide ring. Nonetheless, the properties of the presumed wave in the jet of 3C 66A agree with the CDI scenario.

The second instability of interest is the Kelvin-Helmholtz instability (KHI), which arises at the interface between the jet and the external medium where the velocity shear can lead to turbulence and wave-like disturbances (Hardee 1987; Savolainen et al. 2006; An et al. 2010; Perucho et al. 2012; Nishikawa et al. 2019). It typically occurs over timescales of years but is suppressed by high magnetic fields, typically ruling out the possibility of its development at the jet base. For 3C 66A, the inner jet PA probes the orientation of the inner 10–25 pc, while the downstream jet PA probes the 25–60 pc approximately. At these distances from the central engine, it is possible that the magnetic field is low enough for the KHI to propagate. We draw a comparison with M87, for which Walker et al. (2018) identify a quasi-period of ∼8-10 years in the jet offset with respect to a central line. While this does not correspond to the PA change, it equivalently describes the jet orientation evolution since the central line is fixed along a given PA. Walker et al. (2018) find that the quasi-periodicity in the jet offset propagates down the jet, meaning that its phase is progressively shifted when measuring it at further distances from the core. The quasi-periodicity is observed between 2 and 8 mas from the core, which goes up to approximately 2.5 pc in deprojected distance. They interpret this transverse jet position shift as KHI propagation, showing that a KH helical mode can explain the observed pattern motion for reasonable physical parameters of the jet. Matveyenko & Seleznev (2015) also argue that a hydrodynamical instability in the flow ejection is responsible for precession in the jet of M87. However, an equivalent PA periodicity to that found by Walker et al. (2018), of approximately 11 years, was observed at the base of M87’s jet using a larger dataset by Cui et al. (2023). These authors discarded the KHI hypothesis due to its suppression by strong magnetic fields and argue for LT precession instead. This manifests how the interpretation of periodic signals in jet orientation evolution is still much debated, even for extensively studied sources such as M87. Recently, Kino et al. (2025) also considered the possibility of the presence of a SMBHB in M87, based on the ∼11 year precession period and a ∼0.9 year transverse period.

In the case of 3C 66A, the ∼11 year period observed corresponds to an angular frequency ω_j ≈ 2 × 10⁻⁸ rad/s. Walker et al. (2018) showed that a lower limit to the helical mode resonant frequency can be expressed as

$$\begin{array}{c}\hfill {\omega }_{\ast }\approx \frac{3\pi }{4}\frac{c/R_{\mathrm{j}}}{(c/v_{\mathrm{w}}-1)},\end{array}$$(22)

where v_w represents the speed of sound if we consider the sonic limit (low magnetization) or the Alfvén speed if we consider the Alfvénic limit (high magnetization) for the jet. In order to obtain this expression, the authors assume that the external medium satisfies the no flow limit, implying that its speed u_e = 0. In addition, one must assume that the sound and Alfvén speeds are the same inside and outside of the jet, which can be motivated by the presence of a cocooning medium surrounding it. In this case we take R_j ≈ 0.4 mas, from the average FWHM of the fitted downstream jet circular components, since the condition ω_* > ω must also be satisfied in the downstream jet if we interpret the PA periodicity of this region as a KHI signature. Then, for the resonant frequency ω_* to exceed the observed angular frequency ω, we must have v_w ≳ 0.6c. Below this value, we have ω > ω_* and the instability growth rate is vanishingly small (Hardee 2007). The condition on v_w rules out the sonic limit interpretation since the maximum speed of sound is $c/\sqrt{3}$. The Alfvén speed can theoretically satisfy v_A ≳ 0.6c in the jet, although it may not be this large in 3C 66A for the jet regions considered, due to the distance of more than 10 pc from the core implying a moderate magnetization. If we assume that the jet is in the Alfvénic limit and the magnetic pressure is balanced between the jet and the surrounding medium, the KH helical mode can propagate near the resonant frequency. In this case, the helical mode wave speed is (Hardee 2007):

$$\begin{array}{c}\hfill v_{\mathrm{h}}^{\ast }=\frac{{\mathrm{\Gamma }}_{\mathrm{j}}({\mathrm{\Gamma }}_{\mathrm{we}}{v}_{\mathrm{we}})u_{\mathrm{j}}+{\mathrm{\Gamma }}_{\mathrm{e}}({\mathrm{\Gamma }}_{\mathrm{wj}}{v}_{\mathrm{wj}})u_{\mathrm{e}}}{{\mathrm{\Gamma }}_{\mathrm{j}}({\mathrm{\Gamma }}_{\mathrm{we}}{v}_{\mathrm{we}})+{\mathrm{\Gamma }}_{\mathrm{e}}({\mathrm{\Gamma }}_{\mathrm{wj}}{v}_{\mathrm{wj}})},\end{array}$$(23)

where the j and e subscripts indicate parameters of the jet and of the external medium, respectively, u_j is the jet flow speed and Γ_w = (1 − v_w²/c²)^−1/2 is the sonic or Alfvénic Lorentz factor. Here, we let u_e be positive. The jet is in the Alfvénic limit and we assumed that the Alfvén speed is equal in the jet and external medium, leading to: v_wj = v_we = v_A and Γ_wj = Γ_we = Γ_A. The KH helical mode speed from eq. 23 thus simplifies into:

$$\begin{array}{c}\hfill v_{\mathrm{h}}^{\ast }=\frac{{\mathrm{\Gamma }}_{\mathrm{j}}{u}_{\mathrm{j}}+{\mathrm{\Gamma }}_{\mathrm{e}}{u}_{\mathrm{e}}}{{\mathrm{\Gamma }}_{\mathrm{j}}+{\mathrm{\Gamma }}_{\mathrm{e}}}.\end{array}$$(24)

We estimate a lower and upper limit for this speed using the Lorentz factor and apparent component speeds reported by Weaver et al. (2022). For the jet parameters, we take the average values Γ_j ≈ 20 and β_j^app ≈ 18 from the components studied. Assuming the viewing angle is Θ ≈ 5°, the jet speed is u_j ≈ 0.999c. The lower limit of v_h^* is obtained when assuming u_e = 0 and Γ_e = 1. The upper limit is obtained by considering that the external medium properties correspond to the smallest propagation speed observed in the jet and to the smallest estimated Lorentz factor: u_e ≈ 0.997c and Γ_e ≈ 14. The KH helical mode speed is then in the range 0.951c < v_h^* < 0.998c. We can compare this speed to that of the presumed wave in 3C 66A. The latter is given by β_wave^app ≈ 5, equivalent to the propagation speed v_wave ≈ 0.987c. KHI propagation can therefore explain the observed phase-shift in the PA periodic pattern between the inner and downstream jet. In addition, strong magnetic fields can lead to the damping of the KHI as it propagates along the jet (Mizuno et al. 2009), which is in agreement with the observed amplitude reduction of the periodic pattern further down the jet. Once again, this would imply that the jet is strongly magnetized at tens of pc from the core. In addition, in the case of M87 Walker et al. (2018) note that there is no sign of the quasi-period pattern at arcseconds scales until knot A. A similar property is observed for 3C 66A, since we do not observe the periodic pattern starting at approximately 1.3 mas from the core, equivalent to ∼80 pc deprojected. Therefore, the KHI can explain the observed PA oscillation in 3C 66A. It requires a jet in the Alfvénic limit and an equivalent magnetization of the jet and of the immediate external medium.

Another instability known to occur in jets is the Rayleigh-Taylor instability (RTI), that can arise at the interface between the accelerated jet and the ambient medium due to density gradients, growing into what is often dubbed finger-like structures. Recently, Hu et al. (2025) showed through 3D and 2D RMHD simulations that this instability is stabilized by the presence of a magnetic field. For a helical magnetic field they observe the CDI to develop instead, after passing recollimation shocks, supporting the findings by Jorstad et al. (2022). In the case of 3C 66A, RTI may play a role in the complex jet morphology observed, but unlikely explains any of the periodicities observed. A similar conclusion can be drawn regarding the mushroom instability (MI) which is caused by electron-scale velocity shears and leads to mushroom-like structures in 2D simulations, perpendicular to the jet propagation direction. PIC simulations indicate that this instability grows stronger with stronger magnetic fields (Nishikawa et al. 2019). It has also been associated with particle acceleration and possible magnetic reconnection (Kawashima et al. 2022). It is likely that MI acts in the 3C 66A jet and causes the propagation of disturbances as well as flux enhancement, but not on periodic timescales. Other internal instabilities that may occur in the jet, but on timescales too short (∼ seconds to days) to explain our flux or PA oscillations, are the two-stream instability (Horký & Kulhánek 2013), the Weibel instability (Nishikawa et al. 2019) and the pressure-driven instability (Longaretti 2008; Das & Begelman 2019).

Therefore, many internal instabilities can occur in relativistic jets, with the CDI and KHI able to explain one of the PA periods of 3C 66A. We have shown that the CD kink instability or the KH helical mode can account for the the wave-like propagating PA pattern of ∼11 years. For the latter, this requires that the jet is in the Alfvénic limit and its magnetization is equal to that of the surrounding medium. The instability would propagate while being damped by the high magnetic field and disappear beyond ∼80 pc from the core. This finding leads us to consider an alternative interpretation to the SMBHB system, presented in Section 5. We can indeed imagine a combination of LT precession of the disk, explaining the PA monotonic decrease through secular jet precession, and of either CDI or KHI propagating at the jet base, leading to the PA oscillation on decade-scales.

6.4. Disk instabilities and dynamics

On top of processes occurring in the jet, we investigate disk instabilities that may lead to flux variability and or to precession (Armitage & Reynolds 2003; King et al. 2004). For example, Ferreira et al. (2022) explore instabilities in the disk to explain low frequency quasi-periodic oscillations in black hole X-ray binaries, which could also lead to jet wobbling in accreting sources.

Except for the LT effect discussed in Section 6.1, a relevant disk instability is radiation-induced warping of the accretion disk, which could lead to jet precession and hence explain both PA and flux variability (Pringle 1997; Lai 2003). An accretion disk irradiated by a central source (which might in fact be the central disk regions) is unstable to becoming nonplanar: if the disk is warped, light hits its surface unevenly and, in the optically thick case, it causes the re-emitted radiation to exert uneven torques on different parts of the disk. These lead to further warping over time, effectively acting as a precession. For example, such a warped disk has been considered to explain jet precession at the center of the Seyfert galaxy NGC 1275 (Dunn et al. 2006). Assuming an optically thick disk, Pringle (1997) estimates the characteristic disk-warping timescale to be

$$\begin{array}{c}\hfill P_{\mathrm{warp}}\approx 2\times {10}^6{\alpha }_{\mathrm{vis}}^{-1}\frac{M}{{10}^8{M}_{\odot }}\text{yr}.\end{array}$$(25)

Choosing again the typical value of α_vis = 0.1, we obtain P_warp ≈ 10⁷ years, which is several orders of magnitude larger than both the short- and long-term periods observed for 3C 66A. We therefore rule out disk warping as a mechanism explaining jet precession for our system.

Other well-known disk instabilities that could account for flux variability include thermal-viscous instabilities. These broadly include the ionization instability and the radiation pressure instability. They have been mostly studied in the context of cataclysmic variables or X-ray binaries, but some authors have shown interest in their potential impact on AGN disks (Hameury et al. 2009; Janiuk & Czerny 2011). Ionization instability is caused by hydrogen transitions in the outer disk and may lead to quasi-periodic emission variability (which is well-established for X-ray binaries) by modulating the accretion rate. For a BH of mass M ∼ 10⁸ M_⊙, the ionization instability timescale is expected to be of the order of ∼10⁶ years or more (Janiuk & Czerny 2011). Radiation pressure instability, for which there is little observational background, occurs in the radiation pressure dominated inner disk and may also lead to accretion modulation, with a characteristic timescale of ∼10⁶ years between outbursts for a M ∼ 10⁸ M_⊙ BH. This may be reduced to ∼10³ − 10⁴ years in the case of high accretion rates (Czerny et al. 2009; Janiuk & Czerny 2011), although we do not expect such rates for a BL Lac object like 3C 66A (Ghisellini et al. 2010). The timescales for both of these instabilities remain far too large to explain the (possible) periodicity of 3C 66A in the (radio) optical LCs. While Godfrey et al. (2012) consider that radiation pressure instability may cause the observed periodic structure in the jet of PKS 0637−752 through accretion modulation, we disfavor this mechanism in 3C 66A. We note that the presence of a SMBHB is another strong scenario that Godfrey et al. (2012) propose for PKS 0637−752.

To explain the observed flux variability in 3C 66A, especially the optical period of ∼2 years, we can moreover consider a hot-spot rotating in the accretion disk at some radius, r_hs (Gupta et al. 2009; Prince et al. 2023). Assuming the hot-spot is the source of optical flux with a period of P_obs^optical = 2.38 ± 0.19 years, we can estimate the central BH mass as

$$\begin{array}{c}\hfill M=\frac{3.24\times {10}^4{P}_{\mathrm{obs}}^{\mathrm{optical}}}{(r_{\mathrm{hs}}^{3/2}+a)(1+z)}{M}_{\odot }.\end{array}$$(26)

We consider the Schwarzschild BH scenario, in which r_hs = 6 (in units of GM/c²) and a = 0, and the maximal Kerr BH scenario, in which r_hs = 1.2 and a = 0.9982. In the first case, the resulting BH mass is M = 1.23 × 10¹¹ M_⊙ and in the second, it is M = 7.83 × 10¹¹ M_⊙. Both mass values are too large by several orders of magnitude to be in agreement with the expected system in 3C 66A. We therefore rule out the presence of a hot-spot in the disk to explain the optical periodicity. The same conclusion can be drawn for the radio flux variability. This hot-spot scenario is similarly discarded by Prince et al. (2023) and Mao & Zhang (2024) regarding the periodicity features of ∼100 days and ∼600 days found in the LCs of the blazars PKS 0346-27 and PKS J2156-0037, respectively. Hence, we find that none of the disk phenomena explored can be applied to 3C 66A to explain the observed PA or flux periodicity.

6.5. External factors

Alternative scenarios for the observed bending in jets, such as the twisted structure found for 3C 66A, have been proposed. For instance, bent jet structures might occur because of external causes. Homan et al. (2002) argue that the bending structure of the quasar PKS 1510−089 would possibly be due to a density gradient between the host galaxy and intergalactic medium or a ram pressure from the intracluster medium. However, this scenario is unlikely in the parsec-scale jet bending because it is generally related to much larger-scale jet bending (Miley et al. 1972; Jaffe & Perola 1973; Jones & Owen 1979). Another possible explanation which would arise within the host galaxy is a colliding jet with surrounding clouds (Gómez et al. 2000; Gabuzda et al. 2001; Homan et al. 2003). Homan et al. (2003) find a bent trajectory with increasing flux density and decreasing size for the C4 component in the quasar 3C 279, indicating that C4 has been deflected by the interstellar medium. It is less likely responsible for the whole twisted structure of 3C 66A because the deflecting jet model requires a number of well-distributed clouds near the jet for the multiple bendings observed.

7. Conclusions

In this paper, we present three major features of the 3C 66A jet based on KaVA observations in 2014 and VLBA observations over 29 years. We find a twisted jet structure on the KaVA maps and, from the VLBA archival data, for the first time, we show that the inner jet PA varies periodically. We observe a periodic pattern of ∼11 years in both the inner jet and the downstream jet, with a ∼1.5 year phase-shift and an amplitude reduction further from the core. The PA oscillation is accompanied by a monotonic decrease of about −0.8°/year. Moreover, we find strong variability in the radio core flux with possible periods of ∼10 years, ∼7 years, and ∼4 years. In addition, we detect a long-term decrease of approximately −9.4 mJy/year in the baseline. The core flux presents a possible correlation with the periodic inner jet PA swing.

These findings, together with the well-known ∼2 year period of the optical LC in 3C 66A and with multi-wavelength baseline drops coincident with the jet’s PA decrease, suggest a geometric scenario for both the PA oscillation and its monotonic decrease. We interpret these observations as the signature of a precessing jet over a decadal timescale superimposed on a longer-term secular drift. The short-term nozzle precession may excite MHD waves such as Alfvén transverse waves that propagate down the jet, possibly explaining the downstream phase-shift and amplitude attenuation of the oscillation pattern.

The origin of the precession timescales is a topic subject to some debate on its own. We propose that a SMBHB system is the driver behind both precession timescales, with orbital motion and disk-orbit precession accounting for the short- and long-term trends respectively. Such a system can naturally explain all of the observables reported in this paper. The putative binary system orbits with a rest-frame period of approximately 17 years and the separation between the BHs is 1.65 × 10⁻² pc. The SMBHB is expected to merge in ∼934 Myr, currently emitting GWs below the sensitivity of PTA experiments.

Recently, many SMBHB candidates have been reported and some of them have been deemed controversial. We adopt this scenario as a parsimonious explanation, providing a simple interpretation for the large number of complex properties of 3C 66A, in particular, its different periodic timescales regarding flux and orientation. However, we have investigated a wide range of other possible interpretations. We show that LT precession of the accretion disk can reproduce either the short- or long-term precession timescale, especially for a power-law disk density profile which allows a wider range of spin parameters. We also evaluated the potential of assuming jet instabilities as an origin of the observed variability. The CDI kink mode and the KH helical mode are both possible interpretations for the ∼11 year periodic pattern, assuming that it corresponds to a propagating wave. We further explored disk-related mechanisms, such as radiation-induced disk warping, thermal-viscous instabilities, and a rotating hot-spot in the disk. None of these scenarios can account for the PA or the flux variability observed in 3C 66A. Therefore, excluding the SMBHB model, only one alternative appears viable: LT precession of the accretion disk drives long-term jet precession, while the observed decadal PA oscillation arises from CDI or KHI propagating in a magnetized jet. In the case of KH helical modes, the jet must be remain in the Alfvénic limit out to ∼100 pc and we must assume a magnetized external medium with comparable properties to the jet. In the future, additional kinematic, polarimetric, and astrometric analyses for 3C 66A will allow further investigation of the complex origins of its variability.

Acknowledgments

The authors thank Prof. Eduardo Ros, the MPIfR internal referee, for reading the manuscript and providing constructive comments that have improved the paper. The authors also express their gratitude to Dr. Hyunwook Ro, Dr. Andrei Lobanov and Dr. Wu Jiang for insightful discussions. This research was supported by funding from Korea AeroSpace Administration (KASA)) (grant numbers RS-2024-00509838 and R25TA0065942000). P.T. acknowledges support from France 2030 through the project named Académie Spatiale d’Île-de-France (https://academiespatiale.fr/) managed by the National Research Agency under the reference ANR-23-CMAS-0041. J. K. and B.W.S. are grateful for the support of the Korea Astronomy and Space Science Institute under the R & D program (Project No. 2017-1-840-90) supervised by the Ministry of Science and ICT. G.-Y.Z acknowledges support from the project M2FINDERS that is funded from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 101018682). S.-J.Y. acknowledges (a) support from the Mid-career Researcher Program (RS-2024-00344283) through Korea’s NRF funded by the Ministry of Science and ICT, and (b) support from the Basic Science Research Program (RS-2022-NR070872) through Korea’s NRF funded by the Ministry of Education. This study makes use of VLBA data from the VLBA-BU Blazar Monitoring Program (BEAM-ME and VLBA-BU-BLAZAR; http://www.bu.edu/blazars/BEAM-ME.html), funded by NASA through the Fermi Guest Investigator Program. The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. This study makes use of the following open-source packages: NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), Astropy (Astropy Collaboration 2022), emcee (Foreman-Mackey et al. 2013), libwwz, and psresp (Uttley et al. 2002).

References

Appendix A: PA evolution of three jet components

The evolution of the PA in the jet of 3C 66A is studied for three jet components with two distinct methods (annular bins and circular Gaussian fits). These jet components are found progressively further away from the core. The first two are denominated inner jet and downstream jet. Using the binning method to compute their PA, they correspond to the annuli (0.15 < r < 0.45) mas (annulus A), (0.45 < r < 1.00) mas (annulus B) and (1.50 < r < 3.00) mas (annulus C) respectively, where r is the radius from the center (0, 0). Their PA is also evaluated using circular and elliptical Gaussians from modelfit. In Figure A.3, we show the white circular Gaussian fits to the three jet components (denoted J1, J2 and J3) at six different epochs between 1999–2020, highlighting their different evolution. From Figure A.3, as well as Figures A.1 and A.2, it is clear, for instance, that the inner jet and the downstream jet’s PAs do not oscillate in phase.

In Figure A.1 we show the PA evolution for the annuli A, B and C. The data is visible on the top panel, while the corresponding fits are shown in the lower panel for visibility. Figure A.2 presents equivalent results but obtained with the circular Gaussian fitting approach (components J1, J2 and J3). The fit results of Eq.1 for both methods are given in Table E.1. The figures, as well as the fitted values, indicate an almost identical evolution of the PA for each considered component with the two methods.

Fig. A.1. Evolution of the PA for the annuli A (inner jet: black), B (downstream jet: blue), and C (purple) at 43 GHz, obtained with the annular bin method. The data is shown in the upper panel and the analytical fits in the lower panel: sinusoidal+linear fits for the annuli A and B and linear fit for the annulus C.

Fig. A.2. Evolution of the PA for the Gaussian components J1 (inner jet - maroon), J2 (downstream jet - orange) and J3 (brown) at 43 GHz, obtained with the circular Gaussian fitting method. The data is shown in the upper panel and the analytical fits in the lower panel: sinusoidal+linear fits for the J1 and J2 components and linear fit for the J3 component.

Fig. A.3. VLBA intensity maps at 43 GHz overplotted with three Gaussian fitted jet components (white) and core (black) at six epochs. In the leftmost panel, the labels J1, J2 and J3 denote the three jet components. The blue connected dots indicate the jet ridgeline. The pink contour levels start at three times the rms value, scaling twice. The gray ellipse at the bottom of each panel corresponds to the restoring beam.

Appendix B: Periodicity search methods

We describe the three approaches used to search for periodicity in the inner jet PA (from the annular bin results) and the core flux density. We detail the results obtained in each case and summarize the major periodicities found in Table B.1. All of the reported periods’ standard deviations are obtained from the FWHM of the periodicity peaks.

B.1. Jurkevic analysis

We first implement the Jurkevich method (Jurkevich 1971): the data is divided into m bins based on trial periods and, for each period, the total variance of all bins V_m² is computed. A minimum in the variance indicates a possible periodicity. This method allows to detect a periodicity even if it is not sinusoidal, but it doesn’t provide statistical confidence levels. In order to determine the significance of a minimum in the V_m² plot, we rely instead on the F-test suggested by Kidger et al. (1992):

$$\begin{array}{c}\hfill F=\frac{1-V_{\mathrm{m}}^2}{V_{\mathrm{m}}^2}.\end{array}$$(B.1)

If F ≥ 0.5, the period can be considered strong, while it is considered weak if F < 0.25.

We set the bins number to 15 and the number of trial periods to 500. The plot of V_m² as a function of period is shown in Figure B.1 for the PA, where we explore the period range between 1 and 20 years. There is a single significant drop in the V_m² values that reaches a minimum around 12 years, and gaussian fitting of the variance data leads to a period of 14.4 ± 5.1 years. While this is a very wide drop, it reaches clearly below the F = 0.5 threshold used to determine significance. Figure B.1 shows the results for the core flux density, in the period range between 1 and 15 years. This Jurkevic test appears less conclusive: there are a number of peaks, more or less well defined, that fall below F = 0.25 and two that reach F = 0.5. We report three peaks of interest from this analysis, considering the strength of their minima and the clarity of their peak shape. The corresponding periods are: 3.5 ± 0.3 years, 6.9 ± 0.2 years and 10.5 ± 0.6, all reaching the F = 0.5 level. These periods are also found with additional methods.

Fig. B.1. Jurkevic test plot for the inner jet PA. The black line corresponds to the Vm2 values whose minima indicate possible periodicities and the dotted blue (orange) line to the level F = 0.5 (F = 0.25). The vertical gray dashed line indicates the period at the minimum.

Fig. B.2. Jurkevic test plot for the core flux density. The black line corresponds to the Vm2 values whose minima indicate possible periodicities and the dotted blue (orange) line to the level F = 0.5 (F = 0.25). The vertical gray dashed lines indicate the periods at the minima of interest.

B.2. Lomb-Scargle periodogram

The LS periodogram (Lomb 1976; Scargle 1982) is a method based on sinusoidal least-squares fitting. This approach can only detect sinusoidal periodicities, but it is more robust to noise than the Jurkevich method. We use the LS periodogram from the AstroPy package. The significance levels account for red-noise, using the method described in Section B.4.

In the PA search we explore the frequency range corresponding to periods between 1–20 years and we use 1000 frequencies for the logarithmic grid. The top panel of Figure B.3 shows the LS power at a given period for the inner jet PA. We find the highest peak (power ∼ 0.4) at 10.8 ± 2.3 years, denoted by the gray dashed line in the frequency domain. In order to verify the significance of the peak, we firstly test the effect of uneven-sampling on the LS periodogram by plotting the window function (bottom panel of Figure B.3). While there is a small peak around 7 years, the value of the window functinon is below 0.05 at the PA maximum of ∼11 years, so that uneven sampling does not lead to a falsely detected peak here. In addition, the LS power at this period is larger than the 5σ significance level computed from red noise, indicating that this is a strong detection.

The top panel of Figure B.4 shows the LS result for the core flux density, where the explored periods are 1–15 years with 1000 frequencies. We find three peaks above the 2σ significance level: at 6.2 ± 0.9 years, 3.4 ± 0.3 years and 2.4 ± 0.2 years. We note that the ∼6 years and ∼3 years peaks coincide with the two first significant peaks obtained from the Jurkevic analysis. However, while the peak at ∼6 years attains the 3σ significance level, it also corresponds to a strong peak in the window function of power ∼0.20, which indicates that there are likely non-negligible effects from uneven sampling.

Fig. B.3. LS periodogram for the PA (upper panel) and window function (lower panel). The significance levels are indicated by the (yellow, red, blue) dashed lines (2σ, 3σ, 5σ). The gray dashed line marks the peak at 10.8 years.

Fig. B.4. LS periodogram for the core flux density (upper panel) and window function (lower panel). The gray dashed lines mark the peaks at 6.2, 3.4 and 2.4 years.

B.3. Weighted-wavelet Z transform

Finally, we perform the WWZ transform. The WWZ method relies on wavelet analysis to search for periodicity in both time and frequency, and is able to detect signals with varying shape and amplitude. In addition, it is well-adapted to analyze nonuniformly sampled data. We use the libwwz publicly available package⁴, based on the derivations in Foster (1996), to compute the WWZ power.

For a function of time x(t), which describes in our case either the flux density or the PA, the wavelet transform is:

$$\begin{array}{c}\hfill W(\upsilon ,\tau ;x(t))={\upsilon }^{1/2}{\displaystyle \int x}(t)f^{\ast }(\upsilon (t-\tau ))dt,\end{array}$$(B.2)

where υ is the scale factor or frequency, with υ⁻¹ typically called dilation, τ is the time shift, f is the wavelet kernel and f^* is its complex conjugate. We adopt the abbreviated Morlet wavelet (Grossmann & Morlet 1984),

$$\begin{array}{c}\hfill f(\upsilon (t-\tau ))=e^{i\upsilon (t-\tau )-c{\upsilon }^2{(t-\tau )}^2},\end{array}$$(B.3)

which is acceptable because we set c to be small enough for the additional constant term e^−1/4c in the traditional Morlet wavelet definition to be negligible. The Morlet wavelet describes a Gaussian with harmonic modulation and the abbreviated form is very similar to a windowed Fourier transform with frequency-dependent window e^{−cυ2(t − τ)2}.

Following the libwwz approach, we perform the WWZ transform analysis by splitting the time domain into n_τ = 100. The linear frequency range spans 1 − 24 years for the PA and 1 − 15 years for the core flux density, with a frequency step of 0.001. We set c = 1/32π² = 3.17 × 10⁻³ < 0.02 for the decay constant. As for the LS periodogram, we include significance level estimations that account for red-noise (see Section B.4).

The results for the PA are given in Figure B.5 (bottom panel). The left panel shows the 2D WWZ power as a function of both time and frequency and the right panel shows the τ-averaged power with frequency, together with the significance levels. We find a clear peak at 16.8 ± 4.3 years, once again above the 5σ significance level. Looking at the 2D map, this is a persistent periodicity signal over time. Given the strong signals obtained with the three analyses, we conclude that a PA period of the order of a decade is firmly detected. In addition, we notice from the WWZ transform a new period at 7.7 ± 0.8 years, above the 2σ significance level. The strength of this signal being stronger at earlier times, it is not revealed by the Jurkevic or LS methods.

Figure B.5 shows the WWZ transform results for the core flux density (upper panel). We find a signal reaching the 3σ significance level, which peaks at 5.0 ± 0.4 years. From the left panel of Figure B.5, this ∼5 years signal is persistent over time but merges around 2020 with a period of approximately 3.5 years becoming stronger with time. This explains the larger period detected with the WWZ transform compared to the two previous methods. We also find two other peaks above the 2σ significance levels. The strongest one is detected at 10.3 ± 3.5 years, which is very wide and is the continuation of the weaker peak at 6.9 ± 0.1, which barely reaches the 2σ threshold. While the significance of the results is lower than for the PA, we do detect a ∼7 years and ∼4 years signal with all of our periodicity search methods, warranting consideration. The ∼10 year period is also detected with two of the three methods. Additional data and better sampling would allow us to confirm whether these periodicities are meaningful.

Fig. B.5. Left: WWZ power as a function of frequency and time for the inner jet PA (top) and the core flux density (bottom). Right: τ-averaged power and (2σ, 3σ, 5σ) significance levels in the (yellow, red, blue) dashed curves.

B.4. Significance levels estimation

We expect the 43 GHz data of 3C 66A to be affected by red noise. It is especially known that this noise impacts blazar LCs and must be taken into account when looking for periodicity signals. Red noise is the term used to denote stochastic processes, whose Fourier power spectrum 𝒫 can be modeled simply as a power-law following 𝒫(ν) = ν^−α, where ν is the frequency and α the slope. In order to estimate significance levels in the LS periodogram and WWZ transform search methods that account for red noise, we apply the well-established PSRESP⁵ (power spectral response) approach Uttley et al. (2002). It generates artificial LCs following the method by (Timmer & König 1995) that have a power-law PSD and the same sampling as the observed data, for a range of α values and a range of periodogram and LC bins. The periodograms of the simulated LCs and of the real data are then compared. A success fraction based on chi-squared statistics is obtained. Although the PSRESP model was developed for LCs and red-noise estimation is generally applied to flux density data, we choose to perform the same red-noise estimation also for the PA. Indeed, such red noise may be present in the data and results in larger significance levels compared to white noise. We therefore model it in order to ensure the significance of the periodicity found for the inner jet PA.

We obtained α_PA = 1.7 ± 0.7 for the inner jet PA and α_CFD = 1.5 ± 0.3 for the core flux density. We use these mean values in order to simulate artificial LCs that reproduce the level of red noise of our data when performing the LS periodogram and the WWZ transform. In the first case, we produce 10 000 artificial LCs and measure the LS power to obtain the 1–5σ significance levels. However, we use 1000 LCs in the case of the WWZ method because this analysis is computationally expensive.

Appendix C: Precession model fitting

Fig. C.1. Precession fit to the inner jet PA for 100 random MCMC samples (thin gray lines) and the mean parameters (red line).

We perform a likelihood analysis to model the inner jet PA with a precessing cone. We consider the set of five free parameters ϑ = (Ω, P_p, t₀, ϕ₀, θ₀) and a Gaussian log-likelihood function for the PA:

$$\begin{array}{c}\hfill ln\mathcal{L}(\vartheta )=-\frac{1}{2}{\displaystyle \sum _i}\frac{{(y_i-g(t_i,\vartheta ))}^2}{{\sigma }_i^2}\end{array}$$(C.1)

with y_i the PA data at time t_i, σ_i the corresponding uncertainty from the observations and g(t_i, ϑ) the model prediction for the PA at time t_i and for the parameters ϑ. The precession model g is described by Eq. 2, topped with the linear decrease of −0.83°/yr obtained in Section 3.2. This ensures that the precession model considered here is applied only to the short-term oscillation of the inner jet PA.

We assume a uniform prior on the five parameters searched. For the base case: Ω ∈ [0, 10]°, P_p ∈ [6, 20] yr, t₀ ∈ [2014, 2020] yr, ϕ₀ ∈ [0, 90]° and θ₀ ∈ [170, 200]°. In the relaxed constraint case ϕ₀ ∈ [0, 15]° and in the tight constraint case ϕ₀ ∈ [0, 5]°. The second case is motivated by the fact that 3C 66A is a BL Lac object, and the third is motivated by reported viewing angle values of 3C 66A ranging between 1 and 5 degrees approximately (Jorstad et al. 2017; Weaver et al. 2022). Due to the uncertainty on the viewing angle, we choose not to set a lower limit but only an upper limit on the prior even for the tight constraint case. The viewing angle Θ and the angle between the precession cone axis and the line of sight ϕ₀ are related by Eq. 7. We have checked that, for the parameter ranges considered in this analysis, setting ϕ₀ < 15° is equivalent to Θ < 18° and setting ϕ₀ < 5° is equivalent to Θ < 6°. The log-prior reads,

$$\begin{array}{c}\hfill lnP(\vartheta )=\{\begin{array}{cc}0,\hfill & \text{if parameters are within bounds}\hfill \\ -\infty ,\hfill & \text{otherwise.}\hfill \end{array}\end{array}$$(C.2)

Finally, the posterior probability is given by Bayes’ theorem:

$$\begin{array}{c}\hfill P(\vartheta |y)=\frac{P(y|\vartheta )P(\vartheta )}{P(y)}.\end{array}$$(C.3)

We use the emcee MCMC module to sample from the posterior distribution. We set up 30 walkers and iterate over 20 000 steps, discarding the first 2 000 as burn-in to remove transient behavior. We obtain a robust mean acceptance fraction of 0.5. The corner plot results are shown in Figure 5 for the three cases studied. We take the mean and standard deviation of the samples in the tight constraint case ϕ₀ < 5° for our resulting precession model values (see Table 2). The PA evolution with 100 random MCMC samples is shown in figure C.1, highlighting the successful fit of the model to the data.

Appendix D: BH mass estimate

We estimate the primary BH mass M₁ in 3C 66A based on optical variability timescale measurements, with the corrected redshift value of z = 0.340 compared to the value of 0.444 originally used in the literature (Kaur et al. 2017; Fan et al. 2018). We assume that the minimum variability time-scale Δt_min observed for 3C 66A corresponds to the orbital time-scale of the innermost stable orbit around a central BH. This gives an upper limit that we consider to be the primary BH mass M₁ (e.g., Xie et al. 2002; Kaur et al. 2017):

$$\begin{array}{c}\hfill M_1=1.62\times {10}^4\frac{\delta }{1+z}\mathrm{\Delta }{t}_{\min }{M}_{\odot },\end{array}$$(D.1)

where δ is the Doppler factor of the jet. For this parameter, a large number of values have been reported in the literature on 3C 66A. We take the weighted mean and standard deviation of the following list of factors: δ = 1.99 (Lähteenmäki & Valtaoja 1999), δ = 6.0 (Böttcher et al. 2005), δ = 2.6 (Hovatta et al. 2009), δ = 5.69 (Fan et al. 2014), δ = 4.3 (Liodakis et al. 2017), δ = 16.5 (Jorstad et al. 2017), δ = 11.4 (Weaver et al. 2022), where the last two values are in turn the average Doppler factor measured from 2 and 5 radio jet knots. For the values without reported error, we use an informed estimation from the measurements for which we have a standard deviation. We obtain δ = 4.46 ± 0.19. We note that for all the Doppler factors in the literature the authors assumed a redshift z = 0.444. While it is in disagreement with the value used in this paper, the average Doppler factor remains valid for our study. Indeed, in the method used by the above authors to measure δ, the choice between z = 0.340 and z = 0.444 has a negligible impact compared to the uncertainties on the other parameters and ultimately on the Doppler factor value that we use.

We take Δt_min = 37 min, which is the shortest among the optical intra-day variability timescales found by Kaur et al. (2017), to measure the primary BH mass. This leads to M₁ = (1.20 ± 0.05)×10⁸ M_⊙. We consider the typical range of mass ratios 30:1-3:1 derived by Gergely & Biermann (2009) for SMBH mergers. Taking q = 0.18 ± 0.15, we obtain the central BH mass M = (1.42 ± 0.19)×10⁸ M_⊙.

Appendix E: Tables

In Tables E.1 and E.2 are respectively given the parameters of the KaVA maps and of the Gaussian components of the images at each epoch, both at 22 GHz and 43 GHz, shown in Figure 1.

In Table E.3 are given the PA, core flux and total flux at each epoch for the 43 GHz VLBA archival data, shown in Figure 3.

All Tables

All Figures

Fig. 1. KaVA intensity maps of the 3C 66A radio jet at 22 GHz (top) and 43 GHz (bottom) in 2014. The black contour levels start at three times the rms value, scaling twice. The grey shaded ellipses above each contour show the restoring beam. The detailed imaging parameters are summarized in Table E.1. The horizontal spacing is proportional to the observation time gap. The exact observing time is in the unit of years. The red circles are the model components obtained from modelfit. The continuous curved jet extends downwards to ∼3 mas on both panels.

In the text

	Fig. 2. PA of the Gaussian components of the 3C 66A jet as a function of the core separation obtained with the KaVA observations at 22 (left) and 43 GHz (right) in 2014. Each epoch is marked by green circles, blue triangles, and red squares.
In the text

	Fig. 3. (a) Inner jet PA versus time at 43 GHz. The dashed curve in the upper panel is the analytic sinusoidal+linear fit. The lower panel shows the residuals. The fitting parameters are summarized in Table 1. (b) Time evolution of the total (blue) and core (black) flux density.
In the text

Fig. 4. (a) DCF analysis between the PA and the core flux density. The grey dotted and dash-dotted lines show the 2σ and 3σ significance levels, respectively. (b) Core flux density vs. inner jet PA plot. Each panel displays the data between 1996 and 2010 (left) and between 2010 and 2025 (right). We use changing colors to show the time sequence: grey denotes the values from all observing epochs, purple to dark blue to clear blue denotes the data observed between 1996, 2002, and 2010, while green to orange to red between 2010, 2017, and 2020.

In the text

	Fig. 5. Corner plots of the likelihood analysis for the precession model given three ranges of ϕ0. The parameter values shown are the median and 1σ credible intervals of the posterior samples.
In the text

	Fig. 6. Period of LT disk precession as a function of the spin parameter a, outer radius Rout and density power law slope s (upper panel) or density exponential factor σ (lower panel). The black dashed line corresponds to PLTdisk = 10 years and the gray one to PLTdisk = 1000 years.
In the text

	Fig. 7. 2D plot of the BH spin axis precession period as a function of the spin parameter a and the accretion rate Ṁ for αvis = 0.1. The black dashed line corresponds to PBH = 10 years and the gray one to PBH = 1000 years. The area in between the vertical dotted white lines corresponds to the typical range of accretion rates for BL Lac objects.
In the text

	Fig. A.1. Evolution of the PA for the annuli A (inner jet: black), B (downstream jet: blue), and C (purple) at 43 GHz, obtained with the annular bin method. The data is shown in the upper panel and the analytical fits in the lower panel: sinusoidal+linear fits for the annuli A and B and linear fit for the annulus C.
In the text

	Fig. A.2. Evolution of the PA for the Gaussian components J1 (inner jet - maroon), J2 (downstream jet - orange) and J3 (brown) at 43 GHz, obtained with the circular Gaussian fitting method. The data is shown in the upper panel and the analytical fits in the lower panel: sinusoidal+linear fits for the J1 and J2 components and linear fit for the J3 component.
In the text

Fig. A.3. VLBA intensity maps at 43 GHz overplotted with three Gaussian fitted jet components (white) and core (black) at six epochs. In the leftmost panel, the labels J1, J2 and J3 denote the three jet components. The blue connected dots indicate the jet ridgeline. The pink contour levels start at three times the rms value, scaling twice. The gray ellipse at the bottom of each panel corresponds to the restoring beam.

In the text

	Fig. B.1. Jurkevic test plot for the inner jet PA. The black line corresponds to the Vm2 values whose minima indicate possible periodicities and the dotted blue (orange) line to the level F = 0.5 (F = 0.25). The vertical gray dashed line indicates the period at the minimum.
In the text

	Fig. B.2. Jurkevic test plot for the core flux density. The black line corresponds to the Vm2 values whose minima indicate possible periodicities and the dotted blue (orange) line to the level F = 0.5 (F = 0.25). The vertical gray dashed lines indicate the periods at the minima of interest.
In the text

	Fig. B.3. LS periodogram for the PA (upper panel) and window function (lower panel). The significance levels are indicated by the (yellow, red, blue) dashed lines (2σ, 3σ, 5σ). The gray dashed line marks the peak at 10.8 years.
In the text

	Fig. B.4. LS periodogram for the core flux density (upper panel) and window function (lower panel). The gray dashed lines mark the peaks at 6.2, 3.4 and 2.4 years.
In the text

	Fig. B.5. Left: WWZ power as a function of frequency and time for the inner jet PA (top) and the core flux density (bottom). Right: τ-averaged power and (2σ, 3σ, 5σ) significance levels in the (yellow, red, blue) dashed curves.
In the text

	Fig. C.1. Precession fit to the inner jet PA for 100 random MCMC samples (thin gray lines) and the mean parameters (red line).
In the text