Diagnosing Magnetic Field Geometry in Blazar Jets Using Multi-Frequency, Centimeter-Band Polarimetry and Radiative Transfer Modeling

Abstract

We use multi-frequency linear polarization observations from the University of Michigan blazar program (UMRAO), in combination with radiative transfer simulations of emission from a relativistic jet, to investigate the time-dependent flow conditions, including magnetic field geometry, in an example blazar OT 081. We adopt a scenario incorporating relativistic shocks during flaring, and both ordered axial and helical magnetic field components and magnetic turbulence in the underlying flow; these constituents are consistent with the observed periods of ordered behavior in the polarization intermixed with stochastic variations. The simulations are able to reproduce the global features of the observed light curves, including amplitude and spectral evolution of the linear polarization, during four time periods spanning 25 years. From the simulations, we identify the signature of a weak-to-strong helical magnetic field on the polarization, but conclude that a dominant helical magnetic field is not consistent with the UMRAO polarization data. The modeling identifies time-dependent changes in the ratio of the ordered-to-turbulent magnetic field, and changes in the flow direction and Lorentz factor. These suggest the presence of jet-like structures within a broad envelope seen at different orientations.

1. Introduction

Magnetic fields in blazar jets are important, since they control the dynamics of the formation and evolution of the jet and impact the acceleration of particles, giving rise to non-thermal emission across the spectrum from the radio to TeV bands, yet they cannot be observed directly. While it is theoretically predicted and generally accepted that the magnetic fields are predominantly helical near the central engine of blazars, and that this geometry is expected from the combination of jet outflow and a rotating central black hole and accretion disk, controversy remains about the geometry at parsec and sub-parsec scales. Polarimetry data are a particularly powerful tool for studying the geometries of these magnetic fields indirectly, but such studies require well-sampled data in view of the rapid variability of the polarized emission (as short as nightly in the optical band and weekly at centimeter wavelengths) and multi-year time windows to capture the full range of the variations. Multi-frequency polarization data at radio band, unlike total flux density measurements alone, provide strong constraints for evaluating models exploring a range of scenarios, and additionally they identify Faraday effects which potentially give information on the line-of-sight component of the magnetic field. A way to investigate these geometries is through radiative transfer simulations which attempt to simulate the observed variability of the polarized emission and spectral changes during flaring. Our past studies of individual sources included explorations of parameter space, incorporating propagating shocks in a jet plasma with a turbulent magnetic field and a weak ordered axial component to account for the behavior during quiescent phase [1,2,3], and a preliminary investigation of the effect on the emission of inclusion of a large-scale helical magnetic field component which persists to parsec scales [4]. A similar approach incorporating a standing shock, turbulence, helical magnetic fields, and relativistic effects has been undertaken independently in [5]; but this work, using the TEMZ code [6] is focused on the upstream region of the jet and aims to explain the very rapid polarization variability apparent in the optical regime, and the features in 43 GHz Very Long Baseline Array (VLBA) imaging data. There are many theoretical studies investigating the affect of ordered magnetic fields on the observed emission, e.g., [7], but in general these studies do not include turbulence which we believe is required by the well-documented stochastic character of blazar variability in the optical and radio bands. Our radiative transfer simulations of emission from a relativistic jet flow, therefore, include both turbulent and ordered magnetic field components, one of which may be helical. We investigate the effect on the observed emission of varying the ratio of the contribution to the energy density by the magnetic field components with these geometries. The work presented here is directed at understanding the radio-band polarization arising downstream of the recollimation shock plausibly associated with the 43 GHz very long baseline interferometry (VLBI) core. Ultimately we aim to understand the topology of the magnetic field across the blazar zone.

2. Data: Observed Properties of Radio Band Polarization from UMRAO Monitoring Data

