© 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 1510-089 the counter-clockwise ~720° rotation in 2009 (Marscher et al. 2010) and the counter-clockwise ~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) time-scales 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 non-axisymmetric 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. Non-axisymmetric 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 cross-section, on a helical trajectory that probes different parts of a large-scale 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 multi-wavelength 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 multi-wavelength photometry and polarimetry data analysis of the quasar 3C 279” (Kiehlmann 2015).

2. Data

An intensive VLBI and multi-waveband monitoring campaign targeted on quasar 3C 279, including observatories covering radio, millimetre, infrared, optical, ultraviolet, X-ray, 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 R-band 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 long-term 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 AZT-8 telescope of the Crimean Astrophysical Observatory (CrAO; Nauchnij, Russia); (4) the 1.5 m KANATA telescope of the Higashi-Hiroshima 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 Inter-American Observatory (Chile); (10) the 84 cm telescope of the Observatorio Astronómico Nacional (San Pedro Mártir, Mexico); (11) the 40 cm LX-200 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¹ 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.

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 AZT-8 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², 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 AZT-8 CrAO telescope images were reduced (see Doroshenko et al. 2005) using the star N9 from González-Pérez et al. (2001).

We combined all R-band photometry into one light curve after cross-calibrating all data sets with respect to the SMARTS R-band 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 time-scales. 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 counter-clockwise from North. An increase of the polarization angle refers to a counter-clockwise rotation projected on the sky.

Before analysing the data we averaged data points within 0.5 days in the R-band 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 time-scales 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 R-band 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.

3.1. The nπ-ambiguity of the polarization vector

