Issue 
A&A
Volume 590, June 2016



Article Number  A10  
Number of page(s)  20  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/201527725  
Published online  28 April 2016 
Polarization angle swings in blazars: The case of 3C 279^{⋆}
^{1}
Aalto University Metsähovi Radio Observatory,
Metsähovintie 114, 02540
Kylmälä,
Finland
email:
sebastian.kiehlmann@aalto.fi
^{2}
Aalto University Department of Radio Science and
Engineering, PL
13000, 00076
Aalto,
Finland
^{3}
MaxPlanckInstitut für Radioastronomie,
Auf dem Hügel, 69, 53121
Bonn,
Germany
^{4}
Institute for Astrophysical Research, Boston
University, 725 Commonwealth
Avenue, Boston,
MA
02215,
USA
^{5}
Astronomical Institute, St. Petersburg State University,
Universitetskij Pr. 28, Petrodvorets, 198504
St. Petersburg,
Russia
^{6}
Astro Space Center of Lebedev Physical Institute,
Profsoyuznaya 84/32,
117997
Moscow,
Russia
^{7}
Sternberg Astronomical Institute, M.V.Lomonosov Moscow State
University, Universiteskij prosp. 13, 119991
Moscow,
Russia
^{8}
Institute of Astronomy, Astrophysics, Space Applications and
Remote Sensing, National Observatory of Athens Vas. Pavlou & I.
Metaxa, 15 236
Penteli,
Greece
^{9}
Department of Physics and Astronomy, University of New
Mexico, Albuquerque
NM, 87131, USA
^{10}
Main (Pulkovo) Astronomical Observatory of RAS,
Pulkovskoye shosse, 60,
196140
St. Petersburg,
Russia
^{11}
Instituto de Astrofísica de Andalucía,
CSIC, Apartado 3004,
18080
Granada,
Spain
^{12}
Hiroshima Astrophysical Science Center, Hiroshima
University, HigashiHiroshima,
7398526
Hiroshima,
Japan
^{13}
Instituto de Astronomía,Universidad Nacional Autónoma de
México, 04510
México DF,
México
^{14}
Tuorla Observatory, Department of Physics and Astronomy,
University of Turku, Väisäläntie
20, 21500
Piikkiö,
Finland
^{15}
Foundation for Research and Technology – Hellas, IESL,
Voutes, 71110
Heraklion,
Greece
^{16}
Department of Physics and Institute for Plasma Physics, University
of Crete, 71003
Heraklion,
Greece
^{17}
Crimean Astrophysical Observatory, P/O Nauchny, 98409
Crimea,
Russia
^{18}
Special Astrophysical Observatory of the Russian AS,
Nizhnij Arkhyz, 369167
KarachaevoCherkesia,
Russia
^{19}
Southern Station of the Sternberg Astronomical Institute, Moscow
M.V. Lomonosov State University, P/O Nauchny, 298409
Crimea,
Russia
^{20}
Department of Physical Sciences, Hiroshima
University, HigashiHiroshima,
7398526
Hiroshima,
Japan
^{21}
ZAH, Landessternwarte Heidelberg, Königstuhl 12,
69117
Heidelberg,
Germany
^{22}
Instituto de Astronomía, Universidad Nacional Autónoma de
México, 22860
Ensenada BC,
México
^{23}
Abastumani Observatory, Mt. Kanobili, 0301
Abastumani,
Georgia
^{24} Inst. of SolarTerrestrial Physics, Lermontov st. 126a,
Irkutsk p/o box 291, 664033, Russia
^{25}
Center for Backyard Astrophysics – New Mexico, PO Box 1351
Cloudcroft,
NM
88317,
USA
^{26}
Engelhardt Astronomical Observatory, Kazan Federal
University, Tatarstan, Russia
^{27}
Department of Physics and Astronomy, University of
Southampton, Southampton, SO17
1BJ, UK
^{28}
Kavli Institute for the Physics and Mathematics of the Universe
(Kavli IPMU), The University of Tokyo, 515 KashiwanoHa, 2778583
Kashiwa City Chiba,
Japan
^{29}
Finnish Centre for Astronomy with ESO (FINCA), University of
Turku, Väisäläntie
20, 21500
Piikkiö,
Finland
^{30} Department of Physics, University of Colorado Denver, CO,
USA
^{31}
Steward Observatory, University of Arizona,
Tucson, AZ
85721,
USA
^{32}
Lowell Observatory, Flagstaff, AZ
86001,
USA
Received: 10 November 2015
Accepted: 29 February 2016
Context. Over the past few years, on several occasions, large, continuous rotations of the electric vector position angle (EVPA) of linearly polarized optical emission from blazars have been reported. These events are often coincident with high energy γray flares and they have attracted considerable attention, since they could allow us to probe the magnetic field structure in the γray emitting region of the jet. The flatspectrum radio quasar 3C 279 is one of the most prominent examples showing this behaviour.
Aims. Our goal is to study the observed EVPA rotations and to distinguish between a stochastic and a deterministic origin of the polarization variability.
Methods. We have combined multiple data sets of Rband photometry and optical polarimetry measurements of 3C 279, yielding exceptionally wellsampled flux density and polarization curves that cover a period of 2008–2012. Several large EVPA rotations are identified in the data. We introduce a quantitative measure for the EVPA curve smoothness, which is then used to test a set of simple random walk polarization variability models against the data.
Results. 3C 279 shows different polarization variation characteristics during an optical lowflux state and a flaring state. The polarization variation during the flaring state, especially the smooth ~360° rotation of the EVPA in mid2011, is not consistent with the tested stochastic processes.
Conclusions. We conclude that, during the two different optical flux states, two different processes govern polarization variation, which is possibly a stochastic process during the lowbrightness state and a deterministic process during the flaring activity.
Key words: polarization / galaxies: active / galaxies: jets / quasars: individual: 3C 279
The measured and processed optical polarization and Rband photometry data are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/590/A10
© ESO, 2016
1. Introduction
Blazars, a class of active galactic nuclei (AGNs) containing relativistic jets of magnetized plasma that are pointed close to our line of sight, are known to be highly variable in the optical polarization fraction and electric vector position angle (EVPA) since the early days of polarimetric observations of quasars (e.g. Kinman 1967). Monotonic optical EVPA rotations spanning days to months and with amplitudes of up to several full cycles have been observed, for example by Kikuchi et al. (1988), Larionov et al. (2008), and Marscher et al. (2010). However, the nπambiguity of EVPA makes it difficult to detect large continuous rotations (≥180°) in sparsely sampled polarization curves. Therefore, a significant monitoring effort with several participating observatories is required to find and study these events.
During the past few years several optical EVPA rotations have been reported as coinciding with flaring activity in γray emission, as observed by the Fermi γray Space Telescope at GeV energies and atmospheric imaging Cherenkov telescopes at TeV energies; for instance the ~240° rotation in BL Lac in 2005 (Marscher et al. 2008), in PKS 1510089 the counterclockwise ~720° rotation in 2009 (Marscher et al. 2010) and the counterclockwise ~380° and clockwise ~250° rotations in 2012 (Aleksić et al. 2014b), and the rotation in W Comae (Sorcia et al. 2014). These coincident events include rotations in the bright γray quasar 3C 279 (Abdo et al. 2010; Aleksić et al. 2014a). 3C 279 has exhibited large EVPA rotations on both short (days) and long (months) timescales and with opposite senses of rotation (Larionov et al. 2008; Abdo et al. 2010). The coincidence of EVPA rotations and γray flares has raised considerable interest towards optical polarization monitoring, since it may provide clues about the origin of the high energy emission in blazars.
There are several open questions regarding the observed EVPA swings: What physical processes produce the rotations of the EVPA? Is it the same process for all events, in all objects, or different processes for different events, even in the same object? And is the γray flaring activity physically connected to the rotations of the optical polarization angle or are those events just coincidences (Blinov et al. 2015)? Various models have been proposed to explain the EVPA rotations or, more recently, the potential connection between the EVPA rotation and γray flares.
These models can be divided into two classes: deterministic and stochastic models. The deterministic models include, for example, the superposition of two emission components with different polarization characteristics (Holmes et al. 1984), the shock compression of an ordered helical magnetic field in an axisymmetric, straight jet (Zhang et al. 2014, 2015), and several models that are based on nonaxisymmetric structures in the jet. The latter group of models includes, for example, a global bend in the jet, where a change of the viewing angle Δθ can produce a rotation of the EVPA Δθ<χ< 180° owing to relativistic aberration (e.g. Björnsson 1982; Nalewajko 2010). Abdo et al. (2010) explained the 2009 swing in 3C 279 by a bent jet and Aleksić et al. (2014a) used this same model for the swing in 2011. Nonaxisymmetric magnetic field configuration was invoked to explain EVPA swings in BL Lacertae objects by Konigl & Choudhuri (1985), while Kikuchi et al. (1988) introduced a qualitative model of a shock passing through a helical magnetic field with a changing pitch angle. Recently, Marscher et al. (2008) proposed a model containing an emission feature, which does not fill the entire jet crosssection, on a helical trajectory that probes different parts of a largescale helical magnetic field as it moves along the jet. This model has also been applied to the 2006/2007 EVPA swing in 3C 279 (Larionov et al. 2008).
The second class of models is based on a stochastic variation of polarization parameters of multiple cells (e.g. Jones et al. 1985; Jones 1988; D’Arcangelo et al. 2007; Marscher 2014). These studies have shown that it is easy to produce random EVPA changes that appear as a rotation of several hundred degrees. The two most recent models by Zhang et al. (2014, 2015) – in the class of deterministic models – and Marscher (2014) – in the class of stochastic models – in particular explore the connection between EVPA rotations and γray flares.
One key argument in distinguishing between deterministic and stochastic polarization variation is the smoothness of a continuous EVPA rotation (e.g. Marscher et al. 2008). Deterministic models should produce smooth EVPA variation whereas stochastic models are expected to produce more erratic EVPA curves. Here we develop a method, based on this kind of a quantitative measure of a curve smoothness, to test the stochastic models of polarization variability and apply it to the optical polarization data set of 3C 279 that was collected during our intensive multiwavelength campaign in 2010–2012.
In Sect. 2, we describe our almost three years of polarimetry data. We introduce the quantitative measure of the EVPA curve smoothness and analyse the polarization variation in Sect. 3. In Sects. 4 and 5 we test three simple random walk processes against the observations, and these results are discussed in Sect. 6. A summary and conclusions are presented in Sect. 7. This publication is based on the Ph.D. Thesis “Origin of the γray emission in AGN jets – A multiwavelength photometry and polarimetry data analysis of the quasar 3C 279” (Kiehlmann 2015).
2. Data
An intensive VLBI and multiwaveband monitoring campaign targeted on quasar 3C 279, including observatories covering radio, millimetre, infrared, optical, ultraviolet, Xray, and γray bands, was carried out in 2010–2012 (Quasar Movie Project; Kiehlmann et al., in prep.). For the period covered by the campaign we have accumulated Rband photometry and optical (R and V filters, spectropolarimetry at 5000–7000 Å, and white light) polarimetry data from a number of existing blazar monitoring programs and from an ad hoc campaign specifically targeted on 3C 279. Furthermore, we include in our analysis some data taken before the campaign period by the longterm monitoring programs, so that the data presented in this paper cover a time range of JD 2 454 790–JD 2 456 120. The observations were performed by (1) the 70 cm telescope of Abastumani Astrophysical Observatory (Mount Kanobili, Georgia); (2) the 2.2 m telescope of Calar Alto Observatory (Almería, Spain); (3) the 70 cm AZT8 telescope of the Crimean Astrophysical Observatory (CrAO; Nauchnij, Russia); (4) the 1.5 m KANATA telescope of the HigashiHiroshima Observatory (Hiroshima, Japan); (5) the 35 and 60 cm Kungliga Vetenskapakademien (KVA) telescopes; and (6) the 2 m Liverpool telescope of the Observatorio del Roque de Los Muchachos (La Palma, Canary Islands, Spain); (7) the 1.83 m Perkins telescope of Lowell Observatory (Flagstaff, Arizona, USA); (8) the 1 m telescope of the Special Astrophysical Observatory of the Russian Academy of Sciences (SAO RAS; Nizhny Arkhyz, Russia); (9) the 1.3 m SMARTS telescope of the Cerro Tololo InterAmerican Observatory (Chile); (10) the 84 cm telescope of the Observatorio Astronómico Nacional (San Pedro Mártir, Mexico); (11) the 40 cm LX200 telescope of St. Petersburg State University (St. Petersburg, Russia); and (12) the 1.54 m Kuiper and 2.3 m Bok telescopes of Steward Observatory (Mt. Bigelow and Kitt Peak, Arizona, USA). In addition, data were also gathered in the observing campaign of the American Association of Variable Star Observers (AAVSO), and by individual observers (T. Krajci^{1} using a 36 cm telescope and A. Sadun using the robotic 50 cm New Mexico Skies (NMS) telescope 11). The participating observatories are listed in Table 1.
Observatories contributing optical polarization data to the Quasar Movie Project with telescope diameters, filters, number of polarization data points, and the type of data with references to the data reduction in the table footnotes.
Table 1 also includes, for most of the participating telescopes, references to the descriptions of the programs and the used data reduction procedures. In the following we shortly comment on the reduction of data from the rest of the telescopes. Observations with the 1 m telescope of SAO RAS, the 70 cm AZT8 telescope of CrAO, the 50 cm telescope of NMS, and by T. Krajci were made through colour filters intended to implement the spectral response curve close to the one of Cousins (1976)R_{C} system. Various software implementations of the differential aperture photometry technique were applied to measure the brightness of 3C 279 on CCD images. Based on the seeing conditions an appropriate circular aperture size was chosen by each observer. The observers were advised to use the standardized list of comparison stars maintained by LSW Heidelberg^{2}, although this was not always possible. Some observers used ensemble photometry while others employed a single comparison star. In particular, the 1 m SAO and the 70 cm AZT8 CrAO telescope images were reduced (see Doroshenko et al. 2005) using the star N9 from GonzálezPérez et al. (2001).
We combined all Rband photometry into one light curve after crosscalibrating all data sets with respect to the SMARTS Rband data to compensate for systematic offsets between different instrument and telescope combinations and we converted magnitudes into spectral flux densities (Kiehlmann et al., in prep.). We removed one data point from the KVA polarimetry data, identified as an outlier (>60° intraday jump in EVPA) based on the combined optical EVPA curve. Four data points in San Pedro Mártir (SPM) polarimetry data were identified as outliers based on the simultaneous EVPA jumps in both 3C 279 and 3C 273 which was observed together with 3C 279 in the framework of the Quasar Movie Project. Additionally, we removed five data points from the St. Petersburg, four points from the CrAO, and one point from the Liverpool telescope datasets. These points show large uncertainties and offsets from the combined optical polarization fraction curve, that are not typical of 3C 279 on comparable timescales. We corrected the polarization fraction curves for the statistical (Rician) bias following Pavlidou et al. (2014). In optical the polarization fraction and the EVPA depend only weakly on the frequency (see panels c and d in Fig. 1). Therefore, we combine all optical polarization curves. We follow the convention that the EVPA is measured counterclockwise from North. An increase of the polarization angle refers to a counterclockwise rotation projected on the sky.
Before analysing the data we averaged data points within 0.5 days in the Rband light curve and in the combined optical polarization curves. The averaging intervals are selected iteratively instead of using fixed time bins, taking into account the uneven time sampling of the data. The polarization data is converted to Stokes parameters before averaging and converted back afterwards. The Stokes parameters are averaged weighted with the uncertainties. The variability on timescales smaller than half a day is of the order of the measurement uncertainties of individual data points. The averaging thus reduces the noise without removing significant real variations. The Rband light curve and the combined optical polarization curves are shown in Fig. 1, panels b–d.
3. Polarization analysis
Before the description and analysis of the optical light curve and the polarization curves we discuss a solution to the nπambiguity of the polarization angle and introduce a measure to quantify the smoothness of a curve.
Fig. 1 Optical photometry and polarimetry and γray light curve of 3C 279. FermiLAT γray light curve at > 100 MeV binned into 3 day intervals (panel a)) as published in Hayashida et al. (2015). Combined Rband light curve (panel b)). Measured, optical polarization fraction (panel c)) and EVPA (panel d)); red circles: Calar Alto (R), red squares: CrAO70 cm (R), red diamonds: Perkins (R), orange upsided triangles: SPM (R), orange rightsided triangles: St. Petersburg (R), green downsided triangles: KANATA (V), green leftsided triangles: Steward Obs. (spec. and V), blue circles: Liverpool (V+R), blue squares: KVA (white light). Combined, debiased, and averaged polarization fraction (panel e)). Combined, averaged, and adjusted EVPA (panel f)); open symbols are added from the nonaveraged EVPA curve. Pointwise, local derivative of the adjusted EVPA (panel g)). The grey area highlights the period of γray flaring activity coinciding with a rotation of the optical polarization angle. 
3.1. The nπambiguity of the polarization vector
The EVPA, χ, is defined within a 180°interval^{3}. Thus, the identification of EVPA swings is ambiguous because of the πmodulus: χ = χ ± nπ, n ∈N. This nπambiguity causes problems in analysing EVPA variability, since the difference Δχ = χ_{2} − χ_{1} ± nπ is ambiguous and even the direction of the rotation is ambiguous (e.g. Marscher et al. 2008; Larionov et al. 2008; Abdo et al. 2010). To visualize rotations larger than π, data points are usually shifted by nπ minimizing the difference between adjacent data points (1)where int(·) denotes rounding to the nearest integer. This procedure is based on the assumption of minimal variation between adjacent data points and relies on adequate sampling.
It has been suggested to consider the EVPA uncertainties in this procedure (e.g. Sasada et al. 2011; Sorcia et al. 2013). In that approach an EVPA data point χ_{i} is shifted according to Eq. (1) only if it shows a significant offset from the previous data point χ_{i − 1}, i.e. , where e_{χi}, e_{χi − 1} are the corresponding uncertainties. This method results in EVPA curves that may depend on the choice of the initial two quadrants it is measured in. We consider, for example, two data points χ_{1} = 40° and χ_{2} = 140° in the initial interval , both with uncertainty e_{χ} = 10°. Then, Δχ_{red} = 85.9° and χ_{2} is not shifted, yielding a counterclockwise rotation of χ_{2} − χ_{1} = 100°. Considering the same data in the initial interval , i.e. χ_{1} = 40° and χ_{2} = − 40°, gives Δχ_{red} = 65.9°. Again, χ_{2} is not shifted but this time yields a clockwise rotation of − 80°. This method produces results that depend on the choice of the initial interval. Instead, using Eq. (1) yields a consistent result independently of the choice of the initial interval.
In Appendix A we introduce a generalization of Eq. (1), which considers more than a single preceding data point, and the consistency level as a quality check of the adjusted EVPA curve. Under the assumption that the adjusted EVPA curve accurately represents the intrinsic variation of the EVPA curve, the consistency level allows one to estimate the probability that the curve is correctly reconstructed. A higher consistency level indicates better sampling and a more reliable adjusted EVPA curve.
3.2. A quantitative measure of smoothness
To gauge the smoothness of an EVPA curve, we define a variation estimator, a quantitative measure of the (erratic) variability of a curve. The larger the variation estimator, the less smooth the curve is. First, we define the pointwise, local derivative of the EVPA curve in units of degrees per unit time: (2)The mean of the derivative indicates a secular trend of the data. The deviation of the local derivative from the secular trend is calculated at each data point as: (3)A local derivative of the order of the secular trend indicates a smooth variation with s_{i} ~ 0. Deviation of the derivative from the trend, ± s_{i}> 0, indicates erratic variation, i.e. the curve is more jagged. We define the variation estimator as the mean absolute s_{i}: (4)We use s as an estimator for the smoothness of a curve with respect to a potential linear trend. An EVPA curve with s_{1} is considered smoother than a second curve with s_{2} if s_{2}>s_{1}.
The standard deviation of the pointwise derivative could be used as an alternative approach to quantify the EVPA variability. Averaging over increases the weight of larger deviations from the secular trend compared to averaging over s_{i}, which makes the standard deviation of s_{i} more susceptible to outliers than the variation estimator defined in Eq. (4). The measured variation estimator is biased by measurement errors and curvature of the smooth EVPA variation. Measurement errors do not average out because of the absolute mean in Eq. (4). Thus errors increase the variation estimator. Only first order trends – estimated through ⟨(Δχ/ Δt)⟩ – are considered to contribute to real EVPA rotations in this analysis. This is the simplest assumption about the intrinsic, smooth variability. Higher order variation needs a priori knowledge and modelling of the variation, which is not given. This higher order variation (curvature) increases the variation estimator, i.e. for nonlinear deterministic variation the smoothness of the curve can be underestimated. Both biases are discussed in more detail in Appendix B.
Before analysing the smoothness of the EVPA curve using the variation estimator s, it is crucial to average over data points that occur in short time intervals, in which the measurement errors are expected to be larger than the real variation. EVPA offsets of the order of the errors at small time intervals will lead to large local derivatives and artificially increase the variation estimator s.
3.3. Optical polarization variation in 3C 279
Figure 1 shows the combined Rband light curve (panel b), the observed optical polarization fraction (panel c) and EVPA curves (panel d) from individual telescopes, the combined, debiased, and averaged polarization fraction (panel e), the combined, nπadjusted, and averaged EVPA (panel f), and the pointwise, local derivative of the adjusted EVPA (panel g) of 3C 279 between JD 2 454 790 and JD 2 456 120. The notable yearly gaps around August–November are caused by the Sun proximity and split the data into four observation periods (I, II, III, IV).
During the first observing period (period I; Oct. 2008–Nov. 2009) the Rband light curve exhibits diminishing flaring activity. During period II (Nov. 2009–Aug. 2010) the optical light curve only has a mean flux density of 0.5 mJy with a standard deviation of σ(f_{ν}) = 0.2 mJy. Periods III (Nov. 2010–Aug. 2011) and VI (Nov. 2011–Aug. 2012) show flaring activity with the flux density ranging from 1.1 mJy to 6.6 mJy with a mean of 2.7 mJy. These outbursts will be discussed in more detail in another paper (Kiehlmann et al., in prep.).
The optical polarization fraction and the EVPA exhibit strong variability. The polarization fraction ranges from ~0 to 0.37 with a mean of ~0.12 and a standard deviation of 0.08. The EVPA shows swings in both directions with amplitudes of up to 500°. The EVPA variation looks smoother during the periods I, III and IV than during period II, as it is already evident in the 180°interval, while the variation during period II does not show any distinct structure. Furthermore, the lower scatter in the plot of the EVPA derivative (panel g) during the periods I, III and IV indicates a smoother variation.
In Fig. 1, we also mark times at which the EVPA curve shows an abrupt change in its behaviour. For each period of consistent EVPA change (periods IIa to IVc, interrupted by two observation breaks), we measure the EVPA rotation amplitude, given by the difference between the minimum and maximum EVPA, Δχ = maxχ(t_{i}) − minχ(t_{i}), the corresponding duration of the largest rotation, i.e. the time passing between the minimum and maximum of the EVPA curve, Δt, the mean and standard deviation of the polarization fraction, the variation estimator, s, and the consistency level of nπadjustments, N_{cons} (see Appendix A). Errors on the polarization fraction mean and standard deviation are estimated using a bootstrap method with 10 000 iterations. Furthermore, we estimate the error bias of the EVPA variation estimator for each observing period through simulations and calculate the debiased variation estimator for each period using Eq. (B.1). The results are listed in Table 2.
The mean polarization fraction is low during period IIa as compared to the flaring periods, during which it can be up to ~3 times higher. The standard deviation of the polarization fraction ranges from 0.028 to 0.084. There is no clear connection between the mean polarization fraction and the standard deviation.
During period I the EVPA shows little variation within an interval of 93°. The time of the apparent 208° rotation in 2009, reported by Abdo et al. (2010), is included in period I. We discuss the absence of this apparent rotation in our data in more detail in Sect. 3.4. Period II exhibits erratic EVPA variability with an overall clockwise trend up to an amplitude of at least 494°. The low consistency level N_{cons} = 1 indicates that period IIa is not sampled sufficiently to reliably reconstruct the intrinsic EVPA variation. The estimation of the probability of correct reconstruction ranges from 0% to <76% (cf. Appendix A).
There are six periods of smoothlooking EVPA rotations, periods IIIa–c, and IVa–c. The first smooth EVPA rotation during the period IIIa is an increase of 86° with a slow rate of + 0.9°/ d, which is followed by a sharp decrease of the EVPA by 110° with a rate of − 16.0°/ d. This distinct change in the EVPA at JD 2 455 650 takes place at the beginning of an optical flaring period. During the optical flaring activity in the period IIIc, the EVPA quite smoothly increases by 352° over 98 days, corresponding to a mean rate of + 3.6°/ d. In the period IVa, the EVPA rotates 109° with a rate of + 1.6°/ d and this rotation continues in the period IVb for another 131° at a faster rate of + 4.7°/ d. Thereafter the sense of rotation changes again and the EVPA decreases by 140° at a rate of + 3.6°/ d. Periods IIIb and IIIc have low consistency levels, each owing to a sampling gap. Figure 1 shows the unaveraged data points in these gaps as open symbols. These data increase the consistency level, indicating a higher reliability of the adjusted EVPA curve. The EVPA swings during period IVa and IVb are well sampled N_{cons}> 50, N_{cons} = 23, corresponding to correct reconstruction probabilities of > 99.9% and 63–99%.
The EVPA variation estimator reaches its highest values during period IIIb. We note that this result is strongly affected by the sparse sampling of the rapid rotation with only nine data points. The secular trend of the EVPA is not reliably estimated and a single data point diverging from this trend by ~80°/ d drastically increases the variation estimator which would otherwise be of the order of s_{debiased} ~ 5°/ d.
During period IIa, the EVPA variation estimator is significantly larger than in the subsequent periods, confirming the observation that the EVPA variation is more erratic during period IIa and smoother afterwards. The erratic variations during period IIa can be either intrinsic or, if the (heteroscedastic) measurement uncertainties are underestimated, owing to measurement noise. The latter explanation would require the measurement uncertainties during period IIa to be underestimated by a factor of 3.7, which is not plausible. Therefore, we conclude that the erratic EVPA behaviour during period IIa is sourceintrinsic. The debiased variation estimator is ~5°/ d during periods IIb–IIIc, during period IVa it is of the order of the error bias, and it increases towards the end of the data (period IVc) , when the flaring activity is declining.
3.4. The apparent 208° rotation in 2009
Abdo et al. (2010) report an apparent 208° clockwise rotation of the EVPA in 3C 279 from JD 2 454 880 to JD 2 454 900 based on KANATA data shown in Fig. 2. The larger than 180° rotation is a result of the usual scheme of minimizing differences between adjacent data points (cf. Eq. (1)). The rotation is sparsely sampled, particularly, there is one gap close to 90° between JD 2 454 888 and JD 2 454 892. Within the uncertainty of σ_{Δχ} = 12.8° it is unclear whether this difference corresponds to a clockwise or counterclockwise rotation.
Fig. 2 208° EVPA rotation presented in Abdo et al. (2010) based on data from the KANATA telescope. Open squares show the EVPAs within the initial 0–180°interval, black circles the adjusted EVPA curve. The numbers show the difference between adjacent data points and the corresponding uncertainties in degrees. 
Figure 3 shows additional data from four more telescopes. Additional data from Calar Alto, Perkins, St. Petersburg and Steward Obs. indicate changes of the rotation direction and variability within less than ~90°, instead of a continuous, large, clockwise rotation. The KANATA data point at JD 2 454 888 appears to be a potential outlier. Omitting this data point, the EVPA only varies within ~60°.
Fig. 3 EVPA data as shown in Fig. 2 with additional data from Calar Alto, Perkins, St. Petersburg, and Steward Obs. 
The reported continuous, clockwise rotation of 208° cannot be maintained. This example clearly demonstrates that densely sampled data is necessary to treat the nπambiguity.
3.5. Optical polarization variability and γray activity
For a comparison with the optical light curves, panel a in Fig. 1 shows the FermiLAT γray light curve at 0.1–300 GeV binned into 3 day intervals. This light curve was originally published in Hayashida et al. (2015). 3C 279 exhibits strong γray flaring activity (>10^{6} ph s^{1} cm^{2}) at various occasions. During the first period of high γray activity around JD 2 454 800 the optical EVPA is roughly constant. During the prominent flare around JD 2 454 890 reported by Abdo et al. (2010), we have shown that the EVPA varies within less than 90°, though the polarization fraction drops significantly. Two high activity states coincide with observation gaps in the optical data. The high γray activity state at ~JD 2 455 666–2 455 741, showing multiple peaks, coincides exactly with the 352° rotation of the optical polarization angle during period IIIc and flaring activity at optical bands. Although the EVPA rotation and optical flaring activity continues during periods IVa to IVc, γray activity is low (< 4 × 10^{7} ph s^{1} cm^{2}). In contrast to the 352° rotation of the optical EVPA coinciding with strong γray activity, the large rotation during period IIa only coincides with mild γray activity (≲10^{6} ph s^{1} cm^{2}). We do not find an obvious, consistent correspondence between the optical polarization variability and γray flaring activity. However, the fact that the 352° EVPA rotation during the period IIIc exactly brackets a longterm γray active period, strongly suggests that this particular polarization angle rotation event is connected to the γray emission.
4. Random walk models
The prevailing physical scenario for launching relativistic jets (e.g. Komissarov et al. 2007; Tchekhovskoy et al. 2011) involves a strong, outwardly propagating, helical magnetic field. The flow of plasma in the jet accelerates and becomes increasingly collimated with distance from the black hole, with the magnetic energy density decreasing in favour of a higher kinetic energy density. It is likely that currentdriven instabilities disrupt the helical ordering near end of the jet’s acceleration zone (e.g. Nalewajko & Begelman 2012). This, plus transverse velocity gradients in the flow (e.g. Vlahakis & Königl 2004) probably result in the flow becoming turbulent. Electrons can be energized by secondorder Fermi acceleration (Stawarz & Petrosian 2008) and magnetic reconnections in the turbulent plasma (Sironi et al. 2015). We model this turbulence in terms of cells, each of which contains a uniform magnetic field, although the actual geometry is certainly more complex than this. We choose the field direction randomly from one cell to the next. Although more sophisticated schemes have been used involving nested cells of different sizes (Jones 1988; Marscher 2015), here we limit our analysis to the basic case of independent cells.
Truncated power law distribution parameters modelling the random time steps (Cols. 2–4).
We investigate whether the observed polarization variation characteristics (⟨m_{l}⟩, σ(m_{l}), Δχ, and s) can be produced by stochastic processes. We do so by performing random walks in Stokes parameters I, Q, and U, and comparing the properties of the obtained polarization curves to the observed ones. The model consists of N_{cells} cells. The properties of a subset of cells are changed at each time step. The sampling of the total simulation time T is randomized with time steps Δt following a truncated power law distribution P(Δt) ∝ Δt^{α} with α< − 1 within limits [Δt_{min},Δt_{max}]. The parameters best describing the distribution of time steps in each observation period are shown in Table 3. The distribution limits are directly adopted from the data. We derive α_{t} in the following manner: (i) from the observation time steps we calculate the maximumlikelihood estimator (MLE) of the power law index ; (ii) we simulate time steps testing different values of α_{t} and calculate the MLE power law index of the simulated time steps ; (iii) we use the value of α_{t} for which . The variation rate n_{var} sets the number of cells that change their properties per unit time step (one day). The total number of cells changed between time t_{i} and t_{i − 1} is given by , where int(·) denotes rounding to an integer. If N_{var}(t_{i},t_{i − 1}) >N_{cells} all cells are changed during that time step. In the following we define three simple random walk processes, that differ in the properties of the cells and in the selection of cells be to changed.
4.1. Simple Q,U random walk process
We create N_{cells} initial cells with uniform intensity I_{i}. Each cell has a uniform magnetic field oriented randomly. The EVPA is oriented accordingly. The cell size, thus, corresponds to the largest scale of uniform magnetic field. Uniformity implies that each cell is maximally polarized. The polarization fraction m_{l,max} of synchrotron radiation in a uniform magnetic field is (Longair 2011, p. 217) (5)where p is the power law index of the electron energy distribution. The maximum polarization is m_{l,max} ≈ 0.72 for p = 2.5 (Longair 2011, p. 217). The initial random variables and Û_{i} are drawn from a Gaussian distribution for each cell i = 1...N_{cells}. The Stokes parameters Q_{i}, U_{i} for each cell are obtained through the following normalization which yields the same intensity and the maximum polarization for each cell: The polarization fraction and angle are independent of the total intensity I and we set I = 1. For the simple random walk process the N_{var}(t_{i},t_{i − 1}) cells that change are randomly selected every time step. This implies that a single cell may have changed several times between t_{i − 1} and t_{i} and that the total number of different cells that change is ≤N_{var}(t_{i},t_{i − 1}). The properties of the selected cells are newly set following Eqs. (6) to (10). In the following we abbreviate this simple random walk process siQU.
Jones et al. (1985) note that choosing and Û_{i} from a uniform distribution in the interval [ − 1,1 ] produces almost identical results. We point out that this prior selection leads to preferred directions of the EVPA every 45° and not a uniform EVPA distribution.
4.2. Ordered Q, U random walk process with constant I
This process implements the basic concept of a disturbance passing through a turbulent medium, locally increasing the emissivity and successively highlighting part of the medium. In this case, the model cells are numbered. The first cell represents the disturbance and all the following cells are a trailing afterglow. The initial cells are created exactly as in the simple Q, U random walk process, but the subsequent cell selection differs. Assuming that N(t_{i}) cells change between time t_{i − 1} and t_{i}, the properties of all cells are shifted by N(t_{i}) cells. The first N(t_{i}) cells, representing the region right behind the disturbance, are newly set with random properties. The properties of the last N(t_{i}) trailing cells vanish. We call this process orQUc in the following. The orQUc process is less randomized than the siQU process. Each cell the disturbance passes through has its properties maintained for a given amount of time. For the siQU process, on the other hand, the properties of a single cell are maintained a random amount of time depending on when it is randomly selected to change.
4.3. Ordered Q, U random walk process with decreasing I
For this random walk process we assume the intensity in a cell is increased by a disturbance and gradually degrades as the disturbance moves on. The highest intensity of the emitting region is located at the disturbance and decreases as function of distance from the disturbance. We scale the cell intensities I_{i} linearly between the highest intensity at the disturbance (cell one) and zero intensity at the N_{cells} + 1th trailing cell: (11)The total intensity is again set to I = 1. The cell properties are set by Eqs. (6), (7) and (9), (10), (11). Cell properties are changed following the scheme of the orQUc process (Sect. 4.2). We denote this process orQUd.
The three random walk processes are simplistic implementations of a finite emission region in a turbulent flow that is either randomly changing its magnetic field structure (siQU) or is energized by a disturbance moving through the plasma (orQUc, orQUd). These implementations face several limitations. The siQU process does not include any structure information, the individual cells do not have a specific position. The ordered random walk processes assume the cells are all in line, whereas multiple cells could be next to one another. This would affect the intensity scaling of the orQUd cells. All processes neglect that the magnetic field orientation cannot be entirely disconnected between different cells. Furthermore, a disturbance such as a shock would locally order the magnetic field and leave the field orientation in trailing cells correlated to some extent. Time delays which would arise from the different positions of the cells are not considered and the modelling is performed in the observers frame under the assumption that neither the speed nor the direction of the motion of the emitting region changes during the polarization event. Relativistic aberration, which would make randomly oriented cells appear to be more aligned than they intrinsically are, is also not considered here.
4.4. Integrated polarization
The integrated Stokes parameters and Q, U, accordingly, determine the integrated intensity I, the linear polarization fraction m_{l} and the EVPA χ at each simulation time step: Since the plasma is optically thin at optical frequencies and the Faraday rotation scales with the squared wavelength of the radiation, we can safely ignore the effects of synchrotron selfabsorption and Faraday rotation on the polarized signal at optical frequencies.
4.5. Simulated EVPA errors
Fig. 4 Empirical cumulative distribution function (ECDF) of the combined Stokes q and u uncertainties in the entire data (red curve), period IIa (orange curve), and period IIIc (green curve). 
We have tested two schemes of implementing the simulation of observational uncertainties. In the first scheme, we apply uncertainties to the integrated Stokes parameters q = Q/I, u = Q/I. For each simulated data point the uncertainties σ_{q,i} are randomly drawn from the empirical cumulative distribution function (ECDF) of the combined q and uuncertainties measured in the corresponding period (Fig. 4). In the simulation we set σ_{u,i} = σ_{q,i}. Noise q_{err,i}, u_{err,i} is then drawn from a Gaussian distribution and it is added to each random walk data point. The data points and uncertainties are then transformed into m_{l}, χ and the corresponding uncertainties following King et al. (2014), Eqs. (5), (6). In the observed data the uncertainties in q and u is generally not the same, but they scatter around this linear correlation. As a consequence this method does not exactly reproduce the measured distribution of EVPA uncertainties.
In the second scheme we apply uncertainties directly to the EVPA. The measured EVPA uncertainties can be described by a lognormal distribution. The distribution parameters μ_{unc} and σ_{unc} are estimated for each period individually. To simulate EVPA “measurement” errors we draw a randomnumber σ_{χi} from a lognormal distribution representing the measurement uncertainty of each simulated data point χ_{rw,i} = χ_{rw}(t_{i}). To each data point we then add noise χ_{err,i} drawn from a Gaussian distribution with zero mean and standard deviation σ_{χi}. Finally, each simulated EVPA data point is given by χ_{sim,i} = χ_{rw,i} + χ_{err,i}. This scheme produces EVPA uncertainties more close to the observed distribution than the first scheme. Limitation of this approach are that uncertainties in the polarization fraction are neglected and that the simulated EVPA uncertainties are independent of the simulated polarization fraction.
We ran all simulations presented in the following sections with both schemes to test if the results depend on the specific implementation of the uncertainty simulation. All results presented are based on the first scheme, applying uncertainties to the Stokes parameters. Whereas the absolute numbers, distributions and figures shown in the following slightly differ using the second scheme, the results and conclusions do not depend on the choice of the uncertainty implementation.
4.6. Random walk simulations
We vary the primary input parameters in the ranges N_{cells} ∈ [ 2,600 ], n_{var} ∈ [ 0.1 d^{1},100 d^{1} ]. We run the simulations in two modes. In the first one, the cell number and the cell variation rate are randomly drawn from (discrete/continuous) uniform distributions: and . We run a total of 1 000 000 simulations for each random walk process. Additionally we run simulations with certain chosen input parameters. For each parameter combination and random walk process we run up to 10 000 simulations. The choice of input parameters is described in the following sections.
The fixed secondary input parameters for the random walk simulations are taken from Table 3. The total simulation time T equals the duration of the corresponding period. For each simulation we measure the mean polarization fraction ⟨m_{l}⟩ and its standard deviation σ(m_{l}), the EVPA amplitude Δχ, and the EVPA variation estimator s.
5. Random walk simulation results
In this section we first describe general results of the random walk simulations: expectation values of the mean and standard deviation of the polarization fraction, a general comparison of the three random walk processes, and dependencies between the model parameters and various measured parameters. Then, we test the periods of the two largest rotations in the observed data (period IIa and IIIc) against the models. Finally, we test the rotation amplitude distribution of the entire observed EVPA curve against the random walk models.
5.1. Expectation values of the fractional polarization
For each random walk model we run simulations for selected pairs of cell numbers and variation rates to test the expectation values of the mean and standard deviation of the polarization fraction. The fixed parameters are taken from Table 3, period “total”. The simulation time is T = 260 d, corresponding to the average observation period between sun gaps. We run 10 000 simulations per a pair of parameters and a random walk model.
Fig. 5 Expectation values of the polarization fraction mean (upper panel) and standard deviation (lower panel) for different numbers of cells depending on the cell variation rate. Each data point is based on 10 000 simulations. Solid curves and circles correspond to the siQU model, dashed curves and squares to the orQUc model, and dotted curves and triangles to the orQUd model. 
The mean and standard deviation of the simulated polarization fraction follow lognormal distributions for each input parameter combination. We estimate the distribution parameters via maximum likelihood and from those we calculate the expectation value and variance for ⟨m_{l}⟩ and σ(m_{l}), for each combination of N_{cells} and n_{var}. The results are shown in Fig. 5. The expectation value of the mean of the polarization fraction is independent of n_{var}.
Fig. 6 Distributions of a) the mean and b) the standard deviation of the polarization fraction; c) the EVPA amplitude; and d) the EVPA variation estimator for each random walk process with a range of model parameters. 
The siQU and the orQUc process produce the same amount of polarization. Simulated curves from the orQUd process are on average more polarized by a factor of 1.13–1.15. As anticipated, for all processes the expected mean of the polarization fraction depends on the cell number by . The expectation value of the standard deviation of the polarization fraction increases with increasing n_{var} and decreases with increasing N_{cells}. For high variation rates it saturates at .
Typically, similar relations connecting the cell number and cell variation rate with the expectation values of the polarization fraction mean and standard deviation are used to fix the model parameters based on the observed polarization fraction, when testing stochastic models against the data (e.g. D’Arcangelo et al. 2007). We do not recommend to fix these parameters, as each pair of model parameters results in a broad distribution of ⟨m_{l}⟩ and σ(m_{l}) with standard deviations up to and . The distribution is also righttailed. Thus, the expectation value does not match with the distribution mode. Furthermore, the parameter combination most likely producing the observed polarization fraction may not be the best choice for producing the observed EVPA variability. Instead, we recommend to test a large parameter space and to identify the parameter region that is most likely to produce the observed polarization fraction and the EVPA characteristics.
5.2. Model comparison and parameter dependencies
Figure 6 shows the distribution of the mean and standard deviation of the polarization fraction and the EVPA amplitude and variation estimator for each tested random walk process with primary model parameters in the ranges N_{cells} ∈ [ 2,600 ], n_{var} ∈ [ 0.1 d^{1},100 d^{1} ], secondary parameters as listed in Table 3, period “total”, and a total time of T = 260 d. Each distribution is based on 1 000 000 simulations. The distributions are affected, especially at the lowest values, by the limits of the tested parameter space, and the EVPA amplitude distribution additionally depends on the simulation total time.
The comparison of the three random walk processes based in Figs. 5 and 6 shows, that the siQU process is less variable than the orQU process, resulting on average in smaller EVPA amplitudes but smoother EVPA curves. The reason is, that in the siQU process single cells may change multiple times within one time step, thus the number of cells that stay the same is larger than in the ordered random walk processes. The orQUd process is generally more polarized and more variable, resulting in larger EVPA amplitudes and more erratic EVPA behaviour. Given the intensity scaling of the cells in the orQUd process, this process is comparable to the orQUc process with de facto fewer cells.
In general, increasing the cell number decreases the mean and standard deviation of the polarization fraction, the average EVPA rotation amplitude and the variation estimator, while increasing the cell variation rate increases the distribution means of all measured parameters. Consequentially, there is a correlation between the EVPA rotation amplitude and the variation estimator in the simulations. Larger EVPA rotations are indeed on average less smooth, as one would intuitively expect.
The results furthermore show that in principle very smooth EVPA curves s ~ 1°/ d and very large rotations Δχ ≫ 360° can be produced by a random walk process. However these properties are to some extent mutually exclusive. Therefore, it is mandatory to simultaneously consider both the rotation amplitude and the quantified smoothness in the analysis of the random walk simulations.
5.3. Probabilities of observed polarization variability
The main aim of these simulations is to determine the probability of producing the observed polarization variability characteristics by a random walk process. We use ⟨m_{l}⟩ and σ(m_{l}) to characterize the polarization fraction variation, while Δχ and s are used to characterize the EVPA variation. In the following, we define four test conditions.
First, we define two polarization fraction conditions (PF1, PF2): We search for simulations with ⟨m_{l}⟩ and σ(m_{l}) in ranges [ x − Δx,x + Δx ], where x is the corresponding observed value and Δx is equal to the scatter of the observed ⟨m_{l}⟩ and σ(m_{l}) (cf. Table 2). We note, that the interval range is an arbitrary choice and the final probabilities depend on the accepted range of parameters. Our choice is to set the accepted parameter range based on the data instead of preselecting a fixed value.
Second, we define two polarization angle conditions (PA1, PA2): We search for EVPA rotations with amplitudes as large or larger than observed Δχ ≥ Δχ_{obs} which are at the same time as smooth or smoother than the observed data, s ≤ s_{obs}. Finally, we search for simulations fulfilling all four conditions (PF1+PF2+PA1+PA2).
We run simulations for various input parameter combinations N_{cells,i}, n_{var,j}. For each parameter combination ij we run a set of N_{sim} = 10 000 simulations. From each simulation set ij we find the number of simulations in agreement with our conditions, N_{ij,cond}. The probability of producing a polarization curve in agreement with this condition from a random walk process with input parameters (N_{cells,i}, n_{var,j}) is P_{i,j}(cond) = N_{i,j,cond}/N_{sim}. The highest probability is P_{max}(cond) = maxi,jP_{i,j}(cond) with the optimal input parameters (, ). The accuracy of this probability estimation depends on the sampling of the input parameters.
Large rotations of the polarization angle are of particular interest. Therefore we test the random walk simulations against the two largest rotations in our data, observed during the periods IIa and IIIc. In particular, we test period IIIc because of its contemporaneity with strong γray flaring activity. We compare these two periods as they show very different behaviour in the smoothness of the EVPA rotation, the polarization fraction and the total flux. We discuss the other periods qualitatively after the detailed study of periods IIa and IIIc.
5.3.1. Period IIa
We sample the input parameter space in the ranges N_{cells} = [ 80,90,...,260 ] and n_{var} = [ 4,8,...,60 ] d^{1}. For each parameter pair we run 10 000 simulations. The fixed parameters are taken from Table 3, period II. The simulation time is T = 154 d. The simulation selection conditions are: Figure 7 shows the probability distributions over the input parameter space for the polarization fraction conditions (PF1/2, first row), the EVPA conditions (PA1/2, second row), and for all four conditions (PF1/2+PA1/2, third row) for each of the three random walk processes. The third row of Fig. 7 indicates that the tested parameter space is sufficiently large to capture the region of highest probability given all four conditions. The highest probabilities and the corresponding optimal input parameters are listed in Table 4.
Fig. 7 Probability distributions over the parameter space N_{cells}, n_{var} for three different conditions (rows 1–3) based on three random walk processes (Cols. 1–3) for period IIa. 
Highest probabilities of producing polarization curves fulfilling the conditions on the polarization fraction (PF1/2), on the EVPA variation (PA1/2) and fulfilling all conditions from three different random walk processes with the corresponding, optimal model parameters for period IIa.
The highest probability of producing the observed polarization fraction characteristics (conditions PF1/2) is ~18%. It is relatively low because the observed standard deviation of the polarization fraction of 0.03 is larger than the maximum value of σ(m_{l}) expected from a random walk producing a mean polarization fraction of 0.05, which is 0.025. EVPA rotations of equal or higher amplitude and comparably smooth or smoother than observed (conditions PA1/2) are occurring with probabilities up to 10–18% in the tested parameter space. The probability of producing the observed variability in both the polarization fraction and angle (conditions PF1/2+PA1/2) is relatively low compared to the individual conditions. The reason is that the comparably high standard deviation of the polarization fraction requires higher cell variation rates than the observed EVPA variability. The discussed processes are capable, however, of producing the observed polarization characteristics with a probability up to ~3%.
Figure 8 shows the distribution of the EVPA rotation amplitudes (upper panel) and the EVPA variability estimator (lower panel) for those simulations fulfilling the polarization fraction conditions (PF1/2). The histograms are based on 1 000 000 simulations per random walk process with random primary model parameters. The observed values during period IIa are marked by vertical lines. The distribution of rotation amplitudes shows that rotations as large as observed commonly occur in random walk processes, given the observed polarization fraction. The observed EVPA variability estimator is close to the distribution mode. The erratic behaviour of the EVPA during the period IIa is characteristic of a stochastic process.
5.3.2. Period IIIc
To test the random walk processes against the period IIIc data we sample the cell number in the range N_{cells} = [ 5,10,...,50 ] and the cell variation rate in the range n_{var} = [ 0.1,0.2,...,1.0,2.0,...,20 ] d^{1}. The secondary parameters are listed in Table 3, period III and the simulation time is 98 days. We run 10 000 simulations per pair of parameters and use the following conditions to analyse the simulations: Figure 9 shows the probability distribution over the tested parameter space for all three random walk processes (Cols. 1–3) and different selection conditions. The probability of producing the observed polarization fraction characteristics (conditions PF1/2) during period IIIc is shown in the first row and ranges up to ~7% (cf. Table 5). The observed standard deviation of the polarization fraction of 0.057 is relatively low compared to the value ~0.065 expected from a random walk with a mean polarization of 0.131. Thus, the probability of producing the observed mean and standard deviation is relatively low.
The second row of Fig. 9 shows the probability of producing EVPA rotations with amplitudes at least as large as the observed rotation (condition PA1). Those rotations are most likely produced by a few cells and high cell variation rates. On the other hand, the observed low EVPA variation estimator, which indicates a smooth EVPA curve, only can be produced by random walk models with low cell variation rates and preferentially many cells as shown in the third row of Fig. 9 (condition PA2). Consequentially, the probability of producing an EVPA rotation at least as large and at least as smooth (conditions PA1/2) is low (cf. Table 5).
Fig. 8 Distribution of the EVPA rotation amplitude (upper panel) and distribution of the EVPA variation estimator (lower panel) for all simulations with a mean polarization fraction consistent with the observed value during period IIa. The number of selected simulations is indicated in the legend. The corresponding observed values are indicated by black dotted lines. 
Fig. 9 Probability distributions over the parameter space N_{cells}, n_{var} for three different conditions (rows 1–3) based on three random walk processes (Cols. 1–3) for the period IIIc. 
Highest probabilities of producing polarization curves fulfilling the conditions on the polarization fraction (PF1/2), on the EVPA variation (PA1/2) and fulfilling all conditions from three different random walk processes with the corresponding, optimal model parameters for the period IIIc.
In Fig. 10 we show the distribution of the rotation amplitude (upper panel) and the EVPA variation estimator (lower panel) based on the simulations which show the observed polarization fraction characteristics (conditions PF1/2). On one hand, the observed rotation amplitude is close to the distribution mode. On the other hand, the EVPA variability coming from a stochastic process would be expected to be far more erratic. Comparably smooth or smoother variations than observed do occur, but only rarely and, most importantly, those simulations do not produce large rotation amplitudes. Consequentially, we do not find a single occurrence of variability comparable to the observed period IIIc in both the polarization fraction and angle (conditions PF1/2+PA1/2) at any tested point in the parameter grid. Therefore, we conclude that it is very unlikely that a stochastic process could produce the observed EVPA rotation during the period IIIc.
Fig. 10 Distribution of the EVPA rotation amplitude (upper panel) and distribution of the EVPA variation estimator (lower panel) for all simulations with a mean polarization fraction consistent with the observed value during the period IIIc. The number of selected simulations is indicated in the legend. The corresponding observed values are indicated by black dotted lines. 
5.3.3. Other periods
We have focused on the two most prominent rotations during the periods IIa and IIIc when testing the stochastic models. In this section we discuss the other periods qualitatively. During period I the EVPA shows little variability (≲90°). Though not explicitly tested for the polarization fraction of period I, the tests of periods IIa and IIIc suggest that such small EVPA amplitudes are unlikely to occur from the tested random walk processes (cf. Figs. 8 and 10), given that the variability in the polarization fraction is the highest during period I. Consistent with the likely nonstochastic period IIIc, the polarization fraction and the total flux is mostly high compared to period IIa. The average polarization fraction drops and the EVPA variability becomes more erratic with the fading flare. If both periods are indeed of deterministic origin, the smooth and weak EVPA variability during period I on the one hand, and the smooth and strong variability during period IIIc on the other hand, imply a complex behaviour of 3C 279.
Period IIb shows a lowbrightness state in Rband. Nevertheless, the polarization fraction is exceptionally high in the mean and standard deviation and the EVPA variation estimator is similarly low as during the flaring states. Whether this drastic increase in the polarization fraction marks the beginning of the following flaring state, without showing up yet in the total flux density, cannot be properly tested because of the following observation gap.
The polarization fraction and the EVPA variation estimator during period IIIa are comparable to period IIIc despite the significantly lower rotation rate. Figure 10 suggests that we would expect a larger rotation amplitude from a random walk model, arguing against a stochastic process. Claiming a nonstochastic origin for periods IIIa and IIIc implies that also period IIIb should be produced by another process rather than a random walk. If all periods IIIac originate in the same event, the sharp clockwise rotation opposite to the counterclockwise direction of the enclosing rotations poses a challenge to any nonstochastic model. Although not explicitly tested we suspect period IVa not to be consistent with the random walk processes, given the continuous, smooth rotation.
5.4. EVPA rotation amplitude distribution
In Appendix C we introduce a method to automatically identify continuous rotations in an EVPA curve based on the strict criterion that a continuous rotation begins and ends with a significant change of the rotation direction at a certain significance level ς. We point out that the rotations discussed thus far do not strictly follow this criterion, but were instead defined as the maximum change of the EVPA within a certain time period that shows a general trend in the EVPA. However, in this section we follow the strict definition. At a ς = 3 significance level 109 rotations are identified in the optical EVPA curve. The distribution of the rotation amplitudes is shown in Fig. 11. In Appendix C we discuss characteristics of the rotation amplitude distribution originating from the random walk models.
Fig. 11 Histogram of the rotation amplitudes of the 109 rotations identified in the optical EVPA curve at ς = 3significance and the corresponding cumulative distribution function (ECDF, red curve). 
To test the observed distribution of rotation amplitudes against the random walk models we run long simulations with a total time T = 100 000 d varying the model parameters in the ranges N_{cells,i} ∈ [ 10,20,...,300 ], n_{var,j} ∈ [ 2,4,...,20 ] d^{1} for each random walk process. Secondary model parameters are taken from Table 3, period “total”. For each simulation the rotations are identified at ς = 3 significance level and the distribution of the rotation amplitudes is tested against the rotation amplitude distribution of the entire EVPA curve using a twosample KolmogorovSmirnov (KS) test. We note, that this analysis is not feasible for the individual periods as too few identified rotations do not allow us to estimate the distribution of rotation amplitudes reliably.
Fig. 12 pvalues of the KS test over a grid of model parameters, testing the rotation amplitude distribution of the entire EVPA curve against random walk simulations for three different stochastic processes. 
The pvalues of the KS tests are shown in Fig. 12 over the tested parameter space for each random walk process. Solid curves mark the edge of the parameter space where the hypothesis that the observed distribution of rotation amplitudes originates in the tested random walk process is rejected. For cell numbers lower than 50 the hypotheses that the observed EVPA rotations are produced by the siQu, the orQUc, or the orQUd process are rejected at 1σ, 2σ, and 3σ significance, respectively, regardless of the cell variation rate. In Fig. 9 we have shown that cells fewer than this are needed to produce the high mean polarization fraction for instance during period IIIc. Additionally, cell variation rates >5 cells per day are preferred to produce the observed variability in the polarization fraction. At these variation rates the hypothesis that the observed EVPA rotations are produced by one of the tested processes is rejected at least at 2σ significance. Therefore, we conclude that the distribution of the rotation amplitudes identified in the entire EVPA curve cannot be produced by the tested random walk process with a fixed number of cells.
6. Discussion: two processes?
EVPA rotations of opposite direction in 3C 279 have been reported before in the literature (Larionov et al. 2008; Abdo et al. 2010). In the data presented here, we observe several changes of the rotation direction with two large rotations during the periods IIa and IIIc. These two periods show significantly different variability in the optical flux and optical polarization. During the period IIa the flux and the polarization fraction are lower and less variable, whereas the EVPA is more erratic than during the period IIIc. Testing the polarization data of both periods against random walk models makes the difference even more evident. We have shown that the erratic behaviour of the EVPA during the period IIa has the characteristics of a stochastic process. Although the probability of ≲3% is not high, the discussed models are capable of producing the observed polarization variability. Most likely around 170 or 220 cells (depending on the process) are necessary to produce the observed behaviour. On the other hand, the high polarization fraction during the period IIIc requires significantly fewer cells ~30, implying that this polarized emission is produced in a much smaller region. We have shown that the long and smooth rotation of the EVPA during this period is highly unlikely to originate in a stochastic process and we could not reproduce the observed variability in the polarization fraction and angle with any of the tested stochastic processes. A deterministic origin of this EVPA rotation is furthermore supported by the exact contemporaneity with a strong γray flaring period.
We come to the conclusion that we likely observe two different processes responsible for the polarization variability during the lowbrightness state and the flaring state of 3C 279. During the lowbrightness state the polarization is consistent with a stochastic process. We have argued that the observed variability is unlikely to be an observational artefact owing to low signaltonoise data and underestimated uncertainties, but it is sourceintrinsic. This variability is consistent with a stochastic process as implemented by the siQU, orQUc, and orQUd process, modelling a general turbulent flow (siQU), a disturbance passing through a turbulent medium or, vice versa, a turbulent medium passing through a local disturbance (orQUc, orQUd). These models are not distinguishable, except from requiring different model parameters to produce the observed polarization. We suggest that the polarization during the lowbrightness state is affected by the turbulent flow with a relatively large number of emission cells having individual magnetic field orientations, thus producing a relatively low mean polarization fraction. We suggest that this stochastic variability of the polarization is always present, but only dominant during the lowbrightness state.
In the data we find two phases which might be transitions from a deterministic to a stochastic process. In period I the initial high activity state in optical and γrays dissolves. With the fading flares the polarization fraction drops and the EVPA variability becomes far more erratic (cf. Fig. 1). Again, throughout period IVb and towards the end of IVc, when the flux density is decreasing, the EVPA variability estimator rises. These phases could imply that with the flaring emission fading the turbulent background becomes more dominant again.
6.1. Deterministic process during period IIIc
During the flaring state the polarization fraction is on average higher than during the lowbrightness state, indicating that the polarized flux is dominated by a smaller emission region. The polarization variability cannot be explained by the tested stochastic models. Although we only tested very simplistic random walk processes, this result challenges more sophisticated stochastic models such as the turbulent extreme multizone model (TEMZ; Marscher 2014). In contrast to stochastic models, deterministic models are expected to produce smooth EVPA rotations, some of them following distinct patterns. In the following we discuss several models presented in the literature.
Laing (1980) discusses the shock compression of a tangled magnetic field into an apparently ordered one to explain a high polarization fraction despite a generally random field structure. In principle, this process can produce an EVPA swing if the initially tangled magnetic field produces a low net polarization which is oriented differently than compressed, apparently ordered field. But this effect can only produce swings up to 90° at maximum. The more generic case of two superposed, evolving emission features discussed by Holmes et al. (1984) is also limited to swings <90°. The observed ~360° rotation excludes these models.
Two purely geometric models are based on curved trajectories of the emission region. In the first model, Nalewajko (2010) assumes an axially symmetric magnetic field and a global bend of the flow. Relativistic aberration can drastically change the observed EVPA of the emission region going through even a small bend, if also the viewing angle is small. Yet, the EVPA rotations in this model are limited to <180° for a simple bend on a plane (Nalewajko 2010). Abdo et al. (2010) used this model to explain the apparent ~208° EVPA swing in 3C 279 within the 20 days after JD2 454 880. We have shown, though, that this apparent rotation was an artefact of the sparse sampling and that additional data does not support the smooth, long trend initially reported. Aleksić et al. (2014a) observed the period IIIc rotation only partially with an amplitude of ~140° and used the same bending jet model to explain this event. The full ~360° rotation cannot be explained by a single flaring region on a slightly bent trajectory. This rotation would either require a helical jet structure or at least two successive emission regions passing through the bend. Although we observe several subflares in the Rband light curve during the EVPA rotation, these consecutive flaring regions would likely need significant finetuning in the time separation of the individual events to produce a continuous 360° rotation.
The second geometric model describes an emission feature smaller than the cross section of the jet on a helical trajectory traversing through a helical magnetic field (Kikuchi et al. 1988; Marscher et al. 2008). On its path the feature highlights different parts of the magnetic field structure producing a gradually changing EVPA. In contrast to the previous models this one is not restricted by an upper limit of the rotation amplitude, but can naturally produce in principle arbitrarily long rotations. The deviations from a constant rotation rate can originate, e.g. in changes in the flow speed, varying lighttravel time delays along the helical path, or the superposition of the emission feature and a constantly (or randomly) polarized background (e.g. Marscher et al. 2010; Larionov et al. 2013).
Zhang et al. (2014, 2015) suggest that EVPA swings can occur owing to light travel time effects when a relativistically moving plasma pervading a helical magnetic field encounters disturbance such as a shock. In contrast to previously discussed models this one naturally explains contemporaneous EVPA rotations and flaring activity throughout the entire spectrum from microwaves up to γray emission, which could explain the contemporaneity of the EVPA rotation and the γray flaring during period IIIc (Fig. 1). However, in the model by Zhang et al. (2014) a single event can only produce EVPA rotations up 180°. To explain the observed rotation in the period IIIc, at least two successive flaring events would be needed. There are two sharp drops observed in the polarization fraction during the EVPA rotation in the period IIIc, which qualitatively fits the behaviour expected from the Zhang et al. (2014) model. However, detailed modelling will be needed to test whether the nonstochastic ~360° EVPA rotation during the period IIIc can be produced by the nonaxisymmetric helical motion model or by the axisymmetric model of a shock in a helical magnetic field.
The most critical challenge to deterministic models producing large EVPA rotations (≥ 180°) are the multiple inversions of the rotation direction observed throughout the entire polarization data. Different directions in the nonflaring and in the flaring state can be explained with the twoprocesses interpretation, as a stochastic process during the nonflaring state produces rotations in both directions equally likely. But we observe inversions of the rotation direction even within period IIIc. Deterministic models that are capable of producing twodirectional EVPA swings such as the twocomponent model (Holmes et al. 1984) and some of the models presented in Zhang et al. (2014) typically are limited to amplitudes ≤ 90°, producing a rotation followed by a counterrotation back to the original orientation; whereas a single bend in the jet, a helical path of the emission feature in a helical magnetic field and the models in Zhang et al. (2014) that produce 180° rotations are expected to be geometrically bound to a single EVPA rotation direction. The superposition of these models with a constantly or randomly polarized background and the superposition of multiple events overlapping in time could produce more complex patterns in the EVPA variability. For example, two emission regions, each on its own may produce a clockwise rotation of the polarization angle. Following the twocomponent model of Holmes et al. (1984), the superposition of both regions may temporarily produce a counterclockwise rotation of less than 90°, if these regions change their physical properties as the total intensity or spectral index, for instance when a new emission feature progressively outshines the previous one. This superposition could lead to two inversions of the rotation direction as seen in period IIIc. Considering more than one component, of course, drastically increases the model complexity.
6.2. Comparison with recent RoboPol results
Blinov et al. (2015) investigate the potential connection between γray flares and rotations of the optical polarization angle in a statistical way based on a large sample of blazars monitored with the RoboPol instrument. They show that a stochastic process can in principle produce the rotation amplitudes observed, but it is statistically unlikely that all rotations are produced by a stochastic process. This result is consistent with our conclusion that even within the same object two different processes – stochastic and deterministic – may be responsible for different rotations. Furthermore, Blinov et al. (2015) argue that optical EVPA swings and γray activity are not necessarily physically connected, but it is unlikely that none of the observed events are connected. Particularly, the strongest γray flares had time lags to EVPA rotations close to zero. They conclude there are two processes: one producing strong γray flares and contemporaneous rotations of the optical polarization angle; the other producing moderate γray flares physically not connected to the optical polarization activity, which may be coincident happening owing to a stochastic process. These results are consistent with ours. Period IIIc is an example of a nondeterministic EVPA rotation contemporaneous with strong γray flaring activity, whereas the rotation during period IIa is probably owing to a stochastic process, potentially not connected to the moderate γray variability. Complementary to Blinov et al. (2015), we show that different processes responsible for the polarization variability can occur in the same object. Additionally, we show that strong γray flaring activity can occur without strong variability in the optical EVPA (period I) and vice versa (period IVa–IVc).
6.3. Comparison with BL Lac and PKS 1510089
3C 279, PKS 1510089, and BL Lac are the most prominent sources showing rotations of the optical EVPA and contemporaneous, strong γray activity. The rotations in BL Lac in 2005 and in PKS 1510089 in 2009 were explained by a deterministic process, the motion of an emission feature on a spiral path in a helical magnetic field (Marscher et al. 2008, 2010). The rotations in PKS 1510089 in 2012, on the other hand, showing several inversions of the rotation direction, were interpreted as originating from a stochastic process (Aleksić et al. 2014b). These events show some similarities but also differences in the progression of the total optical flux and the polarization fraction during the rotation of the polarization angle.
The rotations in BL Lac in 2005 and in PKS 1510089 in 2009 both start and end with a strong peak in the polarization fraction (Marscher et al. 2008, 2010). The 380° and 250° rotations in PKS 1510089 in 2012 do not exhibit as pronounced peaks, but indicate maxima in the polarization fraction in the beginning and end of the rotation and a minimum during the event (Aleksić et al. 2014b). The rotation in 3C 279 during period IIIc also starts and ends with pronounced peaks in the polarization fraction, though in contrast it shows an increasing trend throughout the event underlying multiple extrema. This increasing trend of the polarization fraction with a sharp drop at the end of the rotation is seen in all smooth rotation events, periods IIIa, IIIc, IVa, IVb. In contrast, the erratic rotation in period IIa does not show a global trend and no pronounced peak of the polarization fraction at the end of the rotation. The beginning of the rotation is not captured by observations.
All discussed rotation events in BL Lac and PKS 1510089 and the period IIIc rotation in 3C 279 coincide with one or multiple peaks in the optical light curve and γrays, but there is no obviously consistent shape of the curves. Periods IIa, IIIa, and IVa–IVc in the data of 3C 279 show that rotations of the optical polarization angle can happen without coinciding, strong variability in γrays.
7. Summary and conclusions
In connection to an intensive multiwavelength campaign, we have accumulated and combined multiple data sets of Rband photometry and optical polarimetry measurements into wellsampled flux density and polarization curves of the γray loud quasar 3C 279 covering a period from Nov. 2008 to July 2012. These data capture 3C 279 in an optical lowbrightness state between Nov. 2009 and Aug. 2010, followed by an optical flaring state, which coincides with an increased γray activity (Aleksić et al. 2014a). We observe strong optical polarization variation and EVPA rotations during both states. We have critically discussed the nπambiguity and different solutions to treat it and we presented a method to estimate the quality of the data given the time sampling and the observed rotation rates. We have shown that the sparsely sampled EVPA variability reported by Abdo et al. (2010) does not show a continuous, large rotation when additional data is considered. Generally, we find EVPA rotations with both a clockwise and a counterclockwise sense of rotation. These multiple changes of the rotation direction observed in the data eliminate all except perhaps the most geometrically complex models, to explain the entire polarization variability.
Instead, we have tested whether or not the observed EVPA rotations can be of stochastic origin using the smoothness of the EVPA curve as a key indicator. To do this, we introduced the EVPA variation estimator as a quantitative measure of the curve smoothness. We simulated three different processes based on random walks in Stokes Q and U and found that all of them are highly unlikely to produce a smooth ~360° EVPA rotation as observed during the period IIIc (the flaring state), especially coinciding with the high polarization fraction mean and variation. We conclude that the tested class of simplistic stochastic processes based on emission cells that have a random and variable magnetic field orientation cannot produce the observed polarization variation during the flaring state of 3C 279. This result challenges more sophisticated stochastic models such as the TEMZ Marscher (2014).
However, in the lowbrightness state (period IIa) the EVPA variation is much more erratic than during the flaring state and has the characteristics of random walk processes. Hence, we find two different states: the lowbrightness state exhibits a comparably low polarization fraction and erratic EVPA variation, possibly consistent with stochastic variation; the flaring state shows a high polarization fraction and very smooth EVPA variation with a continuous counterclockwise rotation sense. We interpret this result as two different processes contributing to the polarization variation in 3C 279. On the one hand, there probably exists an underlying stochastic variation, which is visible in the lowbrightness, less polarized state. Any stochastic model naturally explains the frequent changes of the rotation direction observed during the corresponding period. On the other hand, during the flaring state a small, highly polarized region of the optical jet dominates the total and the polarized flux. The polarization variation in this small region is not produced by the class of stochastic processes we tested. Yet, deterministic models are challenged by multiple changes of the rotation direction observed even during the flaring period. For period IIIc we can certainly exclude the bending jet scenario, which had been suggested by Aleksić et al. (2014a), who only observed a part of this EVPA rotation, which significantly exceeds 180°.
We have tested three simplistic random walk processes against polarization data. Whereas these toy models neglect various physical effects, we have demonstrated ways to compare also more sophisticated models statistically with polarization data.
Acknowledgments
S.K. was supported for this research through a stipend from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Max Planck Institute for Radio Astronomy in cooperation with the Universities of Bonn and Cologne. T.S. was partly supported by the Academy of Finland project 274477. The research at Boston University was partly funded by NASA Fermi GI grant NNX11AQ03G. K.V.S. is partly supported by the Russian Foundation for Basic Research grants 130212103 and 140231789. N.G.B. was supported by the RFBR grant 120201237a. E.B., M.S. and D.H. thank financial support from UNAM DGAPAPAPIIT through grant IN1162113. I.A. acknowledges support by a Ramon y Cajal grant of the Spanish Ministry of Economy and Competitiveness (MINECO). The research at the IAACSIC and the MAPCAT program are supported by the Spanish Ministry of Economy and Competitiveness and the Regional Government of Andalucía (Spain) through grants AYA201014844, AYA201340825P, and P09FQM4784. The Calar Alto Observatory is jointly operated by the MaxPlanckInstitut für Astronomie and the Instituto de Astrofísica de AndalucíaCSIC. Data from the Steward Observatory spectropolarimetric monitoring project were used. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX14AQ58G. St.Petersburg University team acknowledges support from Russian RFBR grant 150200949 and St.Petersburg University research grant 6.38.335.2015. The Abastumani team acknowledges financial support of the project FR/638/6320/12 by the Shota Rustaveli National Science Foundation under contract 31/77. We acknowledge the photometric observations from the AAVSO International Database contributed by observers worldwide and used in this research. This paper has made use of uptodate SMARTS optical/nearinfrared light curves that are available at http://www.astro.yale.edu/smarts/glast/Bonning et al. (2012). We acknowledge the contributions of Y. Y. Kovalev to the Quasar Movie Project. We thank the anonymous referee for a constructive review that improved this paper.
References
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, Nature, 463, 919 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Agudo, I., Molina, S. N., Gómez, J. L., et al. 2012, Int. J. Mod. Phys. Conf. Ser., 8, 299 [CrossRef] [Google Scholar]
 Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014a, A&A, 567, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2014b, A&A, 569, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Björnsson, C.I. 1982, ApJ, 260, 855 [NASA ADS] [CrossRef] [Google Scholar]
 Blinov, D., Pavlidou, V., Papadakis, I., et al. 2015, MNRAS, 453, 1669 [NASA ADS] [CrossRef] [Google Scholar]
 Bonning, E., Urry, C. M., Bailyn, C., et al. 2012, ApJ, 756, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Cousins, A. W. J. 1976, Monthly Notes of the Astronomical Society of Southern Africa, 35, 70 [Google Scholar]
 D’Arcangelo, F. D., Marscher, A. P., Jorstad, S. G., et al. 2007, ApJ, 659, L107 [NASA ADS] [CrossRef] [Google Scholar]
 Doroshenko, V. T., Sergeev, S. G., Merkulova, N. I., et al. 2005, Astrophysics, 48, 156 [NASA ADS] [CrossRef] [Google Scholar]
 GonzálezPérez, J. N., Kidger, M. R., & MartínLuis, F. 2001, AJ, 122, 2055 [NASA ADS] [CrossRef] [Google Scholar]
 Hayashida, M., Nalewajko, K., Madejski, G. M., et al. 2015, ApJ, 807, 79 [Google Scholar]
 Holmes, P. A., Brand, P. W. J. L., Impey, C. D., et al. 1984, MNRAS, 211, 497 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, T. W. 1988, ApJ, 332, 678 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, T. W., Rudnick, L., Aller, H. D., et al. 1985, ApJ, 290, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Jorstad, S. G., Marscher, A. P., Larionov, V. M., et al. 2010, ApJ, 715, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Kawabata, K. S., Nagae, O., Chiyonobu, S., et al. 2008, in SPIE Conf. Ser., 7014 [Google Scholar]
 Kiehlmann, S. 2015, Ph.D. Thesis, Universität zu Köln, Germany [Google Scholar]
 Kikuchi, S., Mikami, Y., Inoue, M., Tabara, H., & Kato, T. 1988, A&A, 190, L8 [NASA ADS] [Google Scholar]
 King, O. G., Blinov, D., Ramaprakash, A. N., et al. 2014, MNRAS, 442, 1706 [NASA ADS] [CrossRef] [Google Scholar]
 Kinman, T. D. 1967, ApJ, 148, L53 [NASA ADS] [CrossRef] [Google Scholar]
 Komissarov, S. S., Barkov, M. V., Vlahakis, N., & Königl, A. 2007, MNRAS, 380, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Konigl, A., & Choudhuri, A. R. 1985, ApJ, 289, 188 [NASA ADS] [CrossRef] [Google Scholar]
 Kurtanidze, O. M., & Nikolashvili, M. G. 1999, Proc. of the OJ94 Annual Meeting 1999, Blazar Monitoring toward the Third Millennium, 25 [Google Scholar]
 Kurtanidze, O. M., & Nikolashvili, M. G. 2002, Blazar Astrophysics with BeppoSAX and Other Observatories, ASI SP, 189 [Google Scholar]
 Laing, R. A. 1980, MNRAS, 193, 439 [NASA ADS] [Google Scholar]
 Larionov, V. M., Jorstad, S. G., Marscher, A. P., et al. 2008, A&A, 492, 389 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Larionov, V. M., Jorstad, S. G., Marscher, A. P., et al. 2013, ApJ, 768, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Longair, M. S. 2011, High Energy Astrophysics, 3rd edn. (Cambridge University Press) [Google Scholar]
 Marscher, A. P. 2014, ApJ, 780, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Marscher, A. P. 2015, in IAU Symp., 313, eds. F. Massaro, C. C. Cheung, E. Lopez, & A. Siemiginowska, 122, 127 [Google Scholar]
 Marscher, A. P., Jorstad, S. G., D’Arcangelo, F. D., et al. 2008, Nature, 452, 966 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Marscher, A. P., Jorstad, S. G., Larionov, V. M., et al. 2010, ApJ, 710, L126 [NASA ADS] [CrossRef] [Google Scholar]
 Nalewajko, K. 2010, Int. J. Mod. Phys. D, 19, 701 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nalewajko, K., & Begelman, M. C. 2012, MNRAS, 427, 2480 [NASA ADS] [CrossRef] [Google Scholar]
 Pavlidou, V., Angelakis, E., Myserlis, I., et al. 2014, MNRAS, 442, 1693 [NASA ADS] [CrossRef] [Google Scholar]
 Sasada, M., Uemura, M., Fukazawa, Y., et al. 2011, PASJ, 63, 489 [NASA ADS] [Google Scholar]
 Sironi, L., Petropoulou, M., & Giannios, D. 2015, MNRAS, 450, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, P. S., Montiel, E., Rightley, S., et al. 2009, ArXiv eprints [arXiv:0912.3621] [Google Scholar]
 Sorcia, M., Benítez, E., Hiriart, D., et al. 2013, ApJS, 206, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Sorcia, M., Benítez, E., Hiriart, D., et al. 2014, ApJ, 794, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Stawarz, Ł., & Petrosian, V. 2008, ApJ, 681, 1725 [NASA ADS] [CrossRef] [Google Scholar]
 Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Vlahakis, N., & Königl, A. 2004, ApJ, 605, 656 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, H., Chen, X., & Böttcher, M. 2014, ApJ, 789, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, H., Chen, X., Böttcher, M., Guo, F., & Li, H. 2015, ApJ, 804, 58 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: EVPA ambiguity, consistency and reliability
