© ESO, 2011

1. Introduction

The blazar S5 0716+714 at redshift z = 0.31 (Nilsson et al. 2008) is one of the most active sources of its class. It is highly variable on timescales from hours to months across all observed wavelengths. Intra-day variability (IDV) has been detected in this blazar (Wagner et al. 1996; Stalin et al. 2006; Montagni et al. 2006; Gupta et al. 2009). The IDV in radio band is likely to be intrincic to the source because it is either reported to be correlated with the optical band (Wagner et al. 1996) or its spectrum is not consistent with the scintillation (Fuhrmann et al. 2008), which implies the presence of very small emitting regions within the source. At larger time scales this source has frequent flares in both radio and optical bands. Amplitudes of the flares increased abruptly after the year 2003 (e.g. Raiteri et al. 2003; Nesci et al. 2005; Bach et al. 2007; Fuhrmann et al. 2008, and references therein). In 2007 and 2008 it was also detected at very high-energy gamma rays with the MAGIC telescope (Kranich & the MAGIC collaboration 2009; Anderhub et al. 2009).

At the milliarcsecond scale (VLBI observations) 0716+714 has a core-dominated structure with a bright point-like core and a faint diffuse jet pointing at the position angle of  ~20°. At the scale of arcseconds (VLA observations), the source structure is also core-dominated, its jet is double-sided, located along the position angle of  ~100° (misaligned with the milliarcsecond scale jet by  ~80°), and has extended diffuse lobes (Wagner et al. 1996).

The long-term kinematics of this source have been extensively studied, and two mutually exclusive models have been suggested: one implies relatively fast outward motion of the superluminal components in the jet with speeds ranging from 5c to 16c (“fast scenario”, Bach et al. 2005), and the other claims the jet to be stationary, with significant non-radial (“precession”) motion of the jet as a whole (“stationary scenario” Britzen et al. 2009). Both these scenarios are based on the reanalysis of historical VLBI observations that span the years 1992–2001 (“fast”) and 1992–2006 (“stationary”) at frequencies from 5 to 43 GHz.

Our analysis of the inner jet kinematics of 0716+714 at frequencies of 22, 43, and 86 GHz also revealed a fast-moving jet with apparent speeds from 8c to 19c on the spatial scale from 0 to 1.5 mas (Rastorgueva et al. 2009). Here we present results of our short-term kinematic study at 5 GHz. The data set covers six months, and the source was observed nearly every month. With this dense data sampling we are able to study the fast structural changes in the jet of 0716+714.

The standard method of VLBI data reduction, which is used by the vast majority of researchers in the field, includes the Difmap package (Shepherd 1997) as an imaging facility, which became conventional because of its simplicity and accessibility. However, Difmap has its downsides, and the main disadvantage of this package, which uses the standard CLEAN algorithm by Högbom (1974; Shepherd et al. 1994), is a problem with the reconstruction of extended emission (Cornwell 1983). Tests performed on the simulated data have shown that the standard CLEAN, and Difmap in general, are especially inaccurate in the reconstruction of the sources whose structure could be described as bright point-like knots embedded in a region with faint extended emission (e.g. Baikova 2007; Lawson et al. 2004). The morphology of 0716+714 at lower frequencies (5–1.6 GHz) is consistent with this description. Therefore, a conventional approach to this source may have yielded erroneous results. In order to make sure that the structure of 0716+714 is correctly reconstructed, we independently applied both the generalized maximum entropy method (GMEM) and the conventional CLEAN to the same data set.

This paper is arranged as follows: Sect. 2 describes the observations, Sect. 3 contains a description of the a-priori data calibration and an introduction to the imaging methods. In Sect. 4 we present and discuss the kinematic results. The final kinematic scenario, a comparison with previous studies, and possible explanations of the apparent speed differences are presented in Sect. 5.

2. Observations

We observed 0716+714 with the Very Long Baseline Array (VLBA) five times during the year 2004. Each observing session lasted nine hours and covered five frequencies (86, 43, 22, 5 and 1.6 GHz) in dual polarization mode. Here we describe the results of the 5 GHz observations. The five epochs in our experiment were named in alphabetical order: A, Feb. 10; B, May 03; C, Jun. 18; D, Jul. 29 and E, Aug. 29, 2004. Additionally, we used one epoch from the unpublished polarization observation project of 37 northern BL Lac objects from the Kuehr & Schmidt sample (Kuehr & Schmidt 1990) by Gabuzda and Mahmud, in our notation epoch M, Mar. 22, 2004. Epoch M fills the gap between A and B of our experiment, making the time sampling almost uniform with an interval of approximately one month between observations.

3. Data reduction