OT 081 was selected for detailed analysis because of its well-separated, distinct outbursts in both total flux density and linear polarization, and because it has remained highly active for decades, allowing us to examine variability on decadal timescales. While they are more extreme in this source, the variability patterns exhibited in both total flux density and polarization at the centimeter band are representative of those seen in many other UMRAO program sources [8]. The source is a low redshift (z = 0.322) low-frequency-peaked BL Lac object, which is unusual for this class because it has been detected in the TeV spectral band by MAGIC [9]. The source is in the sample of Fermi-monitored bright blazars showing flaring during the mission (https://fermi.gsfc.nasa.gov).

In Figure 1 we show centimeter-band, long-term total flux density and linear polarization light curves based on 30-day averages of the daily observations obtained with the University of Michigan 26-m equatorially-mounted paraboloid (UMRAO) for the example blazar OT 081. A description of the general observing and calibration procedures are given in [10]. The instrumental polarization was calibrated using observations of HII regions which are assumed to be unpolarized during each observing run. The electric vector position angles (EVPAs) were calibrated using a source of polarized emission located at the vertex, and these EVPAs are accurate to 0.5$°$ [10]. As discussed in [11], there are 180$°$ uncertainties in the EVPA determinations, and the EVPA data in this figure were plotted using a criterion which minimizes the jump between successive EVPA measurements. While Figure 1 illustrates the long-term spectral behavior of the flux and linear polarization, the averaging smooths the short-term variations. In Figure 2 and Figure 3 we show daily averages for segments of the data, and here the rapid variations in the linear polarization are apparent.

The temporal behavior of the total flux density (see bottom panel of Figure 1) exhibits several large-amplitude (S ≥ 6 Jy) outbursts in total flux with temporally-associated increases in polarized flux which are ideal for our modeling purposes. While the emission is often blended in both single dish monitoring and very long baseline interferometry (VLBI) images, emission from contributing synchrotron-emitting components can sometimes be deblended using linear polarization observations [12], where resolved structure in the linear polarization light curves is apparent. The fractional polarization remains below 10%, which is significantly lower than the theoretically expected value of 70–75% for an optically-thin, non-relativistic emission region with a uniform magnetic field [13]. The EVPA light curve in the top panel shows a mix of ordered and chaotic variations characteristic of many blazars in the UMRAO program [8]. In [4] we presented preliminary modeling results for three epochs (1985, 2008, and 2010) and noted that a significant change in the flow occurred between 1985 and 2008. Here we describe the behavior in detail for these epochs, and we model the variability during an important additional epoch commencing in 1996 in order to bridge this time gap and to investigate the magnetic field and flow properties during a time window where the data exhibit unusually-strong fractional polarization and markedly different EVPA characteristics from the previously-modeled time windows.

3. Description of the Model and Analysis Method

We adopt the model that is described in detail in [1]. This model assumes that the flow is populated by radiating particles with a power law distribution in Lorentz factor with density given by $n(\gamma$)d$\gamma$ = n${}_o{\gamma }^{-\delta }$d$\gamma$ where $\gamma >{\gamma }_i$ and ${\gamma }_i$ is the Lorentz factor low energy cutoff with index $\delta$ fixed, that the radiating plasma is kinetically-dominated, and that it contains a tangled or turbulent magnetic field which is dynamically unimportant. The rationale for the assumption of a tangled or turbulent magnetic field is summarized in [14] and is based largely on evidence from the UMRAO linear polarization monitoring results that values of only a few percent are typically found in the quiescent phase, and that the fractional linear polarization is rarely higher than 10% during flaring. The turbulent magnetic field is assumed to be static; a single realization was adopted in the simulations, to simplify the computation (see [1]), but the effect on the simulation of using an evolving set of cells was tested by examining the simulated light curves at different azimuthal angles for which the line of sight intersected a different set of cells, and the effect was shown to be unimportant. In addition to this turbulent magnetic field, the model allows for inclusion of ordered magnetic field components. A weak ordered axial magnetic field component, containing 2% of the energy density in the random magnetic field, was initially introduced in order to reproduce the stable EVPAs observed during quiescent time periods, but additional modeling of selected flares subsequently revealed that increasing this axial component improved the agreement between the simulated light curves and the UMRAO data. It was thus included in subsequent simulations as a free parameter.

A general scenario for the production of outbursts in blazar jets based on a shock-in-jet model has been adopted by many researchers since the 1980s when it was proposed to explain the millimeter-to infrared variability observed during a major outburst in 3C 273 [15] and the multi-frequency, centimeter-band flux density and polarization variability observed in BL Lac as part of the UMRAO monitoring program [16,17], and we make the underlying assumption that the outbursts modeled are produced by the passage of one or more propagating shocks. The discrete emission features subsequently detected with very long baseline interferometry (VLBI) have also been associated with shocks; e.g., [18]. In our adopted scenario, a propagating shock compresses the initially-tangled magnetic field and any ordered components associated with the underlying flow. The signature of the shock event in the polarization data is a swing in the EVPA through tens of degrees, typically on a timescale of months, and an increase in the fractional linear polarization. For a transverse shock, this swing is through 90$°$, while for oblique shocks it is through a smaller range. The shocks in the model are allowed to be oriented at a wide range of directions to the flow, but we have selectively picked epochs for modeling which exhibit swings through large ranges of EVPA, indicative of transverse shocks, as these are easiest to identify in the data. The orientation of the shock relative to the flow direction obtained from VLBI observations is specified by two angles [19]. These are the shock obliquity $\eta$, the angle relative to the direction of the upstream flow, and $\psi$, which specifies the azimuthal direction of the shock normal. However, the choice of the value for the azimuthal direction has been shown to have little effect on the simulated light curves [1]. In principle we should allow for the retarded time effect that results from a temporal change in physical parameters as a ray propagates along any line of sight. In the models presented in this paper the change occurs on a time scale quite long compared with the light travel time across the source; indeed, we examined the effect of including retarded time effects in [3] and demonstrated that it produces only a marginal change in the shape of the outburst profile. We do not include it in the analysis presented here because of the required additional computation time.

Before simulations of the light curves can be carried out, the flux contributions to each outburst from the individual, blended flares must be separated in order to determine how many shocks are present: each flare within the outburst is attributed to a single shock. At millimeter wavelengths, separation of the individual flares has been carried out as described by [20]. They adopted a generic profile comprised of an exponential rise, a sharp peak, and an exponential decay approximately 1.3 times longer than the rise, superimposed on a constant baseline. However, that method, based on assuming a generic shape for all flares, does not reproduce the shape of the observed centimeter-band outbursts, where blending is typically more severe. To resolve the individual flares, we have used a combination of the structure apparent in the total flux density and linearly polarization light curves and our knowledge of the flare profile shape obtained from simulations for a single shock [1]. This procedure has the advantage of being both empirical in using the observed structure in the light curves associated with the emission contributions from individual flares within an outburst, and physically motivated by our knowledge of the expected theoretical flare shape for a single shock. The number of shocks is established following this deconvolution (using the minimum number of shocks which can reproduce the data and not including a shock to match every bump in the light curve) and remains fixed throughout the iterative procedure used to establish the model parameters. Calculations were performed for the transfer of polarized emission through a diffuse plasma allowing for emission, absorption, Faraday rotation, mode conversion, and relativistic aberration and boosting [1]. We adopt arbitrary units for lengths, particle densities, and magnetic field strengths. Matching to an observed flux for a source of known distance provides a constraint on an algebraic combination of these quantities, but we do not explore that here, as it is the temporal and spectral form of the light curves that constrain the model. The flow passing through each shock is compressed (including the turbulent pre-shocked magnetic field and any ordered components), thereby increasing the degree of order of the magnetic field and the emissivity. Simplifying assumptions in the modeling are that each shock occupies the full cross section of the flow and propagates rectilinearly with a constant velocity; that shocks contributing to each time window modeled (which may include more than one outburst) have the same orientation; and that shocks do not intersect, so that each flow segment is compressed only once.

Procedure for Determination of the Flow Properties

The fundamental properties defining the jet flow and shock system and the primary observational constraint for each property are given in

Table 2. Free parameters and UMRAO data used as a model constraint.

Type	Parameter	Symbol	Primary Observational Constraint
Internal	low energy cutoff	${\gamma }_i$	EVPA spectral behavior
Internal	axial B field	$B_z$	EVPA and P%
Bulk flow	Lorentz factor of the jet flow	${\gamma }_f$	P%
Bulk flow	viewing angle to the jet	$\theta$	P%
Shock property	shock obliquity	$\eta$	$\Delta$EVPA
Shock Property	shock compression	$\kappa$	$\Delta$S and P%
Shock property	shock sense	F or R	Doppler factor and ${\beta }_{app}$
Shock property	shock length	l	Duration of flare in S
Shock Property	shock onset	$t_o$	Start of flare in S and P

. The parameters describing the jet are the internal state of the quiescent flow, the bulk dynamics, and the orientation of the flow. The shock attributes are specified by the shock obliquity ($\eta$), a measure of the shock strength given by the compression factor $\kappa$ (defined so that unit length is compressed to length $\kappa$); the shock sense (forward or reverse: moving faster than the underlying flow or being overtaken by it); the length of the shocked flow (the extent of the involvement of downstream flow expressed as a percentage of the flow length); and the time of the shock onset chosen to match the onset of flares in the UMRAO data. The onset corresponds to the time at which the shock enters the flow.

Each modeled outburst is described by six key parameters: three for the quiescent flow (Lorentz factor of the flow, viewing angle, and the low energy spectral cutoff), and three describing each shock occurring during the outburst (sense, strength given by the compression factor, and obliquity). As noted above, lengths, densities, and field strengths are in arbitrary units; and as the magnetic field is dynamically insignificant, the field to particle energy density ratio is unconstrained. The axial magnetic field is used as a fine-tuning parameter, while the onset of the shock and the length of the shocked region are set by the start and duration of new variability in the light curves. The value of the optically-thin spectral index $\alpha$, where $\alpha$ is given by $S\left(\nu \right)d\nu =S_o{\nu }^{-\alpha }d\nu$, is set to 0.25 based on the rather flat total flux density spectra found in general in the UMRAO data; this parameter, however, does not have a significant role in the modeling, as the spectral behavior is largely determined by opacity effects below 14.5 GHz. While higher spectral indices have been found for optically-thin jet components in the outer regions of the jet, for example, by the MOJAVE survey [21], which often reveals considerably extended structure well beyond the 15 GHz core, here we are modeling the emission arising in the inner, partially-opaque region of the source which dominates in the single-dish data, corresponding to the core on MOJAVE maps. The fiducial Lorentz factor of the energy spectrum is set to 1000 at 8 GHz. Given a magnetic field strength, Doppler factor, and redshift, the Lorentz factor of particles radiating at the frequency of observation is uniquely determined. However that is not possible here, as the absolute value of the magnetic field strength is not set in the modeling. As this parameter does not play a significant role in the modeling (its purpose is to establish a reference value against which to judge how far the spectrum extends for a given low energy cutoff), we have arbitrarily chosen a plausible value. The cutoff Lorentz factor (${\gamma }_i$) is set to a value of ∼50, to ensure that there are negligible internal Faraday effects unless the data suggest that such effects should be included. The individual shock onset times and lengths are set by the structure under the outburst envelope, which is well-defined by the cadence of the UMRAO measurements. For each flare there are six constraints: the peak total flux density, the percentage polarization, and the EVPA, together with the spectral form of each of these observed properties. Additionally, we require that the same quiescent flow parameters apply to each flare during an outburst, despite their different shock strengths, providing an even tighter constraint on parameters, and ensuring that the flow angle of view in particular is very well-determined.

An initial quiescent flow Lorentz factor and shock sense (always “forward” for the epochs modeled in this paper based on the observed proper motion from complementary VLBA data) is selected, and a typical shock compression is adopted in setting up the initial state of the model. Additionally, an initial shock obliquity is chosen by inspection of the range of EVPA change displayed by the data, using the results based on simulations for a single shock [1] to estimate the obliquity needed to produce a match to the data.

Starting with an arbitrary initial viewing angle, the onset times, lengths, and strengths (compression) of the shocks are then adjusted iteratively, in an attempt to fit the outburst shape and the spectral evolution in total flux density (specifically, the amplitude of the change in the total flux density $\Delta S$, the spectral shape when the emission is most opaque, and the amplitude and position of structure within the light curve) and the observed fractional linear polarization (quantitatively, the peak value during the outburst). If a satisfactory match to the data cannot be achieved, the viewing angle is adjusted, and the process repeated. The quiescent flow Lorentz factor is adjusted to improve the match with the data if no viewing angle is found, which yields a good match to the data. For a given quiescent flow Lorentz factor, the model fractional linear polarization is very sensitive to the viewing angle, especially when $\theta$ is only a few degrees. A change in the viewing angle can be used to refine the simulated fractional linear polarization while leaving the total flux light curves nearly unchanged.

Structure in the observed fractional linear polarization light curves and trends in the EVPA are then examined by eye (in general, the observed EVPA light curves display very complex behavior), and the shock obliquity and the low-energy spectral cutoff are adjusted to reproduce these features. The flow viewing angle is determined primarily by matching the observed fractional linear polarization but uses the observed EVPA and the range of its change as secondary constraints. A decrease in the low energy cutoff (introducing significant internal Faraday effects) can modify the model fractional linear polarization, requiring further iteration of the shock parameters, viewing angle, and the quiescent flow Lorentz factor. The shock speed, ${\beta }_s$ is computed from the assumed quiescent flow speed (${\beta }_f$) and the upstream speed in the shock frame which comes from the compression factor of the shock and the shock obliquity. Details of this procedure are given in [1] following [19].

A major goal is to investigate the impact of including both ordered and turbulent magnetic field components in the simulation to identify the polarization signatures associated with these magnetic field geometries. We examine this question using a series of simulations with a varying ratio of ordered to turbulent magnetic fields to determine for which ratios our model is able to produce the maximum amplitude at all three wavelengths in both total flux and linear polarization and maintain the spectral character of the total and polarized flux as the outburst evolves. We aim to replicate the general character of the polarization variability as a test of the viability of our proposed scenario and not to reproduce the details of the light curves. The polarization is very sensitive to the viewing angle, particularly in sources with jets aligned near to the line of sight, as will be shown for OT 081, and as part of the new work we examine the impact of changing the viewing angle on the polarized emission during the time window T2010 (see

Table 3. Summary of observed emission properties during modeled epochs.

Epoch	Date Range	Length (yrs)	S${}_{\mathit{max}}$(14.5)	P%${}_{\mathit{max}}$(14.5)	S${}_{\mathit{max}}$(4.8)	P%${}_{\mathit{max}}$(4.8)	$\Delta$EVPA (14.5-4.8)
T1985	1985.0–1986.3	1.2	7.4	5.5	3.5	4.5	variable
T1996	1996.8–2002.2	5.4	6.2	13	3.8	9	∼20$°$ when constant
T2008	2008.4–2010.7	2.25	7.2	4	4.3	2-3	generally small
T2010	2010.8–2011.9	1.09	6.5	8	4.2	5	generally <10$°$

).

4. Results

We modeled the spectral variability during four epochs which contain large-amplitude outbursts in both total flux density and polarization; these are labeled in Figure 1, and blow-ups of these time segments are shown in Figure 2 and Figure 3. In the top panel of the lower plots in Figure 2 and Figure 3 showing the EVPA data for each epoch, we have conservatively only plotted polarization measurements with ${\sigma }_{EVPA}\le$14$°$; this criterion corresponds to a >2$\sigma$ detection of the linear polarization [10]. The outbursts included in these time windows were selected because swings in the observed EVPA are through tens of degrees; the temporal and spectral behavior is consistent with the expected ordering of the magnetic field during a compression by a propagating shock.

In Table 3 we summarize the observed emission properties during each modeled epoch that we wish to reproduce in the simulated light curves. Table 3 lists an epoch designation (column 1), the date range of the time window (column 2), the duration of the time window (column 3), the peak total flux density at 14.5 GHz (column 4), the maximum amplitude of the fractional linear polarization (column 5), the maximum total flux density and fractional polarization at 4.8 GHz (columns 6 and 7 respectively), and the typical separation of the 14.5 and 4.8 GHz EVPAs (column 8). The tabulated values are representative values of maintained peaks and not the highest value measured; some single measurements have large error bars, especially in the linear polarization at 4.8 GHz. These error bars are 1$\sigma$ error estimates based on the standard errors of the Stokes parameters and the calibration uncertainties. There are also some deviant total flux density measurements at 4.8 GHz (unexpectedly high or low relative to the neighboring measurements), as can be seen in Figure 2 and Figure 3. Measurements at this frequency were often scheduled during degraded weather conditions, and the best weather was reserved for the 14.5 GHz observations, which were more sensitive to tropospheric conditions.

While the peak flux amplitudes are nearly identical during all of the strong outbursts modeled, the fractional linear polarization is significantly higher during T1996, and there is a persistent, nearly constant, wavelength-dependent separation of the EVPAs consistent with the presence of Faraday rotation in this time window. This is discussed in Section 5.1.

Table 4. Simulated emission properties during modeled epochs.

Epoch	S${}_{\mathit{max}}$(14.5)	P${}_{\mathit{max}}$(14.5)	S${}_{\mathit{max}}$(4.8)	P${}_{\mathit{max}}$(4.8)	$\Delta$EVPA (14.5-4.8)
T1985	7	6	3.5	4.5	generally flat
T1996	6.5	15	4.2	12.5	generally flat
T2008	7	3.5	4.5	2	generally flat
T2010	7	8	3.8	6	generally flat

lists the simulated values of the constraining properties for comparison with the observed values in Table 3. Recall that the simulated total flux density at 14.5 GHz is scaled to match the observed value. Although we list EVPA information in Table 3, the EVPA behavior as a function of time and observing frequency is a mix of chaotic variations interspersed with ordered changes, while the simulated EVPA curves are predominantly flat with little spectral spread; in contrast, the polarized flux in both the data and simulations exhibits well-defined outbursts. The simulated EVPA exhibits a large swing early in each flare, as the leading edge of the first shock competes with the quiescent flow to dominate the emission. This is particularly apparent at the start of the simulation for each time window. Thereafter, the flow is largely filled by a progression of shocks, and variations in the EVPA occur only if and when there is a sufficient time between shocks for the quiescent flow’s contribution to become significant again.

Table 1. Derived flow parameters.

Parameter	T1985	T1996	T2008	T2010
Cutoff Lorentz Factor (${\gamma }_i$)	50	10	50	50
Bulk Lorentz Factor (${\gamma }_f$)	5.0	5.0	10.0	10.0
Viewing angle (${\theta }_{obs}$)	1.7$°$	2.3$°$	1.4$°$	1.1$°$
Number of Shocks	6	6	6	3
Shock Lorentz Factor	14.5	17.6	24.7	27.3
Shock ${\beta }_{app}$ (in units of c)	10.4	16.6	21.8	22.5
Axial B as % of Turbulent B (energy density)	25%	100%	64%	56%

lists the jet parameters identified from matching the general properties of the variability in each epoch and the simulations using the method described. In all four epochs, the fiducial Lorentz Factor (${\gamma }_c$) was assumed to be 1000, the shock sense was forward, and the shock obliquity was taken to be 90$°$ (transverse shocks). The model allows for oblique shocks, but the ranges of change in EVPA are consistent with transverse shocks during these epochs.

Table A1. Shock attributes for T1985.

Shock ID	Start Time	Length	Compression
T1985-1	1985.0	12.0	0.35
T1985-2	1985.22	11.4	0.40
T1985-3	1985.4	11.4	0.30
T1985-4	1985.59	7.44	0.28
T1985-5	1985.77	12.0	0.4
T1985-6	1986.03	8.04	0.35

through

Table A4. Shock attributes for T2010.

Shock ID	Start Time	Length	Compression
T2010-1	2010.83	15.0	0.38
T2010-2	2011.07	15.0	0.40
T2010-3	2011.23	7.5	0.30

list the attributes of the shocks included in the simulations.

In selecting a model, there is a potential degeneracy between the flow speed and the angle of view, but additional features, such as modulation of the total flux density light curve, allow this degeneracy to be lifted in most cases. The Lorentz factor of the underlying flow appears to be rather coarsely determined, but that is because a flow with ${\gamma }_f=10$ is boosted by only $1.26$ from a flow with ${\gamma }_f=5$ (with a similar boost between Lorentz factors of 10 and 20). As seen in §3.4 of [3], the viewing angle is well-determined when the degeneracy between flow speed and angle can be lifted, a change of a few tenths of a degree being sufficient to radically alter the degree and spectrum of the polarized emission for blazars with jets seen close to the flow axis such as OT 081.

While we discuss the degeneracy in the model for T1996 in Section 4.2, the same degeneracy does not occur in the other epochs modeled. In T1985 a lower Lorentz factor of the flow was required to achieve the small modulation of the total flux density where the flares are substantially merged with very little evidence of substructure compared with event T2010, which has a similar value of fractional polarization. The need to achieve comparable levels of P% with different values of the Lorentz factor then fixed the viewing angle. A high value of the Lorentz factor was needed to produce the considerable structure in S during T2008 with the overall (excluding two spikes at 14.5 GHz) very low polarization setting the viewing angle.

Our simulations are able to reproduce the maximum amplitude and spectral behavior of the fractional linear polarization, but they fail to reproduce the complex behavior of the EVPA light curves. In the simulations, these are characterized by a flat spectrum and monotonic changes at each frequency with time. In the data, however, they are highly variable, exhibiting both chaotic and ordered behavior. VLBA images show that in the centimeter band, most of the polarization comes from a single component (occasionally with a contribution from the innermost jet components). The core polarization most likely is the sum of unresolved emission contributions from hot spots with very different EVPA orientations. It is plausible that the observed chaotic behavior, not well-reproduced by the simulations, occurs because of the way these contributions to the Q and U Stokes parameters combine in the upstream end of the centimeter-band emission site.

To illustrate the sensitivity of the simulated light curves to small changes in the input parameters, we consider the effect on the light curves of changing the adopted viewing angle by only a few tenths of a degree (from 1.1$°$ to 1.4$°$) for T2010. We show in Figure 4 the dramatic increase in the amplitude of the fractional polarization in the OT 081 simulations resulting from this small viewing angle change. To asses the effect that this change would have on derived parameters, we determined the effect of changing $\theta$ by this same range (±0.3$°$) on the derived shock speed ${\beta }_{app}$. The latter is computed from the shock speed and $\theta$. For each of the four epochs modeled, the resulting range of values of ${\beta }_{app}$ is shown in parentheses following the values in Table 1: 10.4 (9.1, 11.6); 16.6 (15.7, 17.1); 21.8 (19.1, 23.5); and 22.5 (18.2, 25.2). From this exercise we have estimated that the value of $\theta$ determined from the modeling is accurate to ±0.3$°$ and that the apparent shock speeds are accurate to ≈10%.

Within the accuracy of our parameter determinations, the modeling identifies time-dependent changes in key parameters, including the ratios of the ordered-to-turbulent magnetic field contribution to the energy density, the viewing angle, and the Lorentz factor. These are consistent with the presence of multiple emission regions within a broad jet viewed at different angles by the observer in subsequent time windows.

Because changes in the observer’s viewing angle (geometric effects) are believed to be important in interpreting the variability behavior in some blazars, e.g., [22,23,24], and a viewing angle change is a potential cause of changes in the polarization properties, we examined the published sky-projected innermost jet position angle as a function of time determined from MOJAVE data at 15 GHz for OT 081. The source is one of a dozen in the MOJAVE study of 200 AGN jets in [25] showing an oscillatory trend in plots of inner jet position angle (PA) versus time that have been fit with sine curves. The time period of the well-defined EVPA separation and high values of fractional polarization (1999–2002) in the UMRAO light curves occur during a minimum in the MOJAVE sine curve fitted to the inner jet PA as a function of time for OT 081. We do not identify a clear relationship between the variability in flux and polarization in the UMRAO light curves and variability of the inner jet position angle during T1996. Possibly, VLBA monitoring observations at millimeter wavelengths would be helpful in identifying a relation between jet orientation and the centimeter-band emission, if present in this source.

Comments on the variability during each time window follow.

4.1. T1985

Figure 2 left compares the observations and simulations for T1985. We show time units along the abscissa for ease of comparison, while in [4] these were expressed in arbitrary units. The panels showing fractional linear polarization and total flux density have matched scales for the data and model plots, while the EVPA panel for the simulation has a smaller range. There are six shocks during this time window; the last one was mistakenly omitted in the data figure and table in [4] but included in the simulation. Note that the error bars on the linear polarization data at 8 GHz are often relatively high, that several observations of P% were rejected at this frequency, and that some measurements during this time window were obtained with a receiver measuring total flux density only, such that the 8 GHz EVPA light curves were constructed from data that were less well-sampled than at our other two frequencies. In the data plot we have allowed the range of EVPA values to extend between 45$°$ and 245$°$ and added 180$°$ to measured values below this minimum value; this reduces the number of apparent jumps through 180$°$ in the EVPA light curve. With the cadence of the UMRAO measurements, the ambiguity is generally clearly resolved.

The simulation reproduces the spectral character and maximum amplitude of the total flux density shown in the bottom panel, the detailed behavior of the fractional linear polarization including a small flare which occurred near to the start of the time window, the maximum value of the fractional linear polarization attained at all three frequencies, and the nearly constant value of the EVPAs during portions of the time window modeled. The details of the EVPA light curve, however, are very complex, and they are not well-reproduced by the model. Because this epoch precedes the operation of the VLBA, there are no ancillary imaging data which might be helpful in understanding the physical origin of the EVPA discrepancy in this epoch.

4.2. T1996

The data and simulated light curves are shown in Figure 2 right for T1996. In the data plot, the EVPAs have been rescaled to minimize the jump produced by the 180$°$ ambiguity in the determinations; 180$°$ has been added to EVPAs below 25$°$. This time segment includes several large outbursts in total flux density and spans the longest time window modeled. During these events, the 14.5 and 8 GHz total flux densities often nearly track, while the 4.8 GHz flux is significantly lower, indicative of self-absorption in the emitting region. The most striking features in the observed light curves are the frequency-dependent separation in the EVPAs, which persists from approximately 1999.2 through 2001.6, and is well-defined with the sampling of the UMRAO three-frequency data—the temporally-associated high amplitude of P% which is near 13% at maximum and persistently near 10%, and the substructure in P% at 14.5 GHz. During 1997 to early 1999, the total flux density shows strong opacity effects. The EVPA separation is primarily apparent during the time segment of the total flux density light curve when the spectrum is relatively flat, but it persists as the spectrum begins to steepen circa 2000.8. The simulation includes six shocks over this 5.4-year time period. There is considerable substructure in the light curve, but, as described, we have taken the approach of introducing the minimum number of shocks which can reproduce the primary features in the data, and we have not included a shock to match every sub flare; e.g., the substructure at 1998.7.

While in all other flares modeled the flow parameters were tightly constrained and well-determined, in this time period two acceptable fits were found because of a degeneracy between the Lorentz factor of the quiescent flow and the viewing angle, which resulted in very similar simulated light curves. We initially modeled T1996 with a Lorentz factor of 10 and a viewing angle of 1.1$°$, and all other parameters as shown in Table 1, for this epoch. We subsequently obtained a marginally better simulation of the data features adopting a Lorentz factor of 5 and a viewing angle of 2.3$°$, and these values are listed in Table 1 as the adopted model. While the simulated total flux density and fractional linear polarization light curves using the two models were both able to reproduce the general character of the light curves, the simulated EVPA light curves showed subtle differences, including a marginally-larger separation of the EVPAs with frequency and a better match to the detailed fine structure using the adopted model. Furthermore, the adopted model yielded a shock ${\beta }_{app}$ of 16.6c which is more in agreement with the maximum VLBA component speeds for the same epoch (e.g., ${\beta }_{max}$ is in the range 5-21c from 22 GHz data in [26] and ${\beta }_{max}$ = 6.85 ± 0.72c based on MOJAVE data during this time window [27]). The original model gave a substantially higher shock speed of 32.3c, which is inconsistent with these VLBA determinations, although it is possible that fast components were missed with the imaging observing cadence.

It is not surprising that extreme jet and shock properties are obtained during epoch T1996 in view of the very high fractional linear polarization with an amplitude nearly twice that attained at other epochs and the nearly-constant wavelength-dependent EVPA separation which persists until 2001.6 during the sustained high-amplitude in fractional linear polarization. Our model does not include an external screen, and the assumption that the Faraday rotation is produced internally implies a shift of the cutoff Lorentz factor of the radiating particle distribution to lower energies (10 compared to the value of 50 at the three other epochs). While depolarization is expected for internal Faraday rotation [28], this is offset by the combination of a stronger mean magnetic field during this time window and the series of high-compression (strong) shocks. There is, unfortunately, relatively little VLBA polarization monitoring data available during this time window which might help in understanding the structure and amplitude of the 14.5 GHz polarization and the cause of the very distinct EVPA separation which begins rather abruptly in 1999. The polarization movie on the MOJAVE website shows structural changes occurring after an unfortunate data gap between 1997.1 and 1999.7.

4.3. T2008 and T2010

As shown in Figure 3 left, in T2008 we reproduced the amplitude and spectral variability of the fractional polarization with six shocks. In contrast, in T2010 (Figure 3 right) while we simulated the spectral character of the total flux density and the relatively flat spectrum in EVPA with two long, relatively-strong shocks and one shorter-but-stronger shock (with onsets determined on the basis of the appearance of new flaring in the total flux density light curve), the simulation does not reproduce all the salient features apparent in the data plot. Specifically, it does not reproduce the detailed structure in the fractional linear polarization light curve or the monotonic increase in the EVPA, although the flat EVPA spectrum is reproduced by the model after mid-2011. This discrepancy between the simulation and data suggests that a more sophisticated model is required or that more than one region of the jet (the core and one or more jet components) are contributing significantly to the source-integrated polarized emission. Unfortunately there are only two MOJAVE imaging measurements from T2010, at epochs 2010.98 and 2011.28, and these do not overlap with the period of intense variability in fractional linear polarization (see http://www.cv.nrao.edu/2cmVLBA/data/1749+096), so the possibility of multiple contributing components cannot be verified. The image at 2011.28 is dominated by a weakly-polarized core component.

4.4. Consistency of Model Parameters with Values from Other Methods

As a global check of our overall methodology, we compared our preferred jet flow parameters with those obtained by [29] using VLBA monitoring data at 43 GHz obtained from April 2009 to January 2013 which overlaps with T2008 and T2010. That work finds a Lorentz factor of 11.0 ± 3.6, and a viewing angle of 2.4 ± 1.0$°$ which are in agreement with our values. Source parameters using broadband radio-band observations of flux density variations over the time window 2007 to 2015 to obtain variability Doppler factors [30] identify a Lorentz factor of 7.8 and a viewing angle of 2.3$°$ but a smaller value of ${\beta }_{max}$ = 4.36c than we identify with the modeling; this method assumes that the parameters are not time dependent, and they refer to the full time window of the data.

5. Faraday Rotation or the Signature of a Helical B Field? Wavelength-Dependent EVPA Separation

Here we examine the evidence for Faraday rotation in the single dish and VLBA data and evaluate the effect of a large scale helical magnetic field on the EVPAs using simulations. In principle, evidence for the presence of a toroidal magnetic field is provided by a gradient of the Faraday rotation measure across the jet determined from simultaneous, multi-frequency VLBI images with sufficient resolution to obtain adequate signal-to-noise in the measurement [31]. As discussed in [32], a transverse gradient of the rotation measure indicates the presence of a toroidal component only and does not necessarily mean that the magnetic field is helical since the rotation measure gradient only establishes the toroidal part of the magnetic field; the helical field is a combination of both a toroidal field and a poloidal field, both of which are vector-ordered. The presence of a sign change in the rotation measure gradient, however, does support the presence of a helical magnetic field, but the number of sources with strong evidence for a helical magnetic field based on the identification of this sign change is small [28]. A toroidal magnetic field in OT 081 was identified in [33] using data obtained in 22 March 2004, and this result provided additional motivation for selecting OT 081 for our detailed analysis.

5.1. Faraday Rotation

In our discussion of T1996 we noted a separation in the 14.5–4.8 GHz EVPA light curves that is consistent with Faraday rotation. Faraday rotation is a propagation effect generally attributed to an external screen. Two VLBA studies determine Faraday rotation measures in this jet. These are: (1) a rotation measure study [34] which uses seven frequency VLBA data points between 8.1 and 15.2 GHz obtained on June 20th 2001 (a date within T1996), finding values of 145 ± 24 and 97 ± 25 rad/m${}^2$ for the core; (2) a study of MOJAVE data [28] which finds a median core rotation measure of 136 rad/m${}^2$ based on data obtained on 2006-06-15. Both of these results use VLBA observations obtained between 8 and 15 GHz, and they do not include data at 4.8 GHz where a ${\lambda }^2$ EVPA spread in the light curves would be larger than that based on data in the smaller frequency range used. We verified using UMRAO data that the rotation measure based on our observations at 14.5, 8.0, and 4.8 GHz (2.07, 3.75, and 6.25 cm) in 2001.47 is consistent with the rotation measure determination by [34]. During the time window studied in [28] the polarization was low; hence there is low signal-to-noise in the UMRAO polarization measurements, and additionally, simultaneous (within a week) data were not obtained at the three frequencies, so we could not reliably compare the result based on MOJAVE data with an independent estimation based on UMRAO data which extends to 4.8 GHz.

The VLBA-derived rotation measure values are of the same magnitude during these two epochs, but inspection of the images shows that the EVPAs are oriented nearly perpendicular to the jet flow in the 2001.47 image during T1996 (implying that the magnetic field is expected to be nearly parallel to the flow), while in 2004 the EVPAs corrected for Faraday rotation are parallel to the flow direction with B perpendicular, as expected in the case of transverse shocks. As can be seen in Appendix

Table A2. Shock attributes for T1996.

Shock ID	Start Time	Length	Compression
T1996-1	1996.8	2.5	0.25
T1996-2	1997.48	1.67	0.3
T1996-3	1997.88	1.25	0.2
T1996-4	1998.25	1.25	0.22
T1996-5	2000.31	5.0	0.3
T1996-6	2001.53	3.75	0.275

and the light curves, the 2001.47 image falls near a minimum in the light curve and the start of a strong flare, so the VLBA data during this epoch sample the underlying jet magnetic field.

5.2. Effect of Inclusion of a Helical Magnetic Field on the Simulated Light Curves

In this section we consider the effect on the simulations of adding an ordered helical magnetic field component to the model for T1985. We examine the interplay between the three magnetic field components which we propose are present: the underlying turbulent magnetic field, an ordered axial component associated with the evolution of the flow, and a large-scale helical component, plausibly residual from the formation and collimation of the jet. We incorporate the helical magnetic field component into the model in the form of a magnetic flux rope which permeates the jet plasma. As discussed in [35] this solution for the helical field is parameter-free and depends only on the field intensity: the field strength decreases, and the orientation of the magnetic field becomes more azimuthal as the radius increases.

T1985 is a time window which does not show a prominent and persistent EVPA separation in the UMRAO light curves such as we find in T1996 (see Figure 2 for this comparison). We note that the data in T1985 are more poorly sampled compared to those obtained in T1996, so that we can not exclude the possibility of some Faraday effects.

In Figure 5 we compare simulations for T1985 computed using order multiples, as defined in [1] for a weak helical magnetic field component. The effect of increasing the order multiple from 0.05 (a nearly negligible contribution from a helical magnetic field) to 0.5 in panel b is to suppress the amplitude of the fractional linear polarization and to produce a wavelength-dependent spread in the simulated EVPAs. With the inclusion of this additional magnetic field component, the simulated light curves more closely match the observed amplitude of the linear polarization and the EVPA spectral behavior.

In Figure 6 we show the effect of increasing the helical contribution to be comparable to the contribution to the energy density from the turbulent B field (panel a) and we increase this to be ten times greater in panel b. In these two simulations there is a pronounced spread in EVPA, but it is distinct from the ${\lambda }^2$ separation generally assumed for Faraday rotation, and this EVPA spread between the light curves at the three frequencies is time variable. Note that the structure in the fractional polarization light curves has been dramatically suppressed while the total flux density light curves have been unchanged by the addition of the helical magnetic field. The spectral evolution of the EVPAs for the strong helical magnetic field component does not mimic the spectral patterns in the polarization in any UMRAO sources, including the EVPA pattern seen in T1996. Based on comparison of our simulations with the UMRAO data, we conclude that the presence of a strong helical magnetic field component is not consistent with the data for OT 081. However, inclusion of a modest helical B field component, as shown in Figure 5, in fact, improves the agreement between the simulation and the data. Within the coarse grid investigated here, the best agreement between data and simulation is obtained with an order multiple of 0.5. This result supports a scenario in which a modest large-scale helical B field persists into the parsec-scale region of the jet.

Zavala and Taylor [36] suggest that Faraday rotation could be produced by a helical magnetic field which wraps around the jet. A possible alternative explanation for the unusual character of the linear polarization during T1996 is a increase in the strength of a helical B field component. As discussed in the previous section, the effect of including a helical B field on the light curves is to smooth the structure in the fractional linear polarization and to increase the wavelength-dependent separation of the EVPAs. However, we examined this hypothesis using simulations not shown here, allowing for a helical magnetic field, and we could not account for the spectral variability apparent in the data.

6. Discussion and Future Work

We have demonstrated that our model is able to reproduce the spectral behavior and amplitude range of the centimeter band polarized flux observed in OT 081 which exhibits the variability commonly found in the UMRAO data. We have shown that these simulations are very sensitive to small changes in viewing angle, and illustrated the unique signature in the polarization simulations produced by the inclusion of an ordered helical magnetic field in the radiative transfer simulations, and we have identified the effect of time-dependent changes in the jet parameters on the emission.

We are working to add refinements to the model to improve the agreement between the data and the simulated light curves. As discussed, the simulations do not reproduce the complex behavior in the EVPA light curves. They also do not reproduce the swings in EVPA through more than 90$°$, which we found during some flares in the UMRAO database. In our adopted scenario, we have assumed that shocks are not interacting and that they are moving at a constant velocity in rectilinear motion, while both assumptions are clearly in contradiction to VLBI observations. [37], for example, has identified acceleration in several blazar jet components, including one in OT 081. Furthermore, by assuming that the plasma is not multiply-shocked, we have neglected the impact of the passage of each successive shock on the magnetic field topology in the emitting region, and hence its affect on the evolution of the magnetic field structure.

UMRAO, the facility which provided multifrequency polarization monitoring observations at centimeter band for decades, was closed in 2012, and the data required for the analysis described here depends on well-sampled polarization observations obtained simultaneously at several observing frequencies. No current program provides the spectral coverage required for study of subsequent outbursts. The monthly cadence of the BU Blazar program at 43 GHz [29] and of POLAMI [38] at 1.3 and 3.5 mm are too infrequent relative to the variability time scale to follow in detail the polarization changes in many sources. MOJAVE multifrequency polarimetry observations are currently being carried out as part of a two year program, but this may be too short a time window to sample the range of variability, and the monthly observing cadence may miss structure in the polarization variability which is needed for our modeling procedure. However, a new dedicated millimeter polarization monitoring program using the 14-m radio telescope of the Metsähovi Radio Observatory to obtain polarization observations of a sample of about 90 sources at 22, 43, and 86 GHz is currently being developed which does have the capability to provide the data which are required for quantitative tests of viable models and for investigating the magnetic field in the centimeter-to-millimeter interface. Such measurements have the potential to bridge the parsec to sub-parsec spatial regime, leading to a better understanding of the magnetic field geometry in blazar jets, and we look forward to future results based on the expected millimeter polarization data and further refinement of the method described here.