In the attempt to solve the nπambiguity of the polarization angle, we generally adjust EVPA curves applying Eq. (1), in the following we call this method E1. Here, we introduce a second method (E2) to assess the sampling quality and to test the reliability of solving the phase wraps. E2 determines the median {·} of the N_{ref} previous data points as a reference for an EVPA data point χ_{i}: (A.1)The data point χ_{i} is then shifted according to Eq. (1), where χ_{i,ref} replaces χ_{i − 1}. For N_{ref} = 1 method E2 is identical to E1. For N_{ref}> 1 method E2 considers a longer time frame than E1 as reference for adjusting each data point. It is expected to be less susceptible to outliers but less reliable for sparsely sampled rotations.
The measured EVPA curve χ_{obs} is the 180°modulo of the intrinsic curve χ_{intr}. EVPA adjustment methods E1 and E2 aim to reconstruct the intrinsic EVPA variations from the sampled EVPA curve. We call the probability that the shapes of the adjusted and intrinsic EVPA curves are identical P(χ_{adj} = χ_{intr} ± n·180°) (with n ∈N) the reliability of each method. In the following we test the reliability under various conditions.
Assuming a constant rotation rate , constant time sampling Δt, and no observational errors the EVPA curve is correctly reconstructed if the sampled rotation . Without observational noise N_{ref} = 1 is the optimal choice allowing the sparsest time sampling. For sampled rotations Δχ< 60° methods E1 (N_{ref} = 1) and E2 (N_{ref} ≥ 2) may yield the same adjusted EVPA curve. We call the highest number of reference points N_{cons}, for which all adjusted curves with N_{ref} ≤ N_{cons} are identical, the consistency level: (A.2)Assuming minimal variation between measured data points, the consistency level allows one to estimate the reliability of the EVPA curve reconstruction.
We determine the expected reliability of the EVPA adjustment (E1 and E2) and the expected consistency level for various intrinsic rotation rates given the time sampling and observational errors of our combined polarization data. We simulate linearly increasing EVPA curves up to a total amplitude of 360° with random time sampling and random errors as in the random walk models (cf. Sect. 4) with model parameters given in Table 3, row “total”. The upper panel of Fig. A.1 shows the probability of correctly reconstructing the intrinsic EVPA variation for various intrinsic rotation rates and different numbers of reference points. Method E1 (N_{ref} = 1) is the most reliable for rotation rates , increasing N_{ref} reduces the reliability. For intrinsic rotation rates the reliability is ~98%. In the lower panel of Fig. A.1 we show the expectation value and standard deviation of the consistency level for different intrinsic rotation rates. High consistency rates N_{cons} ≳ 10 indicate low intrinsic variation and a high probability ~98% of correct curve reconstruction.
Fig. A.1 Probability of the correct EVPA reconstruction using different numbers of reference points N_{ref} (upper panel) and average consistency levels (lower panel) over different rotation rates dχ/ dt. 
Given the time sampling and typical uncertainties of our polarization data, method E2 does not improve the reliable reconstruction of the intrinsic EVPA variability. Nevertheless, using both methods E1 and E2 and the consistency level allows us to estimate the quality of the time sampling and the reliability of the reconstructed EVPA curve.
Appendix B: Variation estimator biases
We have defined a quantitative measure of the smoothness of the EVPA curve in Sect. 3.2. This variation estimator measures the average, short timescale (shorter than the analysed data total time), erratic rotation rate of the EVPA corrected for an assumed, underlying, secular trend at the timescale of the analysed data total time. The measured variation estimator is increased (a) by measurement errors, introducing additional variation; and (b) by curvature of the underlying smooth variation, i.e. a nonconstant trend.
Intrinsic EVPA variation and measurement errors independently add to the variation estimator. The debiased variation estimator is given by: (B.1)For constant measurement uncertainties σ_{χ} and even time sampling Δt, the error bias of the variation estimator is s_{bias} ∝ σ_{χ}/ Δt. Equation (B.1) also holds for variable σ_{χ}, but then we need to estimate s_{bias} through simulations, which is how the values in Table 2 were estimated. The uncertainty of the debiased variation estimator is equal to the uncertainty of the error bias.
When calculating the variation estimator s, the local derivative is corrected for (intrinsic) first order rotations, i.e. constant rotation rates, by subtracting a trend . For higher order (i.e. curved) intrinsic rotation curves the variation estimator will be biased depending on the curvature. We test this dependency by simulating EVPA curves (without errors) that follow a power law over time (B.2)with χ_{max} the EVPA at time T equal to the total rotation amplitude during the total time interval T and with α ≥ 1 the power law index. Example curves for α = 1,2,··· ,10 are shown in the upper panel of Fig. B.1.
Fig. B.1 The upper panel shows exemplary power law curves with increasing curvature, i.e. power law index α from 1 (red) to 10 (blue). The lower panel shows the curvature factor (see text for description). 
The trend estimates the total rotation rate given by with an accuracy of 1%. The curvature biased variation estimator is given by the curvature factor f_{s} (plotted in the lower panel of Fig. B.1) and the total rotation rate estimated by the trend : (B.3)Equation (B.3) with the simulation result shown in Fig. B.1 allows us to estimate the impact of a curved, underlying, smooth variation of the EVPA onto the variation estimator. Correcting the variation estimator for higher order variation would require fitting the EVPA curve and subtracting the fit from the local derivatives. To avoid a priori fitting an arbitrary function to the data we stick to the criterion that a linear trend is considered smooth variation and any deviation from that linear trend contributes to the variation estimator.
Appendix C: EVPA rotation identification and distribution
We define an algorithm to automatically identify EVPA rotations in a polarization curve based on strictly fixed criteria. We use this algorithm on polarization curves simulated with our random walk processes to characterize the EVPA rotation amplitude distribution expected from a stochastic process.
C.1. EVPA rotation identification algorithm
We define the start and end point of an EVPA rotation by a change of the rotation direction, i.e. a sign change of the derivative of the EVPA curve. We do not consider a change of the rotation rate. We define the rotation amplitude as the difference of the (local) extrema: Δχ = χ_{max} − χ_{min}. We call a rotation between two data points χ_{i}, χ_{j} significant when the amplitude is larger than the root summed squared errors of the two data points multiplied with a factor ς that characterizes the significance level: (C.1)The algorithm identifies significant rotations by iterating pointwise through the data. The storage contains all identified, significant EVPA rotations. The current rotation contains all data points between the end of the last significant rotation and the data point previous to the current one. The rotation amplitude, direction and significance characterize the current rotation. They are calculated from the minimum and maximum value and the corresponding uncertainties. The data point of the current iteration step and the previous one make the new rotation.
Description of the occurring states and the corresponding operations for the EVPA rotation identification.
Based on the direction and significance of the new, current, and last stored rotation we define eleven states and corresponding operations in Table C.1. Five quantities define each state: the significance of the current and the new rotation (Y/N, Cols. 2 and 3), whether the current rotation changes its direction or not (+/–, Col. 5), whether the current rotation becomes significant by adding the new, current data point(Y/N, Col. 4), and how the newly significant rotation is oriented regarding to the previously stored rotation (+/–, Col. 6). States 6 and 9 may only occur once when no significant rotation has been found yet and the rotation direction cannot be compared to a previous rotation (n/a). This scheme only identifies local maxima and minima, which are significant at a predefined level. Each pair of adjacent identified extrema marks a continuous rotation, which may contain several local, insignificant maxima and minima.
C.2. Random walk EVPA rotation distribution
Fig. C.1 Distribution of the absolute amplitudes Δχ_{rot} of rotations identified at 3σlevel in siQU model simulations with different cell numbers and cell variation rates. Red line: distribution of twodatapoint rotations, blue line: distribution of multipledatapoint rotations, grey filled region: distribution of all rotations. 
We run random walk simulations for various cell numbers and cell variation rates with a total time T = 100 000 d and the other parameters chosen as in Table 3, period “total”. We adjust the EVPA curve with Eq. (1) and identify significant rotations at a significance level ς = 3. In Fig. C.1 we show four examples of EVPA rotation amplitude distributions based on the siQU random walk process with different numbers of cells and cell variation rates. The grey filled region displays the full distribution. Two populations of identified rotations contribute to this distribution: rotations consisting of two data points only, N_{dp} = 2, and rotations containing more than two data points, N_{dp}> 2. The rotation amplitude distributions of the populations are shown as red and blue curves, respectively.
Both distributions peak at a certain value, show a sharp decline towards smaller amplitudes and a long tail towards larger amplitudes. However, because of the nπambiguity the twodatapointpopulation distribution is cut off at ~90°. Slightly larger rotation amplitudes than 90° are possible owing to EVPA errors. The shortest simulation time step and the limits of the tested parameter space (cell number and cell variation rate) affect the lowamplitude tail of the distributions.
When the cell variation rate increases (relative to the number of cells) the EVPA variation becomes more erratic and large, continuous rotations become less likely. The distribution shifts towards smaller rotation amplitudes. When the number of cells is reduced the EVPA variation becomes more erratic and the EVPA uncertainty is lower as the average polarization fraction is increased. As a result, twodatapoint rotations are more likely to be significant and the number of these rotations increases, whereas the number of small rotation consisting of multiple data points decreases.
Though not shown here, the significance level ς used in the rotation identification strongly affects the number of identified rotations and the rotation amplitude distribution. For larger ς twodatapointrotations are less likely to be significant and therefore the number of identified rotation tends to be lower, and rotation amplitudes tend to be larger compared to a lower choice of ς.
All Tables
Observatories contributing optical polarization data to the Quasar Movie Project with telescope diameters, filters, number of polarization data points, and the type of data with references to the data reduction in the table footnotes.
Truncated power law distribution parameters modelling the random time steps (Cols. 2–4).
Highest probabilities of producing polarization curves fulfilling the conditions on the polarization fraction (PF1/2), on the EVPA variation (PA1/2) and fulfilling all conditions from three different random walk processes with the corresponding, optimal model parameters for period IIa.
Highest probabilities of producing polarization curves fulfilling the conditions on the polarization fraction (PF1/2), on the EVPA variation (PA1/2) and fulfilling all conditions from three different random walk processes with the corresponding, optimal model parameters for the period IIIc.
Description of the occurring states and the corresponding operations for the EVPA rotation identification.
All Figures
Fig. 1 Optical photometry and polarimetry and γray light curve of 3C 279. FermiLAT γray light curve at > 100 MeV binned into 3 day intervals (panel a)) as published in Hayashida et al. (2015). Combined Rband light curve (panel b)). Measured, optical polarization fraction (panel c)) and EVPA (panel d)); red circles: Calar Alto (R), red squares: CrAO70 cm (R), red diamonds: Perkins (R), orange upsided triangles: SPM (R), orange rightsided triangles: St. Petersburg (R), green downsided triangles: KANATA (V), green leftsided triangles: Steward Obs. (spec. and V), blue circles: Liverpool (V+R), blue squares: KVA (white light). Combined, debiased, and averaged polarization fraction (panel e)). Combined, averaged, and adjusted EVPA (panel f)); open symbols are added from the nonaveraged EVPA curve. Pointwise, local derivative of the adjusted EVPA (panel g)). The grey area highlights the period of γray flaring activity coinciding with a rotation of the optical polarization angle. 

