© ESO, 2015

1. Introduction

The combination of very long baseline interferometry (VLBI) images and broadband flux density variability is a unique way to probe the emission mechanisms near the base of jets in the blazar class of active galactic nuclei (AGN). The broadband flares are often found to be connected to the ejection of new moving emission components into the jet (e.g., Marscher et al. 2008; Jorstad et al. 2013; Schinzel et al. 2012; Krichbaum et al. 2001, and references therein). Moreover, mm-VLBI observations offer a unique possibility to study the structural evolution in the parsec-scale region, which has been proposed as the site of much of the high-energy emission (e.g., Rani et al. 2013c, 2014; Fuhrmann et al. 2014; Marscher et al. 2008; Schinzel et al. 2012). These observations have therefore provided new constraints on the physical parameters of the emission regions: size, brightness temperature, magnetic field, and motion.

The blazar S5 0716+714 (z ~ 0.3, Nilsson et al. 2008; Danforth et al. 2013) is a BL Lac object with a featureless optical spectrum. It is one of the most intensively studied blazars because of its extreme variability properties across the entire electromagnetic spectrum (e.g., Villata et al. 2008; Fuhrmann et al. 2008; Rani et al. 2010b,a, 2013a,b; Larionov et al. 2013). The broadband flux variability of the source is quite complex with rapid flaring activity (on a timescale of a few hours to days) superimposed on top of a broad and slow variability trend on a timescale of ~1 yr (Rani et al. 2013a; Raiteri et al. 2003). VLBI studies of the source show a core-dominated jet pointing towards the north (Bach et al. 2005; Britzen et al. 2009), while Very Large Array observations show a halo-like jet misaligned by ~90° on kiloparsec scales. Britzen et al. (2009) suggest that there is an apparent stationarity of jet components relative to the core; more recent studies, however, have reported motion as fast as ~40 c (Rastorgueva et al. 2011; Larionov et al. 2013; Lister et al. 2013). Non-radial motion and wiggling component trajectories have often been observed in the inner mas jet region of the source (Britzen et al. 2009; Rastorgueva et al. 2011; Rani et al. 2014).

In Rani et al. (2013a; hereafter Paper I), we presented the densely sampled multifrequency observations of the source between April 2007 and January 2011. These observations allowed us to study the broadband flaring behavior of the source and to probe the physical processes, location, and size of the emission regions. The intense optical and γ-ray monitoring revealed fast repetitive variations (60−70 days) superimposed on a long-term variability trend on a time scale of ~350 days, which propagated down to radio wavelengths with an observed time lag of ~65 days. A detailed investigation of the optical flares found a variability amplitude proportional to the flux level, which can be explained by a variable Doppler factor. We also found that the shock-in-jet model for the evolution of radio flares requires geometrical variations in addition to intrinsic variations of the source. As a possible scenario to explain the observations, we suggested that the geometry significantly affects the long-term flux variations, which could be caused by a relativistic shock tracing a spiral path through the jet. To explain the multifrequency behavior of an optical-γ-ray outburst in 2011, Larionov et al. (2013) also suggest a shock wave propagating along a helical path in the blazar’s jet.

In this paper, we use high-resolution multifrequency VLBI observations to investigate the inner jet kinematics of S5 0716+714, with a focus on the major radio/optical/γ-ray flares over the time period of September 2008 to October 2010. We focus on the morphological evolution of the source to investigate its relation to the broadband flux variations reported in Paper I, in particular the high-energy emission. The paper is structured as follows. Section 2 provides a brief description of observations and data reduction. In Sect. 3, we report and discuss our results. Summary and conclusions are given in Sect. 4.

Fig. 1 43 GHz VLBA images (contours) and Gaussian components of S5 0716+714 at different epochs between September 2008 and October 2010. The date of observation is given at the bottom of each map (shifted Julian day, JD−2 454 000). The fitted model parameters for each epoch are given in Table A.1. The cross at the bottom of each map represents the FWHM of the restoring beam (0.25 × 0.18 mas). Contour levels correspond to −0.003 0.003 0.006 0.012 0.024 0.048 0.096 0.19 0.38 0.76 1.54 3.07 6.14 12.28 24.58 49.15 98.30 mJy/beam.

2. Multifrequency VLBI data: observations and data reduction

To explore the inner jet kinematics of the source, we used the 7 mm (43 GHz) and 3 mm (86 GHz) VLBI data obtained between September 2008 and October 2010, with observations at 26 epochs during this period. The 7 mm data are from the Boston University monthly monitoring program of bright γ-ray blazars with the Very Long Baseline Array (VLBA)¹. The data reduction was performed using standard tasks of the Astronomical Image Processing System (AIPS) and Difmap (Shepherd et al. 1994) software. The imaging of the source (including amplitude and phase self-calibration) was done using the CLEAN algorithm (Högbom 1974) and SELFCAL procedures in Difmap (Shepherd et al. 1994). Several iterations of phase corrections followed by amplitude adjustments were adopted for the self-calibration process. More details of the data reduction can be found in Jorstad et al. (2005).

