CONSTRAINING THE PHYSICAL CONDITIONS IN THE JETS OF γ-RAY FLARING BLAZARS USING CENTIMETER-BAND POLARIMETRY AND RADIATIVE TRANSFER SIMULATIONS. I. DATA AND MODELS FOR 0420−014, OJ 287, AND 1156+295

The UMRAO total flux density and linear polarization observations presented were obtained with the University of Michigan 26 m equatorially mounted paraboloid as part of the University of Michigan AGN monitoring program (Aller et al. 1985a). At 14.5 GHz (2 cm), the primary frequency for the measurements since the launch of Fermi, the polarimeter consisted of dual, rotating, linearly polarized feed horns, symmetrically placed about the paraboloid's prime focus; these fed a broadband, uncooled HEMPT amplifier with a bandwidth of 1.68 GHz. At 8.0 GHz (4 cm) an uncooled, dual feed-horn beam-switching polarimeter system was employed using an on-on observing technique; the bandwidth at this operating frequency was 0.79 GHz. At 4.8 GHz (6 cm) a single feed-horn system was employed with a central frequency of 4.80 GHz and a bandwidth of 0.68 GHz. A change of observing frequency required a horn change, and hence observations were generally carried out for one- to two-day time blocks at each frequency in a rotation. These observations were made within ±2.5 hr of prime meridian passage in order to reduce the effect of position-dependent gain variations on the measurements. The adopted flux density scale is based on Baars et al. (1977) and uses Cas A as the primary standard. Observations of nearby secondary flux density calibrators were interleaved with the observations of the target sources approximately every 1.5–3 hr to verify the stability of the antenna gain and to confirm the accuracy of the telescope pointing. The electric vector position angles (EVPAs) were calibrated using a source of polarized emission mounted at the vertex of the paraboloid; this assembly was surveyed at installation to an accuracy of 012. To verify the calibration and stability of the instrumental polarization, selected Galactic H ii regions were observed several times each day; these objects were assumed to be unpolarized. Each daily averaged UMRAO observation was comprised of a series of 8–16 individual measurements obtained over a 25–45 minute time period. Because of the frequent observations of the calibration sources, the number of targets that could be monitored in a 24 hr run was limited to 20–24 sources.

The initial source sample of 33 sources for the project was selected on the basis of known or expected variability in both the γ-ray and the centimeter-band. In the radio-band the multi-decade, centimeter-band data from the UMRAO monitoring program provided a guide to the expected timescales of the variability, the peak amplitudes reached during flares, and the presence of ordered changes in the linear polarization data. Selecting sources with large-amplitude, resolved variations was crucial to the project goals since time-variability in the flux and linear polarization provide the constraints in the modeling. Sources with highly significant EGRET detections were preferentially chosen initially, but new flaring sources were added based on the Fermi Large Angle Telescope (LAT) measurements. In addition, the sample was restricted to sources in the MOJAVE (Monitoring of Jets in active galactic nuclei with VLBA Experiments⁶), 15.4 GHz imaging program (a frequency which is near to the UMRAO primary frequency and probes the same region of the jet), and priority was given to sources in the Boston University monitoring program at 43 GHz,⁷ a frequency that provides complementary information on the structure and structural changes in the inner jet region. As the project progressed, the sample was reduced to 18–20 blazars in order to observe the current most-variable objects with increased cadence.

The γ-ray photon flux light curves for the three sources studied were produced in the 0.1–200 GeV band using Fermi-LAT data obtained during 2008 August–2012 December. Weekly binned photon fluxes were obtained using ScienceTools–v9r27p1 and P7SOURCE_V6 event selection. The LAT data were extracted within a 10° region of interest centered upon the position of the target. These used an unbinned likelihood analysis (tool gtlike) to determine the photon fluxes by including in the model all of the sources within 15° of the target and by freezing the spectral index of all sources to the values in the 2FGL catalog.

The past variability histories of the three sources based on the UMRAO monitoring data and their relation to EGRET detections are briefly reviewed below. These earlier measurements are useful for placing the variations observed during the operation of Fermi in context, for estimating the duty cycle in the radio band, and for relating the onset of Fermi flares to the radio-band data by identifying opacity effects. Synchrotron self-absorption can produce time delays between the peaks of the flares in the γ-ray and the centimeter bands ranging from zero to several months (Pushkarev et al. 2010).

