Is OJ 287 a Single Supermassive Black Hole?

Light curves for more than century optical photometric observations of the blazar OJ 287 reveals strong flares with a quasi-period of about 12 years. For a long time, this period has been interpreted by processes in a binary black hole system. We propose an alternative explanation for this period, which is based on Doppler factor periodic variations of the emitting region caused by jet helicity. Using multi-epoch very large baseline interferometry (VLBI) observations carried out in a framework of the MOJAVE (Monitoring Of Jets in Active galactic nuclei with VLBA Experiments) program and other VLBA (Very Long Baseline Array) archival experiments at the observing frequency of 15 GHz, we derived geometrical parameters of the jet helix. To reach an agreement between the VLBI and photometric optical observation data, the jet component motion at a small angle to the radial direction is necessary. Such non-radial motion is observed and, together with the jet helical shape, can be naturally explained by the development of the Kelvin-Helmholtz instability in the parsec-scale outflow. In this case, the true precession of the OJ 287 jet may manifest itself in differences between the peak flux values of the 12-year optical flares. A possibility to create this precession due to Lense-Thirring effect of a single supermassive black hole is also discussed.


Introduction
Blazar OJ 287 is the first active galactic nucleus in the center of which it was assumed a system of two supermassive black holes. This assumption was made by [1] to interpret the 100-year optical light curve of OJ 287, which shows powerful flares that repeat approximately every 12 years. These flares were explained by tidal interaction of the secondary black hole with an accretion disc of the primary one at their maximum approach. The idea of the binary black hole in the center of OJ 287 was further developed in [2][3][4][5][6], which suggested that the accretion disc of the primary is not in the orbital plane of the secondary component. In these models, the orbit of the secondary black hole is precessing, thus allowing to explain the observed difference in time intervals between the 12-year flares, which is about a year.
The assumption about the binary black hole system in the center of OJ 287 was made before the unified scheme of the active galactic nuclei (AGN) was proposed [7], which later became widely used in the interpretation of the properties of active nuclei. According to this scheme, blazars are a class of active galactic nuclei with relativistic jets directed close to the line of sight. The blazar OJ 287 shows strong variability in the entire observed spectral range, a high polarization degree of optical (see, e.g., [8] and references therein) and radio emission [9]. Therefore, we suppose that the observed optical emission is synchrotron and comes from the parsec-scale jet. The optical emission of the accretion disc is low enough to be distinguished in the total spectrum of blazars, as it was recently proven [10,11]. Due to relativistic effects, the radiation formed in the approaching blazar jet increases in the observer's reference frame and de-boosts the emission from the receding jet. The relativistic boosting depends on the speed of emitting plasma and the angle of the motion with respect to the line of sight. Thus, the periodicity in the light curve can be created by the motion of the radiating plasma along a helical trajectory, as discussed in [12]. This scenario was applied to OJ 287 [13]. Namely, authors assumed the presence of two approaching relativistic jets, twisted together. These jets are formed in a system of two black holes of the same mass. Flares with a two-peaked profile occur at times when the first jet and later the second jet are at a very small viewing angle [13]. An alternative explanation is that a secondary black hole passes through the accretion disc of the primary component twice during the orbital period [2]. Taking into account the relativistic amplification of the radiation generated in the jet, we think that models explaining the 12-year optical flares by a helical jet are more likely.
Recently, Villforth et al. (2010) [14] argued the existence of a single supermassive black hole in the centre of OJ 287 based on the analysis of optical photo-polarimetric observed data for the 12-year flares. Under this assumption, the 12-year flares arose from avalanche-like accretion of the magnetic field. Britzen et al. (2018) [15] had shown that the observed period of radio flux variability is about 25 years and that the changes of the inner jet position angle can be explained by the jet precession, which can be a result of either the orbital motion in the binary black hole system or the precession of the outer part of the accretion disc of the single supermassive black hole. However, there is no explanation for the difference of the variability periods observed in the radio and optical ranges.
Direct observations of parse-scale jets of active galactic nuclei can be made by very large baseline interferometry (VLBI). For most sources, including OJ 287, VLBI maps show a typical morphology manifesting one-sided jet, the apparent origin of which is the brightest compact feature, called the VLBI core. Due to energy losses on synchrotron radiation and adiabatic expansion, the jet quickly dims and becomes undetectable at a distance of few milliarcseconds from the core. As a result of the analysis of polarimetric observations of two hundred active galactic nuclei performed with the very long baseline array quasi-simultaneously at frequencies from 8.1 to 15.4 GHz, it was showed that the VLBI core observed at a certain frequency is a part of a jet in which the medium becomes optically thick for radiation at this frequency due to opacity caused by synchrotron self-absorption [16,17]. Due to the fact that the magnetic field strength and density of radiating particles decrease with the distance from the true jet base according to the power-law dependence [18,19], the lower the frequency of observations, the further from the true jet base the absolute position of the VLBI core. This implies that the observed optical radiation comes from a region close to the true jet base. Therefore, if we associate the 12-year optical flares with a cyclic change in the angle between the velocity vector of the radiating region and the line of sight, we can expect similar periodicity in the data of radio observations. A comparison of the flux density of AGN measured with VLBI and single-dish observations showed that VLBI core emission significantly contributes and often dominates in a total source radiation [20]. Hence, in measurements of the total radio flux density from the entire source, the periodicity cannot be strongly obscured by radiation from the optically thin part of the jet. Statistical analysis performed in [21,22] based on data from the Owens Valley Radio Observatory (OVRO, 2008−2018) and Radio Astronomy Observatory of Michigan University (UMRAO, 1974(UMRAO, −2011 shows that the period of long-term radio variability is about 25 years. This period washed out the period of about 6 years, obtained previously using shorter data set [23]. We also note that the periods obtained from the total radio flux density (≈25 years [21,22]) and from the motion of the jet ridgeline (≈22 years [15]) are consistent. The periodicity of ≈1 year is found in both optical and radio domains [15,22]. We do not consider this short periodicity here, though it might reflect the rotational motion of the plasma in the jet flow [24] or could be caused by nutation-like motion in the precessing jet [3,15]. We believe that the issue of the absence of a prominent variability in the radio range with a period clearly seen in the optical domain is fundamental for understanding both the processes operating in OJ 287 and for determining the properties of the central object.
The layout of this paper is the following. In Section 2, we derive geometrical and kinematic parameters of the OJ 287 jet assuming that it has a helical shape and using 15 GHz VLBA observations since 1994 through 2019, performed in the framework of the Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE 1 ) project; we also provide interpretation of the difference between the variability periods in the radio and optical ranges. In Section 3, we assume that different values of the maximum flux density of the 12-year flares can be caused by precession of the helical jet. We showed that the period of this precession is consistent with the Lense-Thirring precession of a single supermassive black hole or the inner part of its accretion disc. Discussion of the obtained results and conclusions are presented in Section 4.

Helical Jet Model
In this section, we introduce geometrical parameters used by us to describe a jet helical shape and motion of jet's components. Then, we estimate these parameters from the VLBI observational data. In the framework of the obtained geometrical and kinematic model, we explain the observed difference in periods of flux variability in the radio and optical ranges.

Description of Parameters
Here, we consider the jet as a continuous flow, having a helical shape. According to [17,25], magnetic field strength and number density of emitting particles decrease with the distance from the jet beginning. Due to synchrotron opacity, the emission at different frequencies arises at different parts of this continuous jet. We used the geometrical description of a helical jet proposed in [26]. Namely, we assume that the jet axis forms a helix located on the surface of a notional cone with a half opening angle ξ ( Figure 1). The angle of the cone axis with the line of sight is θ 0 . Note that the cone apex is not necessarily expected to coincide with the supermassive black hole location. The jet swirl is characterized by the parameter ρ, which is the angle between the tangent to some point of the jet axis and the cone generatrix drawn through this point. To describe the jet kinematics, the jet was divided into components by cross-sectional drawing such that the jet axis can be considered a straight line within each component. Then, β is a speed of the jet component and p is an angle between the velocity vector and a cone generatrix passing through the component. The jet component location at the cone surface relative to the observer is characterized by an azimuth angle ϕ. The origin for ϕ is located on the surface farthest from the observer cone and lies in the plane containing the line of sight and the cone axis. The angle ϕ is measured counter-clockwise. We assume that the considered parameters do not change with time and that they are the same for all components located within the considered jet, which is up to a dozen parsecs in the projection on the sky plane. Note that we do not associate the jet components directly with the jet features observed with the VLBA, as their physical nature may be different, but these features reflect the position of the jet flow on the sky plane.

Estimation of Parameters by VLBI Data
A difference of the jet from straight flow appears in changes of jet feature positions onto the plane of the sky. A jet feature position projected onto the sky is characterized by a position angle (PA) measured toward the rise of right ascension between directions from the VLBI core to the north and the given feature. We made use of a set of 145 epochs of VLBA observations of OJ 287 carried out at 15 GHz since 1994 through 2019. Most of the data are taken within a framework of the ongoing MOJAVE program, while the rest are from the VLBA achival experiments. The core component was identified based on results of structure model fitting as an apparent origin of the jet. For each epoch, we constructed a ridgeline of the jet in total intensity using fully calibrated data in the image plane ( Figure 2). Not to be affected by the elliptical restoring beam, which differs from session to session, we convolved the images with the same circular beam, averaged over all epochs. The algorithm for ridgelines starts with making azimuthal slices at progressive separations (with a step of 0.1 mas) from the VLBI core position and finding the points where the integrated intensity across the outflow is equal on both sides of the arc. A cubic spline is then fit to the points and interpolated. The obtained ridgelines are presented in Figure 3. First, the ridgelines are bent. An apparent bend angle is amplified by the projection effects due to a small jet viewing angle. Second, ridgelines at different epochs have noticeable shifts from each other. Third, the direction of the inner jet ridgelines (e.g., within 0.2 mas from the VLBI core) significantly changes, with the magnitude of α obs ≈ 90 • during the first half of 2010. The ridgelines after the rapid change in PA of the inner jet direction (since May 24, 2010) are shown by blue color.  The innermost jet direction can be derived even more accurately by modelfitting the restored brightness distribution performed in the visibility domain. For this, we used the procedure "modelfit" in the Difmap package [27]. The source structure was modeled with an elliptical and circular Gaussian components for the core and the inner jet feature, respectively. Thus, the position angle of the major axis of the core component gives the direction of the innermost jet. This approach provides the highest effective angular resolution as the core is the brightest feature having the highest dynamic range of the restored source structure.
We plotted changes of the innermost jet position angle (PA in ) with time in three ways ( Figure 4). In the first way, the PA is calculated for the position of the ridgeline point closest to the 15 GHz core feature. In the other two ways, the VLBI core was fitted by an elliptical Gaussian component. In the second way, the fitting of the jet was performed with one circular Gaussian component, while in the third way, the fitting was performerd with two circular Gaussian components, in addition to the elliptical Gaussian for the core. The third approach models the rest of the jet better and could be considered the most accurate way of assessing PA in values. In the last two ways, the PA of major axis of the core fitted component was taken as PA in . The median size of the major axis of the core component is 0.19 mas, suggesting that we are sensitive to the innermost jet direction at core separations ≤ 0.1 mas. Figure 4 shows that, despite the available spread of points defined in various ways, there is a general trend of decreasing PA in until 2009-2010, followed by an abrupt increase of PA in by about 100 • and its further gradual decrease. Note that the sudden change of PA in defined by the major axis of the elliptical Gaussian component fitted VLBI core occurs about a year earlier than the sudden change of PA in defined by the innermost part of the jet ridgeline. This suggests that the cause or conditions in the jet which result in the sudden change in PA in values propagate downstream the jet. Noteworthily, PA in values obtained from 43 GHz VLBA observations correspond to those measured at 15 GHz, but the rapid change in the 43 GHz PA in occurred in the middle of 2004 ( Figure 6 in [28]). This can be explained by weaker opacity of synchrotron emission at 43 GHz, which is generated in the jet region located upstream from the 15 GHz VLBI core.
Cohen et al. [28] interpret the abrupt change in PA in as a shift of the jet nozzle in a new direction. For this scenario though, it is problematic to explain the presence of a trend of decrease of PA in with time after the PA in sudden change, the slope of which is similar to that before the PA in rapid rotation. On the other hand, the helical jet considered here can explain the observed change in PA in provided that the ratio θ 0 /ξ is slightly more than unity [26]. According to [26], the change in PA in as a function of the azimuth angle ϕ is given as where PA 0 is a mean position angle (PA of the cone axis). There is a minus sign before the arc tangent because, unlike the considered case in [26], a smooth decrease in PA in occurs before and after the sudden change. That is, in our case, the jet helix twists in the other direction: clockwise. We fitted PA in changes using Equation (1) and expressing the azimuth angle as ϕ = 2πt/T + ϕ 0 , where T is the period of PA in change, t is the time in the observer's reference frame, and ϕ 0 is the initial phase. The curves fitted to the observed change in PA in are shown in Figure 4 (right panel), and the obtained parameters for the fits are listed in Table 1. There is a good agreement between the values of T and θ 0 /ξ obtained from the three different ways of determining of PA in . From Figure 4 and Table 1, it is seen that the approximations of PA in values, which are found from the inner jet ridge lines and from the source structure modeled by two components (an elliptical Gaussian for the VLBI core and a circle Gaussian for the jet), are similar. The difference between them in ϕ 0 and PA 0 can be explained by the fact that, in the second case, PA in is measured by radiation detected from smaller angular scales than in the first case. Additionally, in the second case, the PA in is sensitive to unresolved bright jet parts located upstream of the region where the PA in is measured from the ridgeline. The qualitative difference between the approximation of PA in obtained from PA of an elliptical Gaussian component of the VLBI core under the jet fitting by two circular Gaussian components and the other two approximations is probably caused by a large spread of values.
For further calculations, we use the values T PA = 28.3 years and θ 0 /ξ = 1.5 • . Analysis of the VLBA images stacked over all available epochs at 15 GHz for a given source for a sample of a few hundred AGN showed that jet features propagate within a certain cone [29]. This result is consistent with the assumption considered by us about the jet shape and allows to accept ξ = 1 • to further calculations. Then, θ 0 = 1.5 • , which agrees with other estimations of the jet angle with the line of sight performed by the combined analysis of the jet kinematics and the long-term radio variability [30,31]. The helical jet parameters found by us are listed in Table 2. Table 2. Geometric and kinematic parameters derived for the blazar OJ 287 helical jet on parsec scales.

Parameter
Sign Value half-opening angle of the cone ξ 1 • angle between the cone axis and the line of sight θ 0 1.5 • variability period of PA in T PA 28.3 years speed of the jet component β 0.9979 angle between the jet velocity vector and a radial trajectory p 2.5 • angle between a tangent to the jet flow and a cone generatrix at a given point ρ 10 • We emphasize the following fact. The obtained variability period of PA in , caused by a cyclic change of the azimuth angle of the helical jet components located at the fixed distance from the cone apex, is consistent with the variability period of the radio emission at cm wavelengths [15,21,22]. It can mean that the long-term radio variability is also connected with the jet helical shape causing the periodical variations of the jet angle with the line of sight. As we mentioned above, this scenario seems to be the most probable for the interpretation of the 12-year optical flares as well. The explanation of the difference in the periods of long-term radio and optical flux variability is given in the next section.

Connection between Long-Term Optical and Radio Flux Variability
When the emitting region moves with an ultrarelativistic velocity β at a small angle θ to the line of sight, the spectral flux density of radiation formed in this region increases in the observer's reference frame where α is a spectral index of radiation described by a power-law (F ∝ ν −α ). The Doppler factor is We estimate the component speed from the apparent speed of the jet features β app . Kutkin et al. [32] showed that the jet speed, which is found by time delay of a flare when the observed frequency decreases, corresponds to the maximum β app estimated from the apparent motion of jet components. Therefore, we used the maximum value of β app = 15.31 for OJ 287 [33] to find β = 0.9979 assuming Γ min = (1 + β 2 app ) 1/2 . If β = const, changes in δ occur due to changes in θ = θ (ϕ). In order for the changes in δ to be reflected in the observed flux variability, the radiating region must be compact enough. That is, the angle ϕ should change slightly within this region. Optical emission comes from near to the true jet base region [18]. This region is characterized by high magnetic fields and accelerated particles that have not yet lost a greater part of their energy through radiation. The magnetic field strength and number density of emitting particles decrease according to a power-law with distance from the true jet base [17][18][19]. These gradients together with high radiation energy losses of ultrarelativistic electrons via synchrotron optical emission allow us to assert that the region in which the observed optical radiation is generated is compact. Based on the Blandford-Königl jet model [18], for further calculations, we assume that the optical emitting region is located at a fixed distance from the true jet base. Therefore, we can expect that, at each moment of time, the Doppler factor within the optical radiating region has an approximately constant value. In the radio range, the emission detected by single-dish antennae from flat spectrum radio sources, which includes OJ 287, is mainly generated in a compact VLBI core (e.g., [20]). The distance of the VLBI core observed at a given frequency from the true jet base is almost constant [17][18][19], although changes are possible on average by 0.3 mas for typical offset of the core positions at 2 and 8 GHz of about 0.5 mas [34]. Hence, we supposed that the radio-emitting region is compact enough to have one value of the Doppler factor of all its parts and is located at a fixed distance from the cone apex. Then, the changes in δ caused by variations in the azimuth angle ϕ of the region, from which the observed radiation at a certain frequency comes, should appear in the long-term light curve.
The flux variability period produced by a change in δ depends on a rate of a change in ϕ of components crossing through the emission region. In other words, in their outward motion, the components of the helical jet pass consecutively first through the region in which such parameters are reached, under which optical radiation reaches the observer without a strong absorption. After passing through this region, the magnetic field strength, electron density, and possibly energy spectrum of emitting particles become such that it is impossible to generate a significant flux density of synchrotron optical radiation. Moving further down the jet, the physical parameters in the components change, and at a certain fixed distance from the true jet base, the component becomes less opacity for 15 GHz radiation, forming the 15 GHz VLBI core, and the component becomes optically thin for the 15 GHz radiation further downstream. As shown in [35], when components move along radial (ballistic) trajectories, the change rate of ϕ for components that reach a certain fixed distance from the cone apex does not depend on the distance to the cone apex. Then, the periods of optical and radio flux variability would have to coincide. However, if the jet components move at the angle p = 0 • to the radial trajectory, then the further away from the cone apex, the lower the rate of change in ϕ and, therefore, the longer the period of δ change, which is expressed as [35] T var = 2πr sin ξ cos ρ where r is a distance from the cone apex to a cross-sectional plane passing through the region, in which the observed radiation at a given frequency is generated. The motion of the VLBI jet features of the OJ 287 jet is detected by the VLBA monitoring programs at 15 GHz [33] and 43 GHz [36]. This suggests that p > 0 • . Note that, if p = ρ, the motion of components will occur along the jet flow. The jet will look stationary in time and space, and the periods of variability of both the flux density and the PA in will be absent.
The non-radial motion of components leads to the variations of the Doppler factor of the emitting region caused by ϕ variation being significantly stronger than under the ballistic motion of components [26]. Plot of δ (θ 0 , ξ, p, ϕ) versus ϕ for the adopted above values θ 0 = 1.5 • , ξ = 1 • and a few values of p is displayed in Figure 5. For this, the expression for θ = θ (θ 0 , ξ, p, ϕ) obtained in Appendix cos θ = cos p (cos ξ cos θ 0 − sin ξ sin θ 0 cos ϕ) − sin p sin θ 0 sin ϕ is substituted into Equation (3). Values of θ from Formula (5) are equal to those calculated from Formulae (11)- (13) in [26] under the corresponding parameters. Figure 5 shows that maximum values of δ are reached for small p. For example, if p < 3 • , then δ max > 25, while at p = 5 • , the maximum value is no higher than δ = 16 . Since the 12-year optical flares have large amplitudes, for the interpretation of these flares by a change in δ, it is natural to assume such p at which the following conditions are fulfilled. First, δ (ϕ, p) reaches a high value at the maximum. Second, the ratio of maximum and minimum values of δ is close to the maximum of possible ones. Figure 5 shows that it occurs at p = 2.5 • . To use Formula 4 to find the angle ρ between the tangent to the jet flow at a certain point and the cone generatrix drawn through this point, one needs to know the values of T var and r. From the measurement of astrometric shift of the VLBI core observed at different frequencies, performed in [17], r can be found from r = r core, 15 GHz cos ρ, where r core, 15 GHz = Ω rν / (ν sin θ) is a distance of 15.4 GHz VLBI core from the true jet base, Ω rν is the core shift measure defined in [19], ν is the observed frequency in GHz, and θ is the jet viewing angle, which corresponds to θ 0 in our case. As reported in [17], OJ 287 has Ω rν < 4.18 pc GHz, which gives r core 15 GHz < 10.3 pc. According to this, we will not be much mistaken if we take T var = T PA , since PA in was determined by the inner jet part within about 0.1 mas from the 15.4 GHz core (see Section 2.2). The dependence p (ρ) given by Formula (4) is shown in Figure 6. It can be seen that p > ρ and p ≈ 10 • always for p = 2.5 • . Let us consider whether the formation of the helical jet with such geometric parameters is possible due to the development of the Kelvin-Helmholtz instability, which can be originated due to a velocity difference across the interface between the jet and surrounding medium. According to [37], a helical wave mode can be represented as a wave propagating at some angle to the cone axis. This angle is where R is the jet radius. Based on numerical simulations, Hardee [37] obtained the empirical expression for the wavelength at which the growth rate of perturbation for the helical wave mode is maximum: where M j is the Mach number in the jet and where η is the ratio of densities of a jet and its surrounding medium. The estimation of λ h by Formula (8) gives an accuracy of about 10% [37]. According to [38], the wavelength of the helical wave in the observer's reference frame is where β w = β cos p is the speed of the helical wave along the cone axis. From Formulae (7)-(9), we plotted the dependence η M j for the jet parameters found above ( Figure 7). It is seen that, for creation of the helical jet, which is able to explain the observed period of radio variability, the sets of M j and η parameters lie in the range of adequate values. The difference in the periods of optical and radio (both the flux density and positional angle of the inner jet) variability can be caused by different distances between the jet origin and the regions responsible for the observed variability. The ratio of these distances is equal to the ratio of the corresponding variability periods [35]. Then, the optical emitting region is at the distance of r opt = r 15 GHz /T PA ≈ 4.2 pc (10) from the true jet base. This distance appears to be very large, ≈ 6 · 10 4 gravitational radii for a black hole of 10 9 solar masses. However, the estimation from Formula (10) was made under the assumption that the jet axis lies on the cone surface. In theoretical works [39,40], it is shown that, near its origin, the jet propagates within the paraboloid and starting from a certain distance: within the cone. Based on the VLBA observation for jets of nearby active galaxies, the transition from the parabolic to conical jet shape occurs at distances 10 5 − 10 6 gravitational radii [29,41]. Taking this into account will reduce r opt , but it will not change the interpretation of the difference in variability periods based on a spatial displacement of the regions responsible for dominating optical and radio emission. These jet parts are at different separations from the central machine and thus have different azimuth angle, causing changes in the Doppler factor and PA in and occurring with different rates under other unchanging geometrical and kinematic parameters of the model, listed in Table 2.

Helical Jet Precession
If the 12-year flares are formed only due to the increase in δ, then, according to Formula (2), the maximum value of δ determines the flux peak value. Inserting Formula (5) into (3), we found that the function δ (ϕ) has extremum at For n = 0 and n = 1, the Doppler factor has minimum (δ min ) and maximum (δ max ), respectively. Hence, for constants p and ξ, the values of δ max and δ min are invariable and occur at a fixed ϕ max and ϕ min , correspondingly. If the flux density in the reference frame of the jet is constant, the flux peaks in the observer's reference frame must be equal as well.
Since δ at the maximum of the flux density reaches a fixed value, the difference in peak fluxes of 12-year flares can be caused either by an internal change in the jet or by a change in the angle of the cone axis with the line of sight. In the first case, the ratio of peak fluxes at the maximum Doppler factor (e.g., ≈ 3 for the flares of 1983 and 1995), would have to be maintained at other values of δ, e.g., at δ min . The observed minimum flux density F min between the flares differs only slightly [14]. Therefore, in order to agree with the observed long-term light curve, the increase in the flux density in the reference frame of the jet would have to occur at a time when the Doppler factor for an earth's observer would be close to the maximum value. There is no physical ground for this condition.
In our view, the most possible explanation for the observed difference in the peak fluxes of the 12-year flares is a change in the angle θ 0 , which, for example, may be associated with precession.
Let us check this assumption. The scheme of the helical jet precession is presented in Figure 8. The dependence of θ 0 on the azimuth angle of precession ϕ p is expressed as cos θ 0 ϕ p = cos ξ p cos θ p − sin ξ p sin θ p cos ϕ p , (12) where ξ p is the half-opening angle of the precession cone and θ p is the angle of the precession axis with the line of sight. As it follows from Equation (11), the azimuth angle values, for which the Doppler factor achieves the maximum and minimum values which do not depend on θ 0 . Assuming that the flux density in the reference frame of the emitting plasma is constant, we obtained from Equation (2) δ max (ϕ max , θ 0 , p, ξ) δ min (ϕ min , θ 0 , p, ξ) Inserting into Equation (3) the expression for θ (5), in which θ 0 and ϕ max, min are given by Formulae (12) and (11), correspondingly, and solving Equation (13) with respect to θ 0 , substituting the observed flux densities of the 12-year peaks and minima between them, we obtained where We have taken into account that ϕ min = ϕ max − π for the range of azimuth angle from 0 to 2π. The parameters used in Formula (14) and derived viewing angle θ 0 values for each peak are listed in Table 3.

Epoch (yr)
( deg ) Figure 9. Changes of the angle between the jet helix axis with the line of sight θ 0 due to the precession of the OJ 287 jet.
Fitting the θ 0 values by the expression (12), we obtained θ p = 1.8 • , ξ p = 0.7 • and the precession period in the observer's reference frame T p = 92 ± 8 years. The dependence of θ 0 with time is shown in Figure 9. Table 3. Used observational data [14,42] and obtained values of the angle between the cone axis and the line of sight θ 0 : the columns contain (1)  year F max , mJy F min , mJy θ 0 , deg (1) In order to find the precession period T p in the reference frame of the source, it is necessary to take into account not only the cosmological time change but also that the time interval in the source reference frame is t = δt. Since δ changes over a wide range (see Figure 5), we have divided the whole precession period in the source reference frame into equal time intervals. Assuming that the changes of ϕ p in the jet reference frame occur uniformly and ϕ p changes are 7 times faster than ϕ, we took the time for which ϕ p changes by 1 • as a time interval unit. Then, where δ i is the Doppler factor at the end of the time interval. From (15), we obtained T p ≈ 1.2 · 10 3 years. Thus, by substituting θ 0 from (12) into Formula (5), we can obtain an analytical expression to describe the change in the Doppler factor of the radiating region in the jet. Assuming that the 12-year flares are caused by a change in δ, we modeled the light curve. For this aim, we took the 1983 flare as a reference point. From Formula (11), we found the azimuth angle of the optical radiating region ϕ = 248.2 • . Using it together the data from Table 3, we found the azimuth angle for the precession motion ϕ p = 15.4 • by Formula (12). The Doppler factor at a given time moment depends on ϕ and ϕ p . The rate of ϕ p change was determined from T p , and it is approximately 1.5 • for 5 years in the source reference frame. Figure 9 shows that ϕ changes about 7 times faster than ϕ p . The flux density was found by Formula (2) and normalized to the flux density of the 1983 flare. Figure 10 shows the simulated light curve and observational data from the literature and the archive of the Crimean Astrophysical Observatory. As seen, since 1940, the observational data are in a good agreement with the simulated flux density. Some difference between the observed and modeled flux densities at the maxima may be caused by a small physical change in the flux density in the source reference frame compared to that of the 1983 flare maximum. Here, we do not model the radio light curve because it has no prominent maxima as in the optical range. There is no reliable reference point for which the values ϕ and ϕ p can be confidently found. Since the optical and radio radiating regions are spatially separated, the azimuth angles of the radiating regions are different within the framework of the considered model. Let us check whether the Lense-Thirring precession with such a period can occur in the accretion system of the single supermassive black hole. In this system, the precession of the black hole together with the inner part of the accretion disc can be originated in result of misalignment of the black hole and accretion disc outer part rotation axes, as described in [43]. This precession period can be expressed as where a BH is the dimensionless parameter expressing a black hole's specific angular momentum, α vis is the viscosity parameter, and M BH andṀ BH are the mass and mass accretion rate of the black hole, respectively. We expressed the Formula (16) as whereṁ is the mass accretion rate in Eddington units. Using the obtained period T p in Formula (17), we plotted a dependence of a BH (α vis ) for different masses and mass accretion rates of the black hole ( Figure 11). The range of possible values for the black hole mass is large enough. According to the broad-band observations, black hole masses for BL Lacertae objects were estimated as ∼ (10 8 − 10 9 )M [44,45, and references therein]. Valtonen et al. [42] derived the primary black hole mass of ∼ 10 10 M under the assumption of the binary black hole system in the OJ 287 centre. This mass estimation is model-dependent, but we took it as an upper limit. The mass accretion rate we accepted is in the rangė m = 0.001 − 0.1, according to the values reported in [46,47]. As seen, the dependence on the black hole mass is weak. The precession of a slowly rotating black hole with a high mass accretion rate and a low α vis (as confirmed by observations [49]) can occur with the found period. On the other hand, it is possible that misalignment of the black hole and accretion disc rotation axes can cause precession of the outer part of the accretion disc with a period [48]: where r is a radial distance of the accretion disc precessing part from the black hole and r g is the gravitational radius of the black hole. From Formula (18), it follows that, for a BH = 0.01 − 1, the regions of the accretion disc precessing with the found period are located at distances of (30 − 300)r g from the black hole with M BH = (10 8 − 10 10 )M . These distances seem to be reasonable to form the jet by the Blandford-Znajek mechanism [50]. Note that, here, we investigate the possibility of the formation of the 1200 year period in the accreting system of a single black hole. More precise determination of the parameters of this system, such as the mass, mass accretion rate, and viscosity, should be based on the broad-band photometric and VLBI data.

Discussion and Conclusions
OJ 287 is the first blazar for which periodicity in the optical light curve was detected, revealing prominent flares that occur approximately every 12 years [1]. At that time, the unified scheme of active galactic nuclei [7] had not yet been proposed. Therefore, early investigations of OJ 287 were predominantly focused on the interpretation of optical periodicity. It was supposed that such a short variability period can be formed in a binary black hole system [1]. The secondary component moves along its elliptical orbit and passes the pericenter every 9 years (in the rest frame), invoking a tidal effect on the accretion disc of the primary black hole. As a result, the matter of the accretion disc flows into the black holes, leading to the observed flare on the light curve [1]. This model does not take into account the structure of the 12-year flares, namely, that there are two peaks in each flare. The interpretation of this fact given in [2] is that the secondary black hole passes through the accretion disc of the primary one twice during the orbital period. This model also provided the explanation for some of the observed properties of optical polarization, but it was focused on orbital motion. After the next 12-year flare occurred during 1995-1996, this model has been refined [4][5][6]. The authors of [6] predicted future evolution of the VLBI jet position angle observed at 15 and 43 GHz, but it has not been confirmed by later observations [here and 28,36].
An alternative interpretation of the presence of two peaks in a 12-year flare was proposed in [13] considering two supermassive black holes of the same mass, with accretion discs and jets. These jets have a helical shape and are twisted together. Doppler beaming of a part of each of two jets, responsible for optical emission, can lead to the 2-peaked 12-year flares.
As it is noted in [1], due to the radiation of gravitational waves, a separation between the pericenters of the primary and secondary black holes is expected to decrease, leading to an increase in the intensity of next 12-year flares. The flares observed in 1995-1997 and 2005-2007, on the contrary, had significantly less maximum flux density (≈ 10 − 15mJy in V band) than the previous flares in 1973 and 1984 (55 and 30 mJy in V band, respectively) [14]. Interpretation of the decrease in the maximum flux density of the 12-year flares as well as prediction of the absence of the flares in the future is proposed in [14]. These are explained by something, that occurred before the ≈1970 event, causing an avalanche-like accretion of the disc magnetic field, which repeats every ≈12 years [14]. After each act of this accretion, the poloidal magnetic field near the black hole becomes stronger and more strongly hinders accretion of the poloidal magnetic field next to an act of accretion. This model requires further verification by magnetohydrodynamic simulations. The maximum flux density of the last 12-year flare occurred in 2016 was 17 mJy in B band [42]. It is slightly higher than the ones for two previous flares ( Table 3).
None of the above models can consistently explain the main observable facts: (i) the sudden change in PA in which was first detected at 43 GHz in 2004 [28,36] and about six years later at 15 GHz; (ii) the difference of the long-term variability period in optical (≈12 yr [1]) and radio (≈25 yr [21][22][23]) flux densities.
In this paper, we explain most of the observed long-term properties of the blazar OJ 287, while other properties can also be interpreted within the framework of the proposed model. First of all, we assume that the OJ 287 jet has a helical shape because this shape can be naturally originated by the jet nozzle precession or development of hydrodynamic instabilities. For the certain relation between the angle of the helix axis to the line of sight and the half-opening angle of the cone, the sudden change in the inner jet position angle occurs under the outward motion of the helical jet (θ 0 /ξ 1). As it was shown in [26], for the permanent geometrical parameters of the jet helix, the sudden change in PA in occurs at a certain fixed azimuth angle. For the case of OJ 287, PA in has a minimum and maximum value at ϕ of about 132 • and 228 • , correspondingly. Therefore, the fact that the sudden change in PA in occurred at a higher frequency first and then at a lower one indicates that the 43 GHz emitting region reached this value of ϕ earlier than the 15 GHz emitting region. This means that the regions responsible for dominant radiation at different frequencies are spatially separated. This fact follows from the jet model [18] and has been confirmed by VLBI observations for a large number of AGN jets [17]. The higher the observed frequency, the closer to the true jet base is the location of the region, in which the observed radiation is generated. Hence, the region from which the observed optical radiation comes is located much closer to the true jet base. According to this, the changes in the flux density in the optical and radio ranges, caused by the geometric effects, do not coincide in time, since the Doppler factor of the emitting regions differs. The absence of sharp and prominent flares in radio similar to those present in optical range can be explained by a significantly larger size of the radio-emitting region compared to that in optics. This difference in the size of the regions can naturally be caused by different cooling time scales of ultrarelativistic electrons producing synchrotron radiation in the radio and optical ranges.
Within the framework of the described jet model, the difference between the periods of long-term optical and radio variability occurs when the jet components move along non-radial trajectories. Such motion is observed [28,36] and most simply explained by the development of hydrodynamic instabilities in the jet, for example, Kelvin-Helmholtz instability [37]. Moreover, we have shown that the derived parameters of the helical jet can also explain the observed variability periods. These parameters are attributed to the jet, along which the helical mode of the Kelvin-Helmholtz instability is developed [37].
In the framework of our model, the Doppler factor changes with distance from the true jet base, and for each jet component, it varies with time. Thus, for each jet component, the equal time intervals in the observer's reference frame are different. Also equal time intervals in the observer's reference frame for simultaneously observed jet features are different. These facts lead to different apparent speeds for different jet features, as it is observed for many AGN jets including OJ 287 [33].
With the geometrical and kinematic parameters of the OJ 287 jet unchanged for decades and under the interpretation of the 12-year optical flares by geometry effects, the difference in peak fluxes in these flares may be due to the fact that the helical jet axis precesses. This precession in the our geometrical model of the observed optical flux density variations can additionally explain different flux densities between the different 12-year flares and different time interval in the observer's reference frame between the 12-year flares. We found that the jet precession period in the source reference frame is ∼1200 years and showed that jet precession can be formed by the Lense-Thirring precession of a single supermassive slowly rotating black hole with a high accretion rate and a small viscosity parameter or by precession of the accretion disc at distances of 60-300 gravitational radii from a single supermassive black hole. The observed properties of OJ 287 are explained in the assumption of a single supermassive black hole in the center of the active nucleus, an important result for interpretaion of the blazar long-term flux periodicity. The latter one is often explained by the binary black hole system in AGN, by analogy with OJ 287.
The interpretation of the double-peaked structure of the 12-year optical flares is beyond the scope of this article. However, the discussed above possible explanations for this observed fact do not contradict the proposed model. Namely, the double-peaked flares can occur if some parts of the jet move at a small angle to the general helical trajectory [24], or there are two jet regions radiating in the optical range, located at slightly different distances from the true jet base. The Doppler factor of the more distant region changes as well as the Doppler factor of the region closest to the true jet base, only with a certain time delay. Indications for the existence of extended optical jets have been recently found [10,11,51]. With our model of the helical jet, we can make the following predictions for the blazar OJ 287: (i) an increase in the maximum flux density in the next three-four 12-year optical flares; (ii) a gradual decrease in the position angle of the observed at 15 GHz inner jet until ≈ 2038, when again there will be a sudden increase in the position angle (it will happen a few years earlier at 43 GHz).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A The Angle of the Jet Component Velocity Vector to the Line of Sight
We assume that the jet forms a helix on a surface of a notional cone. Find the angle θ between the velocity vector υ of the jet component and the line of sight. In Figure A1, where R 1 =AB 1 , R 2 =AB 2 , s =B 1 B 2 = υ∆t. Approximate equality is obtained for s R 1 . From Formula A1, it follows that R 1 can be found from ∆OAB 1 taking into account that d R where R = OA, d = OB 1 is a distance of the component at the initial moment of time from the cone apex and where θ 1 is the angle between the line of sight and the cone generatrix drawn through B 1 . By analogy, from ∆OAB 2 , we have where θ 2 is the angle between the line of sight and the cone generatrix drawn through B 2 . The expression for θ 1 (2) was derived in [26] θ 1 (2) = cos ξ cos θ 0 − sin ξ sin θ 0 cos ϕ 1 (2) , where ϕ 1 (2) is the azimuth angle for point B 1 (B 2 ). Inserting Formulae (A3)-(A5) in (A2) and taking into account ϕ 2 = ϕ 1 − ∆ϕ and ∆ϕ 1, we have cos θ = − d s + cos p ∆ϕ sin ξ sin θ 0 sin ϕ 1 + cos p (cos ξ cos θ 0 − sin ξ sin θ 0 cos ϕ 1 ) .