The 3 mm observations were performed using the Global mm VLBI Array (GMVA). The GMVA is currently an array consisting of 14 antennas in Europe and the United States, including the eight VLBA stations equipped with 3 mm receivers, along with the Effelsberg, Onsala, Metsähovi, Pico Valeta, Plateau de Bure, and Yebes antennas. Between October 2008 and May 2010, observations were taken approximately every six months. Amplitude calibration is difficult at 3 mm, with atmospheric fluctuations often having a major effect. We find that amplitudes are frequently too low after the nominal calibration in AIPS. To correct for this, amplitudes are scaled by comparing the flux density on the shortest baselines with single-dish measurements. Details of the 3 mm data reduction can be found in Hodgson et al. (2014).

We modeled the observed brightness distribution (the visibility amplitude and phases of the observed radio brightness of the source) by multiple circular Gaussian components, thereby deriving positions, flux densities, and sizes of the distinct bright features in the jet. The model fits were carried out for all epochs within the DIFMAP package, starting with a point-like model and fixing the position of the brightest component to (0, 0). Jet components were added until adding an extra component did not lead to any significant improvement in the χ² value of the fit to the uv data. The uncertainties of the model component parameters were determined by comparing the parameter ranges obtained after performing model fits with a different number of model components. We used at least four different model fits for each epoch to obtain the parameter uncertainties. However, the uncertainties also depend on the self-calibration and data editing and on the brightness and size of the components. In addition, the uncertainties increase with increasing distance from the core. This is accounted for by following an independent approach for the error estimation proposed by Krichbaum et al. (1998). The fitted model parameters for all of the epochs are listed in Table A.1, where we only list the formal errors obtained from the first method. Figures 1 and 2 show the fitted circular Gaussian components superposed on the clean images.

To investigate the kinematics in the jet of S5 0716+714, the individual model components were identified following the assumption that the changes in the flux density, distance from the VLBI core, position angle, and size should be small for the time period between adjacent epochs. To prevent a potentially large systematic error arising from the incorrect cross-identification of moving features from epoch to epoch, the simplest scheme was adopted while identifying the jet-features. A self-consistent cross-identification is proposed using all available model-fit parameters. These cross-identifications are not necessarily unique, especially if the source evolves in a more complex manner than assumed here.

3. Results and discussion

In this section we present the flux and spatial evolution of the bright radio emission in the jet of S5 0716+714. Special attention is given to the kinematics of the inner jet (<2 mas from the core). We also investigate possible correlations between the jet kinematics and the broadband flux variations reported in Paper I.

Fig. 2 Inner jet structure of S5 0716+714 from GMVA observations at 3 mm. The cross at the bottom of each map represents the FWHM of the restoring beam (0.07 × 0.04 mas). Contour levels correspond to −0.004 0.004 0.008 0.016 0.032 0.064 0.12 0.25 0.51 1.02 2.04 4.09 8.19 16.38 32.77 65.54 mJy/beam.

3.1. Jet kinematics

Figures 1 and 2 display a sequence of VLBI maps convolved with a natural beam. The model fit parameters are given in Table A.1. In the following sections, we discuss the apparent motions of components and brightness temperature gradient in the parsec-scale jet in detail.

3.1.1. Component motion

For the kinematic study of individual components, we chose the VLBI core as a reference point, fixed to coordinates (0, 0). However, it is important to note that the absolute position of the reference point (i.e., core) may change, because of changes, for example, in opacity or instabilities in the jet, but the results presented here are not sensitive to that as we are interested in the relative motion of components with respect to the core. The VLBI data were fitted by using different numbers of components during different activity stages over the time span considered here. During this period, we identified a total of 12 components − C1, K1, A1, A2, A3, and NC1 to NC7 − in addition to the core, C0 (see Fig. 1). Figure 3 plots the evolution of the distance of different knots from the core, and Fig. 4 shows their trajectories in the XY-plane projected on the sky.

Stationary features

The sequence of VLBI images allowed us to investigate the component motion along the jet as a function of time. Most of the components exhibit significant motion down the jet, except A1, A2, and A3, which are comparatively stable in position, although they exhibit significant scatter in position angles that exceed the uncertainties (see Fig. 4). Nevertheless, their radial distance from the core remains at ~0.15 mas. Stationary features are a common characteristic in AGN jets (e.g., Fromm et al. 2013; Jorstad et al. 2001; Britzen et al. 2010). In straight jets, these can be produced by recollimation shocks, by instabilities, such as magnetohydrodynamic (MHD) and/or Kelvin-Helmholtz (KH), or by magnetic pinches (see Marscher 2009; Hardee 2006; Meier et al. 2001, for details). Bends in the jet can also cause quasi-stationary features, either because the jet turns more into the line of sight, thus increasing the Doppler beaming factor or because of the formation of a shock that deflects the flow (Alberdi et al. 1993). Numerical simulations indicate that when a moving knot passes through a standing recollimation shock, the components blend into a single feature, then split up after the collision with no lasting changes in the proper motion of the moving knot (Gomez et al. 1997; Fromm et al. 2012). During the interaction, the feature may move a short distance downstream before returning to its previous position. Therefore, it is quite possible that A1, A2, and A3 are the same standing feature, which becomes disturbed by moving features. We find that the kinematics of three new components, NC2, NC3, and NC6, are significantly different before and after the interaction. The components move more slowly (β_app ≤ 10c) as they approach the stationary feature, then faster (β_app ≥ 20c) after the interaction. This behavior is consistent with the stationary feature that represents a bend in the jet (Alberdi et al. 1993). Indeed, a bend at ~0.15 mas is evident in the 3 mm maps of the source (see Fig. 2). This suggests that the observed stationary features are standing oblique shocks that cause the flow of the jet to bend at ~0.15 mas owing to a transverse pressure gradient or an impact with an interstellar cloud.