In the text 
Fig. 2 208° EVPA rotation presented in Abdo et al. (2010) based on data from the KANATA telescope. Open squares show the EVPAs within the initial 0–180°interval, black circles the adjusted EVPA curve. The numbers show the difference between adjacent data points and the corresponding uncertainties in degrees. 

In the text 
Fig. 3 EVPA data as shown in Fig. 2 with additional data from Calar Alto, Perkins, St. Petersburg, and Steward Obs. 

In the text 
Fig. 4 Empirical cumulative distribution function (ECDF) of the combined Stokes q and u uncertainties in the entire data (red curve), period IIa (orange curve), and period IIIc (green curve). 

In the text 
Fig. 5 Expectation values of the polarization fraction mean (upper panel) and standard deviation (lower panel) for different numbers of cells depending on the cell variation rate. Each data point is based on 10 000 simulations. Solid curves and circles correspond to the siQU model, dashed curves and squares to the orQUc model, and dotted curves and triangles to the orQUd model. 

In the text 
Fig. 6 Distributions of a) the mean and b) the standard deviation of the polarization fraction; c) the EVPA amplitude; and d) the EVPA variation estimator for each random walk process with a range of model parameters. 

In the text 
Fig. 7 Probability distributions over the parameter space N_{cells}, n_{var} for three different conditions (rows 1–3) based on three random walk processes (Cols. 1–3) for period IIa. 