3.1. Calibration

The data were correlated at the Socorro VLBA correlator. A-priori amplitude and phase calibration were made at Tuorla Observatory with NRAO’s Astronomical Image Processing System (AIPS). AIPS’s standard VLBAPROC procedures were used for the amplitude and phase calibrations, scans on 3C 279 were used to calculate instrumental delay and phase residuals, and OJ287 was used as an instrumental polarization (D-term) calibrator. The subsequent procedure for our two methods is described below.

3.2. Image reconstruction: the conventional method

The first method we used was the conventional method, which is a combination of CLEAN deconvolution (Högbom 1974) with self-calibration and uv-plane Gaussian model fitting. This method is used by most of the astronomical community in studies of outflow structure and kinematics.

Fig. 1 Images of 0716+714 obtained using the conventional method: CLEAN maps, convolved with the averaged beam of 1.5  ×  2.5 mas, PA =  −4.5°. Circles with crosses are the Modelfit Gaussian components superimposed on the image. The Gaussian model fitting was performed in the uv-plane. Axes are right ascension and declination in milliarcseconds, measured with regard to the core.

We used the Caltech Difmap package (Shepherd 1997) for imaging, employing a CLEAN algorithm in a self-calibration loop: a model of the source, consisting of point-like clean-components is created with CLEAN, then the antenna gains are self-calibrated with the model, and the procedure is repeated. Self-calibrations were performed in the following way: several phase self-calibration loops were followed by an amplitude one, with the solution interval of amplitude self-calibration gradually decreasing from the time of the whole observation to the data averaging time (0.5 min). Finally, we performed amplitude self-calibration with the solution interval equal to the averaging time. Then the final clean-component model was discarded, and a new source model was created. We used the Difmap procedure, Modelfit, to create a model consisting of circular components with a Gaussian distribution of flux, and fit it to the visibility data. Except for one case, only circular Gaussian components were used to parameterize the data in order to reduce the number of free parameters. For all models, the brightest optically thick component at the southernmost end of the jet was considered to be the core, as in the previous research. Positions of other components were measured with regard to this. The resulting maps with the superimposed model components are presented in Fig. 1, and map parameters for all epochs are listed in Table 1.

The errors of the position and the flux estimates of the components were calculated with the Difwrap package (Lovell 2000), where the parameter to be analyzed (position, flux or FWHM) is varied around the best-fit value in a grid of a given step size and Modelfit is run for each step. The result of this process is plotted as a χ² distribution, and the user can visually inspect the model vs. a uv-data plot or a residual map for each point of the grid. Visual inspection is based on the relative χ² changes. Because the number of degrees of freedom of the model are not exactly known, this method of error estimate based on visual inspection yields more reliable results than the standard χ² analysis. In our case it was not possible to use the model vs. a visibility plot, because we dealt with very faint components of the jet: changes in the model visibilities between different grid points were too small to be detected with the naked eye. Instead, we estimated the goodness of fit based on the appearance of the residual map and its changes from one grid point to another. We applied an intuitive criterion proposed by Savolainen (2006) to determine the maximal acceptable deviation from the best-fit model: the difference between a model and observed data is considered significant if the structure appearing in the residual map would look strong enough that it should be cleaned, if it had appeared during the imaging process. Because the flux of such a feature changes gradually from one grid point to another, it may be difficult to determine at which step the discrepancy between the observed data and the model becomes large enough. In these cases we used an empirical “2σ-rule”, which we worked out by analysis of both the imaging process and the error estimate procedure: the flux is considered ‘excessive’ and should be “cleaned out” by placing a model component in its location if the rms noise in the area under scrutiny is twice as high as the average rms noise at the outskirts of the map. At the grid point where the noise level exceeds twice the rms, we assumed that the parameter value had reached the limit. This deviation from the best-fit value was assumed to be the limiting error, Δ. A more detailed discussion of VLBI data error estimate with Difwrap and a comparison with the method of variance about the best fit can be found in Savolainen (2006).

In most cases the Difwrap package yields asymmetric uncertainties for the model parameters, which leads to a bias in the resulting models. This bias was corrected for with the method suggested by D’Agostini (2003), based on the probabilistic approach to the bias correction: let us assume that the negative limiting error of the component’s separation from the core r is Δ_negr, and the positive limiting error is Δ_pos, and Δ_neg ≠ Δ_pos. We calculated the unbiased limiting error as an arithmetic average of the two limiting errors: Δ_un = (Δ_neg + Δ_pos)/2. The standard deviation of r was calculated as σ_r = Δ_un/3. The unbiased value of the model parameter (in this case, separation from the core) is not equal to the result of the model fitting and is r_un = r + (Δ_pos − Δ_neg)/2.