The EVPA, χ, is defined within a 180°-interval³. 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 Δχ = χ₂ − χ₁ ± 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 ${\mathit{\chi }}_{\mathit{i,}\mathrm{adj}}\mathrm{=}{\mathit{\chi }}_{\mathit{i}}\mathrm{-}\mathit{n\pi }\mathrm{with}\mathit{n}\mathrm{=}\mathrm{int}\left(\frac{{\mathit{\chi }}_{\mathit{i}}\mathrm{-}{\mathit{\chi }}_{\mathit{i}\mathrm{-}\mathrm{1}}}{\mathit{\pi }}\right)\mathit{,}$(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. $\mathrm{\Delta }{\mathit{\chi }}_{\mathrm{red}}\mathrm{=}\left|{\mathit{\chi }}_{\mathit{i}}\mathrm{-}{\mathit{\chi }}_{\mathit{i}\mathrm{-}\mathrm{1}}\right|\mathrm{-}\sqrt{{\mathit{e}}_{{\mathit{\chi }}_{\mathit{i}}}^{\mathrm{2}}\mathrm{+}{\mathit{e}}_{{\mathit{\chi }}_{\mathit{i}\mathrm{-}\mathrm{1}}}^{\mathrm{2}}}\mathit{>}\frac{\mathit{\pi }}{\mathrm{2}}$, 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 χ₁ = 40° and χ₂ = 140° in the initial interval ${\mathrm{[}}^{}\mathrm{0}\mathit{,}{\mathrm{180}}^{\mathrm{◦}}{\mathrm{)}}^{}$, both with uncertainty e_χ = 10°. Then, Δχ_red = 85.9° and χ₂ is not shifted, yielding a counter-clockwise rotation of χ₂ − χ₁ = 100°. Considering the same data in the initial interval ${\mathrm{[}}^{}\mathrm{-}\mathrm{90}\mathit{,}{\mathrm{90}}^{\mathrm{◦}}{\mathrm{)}}^{}$, i.e. χ₁ = 40° and χ₂ = − 40°, gives Δχ_red = 65.9°. Again, χ₂ 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: ${\left(\frac{\mathrm{\Delta }\mathit{\chi }}{\mathrm{\Delta }\mathit{t}}\right)}_{\mathit{i}}\mathrm{=}\frac{{\mathit{\chi }}_{\mathit{i}}\mathrm{-}{\mathit{\chi }}_{\mathit{i}\mathrm{-}\mathrm{1}}}{{\mathit{t}}_{\mathit{i}}\mathrm{-}{\mathit{t}}_{\mathit{i}\mathrm{-}\mathrm{1}}}\mathrm{·}$(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: ${\mathit{s}}_{\mathit{i}}\mathrm{=}{\left(\frac{\mathrm{\Delta }\mathit{\chi }}{\mathrm{\Delta }\mathit{t}}\right)}_{\mathit{i}}\mathrm{-}⟨\left(\frac{\mathrm{\Delta }\mathit{\chi }}{\mathrm{\Delta }\mathit{t}}\right)⟩\mathit{.}$(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: $\mathit{s}\mathrm{=}⟨\left|{\mathit{s}}_{\mathit{i}}\right|⟩\mathrm{=}⟨\left|{\left(\frac{\mathrm{\Delta }\mathit{\chi }}{\mathrm{\Delta }\mathit{t}}\right)}_{\mathit{i}}\mathrm{-}⟨\left(\frac{\mathrm{\Delta }\mathit{\chi }}{\mathrm{\Delta }\mathit{t}}\right)⟩\right|⟩\mathit{.}$(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₁ is considered smoother than a second curve with s₂ if s₂>s₁.

The standard deviation of the pointwise derivative $\sqrt{{}^{\mathrm{⟨}}{\mathit{s}}_{\mathit{i}}^{\mathrm{2}}^{\mathrm{⟩}}}$ could be used as an alternative approach to quantify the EVPA variability. Averaging over ${\mathit{s}}_{\mathit{i}}^{\mathrm{2}}$ 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 non-linear 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 R-band 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 R-band 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 smooth-looking 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 source-intrinsic. 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 counter-clockwise rotation.

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°.

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 Fermi-LAT γ-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 long-term γ-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 current-driven 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 second-order 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.

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 sub-set 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 maximum-likelihood 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 ${\mathit{N}}_{\mathrm{var}}\mathrm{\left(}{\mathit{t}}_{\mathit{i}}\mathit{,}{\mathit{t}}_{\mathit{i}\mathrm{-}\mathrm{1}}\mathrm{\right)}\mathrm{=}\mathrm{int}{}^{\mathrm{(}}\mathit{n}_{\mathrm{var}}{\mathit{t}}_{\mathit{i}}{}^{\mathrm{)}}\mathrm{-}\mathrm{int}{}^{\mathrm{(}}\mathit{n}_{\mathrm{var}}{\mathit{t}}_{\mathit{i}\mathrm{-}\mathrm{1}}{}^{\mathrm{)}}$, 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) ${\mathit{m}}_{\mathit{l,}\max }\mathrm{=}\frac{\mathit{p}\mathrm{-}\mathrm{1}}{\mathit{p}\mathrm{-}\frac{\mathrm{7}}{\mathrm{3}}}\mathit{,}$(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: $\begin{array}{ccc}\mathit{Q̂}\mathit{i}& & \mathrm{~}\mathrm{𝒩}\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{1}\mathrm{\right)}\mathit{,}\\ \mathit{Û}\mathit{i}& & \mathrm{~}\mathrm{𝒩}\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{1}\mathrm{\right)}\mathit{,}\\ {\mathit{I}}_{\mathit{i}}& \mathrm{=}& \frac{\mathit{I}}{{\mathit{N}}_{\mathrm{cells}}}\mathrm{=}\mathrm{const}\mathit{.}\mathit{,}\\ {\mathit{Q}}_{\mathit{i}}& \mathrm{=}& \frac{\mathit{Q̂}\mathit{i}}{\sqrt{\mathit{Q̂}\begin{array}{c}\mathrm{2}\\ \mathit{i}\end{array}\mathrm{+}\mathit{Û}\begin{array}{c}\mathrm{2}\\ \mathit{i}\end{array}}}\mathrm{·}{\mathit{I}}_{\mathit{i}}\mathrm{·}{\mathit{m}}_{\mathit{l,}\max }\mathit{,}\\ {\mathit{U}}_{\mathit{i}}& \mathrm{=}& \frac{\mathit{Û}\mathit{i}}{\sqrt{\mathit{Q̂}\begin{array}{c}\mathrm{2}\\ \mathit{i}\end{array}\mathrm{+}\mathit{Û}\begin{array}{c}\mathrm{2}\\ \mathit{i}\end{array}}}\mathrm{·}{\mathit{I}}_{\mathit{i}}\mathrm{·}{\mathit{m}}_{\mathit{l,}\max }\mathit{.}\end{array}$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: ${\mathit{I}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{2}\mathit{I}}{{\mathit{N}}_{\mathrm{cells}}^{\mathrm{2}}\mathrm{+}{\mathit{N}}_{\mathrm{cells}}}\mathrm{(}{\mathit{N}}_{\mathrm{cells}}\mathrm{-}\mathit{i}\mathrm{)}\mathit{.}$(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 $\mathit{I}\mathrm{=}{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{cells}}}{\mathit{I}}_{\mathit{i}}$ and Q, U, accordingly, determine the integrated intensity I, the linear polarization fraction m_l and the EVPA χ at each simulation time step: $\begin{array}{ccc}{\mathit{m}}_{\mathit{l}}& \mathrm{=}& \frac{\sqrt{{\mathit{Q}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}^{\mathrm{2}}}}{\mathit{I}}\mathit{,}\\ \mathit{\chi }& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\arctan \frac{\mathit{U}}{\mathit{Q}}\mathrm{+}\mathit{n}\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{with}\mathit{n}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}\end{array}$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 self-absorption and Faraday rotation on the polarized signal at optical frequencies.

4.5. Simulated EVPA errors

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 u-uncertainties 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 log-normal distribution. The distribution parameters μ_unc and σ_unc are estimated for each period individually. To simulate EVPA “measurement” errors we draw a random-number σ_χi from a log-normal 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.

The mean and standard deviation of the simulated polarization fraction follow log-normal 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.

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 $\mathit{E}\mathrm{\left[}⟨{\mathit{m}}_{\mathit{l}}⟩\mathrm{\right]}\mathrm{\propto }\mathrm{1}\mathit{/}\sqrt{{\mathit{N}}_{\mathrm{cells}}}$. 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 $\mathit{E}\mathrm{\left[}\mathit{\sigma }\mathrm{\right(}{\mathit{m}}_{\mathit{l}}\mathrm{\left)}\mathrm{\right]}\stackrel{{\mathit{n}}_{\mathrm{var}}\mathrm{\to }{\mathit{N}}_{\mathrm{cells}}\mathit{/}\mathrm{\Delta }\mathit{t}}{\mathrm{-}0\mathrm{\to }}\frac{\mathrm{1}}{\mathrm{2}}\mathit{E}\mathrm{\left[}⟨{\mathit{m}}_{\mathit{l}}⟩\mathrm{\right]}$.

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 $\mathrm{Var}\mathrm{(}⟨{\mathit{m}}_{\mathit{l}}⟩{\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{=}\mathrm{0.058}$ and $\mathrm{Var}\mathrm{\left(}\mathit{\sigma }\mathrm{\right(}{\mathit{m}}_{\mathit{l}}\mathrm{)}{\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{=}\mathrm{0.024}$. The distribution is also right-tailed. 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 (${\mathit{N}}_{\mathrm{cells}}^{\mathrm{opt}}$, ${\mathit{n}}_{\mathrm{var}}^{\mathrm{opt}}$). 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: $\begin{array}{ccc}\begin{array}{c}\\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \end{array}& & \end{array}$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.

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: $\begin{array}{ccc}\begin{array}{c}\\ & & & & & & \\ & & & & & & \\ & & & & & & \\ & & & & & & \end{array}& & \end{array}$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).

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.

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 non-stochastic 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 low-brightness state in R-band. 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 non-stochastic 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 IIIa-c originate in the same event, the sharp clockwise rotation opposite to the counter-clockwise direction of the enclosing rotations poses a challenge to any non-stochastic 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.

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 two-sample Kolmogorov-Smirnov (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.

The p-values 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 low-brightness state and the flaring state of 3C 279. During the low-brightness 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 signal-to-noise data and underestimated uncertainties, but it is source-intrinsic. 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 low-brightness 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 low-brightness 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 low-brightness 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 multi-zone 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 sub-flares in the R-band light curve during the EVPA rotation, these consecutive flaring regions would likely need significant fine-tuning 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 light-travel 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 non-stochastic ~360° EVPA rotation during the period IIIc can be produced by the non-axisymmetric 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 non-flaring and in the flaring state can be explained with the two-processes interpretation, as a stochastic process during the non-flaring 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 two-directional EVPA swings such as the two-component 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 counter-rotation 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 two-component model of Holmes et al. (1984), the superposition of both regions may temporarily produce a counter-clockwise 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 non-deterministic 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 1510-089

3C 279, PKS 1510-089, 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 1510-089 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 1510-089 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 1510-089 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 1510-089 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 1510-089 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 multi-wavelength campaign, we have accumulated and combined multiple data sets of R-band photometry and optical polarimetry measurements into well-sampled 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 low-brightness 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 counter-clockwise 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 low-brightness 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 low-brightness 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 counter-clockwise 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 low-brightness, 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 13-02-12103 and 14-02-31789. N.G.B. was supported by the RFBR grant 12-02-01237a. E.B., M.S. and D.H. thank financial support from UNAM DGAPA-PAPIIT through grant IN116211-3. I.A. acknowledges support by a Ramon y Cajal grant of the Spanish Ministry of Economy and Competitiveness (MINECO). The research at the IAA-CSIC and the MAPCAT program are supported by the Spanish Ministry of Economy and Competitiveness and the Regional Government of Andalucía (Spain) through grants AYA2010-14844, AYA2013-40825-P, and P09-FQM-4784. The Calar Alto Observatory is jointly operated by the Max-Planck-Institut für Astronomie and the Instituto de Astrofísica de Andalucía-CSIC. 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 15-02-00949 and St.Petersburg University research grant 6.38.335.2015. The Abastumani team acknowledges financial support of the project FR/638/6-320/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 up-to-date SMARTS optical/near-infrared 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

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: ${\mathit{\chi }}_{\mathit{i,}\mathrm{ref}}\mathrm{=}\left\{\left[{\mathit{\chi }}_{\mathit{i}\mathrm{-}\mathrm{1}\mathrm{-}{\mathit{N}}_{\mathrm{ref}}}\mathit{,}\mathit{...}\mathit{,}{\mathit{\chi }}_{\mathit{i}\mathrm{-}\mathrm{1}}\right]\right\}\mathit{.}$(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: ${\mathit{N}}_{\mathrm{const}}\mathrm{=}\begin{array}{c}\max \end{array}\mathrm{\{}\mathit{N}\mathrm{\in }{N}_{\mathrm{\ge }\mathrm{1}}\mathrm{:}{\mathit{\chi }}_{\mathrm{adj}\mathit{,}\mathrm{(}{\mathit{N}}_{\mathrm{ref}}\mathrm{=}\mathrm{1}\mathrm{)}}\mathrm{=}{\mathit{\chi }}_{\mathrm{adj}\mathit{,}\mathrm{(}{\mathit{N}}_{\mathrm{ref}}\mathrm{=}\mathit{i}\mathrm{)}}\mathrm{\forall }\mathit{i}\mathrm{\le }{\mathit{N}}^{\mathrm{\}}}\mathit{.}$(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.

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 time-scale (shorter than the analysed data total time), erratic rotation rate of the EVPA corrected for an assumed, underlying, secular trend at the time-scale 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 non-constant trend.

Intrinsic EVPA variation and measurement errors independently add to the variation estimator. The de-biased variation estimator is given by: ${\mathit{s}}_{\mathrm{debiased}}\mathrm{=}\sqrt{{\mathit{s}}_{\mathrm{obs}}^{\mathrm{2}}\mathrm{-}{\mathit{s}}_{\mathrm{bias}}^{\mathrm{2}}}\mathrm{for}{\mathit{s}}_{\mathrm{obs}}\mathit{>}{\mathit{s}}_{\mathrm{bias}}\mathit{.}$(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 de-biased 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 $\mathit{\chi }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{\chi }}_{\max }{\left(\frac{\mathit{t}}{\mathit{T}}\right)}^{\mathit{\alpha }}\mathit{,}$(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.

The trend estimates the total rotation rate given by $\frac{\mathrm{\Delta }\mathit{\chi }}{\mathrm{\Delta }\mathit{t}}\mathrm{=}\frac{{\mathit{\chi }}_{\max }}{\mathit{T}}$ 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 $\frac{\mathrm{\Delta }\mathit{\chi }}{\mathrm{\Delta }\mathit{t}}$ estimated by the trend : ${\mathit{s}}_{\mathrm{curve}}\mathrm{=}{\mathit{f}}_{\mathit{s}}\mathrm{·}\frac{\mathrm{\Delta }\mathit{\chi }}{\mathrm{\Delta }\mathit{t}}\mathrm{\approx }{\mathit{f}}_{\mathit{s}}\mathrm{·}\mathit{\chi ̅}\mathit{t}\mathit{.}$(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: $\begin{array}{ccc}\left|\mathrm{\Delta }{\mathit{\chi }}_{\mathit{ij}}\right|\mathrm{=}\left|{\mathit{\chi }}_{\mathit{i}}\mathrm{-}{\mathit{\chi }}_{\mathit{j}}\right|\mathit{>}\mathit{\varsigma }\mathrm{·}\sqrt{{\mathit{\sigma }}_{{\mathit{\chi }}_{\mathit{i}}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{{\mathit{\chi }}_{\mathit{j}}}^{\mathrm{2}}}\mathit{.}& & \end{array}$(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.

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

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 two-data-point-population 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 low-amplitude 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, two-data-point 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 ς two-data-point-rotations 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

All Figures