The long-term centimeter-band total flux density and linear polarization data obtained by the UMRAO monitoring program are shown in the left plot of Figure 1 for the LSP QSO 0420−014 (PKS 0420−01; 2FGL J0423.2−0120) at redshift z = 0.9161. A precessing binary black hole model has been proposed for this object based on an analysis of nine epochs of archival very long baseline interferometry (VLBI) data at 3.6 cm during the time period 1989–1992 (Britzen et al. 2000); this work includes the UMRAO total flux density variability data from 1982 to 1994. Note that nearly continuous, high-amplitude variability has persisted over decades in both total flux density and in fractional linear polarization. A peak amplitude of S ≈ 14 Jy at 14.5 GHz was attained in late 2003 making this source one of the brightest blazars historically in the centimeter band. Temporally resolved outbursts in fractional linear polarization occurred throughout the monitoring period, reaching a peak value of 5%. This low fractional linear polarization compared with the value of ∼75% expected for a canonical synchrotron source is characteristic of the sources in the UMRAO program and has been used to argue for the presence of a turbulent magnetic field during a quiescent, pre-flare stage. While a source-integrated rotation measure (hereafter RM) of −13 rad m⁻² (Rudnick & Jones 1983) determined from Very Large Array (VLA) observations has been included in the linear polarization data plotted, this small value produces a negligible shift in the observed EVPA from the intrinsic value. To test whether inclusion of Very Long Baseline Array (VLBA) determined values would reduce the spread in the EVPAs in this source and in the other two sources, we generated light curves incorporating a range of rotation measure values taken from the literature, e.g., Hovatta et al. (2012). The EVPA-frequency offset was smallest for 0420−014 and OJ 287 using the VLA-determined values from Rudnick & Jones (1983) and when assuming no correction in the case of 1156+295 for which no VLA values are available from either Rudnick & Jones (1983) or Rusk (1988). The exercise confirmed that significant RM corrections due to external Faraday rotation are not required for the data set presented here.

Zoom In Zoom Out Reset image size

Complementary, long-term spatially resolved VLBI imaging data from both the MOJAVE 15.4 GHz and the 43 GHz Boston University program are available publicly. The MOJAVE data for this source commence in mid-1995, with a typical cadence of a few observations per year; these imaging data track the motion of four components from 1995 through 2011; based on these data a maximum apparent component speed of β_max = 5.74c has been obtained (Lister et al. 2013). A parsec-scale jet orientation of 183° from MOJAVE measurements is presented in Kharb et al. (2010) that aligns with a large-scale VLA position angle of 183°. However, at least at 43 GHz, the projected parsec-scale jet orientation varies with time. The data from the Boston University program identify a significant change in direction from 280° in 1998–2001 to 100° degrees from 2008 to the present time (Troitskiy et al. 2013). These 43 GHz data that resolve the inner jet region identify that eight knots have been ejected since early 2009. The derived values of β_app, the component speed, are considerably faster than those identified from the MOJAVE images, but these components lie within the unresolved core region at 15 GHz. However, the 15 GHz data are able to follow the motions of the individual components to greater distances downstream in the jet.

In the γ-ray band the source was detected by EGRET in the 1990s, with high significance (test statistic TS = 46) during the viewing period centered on late 1992 February (Hartman et al. 1999; 1992.16) which coincides with a plateau in a centimeter-band flare. A maximum flux (E > 100 MeV) of 0.45 × 10⁻⁶ cm⁻² s⁻¹ is reported in von Montigny et al. (1995). Shortly after the launch of Fermi the source flared in the radio-band. The event modeled as part of our program is shown in the right plot of Figure 1. For bins with TS ⩽ 10, a 2σ upper limit is shown. Several low-level γ-ray flares were detected prior to the main flare which started near JD 2455170 (2009.93). Panels 2–4 show daily averages of the UMRAO total flux density and linear polarization monitoring observations at 14.5, 8.0, and 4.8 GHz. Substructure is clearly apparent in the fractional linear polarization light curve indicative of emission from individual flares which are blended in the source-integrated, single-dish, total flux density observations. The centimeter-band linear polarization during the γ-ray flaring exhibits monotonic changes in the EVPAs through tens of degrees, and associated increases in the fractional LP; these are consistent with the signature of the passage of a shock through the emitting region.