3.3. Image reconstruction: the maximum entropy method (MEM), an attempt

Another popular deconvolution algorithm in radio interferometry is the maximum entropy method invented by Jaynes (Jaynes 1957), which was first applied to image reconstruction by Frieden (Frieden 1972). The maximum entropy method is an extremely nonlinear method that leads to maximally smoothed images subject to constraints. The positivity of the sought-for solution forces super-resolution on bright isolated objects.

There are many forms of entropy (Narayan & Nityananda 1986). We used the AIPS task VTESS, which uses the entropy functional in the form (Cornwell 1995) $\begin{array}{ccc}\mathit{E}\mathrm{\left(}\mathit{I}\mathrm{\right)}& \mathrm{=}& \mathrm{-}\mathrm{\int }\mathit{I}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\ln \mathrm{\left[}\mathit{I}\mathrm{\right(}\mathit{x,y}\mathrm{)}\mathit{/}\mathit{M}\mathrm{(}\mathit{x,y}\mathrm{\left)}\mathrm{\right]}\mathrm{d}\mathit{x}\mathrm{d}\mathit{y,}\\ & & \mathit{I}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathit{>}\mathrm{0}\mathit{,}\end{array}$(1)where M(x,y) is a default image that represents the optimum solution in the absence of data. Usually, the default image is taken to be flat or a lower resolution image.

Given a form of entropy, the MEM defines the best image that maximizes E(I) subject to the data constraints: 1) fit the observed visibilities V_k to within the error term ϵ_k: F(I) = V_k + ϵ_k; 2) fit the total flux: ^∫I(x,y)dxdy = F₀; 3) the rms error must be plausible: $\sum \frac{\mathrm{|}{\mathit{ϵ}}_{\mathit{k}}{\mathrm{|}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}}\mathrm{=}\mathrm{\Omega }$.

In AIPS, the Cornwell-Evans algorithm (Cornwell & Evans 1985) is used for the realization of the MEM. This algorithm uses a simple Newton-Raphson approach to optimize the functional (1) subject to the constraints upon the rms and total power enforced by the Lagrange multipliers.

We made an attempt to use the MEM method to recover the faint emission of the 0716+714 jet. We imported visibility data, self-calibrated in Difmap (see Sect. 3.2), to AIPS and used the task VTESS to create the final image.

We used the task IMAGR to clean the central part of the image containing the bright core in order to avoid circular-shaped artifacts around the point-like source. Nevertheless, the image quality of this method was poor and did not allow us to correctly perform the kinematic analysis. Therefore, we conclude that the MEM method is not suitable for recovering the emission of the sources with bright point-like features and faint underlying diffuse emission.

3.4. Image reconstruction: the generalized maximum entropy method (GMEM)

The generalized maximum entropy method has been developed and implemented in the Pulkovo VLBI data reduction software package VLBImager by Baikova (2007). The mapping technique realized in this package utilizes a self-calibration algorithm (Cornwell & Fomalont 1999) combined with the GMEM (Frieden & Bajkova 1994; Baikova 2007), which is used as a deconvolution procedure.

The GMEM is designed for the reconstruction of sign-variable functions, therefore it allows one to obtain unbiased solutions. The bias of the solution is one of the problems of the conventional MEM (Cornwell et al. 1999), which may lead to a substantial nonlinear distortion of the final image if the data contain errors (Baykova 1995).