In the text 
Fig. 8 Distribution of the EVPA rotation amplitude (upper panel) and distribution of the EVPA variation estimator (lower panel) for all simulations with a mean polarization fraction consistent with the observed value during period IIa. The number of selected simulations is indicated in the legend. The corresponding observed values are indicated by black dotted lines. 

In the text 
Fig. 9 Probability distributions over the parameter space N_{cells}, n_{var} for three different conditions (rows 1–3) based on three random walk processes (Cols. 1–3) for the period IIIc. 

In the text 
Fig. 10 Distribution of the EVPA rotation amplitude (upper panel) and distribution of the EVPA variation estimator (lower panel) for all simulations with a mean polarization fraction consistent with the observed value during the period IIIc. The number of selected simulations is indicated in the legend. The corresponding observed values are indicated by black dotted lines. 

In the text 
Fig. 11 Histogram of the rotation amplitudes of the 109 rotations identified in the optical EVPA curve at ς = 3significance and the corresponding cumulative distribution function (ECDF, red curve). 

In the text 
Fig. 12 pvalues of the KS test over a grid of model parameters, testing the rotation amplitude distribution of the entire EVPA curve against random walk simulations for three different stochastic processes. 

In the text 
Fig. A.1 Probability of the correct EVPA reconstruction using different numbers of reference points N_{ref} (upper panel) and average consistency levels (lower panel) over different rotation rates dχ/ dt. 

In the text 
Fig. B.1 The upper panel shows exemplary power law curves with increasing curvature, i.e. power law index α from 1 (red) to 10 (blue). The lower panel shows the curvature factor (see text for description). 

In the text 
Fig. C.1 Distribution of the absolute amplitudes Δχ_{rot} of rotations identified at 3σlevel in siQU model simulations with different cell numbers and cell variation rates. Red line: distribution of twodatapoint rotations, blue line: distribution of multipledatapoint rotations, grey filled region: distribution of all rotations. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.