Fig. 3 Evolution of core separation of the individual components as a function of time.

Fig. 4 Trajectories of the jet knots in the X–Y plane, where X = r cos(θ + 90) and Y = r sin(θ + 90). r is the radial separation from the core and θ is the position angle with respect to an imaginary north-south line drawn through the map center.

Moving features

The other components, C1, K1, and NC1 to NC7, exhibited significant motion in the radial direction. Component C1 represents very faint emission at a core separation of ≥1 mas. It seems to be the remnants of a component ejected in earlier epochs. The bright feature K1 could be a trailing component that forms behind a strong shock (Agudo et al. 2001) because it is identified at four epochs after the ejection of a new component (NC1).

As seen in Fig. 4, all components follow curved trajectories in the XY-plane projected on the sky. In addition, the trajectories differ from one component to the next. The wiggling trajectories in the jet might be a signature of helical motion. Helical jet models (e.g., Gomez et al. 1994; Hardee 2006) represent the jet as an inhomogeneous flow of plasma interacting with ambient matter. Magnetohydrodynamic or Kelvin-Helmholtz instabilities (Hardee 2006; Perucho et al. 2006) can explain the bending and helical structures and, possibly, also the less regular wiggles observed in parsec-scale radio jets.

To investigate the component motions in the jet of S5 0716+714, we used the XY positions of the components relative to the core to fit the trajectories with polynomials of different orders. The resulting fitted radial separation of each component, $\mathit{r}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sqrt{\mathit{x}\mathrm{(}\mathit{t}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathit{y}\mathrm{(}\mathit{t}{\mathrm{)}}^{\mathrm{2}}}$, is shown in Fig. 5. We find that a linear function is sufficient to fit the trajectories of components NC1, NC4, NC5, and NC7. However, a second-order polynomial is required to fit the trajectories of components NC2, NC3, and NC6, which corresponds to an apparent acceleration. Two multiple straight lines that are needed to fit different portions of the trajectories of these components in Fig. 5 clearly demonstrate the inadequacy of simple linear fits. The dashed curves represent the quadratic fits.

For each component, the fits yield an average proper motion and a mean speed. Using the estimated angular speed μ, we computed the kinematic parameters of the jet, such as the apparent speed, β_app, and Doppler factor, δ. We derive the apparent speed of the components, β_app, from the angular speed, μ, using ${\mathit{\beta }}_{\mathrm{app}}\mathrm{=}\frac{\mathit{\mu }{\mathit{d}}_{\mathrm{L}}}{\mathit{c}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\mathit{,}$(1)where d_L is the luminosity distance and z the redshift of the source. The luminosity distance d_L corresponding to z = 0.3 is d_L = 1600 Mpc for a ΛCDM cosmology with Ω_m = 0.27, Ω_λ = 0.73, and H₀ = 71 km s^-1 Mpc^-1 (Spergel et al. 2003). The angular and apparent speed of the individual components are given in Table 1. The calculated β_app values range from 6c to 37c, where the latter is unusually high for a blazar, especially for a BL Lac object.

Fig. 5 Evolution of radial separation of different components from the core. The solid lines represent fits with linear functions, while the dashed curves are quadratic fits.

For the non-ballistic components, NC2, NC3, and NC6, the quadratic fits represent significant acceleration along the jet axis. Formally, we obtain mas yr^-2 for NC2, mas yr^-2 for NC3, and mas yr^-2 for NC6. This apparent acceleration is a commonly observed feature of many objects in the MOJAVE sample (Lister et al. 2013), for which a statistically significant tendency for acceleration in the base of jets has been found. A standing shock slows the flow while deflecting it, after which the flow can accelerate back to the velocity it had upstream of the shock (Gomez et al. 1997). However, since the apparent speed was relatively slow in the upstream region, this does not seem to fit the case of S5 0716+714. The data are more consistent with the scenario that, after a moving component passes through the stationary feature, it bends toward or away from our line of sight. It then appears to move faster because its velocity vector subtends an angle closer to the vector that maximizes apparent motion, arccos(β), where β is the velocity of the component divided by the speed of light. An alternative possibility is that the underlying flow of the jet systematically accelerates outward. Theoretical models involve strong magnetic fields associated with the putative supermassive black hole/accretion disk system that plays a key role in the initial acceleration and collimation of the jet (Blandford & Payne 1982; Blandford & Znajek 1977; Meier et al. 2001). The conversion of Poynting flux to flow energy is gradual and may persist out to parsec scales (e.g., Sikora et al. 2005). With this explanation, however, it is difficult to explain the straight-line kinematic evolution of some of the components.

Variability Doppler factor

A variability Doppler factor using VLBI jet kinematics was defined by Jorstad et al. (2005) as ${\mathit{\delta }}_{\mathrm{VLBI}}\mathrm{=}\frac{\mathit{s}{\mathit{d}}_{\mathrm{L}}}{\mathit{c}\mathrm{\Delta }{\mathit{t}}_{\mathrm{var}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\mathit{,}$(2)where s is the angular size of the component, defined as 1.6 × φ for a Gaussian with FWHM diameter φ measured at the epoch of maximum brightness, and d_L is the luminosity distance. The timescale of variability is defined as Δt_var = dt/ ln(S_max/S_min) (Burbidge et al. 1974), where dt is the time separation in years between maximum (S_max) and minimum (S_min) flux densities. The estimated values of δ_VLBI for the individual components, listed in Table 1, range from 6 to 21. Using the broadband flux density and spectral variability study of the source over the same time period (Paper I, Sect. 3.4), we obtain a range of self-consistent values for the Doppler factor. The different independent approaches suggest that δ ≥ 20, which is consistent with upper range of the values of δ_VLBI obtained here.

3.1.2. Ejection epochs

Back-extrapolation of the components’ motion allowed us to estimate the time of zero separation from the core (i.e., the “ejection” time, T₀). As shown in the previous section, the components follow curved trajectories; therefore, we back-extrapolate the component trajectories in the XY-plane to determine the epoch of (0,0) separation. The error estimates are not very straightforward for calculating for components following non-ballistic trajectories. We used the estimated error on the fitted parameters for a given function to calculate its uncertainty. Back-extrapolation of the two envelopes of the function provided us with the error in T₀. The calculated ejection times for the individual components are listed in Table 1.

3.1.3. Brightness temperature gradient in the jet

The redshift-corrected brightness temperature (T_B,obs) of the bright emission features can be approximated using the following relation (Jorstad et al. 2005): ${\mathit{T}}_{\mathrm{B}\mathit{,}\mathrm{obs}}\mathrm{=}\mathrm{1.22}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}\frac{{\mathit{S}}_{\mathrm{comp}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}{\mathit{d}}_{\mathrm{comp}}^{\mathrm{2}}{\mathit{\nu }}^{\mathrm{2}}}\hspace{0.17em}\mathrm{K}\mathit{,}$(3)where S_comp is the component flux density in Jy, d_comp is the FWHM of the circular Gaussian component in mas, and ν is the observing frequency in GHz. For the calculations here, we have only used 43 GHz data. The calculated brightness temperatures (T_B,obs) of the components generally decline with radial separation from the core (see Fig. 6). This decaying trend can be approximated by a power law, ${\mathit{T}}_{\mathrm{B}\mathit{,}\mathrm{obs}}\mathrm{~}{\mathit{r}}_{\mathrm{jet}}^{\mathrm{-}\mathit{ϵ}}$ with ϵ = 2.43 ± 0.36, where r_jet is the radial distance along the jet axis. The solid line in Fig. 6 represents the fitted power law.

Fig. 6 Component brightness temperature as a function of radial separation from the core. The solid line represents the fitted power law with index ϵ = 2.43 ± 0.36.

We place the observed brightness temperature gradient in the context of the shock-in-jet model (Marscher & Gear 1985). In the common picture of this model, a relativistic shock propagates down a conical jet, slowly expanding adiabatically while maintaining shock conditions during propagation. While the shock propagates down the jet, it undergoes three major evolutionary stages dominated by inverse Compton, synchrotron, and adiabatic energy losses. As a result, the observed brightness temperature declines as a power law, T_B,obs$\mathrm{\propto }{\mathit{r}}_{\mathrm{jet}}^{\mathrm{-}\mathit{ϵ}}$. The value of ϵ can be derived from the spectral evolution of the radio emission (Lobanov & Zensus 1999). Following T_B ∝${\mathit{S}}_{\mathit{\nu }}{\mathit{r}}_{\mathrm{jet}}^{-2}{\mathit{\nu }}^{-2}$ and S(ν) $\mathrm{\propto }{\mathit{r}}_{\mathrm{jet}}^{\mathit{\xi }}{\mathit{\nu }}^{\mathit{\zeta }}{\mathit{\delta }}^{\mathrm{(}\mathit{s}\mathrm{+}\mathrm{3}\mathrm{)}\mathit{/}\mathrm{2}}$ (Lobanov & Zensus 1999), we obtain $\begin{array}{ccc}& & {\mathit{ϵ}}_{\mathrm{Compton}}\mathrm{=}\frac{\mathrm{11}\mathrm{-}\mathit{s}\mathrm{-}\mathit{a}\mathrm{(}\mathit{s}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{8}}\mathrm{+}\frac{\mathit{b}\mathrm{(}\mathit{s}\mathrm{+}\mathrm{3}\mathrm{)}}{\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{,}\\ & & {\mathit{ϵ}}_{\mathrm{synchrotron}}\mathrm{=}\frac{\mathrm{-}\mathrm{\left[}\mathrm{4}\mathrm{\right(}\mathit{s}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{+}\mathrm{3}\mathit{a}\mathrm{(}\mathit{s}\mathrm{+}\mathrm{1}\mathrm{\left)}\mathrm{\right]}}{\mathrm{6}}\mathrm{+}\frac{\mathit{b}\mathrm{(}\mathit{s}\mathrm{+}\mathrm{3}\mathrm{)}}{\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{,}\end{array}$and $\begin{array}{ccc}{\mathit{ϵ}}_{\mathrm{adiabatic}}\mathrm{=}\frac{\mathrm{2}\mathrm{(}\mathrm{5}\mathrm{-}\mathrm{2}\mathit{s}\mathrm{)}\mathrm{-}\mathrm{3}\mathit{a}\mathrm{(}\mathit{s}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{6}}\mathrm{+}\frac{\mathit{b}\mathrm{(}\mathit{s}\mathrm{+}\mathrm{3}\mathrm{)}}{\mathrm{2}}\mathrm{-}\mathrm{2}\mathit{,}& & \end{array}$(6)where a, b, and s parametrize the variations along the jet axis of the Doppler factor, $\mathit{\delta }\mathrm{\propto }{\mathit{r}}_{\mathrm{jet}}^{\mathit{b}}$, magnetic field, $\mathit{B}\mathrm{\propto }{\mathit{r}}_{\mathrm{jet}}^{\mathrm{-}\mathit{a}}$, and power law distribution of energy of the emitting electrons, N(γ) ∝ γ^{− s}. Assuming a constant Doppler factor, i.e., b = 0, and using a typical value of s = 2 (corresponding to optically thin synchrotron spectral index = −0.5), and a = 1 to 2 (for toroidal and poloidal magnetic fields, respectively), we obtain ϵ_Compton = 1.3 to 1.6, ϵ_synchrotron = 4.2 to 5.7, and ϵ_adiabatic = 3.2 to 4.7.

The observed brightness temperature gradient in S5 0716+ 714 (Fig. 8) has a slope equal to 2.43 ± 0.36. The slope of the observed intensity gradient rules out the simple assumptions of a constant Doppler factor. We therefore consider a variable Doppler factor along the jet axis. We find that the observed value of ϵ is consistent with the adiabatic loss phase if δ ∝${\mathit{r}}_{\mathrm{jet}}^{\mathit{b}}$, with b ~ 0.12 to 0.40, which corresponds to a moderate variation in the Doppler factor. However, with relatively larger variations in the Doppler factor b ~ 0.47 to 0.74, the observed value of ϵ would agree with synchrotron loss phase as well.

The evolution of the flux density and frequency (S_m and ν_m) at the spectral turnover from synchrotron self-absorption of the radio flares, discussed in Paper I (Sect. 3.3.2), also suggests a variation in δ along the jet axis. In that study, we found that δ changes as ${\mathit{r}}_{\mathrm{jet}}^{\mathit{b}}$ with b = 0.7 during the rise and b = 2.6 during the decay of the first radio flare. The evolution of the second flare was governed by b = 0.4 during the rising phase and b = −2.0 during the decay (see Paper I for details). Therefore, it is evident that the two flares require a different dependence of δ on distance from the core, which is also inconsistent with what we obtained for the brightness temperature gradient (b = 0.12 to 0.74). This can be explained in two ways. First, the single-dish radio flux correlates well with the VLBI core flux (see Sect. 3.2 for details). This implies that the radio flares originate within the core. The core is more compact and farther upstream than the jet components, with a comparatively higher brightness temperature, and can follow a different dependence of δ on r_jet. The individual components trace the intensity gradient farther downstream of the jet and can therefore individually have a different δ − r_jet dependence. A second possibility is that the brightness temperature gradient plot corresponds to an averaged behavior of the two radio flares.

Fig. 7 Broadband flux light curves of S5 0716+714: a) weekly averaged γ-ray flux light curve at E> 100 MeV; b) X-ray light curve at 0.3−10 keV; c) optical V-band light curve; d) single-dish (43 and 230 GHz) and 7 mm VLBI core flux light curves; and e) flux density curves of moving jet components shown in different colors. The same color is used to shade their respective ejection times, i.e., epochs of zero separation from the core. The green dashed lines indicate the average times when moving components cross the stationary feature at ~0.15 mas. The major outburst is marked with a horizontal arrow; the faster repetitive optical flares are labeled as “0” to “9” and sharp γ-rays flares as “A” to “J” (see text for details).

3.2. Jet kinematics and broadband emission

For the comparison of the observed broadband flaring behavior with the jet kinematics, we used the flux density light curves from Paper I, where the details of observations and data reduction are given. Broadband flux light curves are shown in Fig. 7. Panel a presents the weekly averaged γ-ray photon flux light curve. Panel b shows the X-ray light curve, and panel c displays the optical V-band light curve. Panel d presents the radio (single-dish) and 43 GHz core² flux density curves plotted in the same panel. Panel e exhibits flux density light curves of the individual components shown in different colors. The repetitive optical flares with a typical duration of ~1200−140 days are labeled as “0” to “9”, and the broadband flares are denoted “G” for γ-rays, “X” for X-rays, “O” for optical, and “R” for radio followed by the number adjacent to them; for example, we have flares G0 to G9 at γ-rays, X6 at X-rays, O0 to O9 at optical, and R6 and R9 at radio frequencies. The sharp γ-ray flares last ~10−20 days.

Although VLBI observations of S5 0716+714 are missing at the peak of the two major radio flares (JD′ ~ 1150−1200 and ~1300), there is still an apparent one-to-one correlation between the core flux and single-dish radio flux measurements. This suggests that the single-dish radio flares are produced in the core, which agrees with the high mean core-to-jet flux ratio, which ranges from 1.3 to 21.4 with a mean of 4.9 ± 4.0. The broadband flux cross-correlations (Paper I) also indicate that the flux variations at γ-ray (and optical) and single-dish radio frequencies are correlated such that the peak of the major γ-ray (and optical) outburst leads the peak of the radio outburst by a two-month time period.

Fig. 8 Sketch for the proposed scenario in the BL Lac S5 0716+714 (not to scale). North is to the right in this sketch. A shock propagates down the jet along a helical path, producing a major optical/γ-ray outburst upstream of the 7 mm VLBI core and later a major radio outburst accompanied by ejection of components from the core.

3.2.1. Component ejections and their connection with broadband flares

Because of the superposition of different modes of flaring activity, the broadband flaring behavior is quite complex in the source. The major outburst (Fig. 7) is accompanied by faster repetitive optical/γ-ray flares and sharp γ-rays flares. The jet kinematic study during this period suggests significant changes in the jet morphology of the source; in total, seven new components were ejected. It is therefore plausible that the ejection of new components could be connected to the broadband flares. We have identified a stationary feature at a distance of ~0.15 mas from the core (Sect. 3.1.1). Interaction of a moving component with the stationary feature could amplify the magnetic field and accelerate particles to produce a multiwaveband flare (Marscher 2014; Fromm et al. 2012).

In Fig. 7, the shaded areas represent the ejection times of new jet components, and the figure shows the average epochs when a moving component crosses the stationary feature at ~0.15 mas. In Table 2, we list the new components (Col. 1), along with the flares that coincide with their respective ejection times (Col. 2). Flares observed closer in time to the passage of a moving component across the stationary feature are listed in Col. 3. We note that most of the sharp γ-ray flares (B, C, D, G, I, and J) are observed when a moving component crosses the stationary feature, while most of the repetitive faster flares (0, 2, 3, 5, and 6) coincide with the ejections of new components from the core.

In Paper I, we found that the high-energy (optical and γ-ray) emission is significantly correlated with the radio flares such that the major optical/γ-ray outburst propagates down to radio frequencies (R6 flare) with a time delay of ~65 days, following a power-law dependence on frequency with a slope ~0.3. We also noticed that the evolution of radio flare R6 is in accordance with a generalized shock model (Valtaoja et al. 1992). This suggests that the origin of the major outburst could be related to a disturbance or shock that moves down the jet.

We interpret the broadband flaring events in the following manner (see Fig. 8). A disturbance at the base of the jet could cause a shock wave to form and propagate down the jet (Valtaoja et al. 1992; Marscher & Gear 1985), which could be related to the onset of the major optical/γ-ray outburst (JD′ ~ 840) and disturbs the jet outflow, resulting in the ejection of two new components (NC2 and NC3, Fig. 8a). The core is optically thick at that point, since the shock forms upstream of the core; for this reason, we do not see much activity in the core. Repetitive optical/γ-ray flares and the curved trajectories of the associated components favor the propagation of the shock along a bent trajectory, possibly a helical path. The moving shock reaches its maximum strength at JD′ ~ 1150 (Fig. 8b), close to the peak of the major optical/γ-ray outburst (see Sect. 3.3, Paper I for details), which also coincides with the onset of radio flare R6. At this point, the shock reaches the core, which starts to become optically thin. After this, other disturbances passing through the core correspond to the ejections of three new components (NC3, NC4, and NC5). Major radio outburst R6 is observed ~65 days later, accompanied by the ejection of component NC7 (Fig. 8c). As the moving features cross the the stationary feature at ~0.15 mas, at first, the stationary features are disrupted after the interaction, and later it appears again. The times of passage of moving components through the stationary feature appear to coincide with the sharp γ-ray flares in most cases. Therefore, the sharp γ-ray flares seem to be produced via the interaction of moving features with the stationary one.

X-ray flare X6 coincides with a similar increase in the flux of the 43 GHz core. Optical flare O7 peaks about a month later, which roughly coincides with the expected passage through the stationary feature of the same disturbance that caused flare X6. However, no bright moving knot is apparent in the next VLBA image ~20 days after O7. Also, the dramatic decrease in the millimeter-wavelength flux during this interval implies that the disturbance is strongly quenched following flare O7. There are two anomalies in this sequence of events: there is no optical flare as the disturbance crosses the core, and there is no γ-ray activity during the entire time span covering flares X6 and O7. The first indicates that, for reasons that are unclear (e.g., unfavorable magnetic field orientation relative to a shock; Summerlin & Baring 2012), electrons with energies high enough to produce optical and γ-ray emission were not produced in the core region at this time. The second implies that there was a dearth of seed photons for inverse-Compton scattering ~0.15 mas from the core during flare O7.

3.2.2. Doppler beaming of radio and gamma-ray emission regions

The observed time lag between γ-ray and 43 GHz core flares can be used to estimate the separation between the two emission regions by using the following equation (Fuhrmann et al. 2014): $\mathrm{\Delta }{\mathit{r}}_{\mathit{\gamma ,r}}\mathrm{=}\frac{{\mathit{\beta }}_{\mathrm{app}}\mathit{c}\hspace{0.17em}{\mathit{\tau }}_{\mathit{\gamma ,r}}}{\sin \mathit{\theta }}\mathit{,}$(7)where θ is the viewing angle of the source, β_app the apparent jet speed, and τ_γ,r the observed time lag. Using the observed value τ_γ,r = 82 ± 32 days (Rani et al. 2014), we have Δr_γ,r = (0.069 ± 0.026)β_app(sinθ)^-1 pc.

The measured apparent speeds in the source cover a wide range. The current observations suggest slower motion (≤8 c) on sub-mas scales (~0.3 mas), which later are accelerated to speeds as fast as 40 c on mas scales. The source also has a prominent bent-jet structure at sub-mas scales (see Fig. 2); as a result, the sub-mas and mas regions could have different viewing angles. For β_app = 6−8c, the viewing angle³, θ, is 6−9°, and this gives Δr_γ,r = 2.9−4.4 pc. Using β_app = 40c, we have θ = 1.4° and Δr_γ,r ~ 110 pc. Even for a larger viewing angle, θ = 10°, we obtain Δr_γ,r ~ 16 pc. Pushkarev et al. (2012) have estimated that the 15 GHz core is located at a distance of 6.68 pc from the jet apex. Using the technique to align the optically-thin emission based on model-fitted components and morphological similarities, we derive a shift of 0.042 ± 0.013 mas (0.19 ± 0.06 parsec) between the 15 GHz and 43 GHz core, which puts the location of the 43 GHz core at a distance of ~6.5 pc from the jet apex. Therefore, the distance between γ-ray emission region and 43 GHz core cannot be larger than 6.5 pc. For a typical viewing angle θ ≈ 5°, Δr_γ,r ≤ 6.5 pc gives β_app ≤ 8c. However, we note that for a larger viewing angle (θ ~ 10°), β_app could be as high as 10−12c but not ~40c unless the viewing angle is larger than 30°, which is unlikely to be the case for S5 0716+714. The results therefore suggest that the high-energy and radio emission regions have different apparent speeds and probably different viewing angles as well. This implies that the γ-ray and radio emission regions have different Doppler factors.

4. Summary and conclusion

The monthly sampled mm-VLBI monitoring of S5 0716+714 allows us to study the inner jet kinematics and to investigate its relation to the observed broadband flux variability. Our study reveals significant non-radial motions in the jet outflow, including variations in the orientation of the sub-parsec-scale jet (see Rani et al. 2014, for details) and possibly helical trajectories of moving components.

In radial directions, the individual components exhibit extreme apparent speeds as high as ~37 c. An apparent speed of ~20 c seems to be more typical in S5 0716+714 on parsec scales. Recently, Lister et al. (2013) have reported an apparent speed of 43.6 ± 1.3c. These values would even be extreme for a quasar. Therefore, moving components in the jets of BL Lac objects can be just as fast as in quasars, despite flux-limited surveys at lower resolution (e.g., Lister et al. 2009) finding that the mean apparent velocity of BL Lac objects is lower than for quasars. Further mm-VLBI studies are needed to determine the extent to which angular resolution and the variety of spectral classes of BL Lac objects (low- to high-frequency peaks of the synchrotron spectral energy distributions) skew the apparent velocity distribution. It is possible that the fast BL Lacs are actually quasars in disguise, with their broad emission lines swamped by the non-thermal optical continuum. On the other hand, emission lines have never been observed in S5 0716+714.

The bright components in 0716+714 follow a continuous decay trend in the observed brightness temperature (${\mathit{T}}_{\mathrm{B}\mathit{,}\mathrm{obs}}\mathrm{\propto }{\mathit{r}}_{\mathrm{jet}}^{\mathrm{-}\mathrm{(}\mathrm{2.36}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.41}\mathrm{)}}$), which implies a variation in Doppler factor along the jet axis. The spectral evolution of radio flares (see Paper I for details) also suggests variations in Doppler factor. A significant correlation between γ-ray flux and variations in orientation of the sub-parsec scale jet (Rani et al. 2014) supports the possibility that the observed flux variations in S5 0716+714 are strongly influenced by Doppler factor variations, and therefore may have a geometrical origin.

Our study suggests a connection between jet kinematics and the observed broadband flaring activity. We find that the origin of major optical/γ-ray/radio outbursts in the source seems to be related to a disturbance, such as a shock wave, propagating down the jet. The repetitive, faster optical/γ-ray flares, which are superposed on the major outbursts, and non-ballistic motions of the associated jet components can be interpreted via helical jet models (as discussed by, e.g., Blandford & Payne 1982; Blandford & Znajek 1977; Camenzind & Krockenberger 1992) and/or magnetohydrodynamic instabilities (Hardee 2006). Our analysis favors the major optical/γ-ray flares in S5 0716+714 being produced upstream of the 7 mm VLBI core. The sharp γ-ray flares, however, seem to be produced when the moving components cross a stationary feature 0.15 mas from the core. Finally, our analysis suggests that the upstream γ-ray emission region has a different Doppler factor than observed in the parsec-scale radio jet.

Online material

Appendix A: Model fit component results

Acknowledgments

The Fermi/LAT Collaboration acknowledges the generous support of a number of agencies and institutes that have supported the Fermi/LAT Collaboration. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Énergie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. This research made use of data obtained with the Global Millimeter VLBI Array (GMVA), which consists of telescopes operated by the MPIfR, IRAM, Onsala, Metsahovi, Yebes, and the VLBA. The data were correlated at the correlator of the MPIfR in Bonn, Germany. The VLBA is an instrument of the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This study makes use of 43 GHz VLBA data from the VLBA-BU Blazar Monitoring Program (VLBA-BU-BLAZAR; http://www.bu.edu/blazars/VLBAproject.html), funded by NASA through the Fermi Guest Investigator Program. The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. BR is thankful to Svetlana Jorstad, Benoit Lott, Biagina Boccardi, and Shoko Koyama for useful discussion and comments on the paper. We thank the referee for constructive comments.

References

All Tables

All Figures

Fig. 1 43 GHz VLBA images (contours) and Gaussian components of S5 0716+714 at different epochs between September 2008 and October 2010. The date of observation is given at the bottom of each map (shifted Julian day, JD−2 454 000). The fitted model parameters for each epoch are given in Table A.1. The cross at the bottom of each map represents the FWHM of the restoring beam (0.25 × 0.18 mas). Contour levels correspond to −0.003 0.003 0.006 0.012 0.024 0.048 0.096 0.19 0.38 0.76 1.54 3.07 6.14 12.28 24.58 49.15 98.30 mJy/beam.

In the text

	Fig. 2 Inner jet structure of S5 0716+714 from GMVA observations at 3 mm. The cross at the bottom of each map represents the FWHM of the restoring beam (0.07 × 0.04 mas). Contour levels correspond to −0.004 0.004 0.008 0.016 0.032 0.064 0.12 0.25 0.51 1.02 2.04 4.09 8.19 16.38 32.77 65.54 mJy/beam.
In the text

	Fig. 3 Evolution of core separation of the individual components as a function of time.
In the text

	Fig. 4 Trajectories of the jet knots in the X–Y plane, where X = r cos(θ + 90) and Y = r sin(θ + 90). r is the radial separation from the core and θ is the position angle with respect to an imaginary north-south line drawn through the map center.
In the text

	Fig. 5 Evolution of radial separation of different components from the core. The solid lines represent fits with linear functions, while the dashed curves are quadratic fits.
In the text

	Fig. 6 Component brightness temperature as a function of radial separation from the core. The solid line represents the fitted power law with index ϵ = 2.43 ± 0.36.
In the text

Fig. 7 Broadband flux light curves of S5 0716+714: a) weekly averaged γ-ray flux light curve at E> 100 MeV; b) X-ray light curve at 0.3−10 keV; c) optical V-band light curve; d) single-dish (43 and 230 GHz) and 7 mm VLBI core flux light curves; and e) flux density curves of moving jet components shown in different colors. The same color is used to shade their respective ejection times, i.e., epochs of zero separation from the core. The green dashed lines indicate the average times when moving components cross the stationary feature at ~0.15 mas. The major outburst is marked with a horizontal arrow; the faster repetitive optical flares are labeled as “0” to “9” and sharp γ-rays flares as “A” to “J” (see text for details).

In the text

	Fig. 8 Sketch for the proposed scenario in the BL Lac S5 0716+714 (not to scale). North is to the right in this sketch. A shock propagates down the jet along a helical path, producing a major optical/γ-ray outburst upstream of the 7 mm VLBI core and later a major radio outburst accompanied by ejection of components from the core.
In the text