For the GMEM, the Shannon-entropy functional has the form $\begin{array}{ccc}\mathit{E}\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}& \mathrm{=}& \mathrm{-}\mathrm{\int }\mathrm{\left(}{\mathit{x}}^{\mathrm{p}}\mathrm{\right(}\mathit{t}\mathrm{\left)}\ln \mathrm{\right(}\mathit{\alpha }{\mathit{x}}^{\mathrm{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{)}\mathrm{+}{\mathit{x}}^{\mathrm{n}}\mathrm{(}\mathit{t}\mathrm{\left)}\ln \mathrm{\right(}\mathit{\alpha }{\mathit{x}}^{\mathrm{n}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\left)}\mathrm{\right)} \mathrm{d}\mathit{t,}\\ & & {\mathit{x}}^{\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{>}\mathrm{0}\mathit{,} {\mathit{x}}^{\mathit{n}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{>}\mathrm{0}\mathit{,}\end{array}$(2)where x^p(t) and xⁿ(t) are the positive and negative components of the sought-for image x(t), i.e. the equation x(t) = x^p(t) − xⁿ(t) holds. α > 0 is a parameter responsible for the accuracy of the separation of the negative and positive components of the solution x(t), and therefore critical for the resulting image fidelity. It is easy to show (Baikova 2007) that solutions for x^p(t) and xⁿ(t) obtained with the Lagrange optimization method are connected by the expression x^p(t)·xⁿ(t) = exp(− 2 − 2lnα) = K(α), which depends only on the parameter α. This parameter is responsible for dividing the positive and negative parts of the solution: the larger α is, the more accurate is the discrimination. On the other hand, the value of α is constrained by computational limitations. The main constraint comes from the χ² term in the optimized functional, which depends on the data errors. The larger a standard deviation is, the higher value of α could be set. If data are very accurate, a lower value of α is needed. In practice, α is chosen empirically. In our case we had to compromise between data errors, which determine resolution of the resulting MEM solution, and a need to divide the positive and negative parts of the solution as accurately as possible. It is fair to say that given fixed errors in the data, a maximum possible chosen value of α provides us with the best possible resolution of the MEM-solution. In this work we used α = 100   000. The algorithm takes into account errors of the measurements. It implies a solution of the following two-dimensional problem of conditional optimization, given in the discrete form $\begin{array}{ccc}\begin{array}{}\min \end{array} \{\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{\left[}{\mathit{x}}_{\mathit{kl}}^{\mathrm{p}}\ln \mathrm{\right(}\mathit{\alpha }{\mathit{x}}_{\mathit{kl}}^{\mathrm{p}}\mathrm{)}\mathrm{+}{\mathit{x}}_{\mathit{kl}}^{\mathrm{n}}\ln \mathrm{(}\mathit{\alpha }{\mathit{x}}_{\mathit{kl}}^{\mathrm{n}}\mathrm{\left)}\mathrm{\right]}& & \\ & & \mathrm{+}\sum _{\mathit{m}\mathrm{=}\mathrm{1}}^{\mathit{M}}\frac{\mathrm{(}{\mathit{\eta }}_{\mathit{m}}^{\mathrm{re}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{\eta }}_{\mathit{m}}^{\mathrm{im}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{m}}^{\mathrm{2}}}\}\mathit{,}\end{array}$(3)$\begin{array}{ccc}& & \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{x}}_{\mathit{kl}}^{\mathrm{p}}\mathrm{-}{\mathit{x}}_{\mathit{kl}}^{\mathrm{n}}\mathrm{)}{\mathit{a}}_{\mathit{kl}}^{\mathrm{m}}\mathrm{+}{\mathit{\eta }}_{\mathit{m}}^{\mathrm{re}}\mathrm{=}{\mathit{V}}_{\mathit{m}}^{\mathrm{re}}\mathit{,}\\ & & \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{x}}_{\mathit{kl}}^{\mathrm{p}}\mathrm{-}{\mathit{x}}_{\mathit{kl}}^{\mathrm{n}}\mathrm{)}{\mathit{b}}_{\mathit{kl}}^{\mathrm{m}}\mathrm{+}{\mathit{\eta }}_{\mathit{m}}^{\mathrm{im}}\mathrm{=}{\mathit{V}}_{\mathit{m}}^{\mathrm{im}}\mathit{,}\end{array}$where k and l are the pixel numbers in a source map of the size N × N, m = 1,...,M is a number of the visibility function measurement, ${\mathit{a}}_{\mathit{kl}}^{\mathit{m}}$ and ${\mathit{b}}_{\mathit{kl}}^{\mathit{m}}$ are the constant coefficients corresponding to the Fourier-transform, ${\mathit{V}}_{\mathit{m}}^{\mathrm{re}}$ and ${\mathit{V}}_{\mathit{m}}^{\mathrm{im}}$ are the real and imaginary parts of the complex visibility function (the Fourier spectrum of a source) and ${\mathit{\eta }}_{\mathit{m}}^{\mathrm{re}}$ and ${\mathit{\eta }}_{\mathit{m}}^{\mathrm{im}}$ are the real and imaginary parts of the additive visibility noise with normal distribution over zero value with known dispersion σ_m.

As we can see from Eq. (3), the optimized functional has two parts: a Shannon-entropy functional and a functional that is an estimate of the difference between the reconstructed visibility spectrum and the measured data according to a χ² criterion. This latter functional can be considered to be an additional stabilizing term acting to provide a further regularization of the solution above what is possible with the entropy functional alone. Solution normalization, necessary for the MEM is obtained by limiting the total flux of the source: $\begin{array}{ccc}& & \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{N}}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{(}{\mathit{x}}_{\mathit{kl}}^{\mathrm{p}}\mathrm{-}{\mathit{x}}_{\mathit{kl}}^{\mathrm{n}}\mathrm{)}\mathrm{=}{\mathit{F}}_{\mathrm{0}}\mathit{.}\end{array}$Numerical experiments show that the GMEM combined with the difference mapping technique allows the restoration of the source structures consisting of bright compact features embedded in a weaker extended background with maximal accuracy (Baikova 2007).