Two-week averages of the UMRAO total flux density and linear polarization data are shown in Figure 2 (left) for the ISP BL Lac object OJ 287 (0851+202; 2FGL J 0854.8+2005) at z = 0.306. OJ 287 has been identified as a potential binary black hole based on historical optical data extending back in time to the 1890s; both a 12 yr period attributed to the perturbation of the accretion disk of the primary by the secondary, and a 60 yr period identified with a precession cycle have been proposed, e.g., Valtonen & Wiik (2012). Cyclic variations in the fractional linear polarization in the optical band with a time-scale of 2 days are discussed in D'Arcangelo et al. (2009). A wavelet analysis of more than two decades of the UMRAO total flux density and linear polarization data (1971–1998) identified two persistent signals—a longer timescale period of ≈1.66 yr associated with the quiescent jet, and a shorter periodicity of ≈1.2 yr associated with the series of large flares that occurred during the 1980s (Hughes et al. 1998). The fractional linear polarization exhibits well-defined nulls between flares. A new cycle of intense activity commenced circa 2000 following a relatively quiescent period. While the maximum amplitude of fractional linear polarization is only 5% at 14.5 GHz for 0420−014, similar to values typically found for UMRAO program sources, for OJ 287 we detected unusually high values of fractional linear polarization that are in the range P ≈ 15%–18% at 14.5 GHz; these high levels were reached during relatively quiescent total flux density levels and are not associated with large flares. One of these peaks occurred in fall 1993 during the operation of EGRET; at this time the total flux density was only 2 Jy.

Zoom In Zoom Out Reset image size

