Reversals in the Direction of Polarization Rotation in OJ 287

We have obtained a smooth time series for the Electric Vector Position Angle (EVPA) of the blazar OJ 287 at centimeter wavelengths, by making $\pm n\pi$ adjustments to archival values from 1974 to 2016. The data display rotation reversals in which the EVPA rotates counter-clockwise (CCW) for 180 deg and then rotates clockwise (CW) by a similar amount. The time scale of the rotations is a few weeks to a year, and the scale for a double rotation, including the reversal, is one to three years. We have seen four of these events in 40 years. A model consisting of two successive outbursts in polarized flux density, with EVPAs counter-rotating, superposed on a steady polarized jet, can explain many of the details of the observations. Polarization images support this interpretation. The model can also help to explain similar events seen at optical wavelengths. The outbursts needed for the model can be generated by the super-magnetosonic jet model of Nakamura et al. (2010) and Nakamura and Meier (2014), which requires a strong helical magnetic field. This model produces forward and reverse pairs of fast and slow MHD waves, and the plasma inside the two fast/slow pairs rotates around the jet axis, but in opposite directions.


INTRODUCTION
Many active galactic nuclei (AGN) show a one-sided jet that can be traced inward to a few pc from the massive black hole that powers the system. This one-sidedness is a relativistic effect, in which radiation from the jet, which is composed of plasma flowing relativistically, is strongly boosted when the observer is near the axis, while the counter-jet is strongly de-boosted. The jet may also contain bright features that move superluminally downstream; i.e., their apparent velocity in the plane of the sky is greater than c, the speed of light. In these cases the observed time scale is shrunk, as the emission region follows closely behind its own radiation. This reduced time scale is also partly responsible for the rapid variability that is seen in many AGN.
Many of these jets are highly polarized, and both the fractional linear polarization and the electric vector position an-gle (EVPA) can be variable. The EVPA is measured North through East on the sky, and its variation will be our main concern in this paper. In BL Lac the EVPA tends to point along the jet (O'Sullivan & Gabuzda 2009), and this means that in the jet the transverse component of the magnetic field is dominant. The EVPA can point along the jet even around a bend (O'Sullivan & Gabuzda 2009), and this is taken as a sign that the transverse field is toroidal and that the field configuration is generally helical (Cohen et al. 2015). The jet appears to be a magnetic structure that can support MHD waves. BL Lac has been analyzed with this assumption; the superluminal components were taken as fast or slow magnetosonic waves, and the downstream propagation of the bent structure could be regarded as an Alfvén wave (Cohen et al. 2014(Cohen et al. , 2015. A gradient of the Faraday Rotation Measure (RM) across the jet, especially if there is a sign reversal across the jet, is another indication of toroidal magnetic fields, since the RM is proportional to the component of magnetic field along the line-of-sight. In a recent paper Gabuzda, Nagle, & Roche (2018) provide a list of 52 AGN that have reliable detections of transverse RM gradients, and 5 of these show time variability.
In this paper we are concerned with one particular AGN, the BL Lacertae object OJ 287, which is highly active at all wavelengths. We have made images of its jet with the VLBA 13 , a high-resolution radio instrument with EW resolution ∼ 0.6 milliarcsec (mas) at λ ≈ 2cm. OJ 287 has redshift z = 0.306 giving a linear scale of 4.48 pc mas −1 ; thus we can probe OJ 287 at scales of about one pc. OJ 287 is not in the RM gradient list of Gabuzda, Nagle, & Roche (2018), but Motter & Gabuzda (2017) have tentatively identified it as having a transverse RM gradient.
OJ 287 has provided another reason to think that the jets of AGN are threaded by helical magnetic fields. Cohen (2017) 13 The Long Baseline Observatory and the National Radio Astronomy Observatory are facilities of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. has studied the evolution of the ridge lines of OJ 287, and has shown that they are twisted, and can be interpreted as sections of a rotating helix. In the present paper the model we use contains a rotating helix, and the observations show that it has positive (right-hand) helicity.
At optical wavelengths OJ 287 shows flares, roughly 12 years apart, whose timimg can be fit to a model consisting of a binary black hole system, including spin and gravitational radiation in addition to the orbital parameters. This model has successfully predicted the appearance of flares in 2006-2010(Valtonen et al. 2011. In terms of kpc-scale radio morphology and power, OJ 287 exhibits both Fanaroff-Riley Type I and Type II characteristics; i.e., FR-I morphology and FR-II radio power. It is an exception to the simple Unified Scheme which proposes that BL Lac objects are pole-on counterparts of FRI radio galaxies (Kharb et al. 2015;Stanley et al. 2015).
In this paper we concentrate on the EVPA of OJ 287 at radio wavelengths, and report the observation of rotation reversals. One of these consists of a large counter-clockwise (CCW) swing in EVPA, followed closely by a similar but clockwise (CW) swing. Variations in the EVPA of AGN, including OJ 287, have a long history of study. Holmes et al. (1984) measured the optical polarization of OJ 287 over a 4day period and found rotations in time and also variations in frequency. Roberts, Gabuzda & Wardle (1987) made early VLBI observations of OJ 287 that separated the core and the jet components, and showed that their polarizations changed over a one-year interval. Kikuchi et al. (1988) observed a steady swing of 80 • in the EVPA in 5 days at radio wavelengths, and a nearly-simultaneous swing of 120 • in 7 days at optical wavelengths. A close correlation of radio and optical EVPA rotations has also been reported by Gabuzda et al. (2006) and by D' Arcangelo et al. (2009). Villforth et al. (2010 have made extensive optical observations of the EVPA of OJ 287. They showed that the EVPA has a long-term preferred value, 170 • , although it often appears to be chaotic. Currently, the RoboPol program (Blinov et al. 2015(Blinov et al. , 2016 is making optical polarization measurements for many AGN, including OJ 287. This paper is organized as follows. In Section 2 we discuss the observations, and first show the EVPA data from the archives. These data are erratic in some time periods, and in Figure 2a we smooth the EVPA by adding ±nπ as appropriate. This smoothing allows us to see four rotation reversals. In Section 3 we briefly consider the possibility that the rotations and reversals are spurious, and conclude that they are not. The reversals themselves are described in detail in Section 4. In Section 5 we propose a two-component model to explain an EVPA rotation as a flux density outburst with a rotating EVPA, superposed on a steady jet component. Two of these outbursts in succession, with counter-rotating EVPAs, generate the reversal. We describe a simple geometry with a relativistic jet containing a helical magnetic field that can make counter-rotating bursts in Section 6, and in Section 7 suggest that the super-magnetosonic jet model of Nakamura et al. (2010) and Nakamura & Meier (2014) can help to explain the observations. Section 8 briefly describes some aspects of the optical observations of OJ 287. Section 9, the Discussion, comments on the time scales for the rotation reversals, the many outbursts without an EVPA rotation, and on how our reversals contain 12-year separations, the same separation that is found for the repeating optical flares. Section 10 contains a Summary and Conclusions.

OBSERVATIONS
OJ 287 is a rapidly varying source, and at radio wavelengths the EVPA can change on a time-scale of days. On the other hand, the rapid EVPA changes occur episodically and are unpredictable; so that, to capture the full story of the EVPA, observations need to be made every few days and the series has to last for many years. The archives of the University of Michigan Radio Astronomy Observatory (UMRAO) provide data that meet this need . They comprise measurements of flux density (F) and polarized flux density that were made every few days with a 26-m dish, and span the years 1975-2012. However, OJ 287 passes close to the Sun every year, and 1 or 2-month gaps in the data do occur regularly, as seen in the graphs below. Only points with σ(EVPA) < 14.3 • are used here; this is equivalent to limiting the signal-to-noise ratio of the linearly-polarized flux density, P × F (hereafter simply called PF) to SNR(PF) > 2. Here P is the fractional linear polarization. Each UMRAO point is a one-day average.
We also use data from the MOJAVE program 14 (Monitoring of Jets in Active Galactic Nuclei with VLBA Experiments), which includes archival data back to 1995. This is a continuing program and for this paper we stop at 2016.0. MOJAVE uses the VLBA at 15.3 GHz. An abbreviated version of the data analysis is as follows; see Lister & Homan (2005) for details. At each epoch we make images of Stokes I, Q, and U, with pixel size 0.1 mas, and fit elliptical Gaussians (circular if possible) to the I image, to find a set of 'components'. There typically is a bright component in Stokes I at the NE end of the jet, and the center of this component is defined as the location of the 'core'. We cannot find similar components in the Q and U images because polarization cancellation in close components can result in non-Gaussian structures. Hence, we treat all Stokes parameters the same and find I, Q, and U for the core by averaging over 9 pixels centered on the core. The unit we use for I, Q, and U is Jy beam −1 . The fractional linear polarization is defined as m = Q 2 + U 2 /I, and the EVPA is calculated as EVPA = ξ/2 = (1/2) tan −1 (U/Q). In the following we mix the flux densities (in Jy) from UMRAO with the specific intensities (in Jy beam −1 ) from the VLBA and use the symbol F for all of them; the fractional linear polarization is called P, and the product PF is the linearly-polarized flux density. For the VLBA data, F and PF are the flux densities of the compact core.
We have also used results obtained by other VLBA observers at 15.3 or 15.4 GHz, and placed in the VLBA archive. In these cases the data have been reprocessed by the MOJAVE team, to make a homogeneous data set. The combined points are typically a month apart, and by themselves would be too infrequent for the rotation reversals we study in this paper, but they are useful as a check on the UMRAO points.
In addition to the UMRAO and VLBA data we use the results of Kikuchi et al. (1988), who made polarization observations of OJ 287 at several frequencies ranging from 9.0 to 10.5 GHz, for 6 months in 1986. One of our EVPA rotation reversals (Event A) occurred during their observing period, and we include part of their data in our analysis. In this period they observed on a daily basis, and this is important in reducing ambiguity in Event A. The Kikuchi et al. (1988) data were taken with the 45-m dish at Nobeyama. We use numerical values that were found by digitizing the points in Kikuchi et al. (1988), using the Dexter tool. 15 In Figure 1 we show the EVPA from the 5 archival data sets: 4.8, 8.0 and 14.5 GHz from UMRAO, 15.3/15.4 GHz from MOJAVE, and 9.0-10.5 GHz from Kikuchi et al. (1988). In the archives the data are listed in the range 0 • − 180 • , but for Figure 1 we changed the range to 50 • − 230 • , to better show the continuity of the points. We have ignored the Galactic Faraday Rotation towards OJ 287 because it is only of the order of 30 rad/m 2 (e.g., Rudnick & Jones 1983), which rotates the 14.5 GHz EVPA by less than 1 • .
In some regions of Figure 1 the EVPA varies smoothly, but in others it is highly erratic. Therefore, we sought a smooth EVPA curve by adding ±nπ as appropriate. Kiehlmann et al. (2016) have derived some procedures for this, based on a smoothness criterion, but we followed the common practice of adding ±nπ so that adjacent points differ by less than 90 • . However, we relaxed this rule when there was a substantial time gap in the observations.  have made a statistical study of how such gaps can affect the interpretation of polarization data. We also had a second criterion: make the curve fit all frequencies as closely as possible. This is important in reducing ambiguity when one frequency has a data gap that can be filled by another. Figure 2a shows the result we obtained for the smoothed EVPA when we followed both criteria. In this Figure we identify three major events and one minor event, labeled A, C, D, and B, respectively. Event D is a smooth reversal; the EVPA swings CCW by about 200 • , is stable for roughly 1.5 y and then swings CW by about 160 • . Event C is a similar reversal with the same sign (CCW then CW) and similar amplitude, but is narrower and appears to have a low-amplitude precursor. Event A includes a sharp rotation reversal with the same sign as the others, but with a larger CW swing. Event B has low amplitude and a different shape. All these events are discussed in Section 4.
We have two immediate results for OJ 287 from Figure 2. The EVPA values from UMRAO and MOJAVE generally lie close together, and so the EVPA data obtained with a 26-m 15 https://dexter.edpsciences.org/Dexterhelp.html dish are usually a good proxy for VLBA measurements for the core alone. This assertion can be tested by examining the MOJAVE polarization images (Lister et al. 2018). In most of them (51/59) the core is clearly the strongest component in PF, and so the polarization of the total source is similar to that of the core. In 8/59 images a secondary component is stronger. However, they are not distributed uniformly in time, but all occur during Events C and D. Figure 3 shows four examples of the images. In each panel the left-hand image shows the contours of Stokes I, and the linear polarization fraction is in color. The right-hand image shows the contours of PF, the linearly polarized flux density, with an additional contour that is the same as the lowest contour of the Stokes I image. In Figures 3c and 3d the cores are stronger than the secondary components in PF, but in Figure 3b the core is weaker and in Figure 3a the core and the secondary component have similar strength. We discuss this further in Sections 4.2 and 4.3.
The second result from Figure 2a is that we can assume that the EVPA is largely frequency-independent over the range 4. 8 -15.4 GHz. This is consistent with most of the data. However, the frequency-independence is violated in Event A, from 1985.9 to 1986.5, when the points at 4.8 GHz are separated from those at the other frequencies, as discussed in Section 4.1.
The nπ adjustments that convert Figure 1 into Figure 2a were made by hand. We also used an algorithm similar to that of Kiehlmann et al. (2016) that tests every point against the previous one, and adds ±nπ as necessary to keep the difference below 90 • . This is an automatic procedure that does not allow for any special considerations at a data gap. We did this for points at the different frequencies being treated separately, and also for all the points being used together. For the latter case, the results were similar to the non-automatic solution shown in Figure 2a. Figure 2b shows the flux density F of OJ 287 at the five frequencies. The MOJAVE values are for the core, but the others are total-flux measurements made with a single large dish. At most epochs OJ 287 has an "inverted" spectrum, with S 14.5 > S 4.8 , like many AGN (Kovalev et al. 1999;Fuhrmann et al. 2016). The 15.3 GHz flux densities for the core are usually well below the 14.5 GHz values for the total source, and show that the jet makes a substantial contribution to the total flux density. This is especially noticeable after 2010. Figure 2c shows the polarized flux density, PF, which will be important in the discussion of models for the EVPA rotations.
Tables 1, 2 and 3 contain all the points in the adjusted data sets, as shown in Figure 2. These tables contain 856 points at 4.8 GHz, 917 at 8.0 GHz, 1,207 at 14.5 GHz, and 93 at 15.3 GHz. The archival data can be reconstructed from Tables 1 and 2 by constraining each EVPA point to lie in the range 0 • to 180 • , by adding nπ as needed. Table 3 contains the 19 Kikuchi points at 9.0-10.5 GHz, found by digitizing the plots in Kikuchi et al. (1988). In this process the epochs differ slightly among the points for the EVPA, F, and P; and the mean epoch is shown in column 1 of Table 3.      Note. -Columns are as follows: (1) observation epoch, (2) Stokes I intensity in Jy/beam, (3) fractional linear polarization in per cent, (4) adjusted electric vector position angle in degrees. This table is published in its entirety in the electronic edition of the Jornal; a portion is shown here for guidance regarding form and content.  EVPA rotations can be spurious for two reasons: they can be generated by a random walk process, and they can be both generated and destroyed by statistical noise. This has also been discussed by e.g., Jones et al. (1985); Marscher (2014); Kiehlmann et al. (2016Kiehlmann et al. ( , 2017. In this Section we ask if these effects can be at work in our observations. We believe that they are not, because the probability of a random large double rotation with a reversal must be much smaller than the probability of a single rotation, but in OJ 287 we see three large reversals in 42 years, with no similar single rotations. In addition, the four reversals that we see are all in the same direction: CCW then CW. This alone reduces the probability that the rotations are random by an order of magnitude. Jones et al. (1985) first estimated the probability that a large EVPA rotation could be due to a random process. They considered a source that consisted of turbulent cells with random polarizations, and evolved the system by changing one cell per time step. With Monte-Carlo calculations, they found a rather high probability of a large rotation; with appropriate assumptions the probability of a rotation of 180 • or greater was as much as 0.3. For our purposes we need to multiply this estimate by the probability that the next rotation has a similar amplitude and the opposite sign, occurs shortly after the first one, and is isolated; i.e., there is no third rotation for a substantial period. This appears to call for a Monte-Carlo calculation, which is beyond the scope of this paper. However, it is clear that each of these factors will appreciably reduce the overall probability for the observed double rotations to arise by chance, compared to the probability for a single large rotation.

ARE THE ROTATIONS WITH
Another factor affecting the probability is that the rotations occur at the same time with independent observations at three frequencies. We have coincident events, and this greatly reduces the probability that they are due to random noise. But it may not reduce the probability that they are due to random walks, since the turbulent cells may be frequencyindependent. For this to be the case, however, opacity effects . EVPA for Event A. The step at 1986.52 is cosmetic, and the plot shows that there is a stable EVPA near -70 • that exists before and after the event. See text. must be negligible.

Event A
In this Section we describe the principal rotation events that are seen in Figure 2. We first discuss Event A, because it is bracketed by outbursts in total and polarized flux density, and this motivates the model we describe later. We then discuss, in order of complexity, events D, C, and B.
Our result for the EVPA of Event A is shown in Figure 4 and was obtained by following our two connection criteria: generally keep adjacent points less than 90 • apart except where there is a large time gap, and keep all frequencies on the same curve to the extent possible. To make the curve we first noted that the first five Kikuchi points, marked 17 • /day, form a steep line that is unambiguous, as are the 14 8.0 GHz points that are indicated with the line marked 1.8 • /day. The two lines fit to-gether well and define the main structure of the EVPA curve. The other points for 8.0 and 10 GHz then connect as shown. The points for 4.8 GHz show no evidence for the steep CW rotation seen at the other frequencies, and we dropped the requirement that the 4.8 GHz points had to fit in with the others. The 14.5 GHz points from 1986.2 to 1986.5 do not fit well with the others, and we placed them close to the 10 GHz line, since 10 GHz is the nearest frequency. This is arbitrary, and raising them by 180 • would place them close to the 4.8 GHz points. As we discuss later, in Section 5 in connection with the two-component model, these differences might result from the different behaviors of the polarized flux, at the different frequencies. Figure 4 shows that Event A had a CCW EVPA rotation of about 180 • , followed by a CW rotation of roughly 360 • . The EVPA before Event A was about -60 • , and roughly -260 • after it. But -260 • is the same as -80 • , and so we inserted a step of +180 • at 1986.52 for cosmetic purposes, to make it easy to see that the EVPA was approximately the same before and after the event. Figure 5 gives an extended view of Event A. with the three panels showing the EVPA, the flux density, and the linearlypolarized flux density. Two large outbursts in flux density, A1 and A2, bracket the EVPA event. They show the normal evolution of emission from an expanding synchrotron cloud, with lower frequencies delayed and reduced in intensity. A weak double rotation in EVPA, with a reversal, occurs at the same time as the peak of outburst A1. The strong EVPA Event A occurs during the tail of A1 and the rise of A2, and the reversal itself occurs at the time of the F and PF minimum between A1 and A2.
The polarized flux PF has a complex pattern in this interval. We only discuss the highest frequency, 14.5 GHz. PF has a strong peak, PF1, at 1985.4, at the time of the weak EVPA reversal, and it has a deep minimum at 1985.6, when the EVPA has almost returned to its baseline value. The PF has peaks PF2 and PF3 bracketing Event A, and is in a deep minimum through much of Event A. The CCW swing in EVPA from 1985.8 to 1986.1 occurs during the tail of A1, and the CW swing from 1986.1 to 1986.4 occurs during the rise of A2. Thus A1 itself, or at least its tail, is polarized with CCW rotation, and similarly the rise of A2 has CW rotation. The observed reversal in rotation occurs when A2 begins to dominate the total PF, at 1986.1. The deep minimum in PF at that time implies that the two components have EVPAs that are roughly 90 • apart.
The EVPA swings in Event A are of order 180 • or more, and they cannot be due to two variable sources with fixed EVPA, nor to the evolution of optical depth, since both of these give a maximum swing of 90 • . In Section 5 we present a model consisting of a steady polarized component combined with a variable component that has a rotating EVPA. If the components have similar amplitudes then when the EVPAs are nearly perpendicular the net PF will have a minimum, and the EVPA can have a rapid swing. This model can explain many of the observed features of the EVPA reversals in OJ 287.  Figure 6. It is bracketed by modest outbursts D1 and D2, similar to the way in which Event A is bracketed by A1 and A2. The EVPA has a rapid CCW swing in late 2000 at a rate ∼ 1.7 • d −1 ; the rate is not uniform. After the CCW swing the EVPA is nearly steady for about 1.5 years and then has another rapid swing, this time CW, at the rate ∼ −0.8 • d −1 . These rotations in EVPA occur during the rise of D1 and D2. As with Event A they are greater than 90 • and cannot be solely due to evolution of optical depth, or to a combination of two variable sources with fixed EVPA.

Event D Event D is shown in
The PF for Event D has a peak near 2001.75, in the middle of the steady period for the EVPA, and the PF has minima during the rapid EVPA swings. This is different from the behavior in Event A, where the EVPA reversal at 1986.1 occurs during a PF minimum. (See Figure 5). This will be discussed in terms of the two-component model in Section 8.
The arrows on the abscissa of Figure 6a correspond to the epochs for the images in Figure 3. The PF images show the core and one or two secondary jet components to the W or SW. These jet components are Nos. 1, 5, 4 and 11 in the MOJAVE list (Lister et al. 2013(Lister et al. , 2016 and are labeled C1, C5, C4 and C11 in Figure 3. 16 These four components are all superluminal, and C4 is the fastest one, with β app = 15. We are especially interested in C5, because it is intimately connected to Event D. C5 is moving at the rate µ = 0.54 ± 0.07 mas y −1 in the direction PA = −103 • . It is not moving radially but projects back close to the core, and was near the location of the core around 1999-2000, assuming that it was in uniform motion (Lister et al. 2016).
In Figure 6 both F and PF, at 14.5 GHz (green crosses in panels b and c), begin to increase around 2000.5. F rises to nearly 4 Jy by 2001.4, while PF continues to rise until nearly 2002. The PF image in Figure 3b shows that the PF rise is due to component C5, which dominates the image, and presumably first became visible around 2000.5 when the total  Figure 2 , 1988-1991. flux density started to increase. In this case the simultaneous increase in F is also due to C5, although this is not so obvious in the total flux image in Figure 3b.
The rise of the outburst D1 in Figure 6b has CCW EVPA rotation, but the subsequent decline has a steady EVPA, as seen in Figure 6a. D2 has CW rotation and the combination of the core and the two bursts starts to rotate CW when D2 starts to dominate the flux density. This happens around 2002.3. PF has a minimum then because, according to the model in Section 5.1, the sum of the Stokes vectors for the three components becomes small. At that time the phase of the sum, ξ, can sweep rapidly, and so the EVPA = ξ/2 also sweeps rapidly.

Event C
The events in Figure 6 start at 1995.5 with a CCW swing in EVPA of about 100 • , coincident with small bursts in F and PF. They are presumably due to the emergence of C1 from the core; C1 is seen 1.3 years later in Figure 3a. The epoch of this image is shown in Figure 6a with the arrow marked "a". The left-hand image in Figure 3a shows that C1 is highly polarized, but the right-hand image shows that the core, while weakly polarized, has more polarized flux density. The flux burst at 1996.0 does not show the common high-to-low frequency evolution, and we can ignore the possibility of EVPA changes due to optical depth effects. The 100 • EVPA swing could be due to a combination of variable sources that have fixed EVPA, but it could also be due to sources with a variable EVPA. In the scenario presented in Section 5 the new component C1, responsible for the flux burst at 1996.0, would have an EVPA with CCW rotation.
In 1997.2-1998.5 Event C has a CCW EVPA swing of about 160 • , followed by a CW swing of about 200 • . F and FP change little during the event, The large EVPA swing is about the same in the UMRAO and MOJAVE points, however, suggesting that the EVPA rotation is in the core. Note that the PF values from MOJAVE are very low, implying that the errors in EVPA are high, and so the individual MOJAVE EVPAs should be treated with caution.
The MOJAVE EVPA points in Event C, starting near 1996.0, were plotted in a different way by Cohen (2017), who did not show the reversal at 1998.5. This resulted from large gaps in the MOJAVE data; and without the closely-spaced points from UMRAO it is difficult to obtain the correct curve. The earlier work by Homan et al. (2002) also shows a different curve, and may have similarly suffered from the lack of closely-spaced points. Figure 7 shows an expanded view of the period 1988-1991. Again there is a sun gap, at 1989.6, that interferes with the interpretations. Event B includes two EVPA reversals, one near 1989.5 that is preceded by a shallow CCW rotation of about 100 • and followed by what is probably a steep CW rotation of at least 60 • . Unfortunately, the data are missing for this last rotation. The second reversal, near 1990.0, is symmetric. These two reversals have the sign that is common to all the reversals in OJ 287, CCW then CW. Detailed modeling is needed to investigate event B.

TWO-COMPONENT MODEL AND ROTATING STOKES VECTORS
In this Section we present a two-component model that can reproduce many of the polarization features seen in the preceding Sections. Two-component models have frequently been used to describe polarization events. Björnsson (1982) analyzed polarization changes due to relativistic aberration, and compared them to changes that can be produced by a non-relativistic two-component model. Holmes et al. (1984)  Our model has a steady component that we call the jet, and a time-dependent component that we refer to as the outburst. The outburst has a Gaussian shape with truncated tails, and an EVPA that rotates uniformly. The amplitude ratio of the two components, and the EVPA rotation rate and phase, are picked so that the results mimic some of the observations. A single rotation is modeled with one outburst. A double rotation, with a reversal, can be modeled with two successive outbursts that have opposite senses of EVPA rotation. The direct observation of outbursts A1 and A2, seen in Figure 5, motivates this model. Figure 8a shows our model for the case where the peak of the Gaussian outburst is weaker than the jet (Gaussian/jet = 0.8); and the other case, with the Gaussian stronger than the jet (Gaussian/jet = 1.2), is shown in Figure 8c. The relative size of the two components is important, for it controls the details of the EVPA of the combination. The Gaussians are truncated at t = −32 and t = 32. The rotation rate for the EVPA of the outbursts is +7.5 • per unit time step. The EVPA of the jet is 90 • from that of the outburst at its maximum, at t = 0. Figures 8b and 8d show the results of combining the two components. In both cases PF has a deep minimum where the Gaussian and the jet have similar amplitudes, and where their EVPAs are nearly perpendicular. But the resultant EV-PAs behave differently. In Figure 8b the EVPA curve of the combination has 6 extrema, or reversal points, at epochs indicated by the arrows on the abscissa. There is a weak maximum at point a at t = −27.2, then a shallow CW swing to point b at t = −17.1, where the rotation direction reverses, a CCW swing to point c at t = −2.3, then another reversal and a rapid CW swing to point d at t = +2.3, where the process repeats in reverse. In Figure 8d the EVPA is similar to that in Figure 8b at early and late times, but is continuously CCW for −17 < t < +17. The total EVPA rotation in Figure 8d to h, is 199 • . These changes are most easily understood with the Stokes Parameters. Figure 9 shows the Stokes plane for Figure 8. Here we adopt the IAU recommendations for the sign of the Stokes parameters (Hamaker & Bregman 1996) with the Stokes plane overlaid on the sky plane; Q increases to the North, U increases to the East, and the Stokes angle ξ = tan −1 (U/Q) = 2 × EVPA.
In Figure 9 the vertical arrow labeled "jet" represents the steady jet, which is the same in Figures 8a and 8c. The Stokes vectors representing the outbursts are added to the jet vector, to form the sum vectors, as shown at time t=2 for the weaker outburst in Figure 8a. As time advances, the EVPA of the outburst rotates CCW, and the sum vector traces out the inner loop. The loop is parametric in time, and the times a-f on the loop are EVPA reversal points that can be seen in Figure 8b. At the reversal points the vectors are tangent to the loop. The total excursion of ξ (between points c and d) is 102.6 • ; the EVPA excursion, seen in Figure 8b, is 51.3 • .
When the peak of the outburst is stronger than the jet, as in Figure 8c, the loop encloses the origin, as shown by the outer loop in Figure 9. In this case the sum vector rotates continuously CCW between points g and h. The full excursion of ξ is 398 • , and the corresponding EVPA rotation in Figure 8d is 199 • . This striking difference in EVPA rotation, caused by the relative size of the jet and the outburst, could be responsible for the differences in EVPA behavior seen in Figure 5.  GHz, with a small rotation at 4.8 GHz, as seen in Figure 5a.

Double Rotation with a Reversal
A double rotation with a reversal can be obtained with two successive outbursts, with opposite senses of rotation. Figure 10 shows an example where both outbursts are stronger than the jet. In panel (a) the two outbursts are the same as in Figure 8c but with opposite rotations, and they are separated by ∆(t) = 32. The sum in (b) has some similaritites to Event D, seen in Figure 6. In both Event D and in the model (Figure 10b) the EVPA has rapid swings of order 180 • , and PF has a smooth top and deep minima centered near the EVPA swings.
The Stokes plane plot for the model in Figure 10b is in Figure 10c. The early part of this diagram is the same as the corresponding part of Figure 9. When the second outburst becomes appreciable, at t ∼ 2, the loop opens out and, at the star, where the second outburst begins to dominate the amplitude, the loop reverses and goes CW back along the same track. This motion gives the flat-top amplitude in Figure 10b, and the steep-sided EVPA curve.
Note that in Figure 10b the central parts of the EVPA swings are much steeper than the linear EVPA curves for the two outbursts. In the context of models where the synchrotron source rotates around the jet axis (Section 7), this means that the physical rotation rate can be much less than the apparent rate, seen as the rapid change in EVPA. It is likely that relativistic effects also affect the apparent rotation (Björnsson 1982). In Figure 11 we show the two sides of Event D on the Stokes plane, for 14.5 GHz. Both of these loops are like the outer loop in Figure 9 in that they enclose the origin. Hence the EVPA swing for each is of order 180 • . In Figure 11a the jumble of points near the star contains both the beginning and the end of the swing. The circled point is at 2000.94 and has ξ = 97 • or EVPA = 48.5 • . The polarization is exceptionally high for this point, but there is no reason to exclude it as an outlier. Because of it we can claim that the loop surrounds the origin, and that the swing at 2001 is of order 180 • . On the CW side, shown in Figure 11b, there are more points defining the loop, and, crucially, we see that in Figure 6a the swing is well-defined when data from all the frequencies are included.

Stokes Plot for Event D
It is easier to study the EVPA with a time series like that in Figure 6a rather than with loops on the Stokes plane, because time is not uniform on the loop. Further, it would be difficult to plot all the frequencies together on the Stokes plane, because they would need to be normalized in some way to the frequency spectrum of the polarization. Still, the Stokes plot is useful in visualizing how polarized radiations combine, and in arguing that, because the loop encloses the origin, the EVPA rotation really is 180 • or more (Villforth et al. 2010).

A SIMPLE GEOMETRY
In the preceding Section we modeled the EVPA rotation reversals with a pair of outbursts whose EVPAs are counterrotating. We now present a geometric model that can generate these counter-rotating outbursts, in a simple and intuitive way.
Consider a plasma jet with a relativistic flow, with a helical magnetic field. Let a disturbance generate a sub-relativistic shock pair, a forward shock traveling downstream and a reverse shock traveling upstream, each with β jet sh = 0.1 in the jet frame. Let the Lorentz factor of the jet be Γ gal jet = 10 in the frame of the host galaxy. Then in the galaxy frame both shocks are moving forward relativistically, with Lorentz factors Γ gal fwd = 11.05 and Γ gal rev = 9.05. An observer on axis sees Doppler shifts of 22.1 and 18.1 for radiation from these shocks, and if they have similar synchrotron sources then their flux density ratio is about 1.5. The observed radiation from the reverse shock is not substantially weaker than that from the forward shock.
Let the magnetic field lines have the structure of a righthand helix. If the shocks travel along this helix then the forward shock is seen to rotate CCW and the reverse shock (moving upstream in the jet frame) is seen to rotate CW. With appropriate synchrotron sources whose aspect to the axis is fixed, the EVPA rotation will follow that of the shocks (Marscher et al. 2008). In this geometry the observed rotation automatically reverses when the second shock becomes the dominant source for the polarized flux density. The righthandedness is required by the observed sense of the EVPA reversal, CCW then CW.
In this model we have implicitly assumed that the plasma jet is not rotating, so that the upstream shock is not carried into CCW rotation as seen by the observer. But the jet is moving forward and cannot cross the magnetic field. Hence, the helical field must be rotating at a rate such that the screw action drives the plasma straight forward. The non-relativistic condition for this is β = ΩR cosα where β is the longitudinal jet velocity in units of c, Ω is the angular rotation rate of the magnetic field in y −1 , R is the radius of curvature of the field in ly, and α is the pitch angle. But we have a relativistic flow and must be concerned with the velocity-of-light cylinder around the axis. The situation here is similar to that in a pulsar atmosphere, where the field lines bend backward and the toroidal component of the field slips through the plasma (Meier 2012). In this way the helical field continues across the light cylinder while the plasma velocity stays below c.
We emphasize that we have proposed here a purely geometric model, and the nature of the shock waves is not specified, nor is the mechanism by which the source is guided by the helical field, and why the EVPA itself stays fixed with respect to the helix. In the next Section we describe a physical model that has many of the features of the geometric model, and sug-gest that it may explain the observations. 7. MODELS USING HELICALLY MAGNETIZED JETS 7.1. Sub-fast, Super-slow Magnetosonic Jet Models Nakamura (2001) and Nakamura & Meier (2004) simulated 1.5-D and 3-D helically magnetized jets whose flow speed was slower than the jet's internal MHD fast-mode magnetosonic wave speed, V jet < V fast ≈ (V 2 A + V 2 s ) 1/2 , where V A is the internal jet Alfvén speed and V s is the internal sound speed (with V A > V s ), but the jet was faster than the internal slowmode wave speed (V jet > V slow ≈ V s ). Furthermore, the jet flow speed also was greater than the fast-mode magnetosonic wave speed in the material into which the jet was flowing. 17 Jets like this, with a sub-magnetosonic internal Mach number but super-magnetosonic external Mach number, develop three shocks in the flow: a forward fast-mode shock (FF), a forward slow-mode shock (FS), and a reverse slow-mode shock (RS). Furthermore, because of conservation of the combined plasma and magnetic field angular momentum at the FF shock, the material between the FF and FS shock has an enhanced (compressed) helical magnetic field strength and a rotation velocity significantly greater than the rotation rate of the main jet body near the contact discontinuity and which also exceeds V slow . Therefore, if the supersonically-rotating plasma in the FS/FF region develops a non-axisymmetric shock feature (e.g., near the FS shock itself), then an observer viewing this jet end-on would observe synchrotron emission from that feature that exhibited a physical rotation about the line of sight of perhaps several radians.
Thus, a sub-fast, super-slow helically magnetized jet could be a promising model for sources that exhibit a single, onedirectional rotation of the EVPA. However, such jets do not produce the double rotations, with reversals, seen in OJ 287.

Super-fast Magnetosonic Jet Models
On the other hand, Nakamura et al. (2010) and Nakamura & Meier (2014) performed similar 1.5-D simulations of helically magnetized jets, but whose flow speed this time was greater than the jet's internal fast-mode magnetosonic wave speed. These jets developed four shocks in the flow: FF, FS, RS, and also a reverse fast shock (RF). Figure 3d of Nakamura et al. (2010) shows that, in the galaxy frame, the toroidal component of magnetic field is substantially enhanced between the FF and FS shocks, and also between the RS and RF shocks. This leads to two moving synchrotron sources. Further, Figure 3e of this paper shows that an azimuthal motion of the plasma is established between the FF and FS shocks, and between the RS and RF shocks, but that the sense of rotation is opposite in the two regions. Thus, two oppositely-rotating synchrotron-emitting regions are established, moving relativistically downstream because V jet > V fast . If the emission regions are not axisymmetric, then an observer near the axis will see the EVPAs of the two outbursts rotate in opposite directions.
This super-magnetosonic jet model is a good candidate to explain the observations of the EVPA reversals seen in OJ 287. It produces the main feature used in constructing Figure 10, namely, the two oppositely-rotating emission regions  (Kikuchi et al. 1988;D'Arcangelo et al. 2009;Villforth et al. 2010;Blinov et al. 2015) have reported optical and IR observations of OJ 287 that are closely-enough spaced in time to be useful for studying EVPA rotations. Two of them (Kikuchi et al. 1988;D'Arcangelo et al. 2009) also show that the optical and radio variations of polarization are nearly synchronous. The observations of Kikuchi et al. (1988) Figure V1 shows the R band flux density, the polarized flux density PF, and the EVPA for OJ 287, during 2004OJ 287, during .9-2009. The EVPA has points restricted to 0 • -180 • , and the figure is analogous to our Figure 1. It shows a number of rotations, several of which are described in detail and are shown with Stokes vector plots like the one in Figure 11 for Event D. The event in April 2006 is shown expanded in Figure V2 and the Stokes plot is in Figure V16. 18 The loop in Figure V16 does not enclose the origin, and the EVPA swing as seen in Figure V2 appears to be about 50 • . A better view of this event is in Figure D5, where the swing is seen to be ≈ 45 • CW. In this Figure the peak of PF is in the middle of the EVPA swing. This is the opposite of what happens in Event D, seen in Figure 6, where PF has minima in the middle of the EVPA swings. An easy way to accommodate this difference, and keep the optical result within the two-component model, is to shift the phase of the EVPA. Figure 12 shows the Stokes plane for the model in Figures 8c and 9 when the EVPA of the jet is rotated by 90 • ; i.e. the Stokes vector for the jet is reversed. Although the outburst is stronger than the jet, the loop does not enclose the origin, and the CCW swing of the sum, from j to k, is 116.7 • . Figure 12 is analogous to Figure V16. Both have the maximum PF in the middle of the EVPA swing. However, any physical significance attached to the shift of the jet EVPA relative to the outburst will depend on the specific model used to describe the event.
From their Stokes plane plots like that in V16, Villforth et al. (2010) suggest that OJ 287 has two components of emission that generate the EVPA rotation, the "optically polarized core" (OPC) and the "chaotic jet emmision". Our model corresponds closely to this; the jet corresponds to the OPC, and their chaotic jet emission corresponds to our outbursts that are superposed on the jet. When, as in Figures V16 and V21, they show a loop on the Stokes plane, it is their chaotic jet emission that has a systematic CCW swing in EVPA. In Figure V16 the OPC is roughly lined up with the maximum PF as, in Figure 12, the vector for the jet is aligned with the maximum PF. This relationship also holds for Figure  V21 and the other Stokes plane plots in Villforth et al. (2010). Thus, we see that the two-component model that we use for the radio observations may be useful for the optical data also.
Villforth et al. (2010) describe rotations in the EVPA of OJ 287 at optical wavelengths, but these are all single rotations, not double with a reversal like those we have seen at radio wavelengths. However, the EVPA data in Figure V1 are restricted to 0 • to 180 • , and a double rotation of order 180 • might not be recognized if it did exist there. It would be useful to smooth the data in Figure V1 (Blinov et al. 2015) and J1512-0905 (PKS 1510-089) (Blinov et al. 2016). J1512-0905 was also studied at R-band by Beaklini et al. (2017). There was no overlap in these observations, but the Blinov et al observations ended with a strong CW rotation and the Beakilini et al observations started about 20 days later with a strong CCW rotation. This appears to show a double rotation of about 200 • , with a reversal. Figure 10 of Beaklini et al. (2017) shows a plot of the EVPA of J1512-0905 that combines the data from a number of observers. 9. DISCUSSION 9.1. Time Scales for EVPA Rotation We have found various rotation rates in OJ 287, from the fastest, 17 • d −1 in the CW swing in Event A (Figure 4), to the over-all long-term trend of roughly 90 • in 30 y (Figure 2). The reciprocal of a rate is a time scale, which we take here to be the time to rotate by one radian. Thus our time scales run from 3.3 days to 30 years. On the short end, measureable time scales are limited by the sampling interval, which for Kikuchi et al. (1988) is one day.  have argued that for reliable recovery of the intrinsic time scale, the sampling interval should not exceed ∼ 30% of the intrinsic scale; i.e. we need at least three samples per intrinsic time scale. Thus 17 • d −1 is about the fastest rate that Kikuchi et al. (1988) could have reliably determined. The UMRAO points are typically 3 days apart and so 5 • d −1 is about the fastest that can be found in the UMRAO data. The MOJAVE points are a few weeks to several months apart and the most rapid CCW and CW swings in the rotation reversal events cannot be determined reliably from the MOJAVE data alone. This has already been noted in Section 2. The CCW rate in Event D is about 1.8 • d −1 , and its detection requires a sampling interval of 10 days or less.
To find a time scale in the frame of the jet, τ jet , we multiply the observed time scale τ o by δ/(1 + z), where δ is the Doppler factor of the radio source. There are two different values for δ in the literature, and we refer to them as "early" and "recent". Two early values are δ = 17.0 (Hovatta et al. 2009) and δ = 18.9 ± 6 (Jorstad et al. 2005), and they are derived from 43-GHz flares in 2003 and 1998-2000, respectively. Two recent measurements both give δ = 8.7 (Jorstad et al. 2017;, and they are derived from mm-wave flares after 2007. Apparently, δ changed around 2005; why? We suggest that there was a change in the direction of the inner jet. Agudo et al. (2012) showed that the PA of the inner jet jumped in 2004 at 43 GHz; and at 15 GHz there was a similar change in 2006 (Cohen 2017). Presumably, this PA jump reflects an increase in the viewing angle to the jet, and as a consequence the Doppler factor changed by a factor of about 2. Since all the major rotation/reversal events at 15 GHz that we see in Figure 2 happened prior to 2006, we use the early value, δ ≈ 17 in the following discussion. However, the use of the lower value would not make a substantive change.
With δ ≈ 17, τ jet /τ o ≈ 13. The fastest swing, 17 • d −1 now becomes the shortest time scale in the jet, τ jet,min ≈ 44 d. This time scale may be too short to represent a shock circulating around a helix, as it would make the radius of curvature of the field line much less than 1 ly. Our preferred explanation, using the model in Section 5, has no analogous velocity-of-light limit. As seen in Figure 9, the rotation speed of the resultant Stokes vector can be very high, when the amplitudes of the two components are nearly the same.
Longer timescales are seen in the slowly-changing EVPA baseline in Figure 2a; e.g., in 1993 where the apparent rate is about 20 • y −1 , or τ jet ∼ 50 y. The changes in 2005 -2012 are coincident with changes in the orientation of the inner jet, as defined by the appearance of a new superluminal component (Cohen 2017).

Outbursts Without Rotations
In Figure 5 Event A appears to be associated with the strong outbursts A1 and A2 in flux density, but when we look at In the two-component Gaussian model, a large rapid rotation is only seen when the conditions are right; the outburst must be stronger than the jet, and the EVPA phase and rotation rate must be appropriate. On the other hand, outbursts without large rotations may simply show that there is more than one cause for the outbursts. The intervals A-C and B-D are both roughly 12 years. This is interesting, because the optical flares that match a binary black-hole model have a period of about 12 years (Valtonen et al. 2011). The top axis of Figure 2 has six bars that indicate the epochs of the flares. The epochs are 1983.0, 1984.2, 1994.8, 1996.0, 2005.8, and 2007.7 (Sillanpää et al. 1988Valtonen et al. 2006;Sillanpää et al. 1996a,b;Valtonen et al. 2008a,b). These are the original references reporting the flares, except for 1984.2, where the reference is to the compilation in Valtonen et al. (2006). Additionally, a strong flare was seen at 2015.9 after having been predicted (Valtonen et al. 2016); showing that the model closely matches the observations. We also note that light curves for OJ 287 are highly variable, and that an analysis of 9.2 years of well-sampled optical data yielded evidence for quasi-periodic oscillations of periods ∼ 400 and ∼ 800 days (Bhatta et al. 2017).
In Figure 2 the optical pairs in 1983-84 and 1994-96 precede the radio rotation pairs A-B and C-D by about three years. However, there are no radio events corresponding to the optical pair in 2005-07, and this suggests that the radiooptical 12-year similarity is a coincidence.

SUMMARY AND CONCLUSIONS
We report what appears to be a new phenomenon, rotation reversals of the EVPA at radio frequencies, in the BL Lac object OJ 287. These consist of a large ∼180 • CCW rotation followed by a similar CW rotation. Three of these events were seen in 40 years, and a fourth, smaller one, was also seen. They were all in the same direction, CCW followed by CW. We suggest that a rotation can be explained with a twocomponent model consisting of an outburst superposed on a steady jet, with the EVPA of the outburst rotating steadily in time. This reproduces many of the observed features of a rotation. A three-component model consisting of two successive outbursts with oppositely rotating EVPAs, together with the jet component, explains the reversals. This model is also applicable to polarization rotations seen at optical wavelengths.
In a more physical model, we consider that the reversals take place in a super-magnetosonic jet; i.e., one in which the bulk speed of the plasma is greater than the speed of the fast magnetosonic wave. The jet is threaded by a helical magnetic field. We use the mechanism analyzed by Nakamura et al. (2010) and Nakamura & Meier (2014) that produces four MHD waves, forward and reverse fast and slow magnetosonic waves. Between the forward fast and slow waves the toroidal component of magnetic field is compressed; this increases the angular momentum of the field, and to conserve angular momentum the plasma rotates around the axis in the opposite direction. This happens also to the reverse fast and slow pair of magnetosonic waves, but the rotation is in the opposite sense. This forms two regions of enhanced plasma density and magnetic field, rotating in opposite directions. Both regions move relativistically downstream because V jet > V fast . The resulting synchrotron radiation, as seen by an observer near the axis, consists of two outbursts that have oppositely rotating EVPAs. The observed rotation sense, CCW followed by CW, requires a right-hand helix.
We conclude that our observations of EVPA reversals provide evidence for a strong helical magnetic field in OJ 287. This is consistent with the observations and conclusions of many others e.g. Cohen et al. (2015); Gomez et al. (2016); Motter & Gabuzda (2017). The observations also provide evidence that the jet of OJ 287 is super-magnetosonic, and this can provide a constraint on B 2 /n, where B is the strength of the magnetic field and n is the particle density.