Fig. 2 GMEM solutions for each epoch, convolved with the averaged CLEAN beam: 1.5  ×  2.5 mas, PA  =  −4.5°. Gray ellipses are the Gaussian components superimposed on the image. The Gaussian model-fitting was performed in the image plane for the original GMEM solutions (Fig. A.1) convolved with the beam of 0.5  ×  0.5 mas (Fig. A.2). Axes are right ascension vs. declination in milliarcseconds, measured with regard to the core.

After the a priori calibration, the raw visibility data were averaged in frequency (averaged IFs) in AIPS, then exported to Difmap, where they were rebinned, phase self-calibrated with a point source model with a 0.5 min solution interval, and finally averaged in time (over a period of four minutes). After that the data was exported to VLBImager, where the automatic GMEM deconvolution – self-calibration loop iterations were performed until convergence. The size for the maps is chosen 256    ×    256 pixels with a pixel size of 0.15 mas. The GMEM solution convolved with a clean beam of  ~ 2.4    ×    1.5 mas, PA =  −4.5° is presented in Fig. 2, and parameters for those maps are listed in Table 2.

In the Appendix we present two additional sets of maps to demonstrate the capabilities of this deconvolution method: the original GMEM solutions (Fig. A.1) and the GMEM solutions convolved with a circular clean beam of the size 0.5 × 0.5 mas, which were used for the modeling of the source structure (Fig. A.2). The model-fitting of the maps was made in the image plane. The source structure was modeled using Gaussian elliptical components. Each of the components was described by the following six parameters: r, distance from the core (the core was determined like in the conventional method and was assumed to be stationary); φ, position angle measured counter-clockwise between the ordinate axis and the radius-vector; a_maj and a_min, sizes of the major and minor axes of a Gaussian ellipsoid at half maximum amplitude and θ, the angle describing the orientation of the component ellipse, an angle between the ordinate axis and the major axis of the ellipse, measured counter-clockwise. Cross-identification of the components was based on the assumption that the core separation, the position angle and the flux of the component change smoothly from epoch to epoch. In order to obtain a model of the inner jet (r < 1.5 mas) with higher accuracy, we interpolated the brightness distribution to a finer grid: the pixel size in thit case was 0.0375 mas, which was four times smaller than that for the large-scale jet. The errors of the component core separation determination were estimated by the formula (Fomalont 1999) $\begin{array}{ccc}& & \mathrm{\Delta }\mathit{r}\mathrm{=}\frac{\mathit{\sigma }\mathrm{\Theta }}{\mathrm{2}{\mathit{S}}_{\mathrm{peak}}}\mathit{,}\end{array}$where σ is the rms after the subtraction of the model from the source map, S_peak is the peak flux of the component, Θ is the size of the component at the amplitude half maximum (in our case we assume Θ = a_maj). We point out that that this formula likely underestimates position errors, which is clear from the comparison with the Difwrap errors (see Fig. 3).

4. Results

Our aim was to study the structure and kinematics of the large-scale (from 1 to 12 mas) jet, and compare it with the results for the inner (< 1 mas) jet, which were published in Rastorgueva et al. (2009).