The MOJAVE data at 15.4 GHz identify five primary moving components and a maximum apparent component speed, β_max, of 15.14c based on seven moving features⁸; from an analysis of the variation of the projected innermost jet position angle on the sky with time, a monotonic swing in jet position angle of a few degrees per year from 1995 to the end of 2010 is discussed in Lister et al. (2013). Campaign VLBA data at 43 GHz obtained in 2005 October–November and 2006 March–April have been modeled with a fast spine and a slow sheath (D'Arcangelo et al. 2009) with a maximum Lorentz factor (LF) of 16.5 for the spine, and an LF of only 5 for the sheath. The presence of non-periodic jet wobble has recently been suggested based on longer-term 43 GHz VLBA data (Agudo et al. 2012). Differences in the behavior at 43 and 15.4 GHz are attributed to the different angular resolutions of the data sets (Lister et al. 2013). As discussed in a later section of the paper, there is also significant opacity in the source in the centimeter-band emission region which may contribute to these differences.

In the γ-ray band, this source was detected by EGRET with 3σ detections (TS = 9) during Compton Gamma Ray Observatory pointings centered on late 1992 September and in 1994 mid-November (Hartman et al. 1999), but no strong EGRET detections were reported during the epoch of the high fractional linear polarization in 1993 apparent in the UMRAO data. During the operation of EGRET the radio-band total flux density was on the decline. In contrast, the source has been unusually bright in the past few years, reaching an amplitude of S = 10 Jy, which is just below the highest amplitude measured historically by the UMRAO program data. Only one strong γ-ray flaring period (photon flux >5 × 10⁻⁷ photons s⁻¹ cm⁻²) was detected by Fermi during the time window shown in Figure 2 (right).

The LSP, FSRQ 1156+295 (TON 599; 4C +29.45; 2FGL 1159.5+2914) at z = 0.725 has exhibited nearly continuous, centimeter-band flaring over the >30 yr time period shown in Figure 3 (left). Until 1998.0 the amplitude of the total flux density variations was relatively modest. Subsequently the variability changed in character, and several well-separated flares in total flux density occurred. During these events the total flux density is characterized by a flat spectrum, and the events exhibit an unusual shape characterized by nearly equal rise and fall times, reminiscent of strong, centimeter-band events in 1510-089 (Marscher et al. 2003). The time segment since the launch of Fermi includes two of these large and distinctive flares. Overall, the fractional linear polarization has exhibited a long-term base level near 2% and reached a maximum fractional linear polarization of about 10%. The EVPAs show no preferred orientation, and they exhibit a large range of values in the light curves that show complex structure.

Zoom In Zoom Out Reset image size

A maximum apparent component speed of β_max = 24.59c has been found from MOJAVE data, making the flow in this source the fastest of the three sources studied here, assuming that this pattern speed is an indicator of the bulk flow speed. Bending of the relativistic jet on VLBI and VLA scales away from the observer was invoked to explain earlier X-ray variability (McHardy et al. 1990).

In the γ-ray band, EGRET obtained a highly significant detection with TS ⩾ 36 during early 1993 January, and an additional detection with TS > 36 with lower photon flux was obtained in late April–early 1995 May (Hartman et al. 1999). The Fermi time window shown in the lower panel of the right plot in Figure 3 captures two similar radio-band flares. The Fermi data train commences just after the start of the first outburst, and there is no indication of an enhanced flux density level in the early measurements indicative of a prior flare. The first radio-band flare does not have a clearly associated γ-ray counterpart (an example of an orphan flare), while the second, slightly flatter-spectrum event does. While changes in the flow may have occurred between the times of these radio-band events upstream of the centimeter-band region, no change in opacity within the centimeter-band emission region is indicated by the UMRAO data. Further, only modest changes in the projected inner jet position angle between 2008 and 2010 are identified by the MOJAVE data (Lister et al. 2013), and hence geometrical differences do not appear to account for the occurrence of the radio-band flare with no γ-ray counterpart.

The formulation and assumptions adopted in the radiative transfer models incorporating a propagating shock scenario and used in the analysis presented here are described in detail in Hughes et al. (2011). We summarize some of the salient features here.

The models assume that a power law distribution of radiating particles permeates the quiescent flow and that the density distribution of these particles follows a power law of the form n(γ)dγ = n₀γ^−δdγ where γ > γ_i. The low-energy cutoff of this distribution is specified in terms of a thermal⁹ Lorentz factor (LF), γ_i, and a fiducial value, γ_c, indicates from where in the particle spectrum the observed emission arises. The power law index, δ, (δ = 2α_thin+1; α_thin is the optically thin spectral index) is fixed, and the density constant, n₀, is assumed to fall off due to the adiabatic expansion of the flow.

The models assume the presence of a passive, turbulent magnetic field in the synchrotron-emitting region of the jet before the passage of the shock. The magnetic field structure in this region is represented by cells with randomly ordered field directions. The passage of a shock through this region produces a compression that increases the particle density and the magnetic field energy density (and hence the emissivity). The degree of order of the magnetic field is increased by the compression associated with the passage of the shock, and this leads to an increase in the fractional linear polarization.

The formulation of the models allows for shocks to be oriented at any angle to the flow direction (Hughes et al. 2011) and builds on earlier work restricted to the special case of transverse shocks (Hughes et al. 1985, 1991). Each shock's orientation is specified by two angles; these are its obliquity, η, measured with respect to the upstream flow direction, and the azimuthal direction of the shock normal, ψ. Earlier work showed that the simulations are relatively insensitive to changes in the azimuthal direction, however. The shocked flow is characterized by a compression factor, (κ), a length (l) defined as the evolved extent of the shocked flow, and the shock sense (forward or reverse: F or R). The outburst seen in the variability data is associated with the propagation of the region bounded by the limits of the shocked flow. To simplify the computations, each shock is assumed to span the entire cross section of the flow, and it is assumed to propagate at a constant speed. Additionally, multiple shocks contributing to a single outburst are assumed to have the same orientation, and subsequent shocks contributing to an outburst do not shock pre-shocked plasma. Inclusion of retarded-time effects for test cases revealed that only small differences in the simulated flare shapes and spectral behavior result (Aller et al. 2013a), and hence this effect is neglected in the simulations presented here.

To reproduce a well-defined EVPA in the quiescent state a small fraction of the magnetic energy (2%) is initially assumed to be in an ordered component of the magnetic field; support for this assumption comes from the low levels of fractional linear polarization identified during relatively quiescent states in the UMRAO data (Hughes et al. 2011). An important factor in the modeling is the observer's viewing angle. Initial estimates for this parameter for each source were taken from VLBI studies (Table 1) limiting our need to explore parameter space. Table 1 additionally includes the maximum apparent component speed, β_max, determined from 15.4 GHz MOJAVE data, (Lister et al. 2013) in Column 2; the variability Doppler factor, δ_var, in Column 3, derived from fits to archival Metsähovi single-dish monitoring data at 22 and 37 GHz (Hovatta et al. 2009); a flow LF, Γ, and viewing angle, θ_var, in Columns 4 and 5, respectively, computed from these data following the prescription in Hovatta et al. (2009) and using the revised apparent speed determined from longer-term MOJAVE data; and, in columns 6 through 8 the Doppler factor, LF, and viewing angle determined by Jorstad et al. (2005) from 43 GHz VLBI data which probe the flow upstream of the centimeter-band emission region. These results provide useful comparisons with our derived parameters.

Table 1. Jet Flow Parameters from Centimeter and Millimeter Band Variability

Name	βmax	δvar	Γ	$\theta _{{\rm var}}^{\circ }$	δ43	Γ43	$\theta _{43}^{\circ }$
0420−014	5.74c	19.9	10.8	1.5	16.2 ± 1.2	11.0 ± 0.5	3.0 ± 0.6
OJ 287	15.14c	17.0	15.3	3.4	13.6 ± 2.4	16.3 ± 3.8	4.0 ± 0.6
1156+295	24.58c	28.5	24.9	2.0	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅

Notes. β_max from the MOJAVE Website; variability Doppler factor and viewing angle as described in the text; 43 GHz Doppler factor, Lorentz factor, and viewing angle from Jorstad et al. (2005).

Download table as:  ASCIITypeset image

A complication in analyzing single-dish, centimeter-band light curves is the blending of emission contributions from individual, evolving flares. In general, the UMRAO light curves trace out an envelope over activity that may include contributions from several flares/shocks. Past work on separating blended flares, particularly in the analysis of mm-wave data, has often assumed a generic flare shape with an exponential decay (Valtaoja et al. 1999; León-Tavares et al. 2011). This exponential shape, however, is not a good match to the centimeter-band events, and, further, the method requires knowledge of the baseline flux level for removal of this emission contribution prior to the fitting process. Instead, we have resolved the blended events into individual flares by using a combination of the structure apparent in the total flux density and linear polarization light curves at our highest frequency, 14.5 GHz, where they are best resolved, and the expected flare shape for a single shock based on a library of simulated light curves generated as part of this work. Example simulated light curves for a single shock are shown in Hughes et al. (2011).

The number of shocks contributing to a modeled event, established by this deconvolution procedure, is fixed during subsequent modeling. A quiescent flow LF and shock sense (always "forward" for the sources discussed in this paper for the reasons discussed in Section 5) are chosen based on the expected viewing angle and the observed proper motions as discussed in Section 3.1. The component proper motions predicted by the model depend on the quiescent flow LF and shock strength, so a typical shock compression is adopted in setting up the initial state of the model. An initial shock obliquity is chosen by inspection of the range of EVPA change displayed by the data, using the results of Hughes et al. (2011) to estimate the obliquity needed to produce a match to the data.

The fiducial LF (γ_c) is arbitrarily set to 1000, and the cutoff LF (γ_i) is set to a high value (typically 50 or 100), to ensure negligible internal Faraday effects. The optically thin (frequency) spectral index, α, is fixed at a value of 0.25 based on the rather flat total flux density spectra exhibited by most UMRAO sources even when there is no evidence that opacity is important.

The onset time, length, and strength (compression) of the shocks are then adjusted iteratively, in an attempt to fit the event shape in total flux (quantitatively, the flux increase at peak, the amplitude of the change in the total flux density Δ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 event). If a satisfactory match to the data cannot be achieved, the viewing angle is adjusted, and the process repeated. The quiescent flow LF is adjusted similarly if no viewing angle is found that yields a good match to the data. For a given quiescent flow LF the model fractional linear polarization is very sensitive to the viewing angle, so that a change in the viewing angle can be used to refine the model fractional linear polarization while leaving the total flux light curves largely unchanged; importantly, the value of the viewing angle is very well-constrained by the modeling.

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. A change in the 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 LF.

We emphasize that the intent of the fitting procedure is to ascertain whether the model can reproduce the general features of the data. We do not provide proof of uniqueness of the models; nor can we quantify a goodness of fit, as the final phase of modeling attempts to refine the model parameters using trends in the EVPA variability established by inspecting the data.

The UMRAO data for the event modeled are shown in Figure 4 (left) which spans the two year period 2009.0–2011.0. (circa JD 2454832–2455562). These data exhibit the characteristic behavior expected for shocks propagating through the emitting region: an outburst in total flux density, time-associated outbursts in polarized flux, and ordered changes in the EVPA light curve. A characteristic behavior in the total flux density light curves is blending of the contributions from individual flares; these individual flares are better-resolved in the linearly polarization light curves as is illustrated by the structure in this fractional polarization light curve. There is a 180° ambiguity in each EVPA determination that can affect the interpretation of the temporal changes exhibited by the EVPA data. In general, the choice of whether to add or subtract a multiple of 180° when generating UMRAO EVPA light curves is based on an algorithm which makes this decision by requiring the smallest point-to-point jump in the time sequence of the EVPAs using the data at all three frequencies, and additionally requiring that differences between consecutive data points be less than 90°. Further, the range of EVPA values is restricted to 180°. However, for ease in comparison of the data with the simulated light curves for this source, we relaxed these restrictions on the EVPA range and added multiples of 180° to reproduce the character of the simulated light curve shown in Figure 4 (right). Following Aller et al. (1985b), we have only included data in this plot, and in those for the other two sources studied, for which the differential errors in EVPA derived from the standard errors of the Stokes parameters are less than 14° (P/σ_P >2). Several observations were rejected at 4.8 GHz because they did not meet this acceptance criterion: the observations at this frequency were often made during degraded weather conditions, while the best weather was generally reserved for observations at 14.5 GHz. As a consequence, the observed 4.8 GHz data is more poorly sampled than at the two higher frequencies. The EVPAs during the outburst show monotonic swings in EVPA through ∼90° as expected for the passage of a transverse shock through the emission region. The best-sampled case is the systematic change in EVPA from late 2009 through early 2010. These changes occur during the rise portion or plateau of the individual flares within the outburst envelope, as expected for the shock scenario. While evidence for large, systematic rotations of the EVPA have been found for a few sources based on the UMRAO data (Ledden & Aller 1979; Aller et al. 1981), these rotations have occurred when the total flux density had reached a maximum during an outburst or was declining and are characteristically different from the linear polarization changes produced by the passage of a shock. Shocks produce an increased degree of order of the initially turbulent magnetic field and an increase in the fractional linear polarization apparent in the light curves as a polarization outburst. Polarization rotations occur when the fractional linear polarization is relatively small (<2%; they can be produced via a random walk resulting from the evolution of the turbulent magnetic field within the emission region (Jones et al. 1985).

Zoom In Zoom Out Reset image size

Based on the structure in the total flux density and linear polarization light curves, combined with the burst profile shape expected for a single shock, we introduced three shocks oriented transversely to the flow direction into the simulations. The onset times of the shocks included in the simulation are indicated by upward arrows along the abscissa of the lower panel, while the γ-ray peak photon flux is marked by a downward arrow at the top of this panel. The parameters specifying each shock start time, the length as a percentage of the flow, the compression factor, and the location of the peak flux in dimensionless time units are summarized in Table 2. The first shock is weak, the second shock which starts at 2009.6 is the main shock, and the third shock which commences around 2010.0 is the shock that coincides temporally with the brightest γ-ray flare. Figure 4 (right) shows the simulated light curve. The amplitude of the fluxes presented here and in the other simulations shown are scaled to match the peak total flux density in the observed light curve. Twenty time steps were used in all of the simulations presented in order to resolve the structure in the light curves during the event modeled.

Comparison of the data with the simulations reveals that the simulations are able to reproduce the global features of the light curves; these include the range in the total flux density at the highest frequency (14.5 GHz), the evolution in the spectral index of the total flux density during the outburst based on the amplitudes of the total flux density at 14.5 and 4.8 GHz, the amplitude of the change (maximum to minimum) of the fractional linear polarization, and the change in the spectrum of the EVPAs during the evolution of the outburst; the latter is complex. The flow properties adopted in the simulation shown are: an optically thin spectral index α_thin = 0.25 (fixed for all three sources), a fiducial "thermal" LF (associated with the random motions of the emitting particles) γ_c = 1000, a low-energy cutoff "thermal" LF of γ_i = 50, a bulk LF of the flow γ_f = 5.0, and an observer's viewing angle of θ_obs = 4°. From the modeling, the derived LF of the shocked flow is 8. At a viewing angle of 4° the shock transition LF (which is the emission feature that one would see moving in the VLBA images) implies β_app ∼11c. While the overall spectral behavior is very well reproduced by the simulation, note that there are some differences, including the range of the EVPA swing. The extent of the swing is best judged from the behavior while the trailing shock (shock 3) dominates. In the data the range of the swing is through approximate 135°. A comparable EVPA swing occurred during the first shock; however, with the sampling of the data, the time at which the shock enters the flow and becomes manifest is not as well defined in the light curve. Additionally, the null in the observed fractional linear polarization light curve and an associated sharp variation in the EVPA light curve near 2010.6 are not reproduced by the simulation. These differences between the data and model may result from the adopted simplifying assumptions in the model formulation.

This source is included in both the MOJAVE and the 43 GHz BU VLBA programs, and it is instructive to compare the shock start times with those of the propagating jet components apparent in these VLBI images. A detailed analysis of the 43 GHz imaging data during the γ-ray flaring is presented in Troitskiy et al. (2013) based on monthly sampled data. A comparison between onset times and apparent speeds of the shocks and components is presented in Table 3. Recall that the shock onset times correspond to the time when the leading edge of the shock enters the flow; this is not the same parameter as the time at which the separation of the component from the radio-band "core" occurs. At 43 GHz, five components have been ejected during the time period 2009.0 through 2011.5, and four of these occur during the time window included in the simulation. Within the error bars, the injection times for components K2-K4 are consistent with the onset of our three shocks. The fourth component is very fast, and there is no corresponding shock required from our analysis. Overall, however, this comparison reveals an excellent correspondence between the moving emission features and the model shocks. Note that the model β_app and the observed component speed are nearly equal for shocks 1 and 3.

The radio-band event that was selected for modeling occurred during the 10-month time period from 2009 September to 2010 July and commenced just after the source exited from the summer 2009 Sun gap (a time at which it was located too close to the Sun for observation by UMRAO). The observed light curve is shown in the left plot in Figure 5. A large, self-absorbed total flux density flare took place during this time window which temporally coincides with a large γ-ray flare as shown in Figure 2 (right). A large flare in fractional linear polarization is also apparent in the UMRAO monitoring data which reached a maximum fractional linear polarization near 8% during the rise portion of the outburst. The EVPAs show systematic changes, most apparent at 14.5 GHz; these changes are through a smaller range of EVPA values than observed for 0420−014, and the pattern is consistent with the presence of oblique rather than transverse shocks. The simulation shown in Figure 5 (right) incorporates three forward-moving shocks; the parameters specifying these shocks are tabulated in Table 4. Comparison of the model light curves and the data shows that temporally resolved swings in EVPA are associated with the entry of the shocks into the flow. As in the case of 0420−014, these EVPA swings occurred during the rise or plateau phases of the blended, total flux density light curve. Note that the EVPA swings at 14.5 GHz in the simulation and in the data are through a comparable range of about 30°. The range of the swing at 4.8 GHz, however, is through a larger range in both the data and in the simulation. There is an indication of a swing through a similar range at 8 GHz for the observed changes associated with shocks 1 and 2, but unfortunately, the data are not sufficiently well-sampled at either 4.8 or 8.0 GHz to capture the full range of the observed swing. Additionally, the presence of larger swings at 4.8 and 8.0 GHz is masked by the lower signal-to-noise of the data. The larger swings in the simulated light curves at the two lower frequencies through approximately 90° are the result of Faraday effects; at 14.5 GHz these effects are reduced, and the EVPA swing reflects the underlying oblique shock. While the modeling identifies Faraday effects within the source during the time period modeled, note that the separation of the EVPAs in the data plot does not follow a simple λ² relation. We attribute the observed spectral evolution of the EVPAs to the complex effects of shocks within the flow, Faraday effects, and opacity. In this source, there is evidence for substantial self-absorption in the total flux density light curve, and this important feature of the total flux density is well-reproduced by the simulated light curve shown in the right plot. This model additionally reproduces the range of change in total flux density (2–9 Jy), the spectral evolution of the total flux density (inverted near outburst peak, flattening in the decline portion of the event), the approximate time delay of the peak flux at each frequency; the range of the change in fractional linear polarization (0%–6%); and a swing of the EVPAs through a smaller range of values than expected for a transverse shock. In the simulation shown we adopted a fiducial "thermal" LF (associated with the random motions of the emitting particles) γ_c = 1000, a low-energy cutoff LF γ_i = 20, a bulk LF of the flow γ_f = 5.0, and an observer's viewing angle of θ_obs = 15.

Zoom In Zoom Out Reset image size

As in the case of 0420−014, it is useful to compare these modeling results with the parameters identified from VLBI imaging studies. The long-term MOJAVE data identify a range of β_app from 0.438c (a nearly stationary, long-lived component) to 15.13c (Lister et al. 2013); component 11, ejected nearest to the launch of Fermi, has a speed of 6.58 ± 0.56c. The LFs of the modeled shocks, denoting the speed of the shock transition, are near 20c; these show relatively little spread since the compressions are similar, at least for shocks 2 and 3. The derived values of β_app from the modeling are 17c. We note that there is evidence for a change in the flow orientation with time (Lister et al. 2013) and that the modeled flow seen at angles between 0.5 and 25 would produce β_max values of 5c–20c consistent with the range of MOJAVE component speeds. This result suggests that geometric effects may be important in interpreting the data for this source. As noted earlier, there are differences between the flows at 43 and 15 GHz based on the VLBI monitoring data (Lister et al. 2013).

In Figure 6, we compare the data (left) and the simulation (right) for the outburst modeled in 1156+295. The shock parameters are given in Table 5. In this source, four forward shocks were required. The shock compressions range from 0.5 to 0.8, while the lengths of the shocked flows are identical. The simulation shown reproduces the relative amplitudes of the total flux density at the three frequencies and the spectral evolution of the total flux density, including the cross-over of the 14.5 GHz total flux density relative to the values at the other two frequencies during the outburst decline, the peak amplitude (6%) of the fractional linear polarization, and the monotonic increase in the linear polarization during the decline portion of the total flux density outburst. The simulated light curve shows two swings through approximately 90°. The first is associated with the entry of shock 1 into the flow. The data during this portion of the outburst are not sufficiently sampled to define the range of the swing using the 14.5 GHz alone. However, the combined 14.5 and 8 GHz data are consistent with a swing through a comparable range. The second swing terminates as the amplitude of the fractional linear polarization begins a new rise. This characteristic behavior is well-matched by the simulation. However, the swing is through approximately 120° exceeding the theoretical prediction of 90°. In addition to the proposed explanations for 0420−014, here shock interactions may contribute to the complexity of the flow (Aller et al. 2014).

Zoom In Zoom Out Reset image size

To characterize the flow, we adopted a fiducial "thermal" LF γ_c = 1000, the low-energy cutoff LF γ_i = 50, a bulk LF of the flow γ_f = 10.0, and an observer's viewing angle of θ_obs = 2°. The bulk LF of the flow is the highest of the three sources modeled, a result constant with the values listed in Table 1. This source is included in both the BU and the MOJAVE VLBA monitoring programs. Component speeds are presented for five jet components in Lister et al. (2013). These range from 1.65 ± 0.30c to 24.6 ± 1.9c. Note that the model β_app of 22c is in very good agreement with this VLBI result based on contemporaneous data.

In P. A. Hughes et al. (2014, in preparation) we present a detailed discussion of the effect of changes in the input parameters on the simulated light curves. However, to aid the reader in assessing the sensitivity of the simulations to the choice of these parameters and hence to judge the quality of the fits, we illustrate the effect of independently changing the fraction of energy in an ordered axial component of the magnetic field for all three sources, and the effect of increasing the cutoff LF for OJ 287 to match the value of 50 found for the other two sources.

During the analysis of the third source in the trio presented, 1156+295, we found that an increased contribution from an ordered axial magnetic field component was required to reproduce the features of the observed light curves. Motivated by this result, an increased ordered component of the magnetic field was included in the simulations for OJ 287 and 0420−014, and these results are presented in this paper. These new light curves are better able to reproduce the structure within the observed light curves and the amplitude of the linear polarization compared with our earlier results in Aller et al. (2013a). To illustrate the affect of this modification, we show in Figure 7 (left and right) and Figure 8 (left) light curves generated assuming the weakest axial magnetic field able to yield a well-defined EVPA in the quiescent state (2%). For 0420−014 comparison of the observed light curves with the simulations based on the adopted model (Figure 4) reveals that the fractional linear polarization is 50% too high at all three frequencies (nearly 6% at 14.5 GHz compared with the observed value of 4%), and that there is no significant structure in the fractional linear polarization contrary to the data and the adopted model; additionally the EVPA is constant from early in the rise of the total flux, again contrary to both the data and adopted model for this event. For OJ 287 as shown in Figures 8 (left) adopting the minimal mean field yields a peak fractional linear polarization that is similar in amplitude to the data and to the adopted model; however, the value remains high until late in the outburst, while in the data and adopted model it drops by almost a factor of two well before the peak of the total flux. Again, the alternate model light curve is smoother with less "structure" than is exhibited by the data. Further, subsequent to the initial rise in total flux density the EVPA has a smooth variation spanning only tens of degrees, not reproducing the rapid swings apparent in the observed light curve. In the case of 1156+295, the simulated fractional linear polarization exceeds 10% and the time variation in the EVPAs is too flat when compared with the observed features. Hence adopting a value for the axial magnetic field above the minimum level improves the ability of the simulations to reproduce the detailed structure apparent in the linear polarization light curves in this source.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For OJ 287, we show in Figure 8 (right) the simulated light curve produced by increasing the low-energy cutoff from 20 to 50 to match the values adopted for the other two sources, with all other parameters as tabulated for the adopted model. Here the peak fractional linear polarization is about 25% too high compared with the data and the adopted model, and it remains high until late in the outburst, while in the data and adopted model it drops by almost a factor of two well before the peak of the total flux. Again, the light curve is smoother, exhibiting less structure. Beyond the initial rise in total flux the EVPA exhibits a smooth variation spanning only tens of degrees, not reproducing the range in EVPA during the rapid swings.