We applied a conventional method and a GMEM to the same data set to compare their kinematic results. For each of them, we calculated the core separation of the components as a function of time, the proper motion μ, and the apparent speeds β_app according to the formula $\begin{array}{ccc}& & {\mathit{\beta }}_{\mathrm{app}}\mathrm{=}\frac{\mathit{\mu }{\mathit{d}}_{\mathrm{L}}}{\mathit{c}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\mathit{,}\end{array}$where the luminosity distance d_L is calculated with the analytical formula derived by Pen (1999) for a flat cosmology. The same analytical expression was used by Bach et al. (2005) for their kinematic analysis of this source, therefore our results could be compared directly without conversion. For this cosmological model, a proper motion of 1 mas/yr is equivalent to the apparent speed of 19.26c. The flat Universe model with the Hubble constant of H₀ = 71 km s^-1 Mpc^-1 and with Ω_m = 0.3 and Ω_Λ = 0.7 was used. We are aware that Wickramasinghe & Ukwatta (2010) recently obtained a new and more precise analytical expression for the luminosity distance, but for z = 0.3 the difference with Pen (1999) is negligible (1 mas/yr corresponds to the apparent speed of 19.22c).

The results of the Gaussian model-fitting and component identification are presented in Table 5 (conventional method) and in Table 6 (GMEM). The original GMEM-solutions are shown in Fig. A.1.

4.1. Component cross-identification between different methods.

4.1.1. Positions and apparent speeds of components

Based on core separation and fluxes, we identified four nearly stationary components in the large-scale jet of 0716+714: one at a separation of about  ~1.5 mas (C3–A3), one at  ~3.0 mas (C4–A4), one at  ~4.6 mas (C5–A5) and one at  ~12 mas (C6–A6). Their proper motions and apparent speeds can be found in Tables 3 and 4. A comparison of their positions, obtained with different methods, and linear fits are shown in Fig. 3. The components C1–C6, found by the conventional method, have a larger scatter around the mean position and larger positional errors than those of the GMEM components (compare Tables 3 and 4), which complicated the component cross-identification and made a choice of the kinematic scenario ambiguous. However, the GMEM method confirmed the stationary scenario.

However, the component C1–2, determined by the conventional method, and the component A1 found by the GMEM in the inner jet have a clear outward motion:  ~10c and  ~2c, respectively (Tables 3 and 4). The difference in apparent speed could be explained by the fact that the GMEM determined an extra stationary component at the separation of about 1 mas from the core (A2), while the conventional method did not. On the other hand, A2 is only found at the last four epochs, and is separated from the component A1 by only  ~0.5 mas, and the component C1–2, determined by the conventional method, is located between them (see Fig. 3). Therefore, we believe that the component C1–2 is resolved by GMEM in the four last epochs, making up two components A1 and A2, and this effect causes the difference in the apparent speed.

Note that the kinematics of the outer jet are consistent with the results of Britzen et al. (2009).

4.1.2. Core flux

In this section we compare the ability of the two methods to restore the flux of the source. In the optically thin jet, the fluxes of the components found by the conventional method and by the GMEM agree within two sigma. At the same time, the performance of the two methods in the core region is significantly different: the conventional method places one very bright component C1–2 (the flux density is about 0.1 Jy, so it is bright compared to other jet components) next to the core, while the GMEM method finds two components, A1 and A2, with the flux density comparable to the other components in the jet (on the order of 0.01 Jy). However, if one adds up the flux densities of the core and nearby components (C0 + C1–2 and A0 + A1 + A2), the results of the two methods would agree within two sigma for all six epochs. We attribute this difference to the fact that the GMEM restores the bright point-like core accurately, whereas the conventional method, because of significantly lower resolution at 5 GHz, is likely to mix up emission of the core and of the nearby components.

5. Discussion

5.1. Inner jet kinematics: comparison of the 22 and 5 GHz data

The inner jet of 0716+714 is well resolved at higher (> 5 GHz) frequencies and reveals a clear outward motion, although the emission fades away at about 1 mas as the spectral peak moves to lower frequencies (see Rastorgueva et al. 2009). At 5 GHz, only one or two components are resolved in the inner jet region, depending on the method, and they also have a clear outward motion. Because the data at higher frequencies were reduced by the conventional method, we discussin in this section only the component C1–2. Its apparent speed can be found in Table 3. In Fig. 4 we show the motion of the component C1–2 at 5 GHz together with components C5 and C6 at 22 GHz. They are obviously located in the same area in the inner jet and have speeds of a similar values: the apparent speeds of the components C5 (22 GHz) and C6 (22 GHz) are 21.3 ± 2.5   c and 18.9 ± 0.8   c respectively, which agrees with the speed of C1–2 (5 GHz) within 3σ. The difference in the apparent speeds of the inner jet components, detected at different frequencies, is likely caused by a resolution effect.

Fig. 3 All components obtained with the conventional and GMEM methods are shown at the same plot.

Fig. 4 22 GHz inner jet components C5 and C6 from Rastorgueva et al. (2009) and the component C1–2 at 5 GHz coincide spatially and have similar motion.

Fig. 5 Components of the inner and outer jet on the sky plane. The change of the orientation occurs at the distance of about 1 mas from the core. The components of the jet are well aligned, which is an indication of the collimated jet.

5.2. Jet bending and beaming

The jet components that were found beyond 1 mas from the core at both 5 and 22 GHz have a very high dispersion in the apparent speeds, and in most cases move slower than the inner ones (see Table 3). Therefore, we observe two kinematically different regions in the jet of 0716+714: a fast moving inner jet (< 1 mas) and a slowly moving outer jet with nearly stationary components (> 1 mas). The orientation of the inner and outer parts of the jet on the sky plane are also different: the jet changes its position angle from  ~ 25° (inner region) to  ~ 15° (outer region). The bend occurs at the distance of about 1 mas from the core, see Fig. 5. In addition to that, we found that the flux of the jet components in the inner jet is systematically higher than the flux of the outer jet components. Based on this we propose that those effects could be caused by the jet bending away from the line of sight at the apparent distance of 1 mas from the core. We tried to find out whether it is possible to explain the observed effects by the jet bending, without additional assumptions of jet parameters. If the jet bends away from the line of sight, the beaming factor δ decreases, which causes a decrease in the apparent luminosity of the jet, because it is proportional to δ^2 + α for a smooth jet (Cohen et al. 2007), where α is a spectral index (S ~ ν^α). The apparent speed of the jet also changes with the change of the angle to the line of sight, for 0716+714 it decreases. We estimated how much the flux of the jet components change from the inner jet to the outer, and compared these results with the theoretical predictions for jet bending. We assumed the jet to be a smooth collimated flow with a constant plasma speed (and therefore constant Lorentz factor Γ). This assumption is indirectly supported by the fact that the jet components are well aligned (see Fig. 5), and the structure of the jet on VLBI maps is smooth and continuous. We used the results of the conventional method for this analysis.

For each epoch we calculated the ratio of the average component flux at 5 GHz in the outer part of the jet to the average flux at 22 GHz of the inner ones: r_I =  ⟨ I_o ⟩ / ⟨ I_i ⟩ . Before calculating the flux ratio r_I for each epoch, we corrected the spectral slope of the 5 GHz fluxes in order to compare them directly with the 22 GHz data. The 5–22 GHz spectral index α was estimated for the inner jet components C1–2 (5 GHz) and C5 and C6 (22 GHz) for each epoch separately. Individual values of the flux ratios were averaged, which gave  ⟨ r_I ⟩  = 0.11 ± 0.08. Assuming that the flux change is caused by jet bending, we can calculate the ratio of the beaming factors of the jet after and before the bent using the flux ratio we have just obtained: $\mathrm{⟨}{\mathit{r}}_{\mathit{\delta }}\mathrm{⟩}\mathrm{=}{\mathit{r}}_{\mathit{I}}^{\mathrm{2}\mathrm{+}\mathit{\alpha }}\mathrm{=}\mathrm{0.24}\mathrm{\pm }\mathrm{0.12}$. These values were also calculated separately for each epoch, using the individual values of the spectral index, and then averaged. Errors of the averaged value consist of the standard deviation of averaging the ratios between epochs and the standard deviation of averaging the component fluxes for each epoch, added in quadrature.

The apparent kinematics of the relativistic blazar jet are determined by four parameters: the apparent speed β_app, the Lorentz factor Γ, the beaming factor δ and the angle of the jet to the line of sight Θ, any two of which could be used to find two others (e.g. Cohen et al. 2007). In the inner and outer parts of the jet, the Lorentz factor Γ is the same, but the other three parameters differ: Θ₂ > Θ₁, and therefore β_app2 < β_app1 and δ₂ < δ₁. From the 22 and 43 GHz data analysis (Rastorgueva et al. 2009) we estimated the kinematic parameters for the inner jet: the maximum apparent speed is β_appmax = 19.4 and the corresponding minimum Lorentz factor is Γ_min = 19.5. The corresponding minimum angle of the jet to the line of sight is Θ_mins = 3.0°. In the outer jet, the maximum apparent speed was β_app ≃ 10. Taking into account the above mentioned observed quantities and assuming that the jet has a constant speed, we can calculate its angle to the line of sight before (Θ₁) and after (Θ₂) the bending as well as other kinematic parameters of the jet.

We can consider the beaming factor of the jet as a function of the initial (inner jet) angle to the line of sight Θ₁, which is larger than the Θ_min = 3.0°, but smaller than the maximal value of ${\mathrm{\Theta }}_{\max }\mathrm{=}\arcsin \mathrm{(}\mathrm{2}{\mathit{\beta }}_{\mathrm{app}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\beta }}_{\mathrm{app}}^{\mathrm{2}}\mathrm{)}\mathrm{=}{\mathrm{5.9}}^{\mathrm{◦}}$ (all possible values lie in the lower part of a β_app − δ diagram). The ratio of the beaming factors r_δ is a monotonic function of Θ₁ and changes from 0.15 to 0.26, while Θ₁ changes within its limits, and reaches the value of 0.24 at Θ₁ = 5.2°. Therefore, the observed change of kinematics is possible if the Lorentz factor of the jet is Γ ≃ 30.0, and at the apparent distance of 1 mas from the core the jet changes its direction, with the angle to to the line of sight increasing from  ~ 5° to  ~ 11°. At the same time the beaming factor is changing from  ~ 7 to  ~ 2, respectively.

The errors of the observed flux factor ratio determination, which give us a beaming factor ratio, are fairly large, and precise calculations of the angle to the line of sight and beaming factor are not possible. The values we obtained are only estimates, and we can only say that the observed change of the jet apparent speed and brightness is possibly caused by jet bending. Also, the proposed beaming factor for the inner jet is consistent with the value found by Hovatta et al. (2009) from the radio-frequency variability.

6. Conclusions

We studied the large-scale jet kinematics of the blazar S5 0716+714 by analyzing six epochs of VLBA observations at 5 GHz. Two imaging methods were applied to this data set: a conventional method (CLEAN self-calibration loop in Difmap combined with Modelfit) and GMEM (Pulkovo VLBImager package combined with the image plane model-fitting).

6.1. Performance of the methods

S5 0716+714 has a bright point-like core and a faint diffuse jet. The emission fades down very quickly in the vicinity of the core, and the CLEAN algorithm, producing the model that is composed of the discrete point sources, is not adequately restoring this emission region. In the Gaussian model, some fraction of the core emission was added to the inner jet component, and the kinematic results were also ambiguous because of the large uncertainties of the positions of jet components.

The GMEM method delivered a smooth image of the source, which restored both the bright core and the diffuse jet. However, the jet of 0716+714 is very faint, therefore a convolution of the original solutions with the small (0.5  ×  0.5 mas) beam was needed in order to obtain accurate Gaussian models of the jet structure. The original GMEM solutions convolved with the CLEAN beam in general resemble the CLEAN images (Figs. 2 and 1).

We recommend GMEM in combination with the difference mapping technique for restoring the structure of compact AGN and other objects that have a diffuse structure with bright point-like features. In addition, the original GMEM solutions have a higher resolution than other methods, and in principle they could be used directly to derive kinematics of the source, provided that the jet features are bright enough (e.g., Bajkova & Pushkarev 2008).

6.2. Structure and kinematics of the jet

The large-scale jet of 0716+714 is faint, diffuse and smooth, without prominent brightness enhancements.

Both the conventional method and GMEM deliver similar kinematic results: the large-scale jet (1–12 mas) is mostly stationary, with a large scatter of the component positions.

The position angle of the jet on the sky plane changes by  ~ 10° at about 1 mas from the core.

Large-scale jet components are systematically slower and fainter than those of the inner (< 1 mas) jet components.

We attribute these last effects to the bending of the jet at the apparent distance  ~ 1 mas from the core from  ~ 5° to  ~ 11°.

Acknowledgments

We thank the anonymous referee for the useful comments which have helped us to revise and improve the paper.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

This study was supported in part by the “Origin and Evolution of Stars and Galaxies” – “Program of the Presidium of the Russian Academy of Sciences and the Program for State Support of Leading Scientific Schools of Russia” (NSh–3645.2010.2).

References

Appendix Appendix

Fig. A.1 Original GMEM solutions.

Fig. A.2 GMEM solutions, convolved with the beam of size 0.5 × 0.5 mas.

All Tables

All Figures

	Fig. 1 Images of 0716+714 obtained using the conventional method: CLEAN maps, convolved with the averaged beam of 1.5  ×  2.5 mas, PA =  −4.5°. Circles with crosses are the Modelfit Gaussian components superimposed on the image. The Gaussian model fitting was performed in the uv-plane. Axes are right ascension and declination in milliarcseconds, measured with regard to the core.
In the text

Fig. 2 GMEM solutions for each epoch, convolved with the averaged CLEAN beam: 1.5  ×  2.5 mas, PA  =  −4.5°. Gray ellipses are the Gaussian components superimposed on the image. The Gaussian model-fitting was performed in the image plane for the original GMEM solutions (Fig. A.1) convolved with the beam of 0.5  ×  0.5 mas (Fig. A.2). Axes are right ascension vs. declination in milliarcseconds, measured with regard to the core.

In the text

	Fig. 3 All components obtained with the conventional and GMEM methods are shown at the same plot.
In the text

	Fig. 4 22 GHz inner jet components C5 and C6 from Rastorgueva et al. (2009) and the component C1–2 at 5 GHz coincide spatially and have similar motion.
In the text

	Fig. 5 Components of the inner and outer jet on the sky plane. The change of the orientation occurs at the distance of about 1 mas from the core. The components of the jet are well aligned, which is an indication of the collimated jet.
In the text

	Fig. A.1 Original GMEM solutions.
In the text

	Fig. A.2 GMEM solutions, convolved with the beam of size 0.5 × 0.5 mas.
In the text