© ESO, 2016

1. Introduction

M 87 (Virgo A, NGC 4486, 3C 274) is a giant elliptical galaxy. Its proximity (D = 16.7 Mpc; Mei et al.2007) combined with the large mass of its central black hole (M_BH ≃ 6.1 × 10⁹ M_⊙; Gebhardt et al. 2011) make M 87 one of the primary sources to probe jet formation and acceleration down to smallest linear scales (with 1 mas ≈ 0.08 pc ~ 140 Schwarzschildradii (R_s)).

The kiloparsec-scale structure of the jet in M 87 has been extensively studied in the radio regime using the Very Large Array (VLA; Owen et al. 1989), the optical regime using the Hubble Space Telescope (HST; Biretta et al. 1995), and the X-ray regime using the Chandra space telescope (Marshall et al. 2002). These observations have revealed an edge brightened conical jet with an apparent opening angle Θ_obs ~ 3.2° (in this paper, “opening angle” refers to half the full opening angle). The jet manifests several bright knots, with a particularly strong feature HST-1 located at a projected distance of z_obs ~ 0.84″. The jet expends uniformly up to knot A at a distance z_obs ~ 12″. Analysis of multi-epoch observations has yielded detections of superluminal speeds of up to β_app ~ 6c in the optical (Biretta et al. 1999) and β_app ~ 4c in the radio (Biretta et al. 1995; Cheung et al. 2007; Giroletti et al. 2012) at the location of HST-1, with evidence of deceleration downstream from this point. Several components with subluminal speed (β_app ~ 0.5c) were also detected all along the jet (Biretta et al. 1999; Meyer et al. 2013), concurrently with the detections of superluminal features made at similar distances.

Both the radio and optical maps exhibit bright filaments, twisted into an apparent double-helical pattern. This has been successfully modeled as a consequence of Kelvin Helmholtz (K-H) instability developing in the jet (Lobanov et al. 2003; Hardee & Eilek 2011), with the resulting expected pattern speed of ≈0.5c that can explain the subluminal features found in the jet. In the optical regime, jet structure is dominated by knots, while the filaments are more pronounced in the radio images. This difference has been interpreted in the framework of a stratified jet by Perlman et al. (1999), in which the high energy plasma particles emitting optical synchrotron radiation originate from regions close to the spine of the flow, while the lower energy radio emitting particles are concentrated in the outer layers of the flow.

The physical nature of the inner jet in M 87 has been the subject of a number of studies (cf., Biretta & Junor 1995, for a review). Recently, multi-frequency phase-referencing VLBI observations (Hada et al. 2011) were used to locate the central engine at a projected distance of only about 40R_s upstream from the base (or the core) of the radio jet observed at 43 GHz. Jet expansion can be described with a parabolic profile (Asada & Nakamura 2012), indicating that the jet is collimated by magnetohydrodynamic (MHD) processes (Meier et al. 2001). The transition to a conical jet was found to occur at a projected distance of ~350 mas from the core.

Proper motions detected in the jet in M 87 at parsec scales have been the subject of an intensive debate (cf., Biretta & Junor 1995; Kovalev et al. 2007; Ly et al. 2007; Walker et al. 2009). The jet has been observed at 2 cm wavelengths as part of the MOJAVE project, with a reported slow subluminal speed (β_app ~ 0.01c; Kovalev et al. 2007) which seems to be in contradiction with proper motions of up to β_app ~ 4−6c observed at HST-1 (Biretta et al. 1999). Recently, VLBI observations at 1.6 GHz (Asada et al. 2014) have been used to claim a detection of gradual acceleration of the jet between 100 mas and 900 mas, suggesting a link between the subluminal acceleration found in the first 20 mas of the jet, and the relativistic speed at HST-1. A major issue persists in this picture: assuming that the jet is intrinsically bidirectional and symmetric, detection of a subluminal speed would imply a jet to counter-jet intensity ratio near unity, while the jet in M 87 is essentially one sided. The presence of a counter-jet is suggested from 15 GHz, 43 GHz, and 86 GHz VLBI observations (Kovalev et al. 2007; Ly et al. 2007; Hada et al. 2016). However, the counter-jet emission is only about 2 mas in extent, and jet to counter-jet ratios of ~5–20 are reported. These findings support an alternative scenario in which the ultra-slow speed of ≈0.01c is either a pattern speed or a slow wind from an outer sheath (Kovalev et al. 2007). This conclusion is further strengthened by the observed limb-brightened morphology of the jet. Sparsity of the MOJAVE time sampling is a likely reason for the non-detection of superluminal speeds in M 87. The MOJAVE observations of M 87 are sparse, with an average time interval of about 200 days between successive epochs, resulting in a maximum reliably detectable speed of about 1c (taken as a displacement of about the beam size between consecutive epochs). While Kovalev et al. (2007) could not identify any fast component, Ly et al. (2007) reported possible features with speeds between 0.25c and 0.4c at a distance of 3 mas from the core, with the detection obtained from only two epochs.

In an attempt to finally assess the true bulk speed of the jet at sub-parsec scales, systematic VLBA¹ observations were performed by Walker et al. (2008) at 43 GHz (λ = 7 mm) with a high cadence in time. The first set of pilot observations indicated that a fast moving component could indeed be traced. This led to a program that produced a batch of eleven observations between 27 January 2007 and 26 August 2007, with an average time interval of 21 days, which allowed detection of apparent speeds as high as 3–4c. The result is a set of excellent-quality, high-resolution VLBI maps with nearly homogeneous image rms noise. A more complete presentation of the 2007 data used here along with more densely sampled, but lower quality, data from 2008 plus roughly annual observations between 1990 and 2016 will be made separately (Walker et al., in prep.)

In the images of M 87 obtained from the VLBA movie project, nearly the entire jet is well resolved, with up to 3–4 beamwidths across the flow. The limb-brightened morphology is striking, especially if one looks at the stacked image that comprises the eleven observations into a single, high-sensitivity map. Determination of the velocity of the flow is however difficult from these data. The standard procedure for making proper motion measurements relies on modeling the source with a set of Gaussian components and cross-identifying those components between epochs. In this case, the complexity of the transversely resolved structure of the jet precludes application of this approach. Visual inspection was done in an attempt to identify apparently related features and a speed of about 2c was consistently found (Walker et al. 2008). This method obviously lacks objectivity and robustness, and it is difficult to use it for obtaining a complete velocity field map of the jet.

A robust analysis of the M 87 VLBA movie dataset is required in order to make a quantitative assessment of the jet morphology and kinematics. Such an assessment can be made using the wavelet-based image segmentation and evaluation (WISE) method (Mertens & Lobanov 2015) and the stacked cross-correlation algorithm (Mertens & Lobanov 2016) developed for determining two-dimensional velocity fields. In the following a complete WISE analysis of the internal structure and dynamics of the jet in M 87 at sub-parsec scale is presented. In Sect. 2, the WISE method is used to determine the two-dimensional velocity field, and the collimation profile of the jet from the 43 GHz VLBA observations of Walker et al. (2008). The complementary WISE analysis of the 15 GHz VLBA observations from the MOJAVE project is presented in Sect. 3. The physical properties of the jet are discussed in Sect. 4. In Sect. 4.1, we present a physical model for the flow stratification observed in M 87. In Sect. 4.2, the viewing angle is derived from the velocity observed on the counter-jet side. In Sect. 4.3, evidence for jet rotation is discussed, as inferred from the differential velocity between the northern and southern limbs of the jet. Based on these findings, the acceleration and collimation of the jet are discussed in Sect. 4.4 in the framework of MHD jet acceleration. The physical nature of the jet spine is discussed in Sect. 4.5, and a plausible jet launching mechanism is proposed in Sect. 4.6.

Throughout the paper, a cylindrical coordinate system (z, r, φ) is adopted for the jet description. In this coordinate system, the coordinates z and r correspond to axial and radial directions in the jet. We set the jet axis (z) at a position angle of 20° with respect to the RA axis, and we denote the projected axis as z_obs. The black hole mass estimates reported in the literature for M 87 are typically in the range of ~3–7 billion solar masses (cf. Walsh et al. 2013, and references therein). We adopt the value of M_BH = 6.1 × 10⁹ M_⊙ in this paper.

2. WISE analysis of 7 mm VLBA observations

We applied WISE analysis (Mertens & Lobanov 2015) to 11 images of the jet in M 87 obtained from 43 GHz VLBA observations made between 27 January 2007 and 26 August 2007, with an average time interval of 21 days between the successive epochs and a common restoring beam of 0.43 × 0.21 mas in PA −16°. We refer to Walker et al. (2008; in prep.) for a complete description of the data reduction. At the cadence of the 43 GHz VLBA measurements, a displacement comparable to one fifth of the beam corresponds to a minimum detectable proper motion of 0.74 mas/yr, or 0.19c. The WISE algorithm has been shown to provide robust detections of up to ≈4 FWHM between pairs of successive epochs (Mertens & Lobanov 2015), which translates into a maximum apparent speed of ~5c that can be detected from the VLBA images of M 87. Hence, the WISE analysis should be able to detect essentially the full range of speeds reported in the compact jet of M 87.

2.1. Velocity field analysis

Analysis of the M 87 images is performed as described in Mertens & Lobanov (2015). Each map is decomposed and segmented using the segmented wavelet decomposition (SWD) method, which provides a description of the two-dimensional jet structure based on a set of significant structural patterns (SSP). These features are then cross-identified in pairs of adjacent epochs using the multiscale cross-correlation (MCC) method. This procedure yields spatial displacements and proper motions for each of the SSPs identified.

To ensure robust detectability of displacements larger than one beamwidth (corresponding to apparent speeds of ≳2c), we applied the SWD on four spatial scales (ranged as powers of two) of 0.07, 0.14, 0.28, and 0.56 mas and amended them with the intermediate wavelet decomposition (IWD) performed on scales of 0.105, 0.21, 0.42, and 0.84 mas. The IWD implementation, described in Appendix A, is made to improve the efficiency of the MCC algorithm for recovering structural changes in optically thin, stratified flows. With these settings, maximum displacements of about 0.84 mas (corresponding to proper motion of ~18 mas/yr and β_app ~ 4.5c) can be robustly detected in each successive pair of images even for overlapping optically thin regions with different apparent speeds. The noise was estimated by computing σ_j at each wavelet scale, as described in Mertens & Lobanov (2015) and a 3 σ_j thresholding was subsequently applied at each scale. The typical S/N found in different parts of the jet and for different Scales is summarized in Table 1.

Figure 1 shows an example of the SWD/IWD decomposition of the VLBA image of M 87 from the epoch 2007/02/17. The internal structure of the flow is well sampled at Scales 1 and 2 of the SWD, while the SWD Scales 3 and 4 are essential for providing an effective reference frameset for the successful determination of the SSP displacement made using the MCC algorithm.

The MCC algorithm has been applied to each consecutive pair of the jet images. Based on the previously reported measurements made at these distances in the jet (Reid et al. 1989; Biretta et al. 1995; Dodson et al. 2006; Kovalev et al. 2007), one can expect the jet to manifest a mixture of large and small displacements. To ensure robust detectability of the entire range of potential displacements, we use a tolerance factor of 1.5 and a correlation threshold of 0.6 in our MCC analysis. Finally, velocity constraints are set according to the expected largest displacement ranges: on the jet side, the longitudinal velocities are allowed to range between −2 mas/yr and +18 mas/yr, transversally they are allowed to range between −5 mas/yr and +5 mas/yr. On the counter-jet side, velocities ranging between −2 mas/yr and 2 mas/yr longitudinally and transversally are allowed. For the analysis of the fine two-dimensional velocity structure of the jet, only the displacements determined at the lowest SWD scale (Scale 1: 0.07 mas) are considered. The displacements determined at coarser scales were used in the MCC and overall are consistent with the ones obtained at the 0.07 mas scale.

The SWD yields a total of 68 SSPs that are detected on average in each single image of the jet. At the smallest spatial scale of the SWD, SSP matching is successfully made for 51% of the SSPs, resulting in detection of a total of 391 displacement vectors. The normalized cross correlation between matching sets of SSPs is larger than 0.8, on average. The errors of the displacements and velocity vectors are obtained from the uncertainties of the SSP positions which are computed using the stacked cross-correlation of the SSP (Mertens & Lobanov 2016), which yields: (1)where b_x and b_y describe the beam size along the x and y coordinates, respectively, and S/N is the signal-to-noise ratio of the SSP detection.

2.2. SSP displacements and jet kinematics

The SSP displacements obtained using this approach are shown in Fig. 2, and the transverse and longitudinal components of the corresponding velocities are plotted in Fig. 3.

Figures 2, 3 indicate that the jet exhibits three kinematically distinct regions of organized stream-like features with predominantly longitudinal displacement vectors. These features will hereafter be referred to as the northern and the southern limbs and the central stream (we prefer not to use the term “spine” for this feature as it is not clear whether it is indeed related to the physical relativistic spine of the stratified outflow in M 87).

The kinematics is more complex in the central stream region which manifests both outward motions and substantial transverse displacements. The displacements can be reliably detected in the central stream at z_obs ≳ 1.75 mas, when the width of the jet exceeds three beam sizes.

The top panel of Fig. 3 shows the presence of both slow (subluminal) and fast (superluminal) speeds in the central stream and each of the limbs. This spread of apparent speeds might indicate jet stratification. In each of the three regions, we also detect several outwards motions which are most of the time compatible with no motion (taking into account position uncertainty) but might also be related to instability pattern motion.

An apparent acceleration can be seen in both limbs over the innermost 2 mas of the jet. At a distance z_obs ~ 0.25 mas, the maximum apparent velocity is β_app ~ 0.2c, which then linearly increases up to β_app ~ 2.5c at z_obs ~ 2 mas. Between z_obs ~ 2 mas and z_obs ~ 6 mas, we measure either a slow acceleration or no acceleration at all. In this region, a maximum velocity of β_app ~ 3c is found at z_obs ~ 5 mas. We note that because this acceleration is observed at the limbs of the jet, the measured apparent velocities are only marginally affected by the jet opening angle (see Eq. (6)), and a fast decrease in opening angle between 0 and 2 mas would not be able to explain the observed increase in apparent velocity as this decrease would strongly affect only the velocities measured in the central stream, which is not detected at distances smaller than ≈2 mas.

We also observe an oscillation of the speed measured along the jet, with minima in the jet speed found at z_obs ~ 1.2 mas, 2.25 mas and 4.25 mas.

On the counter-jet side, there is only one SSP component that shows consistently outward motion over five subsequent epochs (orange arrows in Fig. 2). At this distance from the core, the counter-jet also appears limb-brightened, and the apparent position of this feature is consistent with it being located in the northern limb. We can use it’s speed and brightness relative to the main jet to obtain estimates of the jet angle to the line of sight and intrinsic speed.

Several other features in the outer part of the counter-jet exhibit slower motions tending toward the northward direction, which may be viewed as a hint for the clockwise rotation, as viewed by the observer, of the emitting material. The robustness of detection of the displacements of these features and potential physical interpretation of these displacements should be a matter of further, more detailed investigation.

In order to assess the robustness of the SSP and displacement identification, we perform a bootstrapping analysis by randomly shuffling the individual images and carrying out the WISE analysis on a number of such random sets. We found that on average only 35% of the SSPs can be matched in the pairs of shuffled images and only 13% of those can be traced over five or more epochs, compared to 24% for the non-shuffled set of images. The bootstrapping datasets also do not produce regular velocity patterns and do not show the correlation between z_obs and β_long observed in the original sequence of images. This demonstrates that the solution found by WISE does not result from a spurious correlation.

2.3. Flow stratification analysis

The large variations of velocity measured along and across the jet indicate the complex physical nature of the flow, with both acceleration and stratification of the flow likely to play important roles. To investigate the potential effect of flow stratification on the observed velocity field we combine the WISE velocity field analysis with the stacked cross correlation (SCC) analysis (Mertens & Lobanov 2016). The SCC method combines the cross-correlation results for all SSPs identified at different scales of the SWD in all different pairs of epochs. This procedure yields a two-dimensional distribution of the cumulative correlation coefficients, with peaks corresponding to different velocity components. The significance and uncertainty of each of these peaks is evaluated using Monte Carlo simulations.

For the present analysis, the SSPs detected at the three lowest SWD scales are used. The jet is divided in two different regions (region “A” at 0.5–1.0 mas distance from the core and region “B” at 1–4 mas distance from the core), and the SCC analysis is performed independently for each of these regions. The extent of these two regions is constrained by the requirement for the detection of statistically significant velocity components in the SCC analysis of both the limbs and the central stream of the flow.

2.3.1. Region “A”

In the northern part of the jet, the SCC is computed from a total of 42 SSPs. The resulting velocity cross correlation maps and the corresponding significance estimates are plotted in Fig. 4. For the northern limb of region “A”, the SCC results imply a single velocity component, with an apparent speed of and a significance of 5.9σ. In the southern limb, the SCC has been applied to a total of 47 SSPs, also revealing a single velocity component with and a significance of 8.7σ. The results are further summarized in Table 2.

The SCC analysis did not unveil any stratification close to the core. This result confirms however the speeds measured in Sect. 2.1, further strengthening the conclusion that subluminal speeds are generally found in this region. If any stratification would be present, the respective velocity difference must have been too small so that it would remain undetected with the SCC. In this case, the apparent velocity measured by this method can be interpreted as a weighted mean velocity of the different jet layers along a given line of sight. This may also explain the difference between the speeds measured in the northern and southern limb (an alternative explanation may involve rotation of the flow, which will be examined later).

2.3.2. Region “B”

For region “B” (1–4 mas from the core), the SCC analysis was performed independently for each of the limbs and for the central stream of the jet, using 245, 267, and 217 SSPs, respectively. The resulting cross-correlation maps are plotted in Fig. 5. The measured correlations reveal the presence of two distinct velocity components, which suggest a strong stratification of the flow with a slow, mildly relativistic speed (β_app ~ 0.35c) and a faster, relativistic speed (β_app ~ 2.3c). The velocity parameters derived for these components are summarized in Table 3.

The statistical significance of most of the velocity components found in region “B” is better than ≈4σ. The overall consistency of the speeds found in the limbs and the central stream lends further support for stratification of the flow. A somewhat low (2.5σ) formal significance of the correlation is found for the fast velocity component in the southern limb of the jet. This may result from a broader range of measured displacement velocities registered there as compared to the northern limb.

2.4. Jet collimation

In the magnetically accelerated jet model, the conversion from Poynting flux into kinetic energy depends on the configuration of the magnetic field lines (Begelman & Li 1994; Lyubarsky 2010). The jet collimation profile is therefore an important indicator of the jet acceleration mechanism. The segmented wavelet decomposition (SWD) method introduced in Mertens & Lobanov (2015) provides a reliable way to analyze the evolution of the jet width (jet radius) with distance from the core.

To determine the jet collimation profile, three stacked-epoch wavelet scale (SWS) images of M 87 are prepared. Each of these images comprises all SSPs detected at a given scale j = 1,2,3 at all 11 epochs. Near to the jet base, the jet limb is well described by the SWS at the smallest wavelet scale. As the jet expands, the width of each of its limbs is best described by progressively larger SWS scales. Scale 1 (0.07 mas) of the SWS is used between z_obs = 0.4 mas and z_obs = 0.8 mas, scale 2 (0.14 mas) is then used up to z_obs = 1.6 mas, and finally scale 3 (0.21 mas) is selected at larger core separations.

The jet radius is estimated by taking transverse profiles of the stacked-epoch SWS images of the jet at regular intervals (δz = 1 px = 0.035 mas) and measuring the distance between the two limbs of the jet. A Gaussian profile is fit to each of the limbs to further improve this calculation. It is particularly difficult to obtain a robust uncertainty estimate for these measurements. Our conservative approach is to assume that the position error of limb determination is 1/4 of the beam (0.07 mas) at the core of the jet, and linearly increases to reach 1/2 of the beam at the last profile, located at 6.6 mas from the core. The result is plotted in Fig. 6. The instantaneous apparent jet opening angle, Θ_app, is computed from the measured jet radii using the relation Θ_app = arctan(dr/ dz_obs). The local curvature dr/ dz_obs is approximated from a linear fit of the jet radius 0.2 mas before and after each point of the calculation.

The jet collimation profile is well described by a power law r ∝ z^k with k = 0.60 ± 0.02 (reduced χ² = 0.55). This result is consistent with previous measurements made by Asada & Nakamura (2012; k = 0.58 ± 0.02) and Hada et al. (2013; k = 0.56 ± 0.03). As seen from Fig. 6, expansion of the jet reveals an oscillatory pattern, most likely reflecting repeated over-collimation and over-expansion of the flow. The over-collimation is visible at z_obs ~ 1, 1.9 and 4 mas, and the over-expansion manifests itself at z_obs ~ 0.75, 1.4 and 2.6 mas, suggesting that the spatial period of these oscillations increases with the distance as ~z_obs. Oscillations of the jet width seen at distances ≳5 mas have a shorter period of ≈0.5 mas.

It is interesting to note that oscillations of the jet width correlate with the changes in the apparent jet speed reported in Sect. 2.2. The locations of maximum over-collimation of the jet correspond well to the observed minima in the apparent speed of the jet. If these oscillations reflect the hydrodynamics of the flow, the pressure ratio between the jet and the external medium might be close to unity (Daly & Marscher 1988). Presence of oscillating patterns in an expanding flow is also expected for an expanding magnetized jet (Lyubarsky 2009; Komissarov et al. 2015). Alternatively, these oscillations might reflect the evolution of the pinch mode or the elliptical surface mode of K-H instability (Hardee 2000; Lobanov & Zensus 2001).

3. WISE analysis of 2 cm VLBA observations

Between 1995 and 2010, the jet in M 87 was regularly observed as part of the MOJAVE project (Lister et al. 2009, 2013). The MOJAVE database contains 30 images of this jet, providing, on average, one observation every seven months. With this cadence, WISE can reliably detect jet speeds of up to β_app ~ 2c, assuming a maximum detectable displacement of four times the beam size. In an analysis of this data using the Gaussian model-fitting technique, Kovalev et al. (2007) reported only a slow, subluminal motion with a maximum speed β_app ~ 0.05c. Finding this result to contradict several physical aspects of the jet, including the counter-jet flux density ratio, Kovalev et al. (2007) suggested that the measured subluminal speed might represent pattern motions of either shocks or MHD instabilities. Using the WISE analysis, we can attempt to detect a faster speed from the MOJAVE data. The application of the SWD made for this purpose is performed using four scales, ranging from 0.4 mas (scale 1) to 3.2 mas (scale 4).

The interval between two epochs is generally too long to warrant detection of the fast velocity β_app ~ 2.5c that has been measured in the 43 GHz VLBA data. We concentrate therefore on the pair of observations separated by the shortest time interval. The complete MOJAVE database has two epochs (2000/04/07 and 2000/06/27) that matched this criteria, and we complement them with two full track 15 GHz VLBA observations obtained at epochs 2000/01/22 and 2000/05/08 (Kovalev et al. 2007).

The result of applying the WISE algorithm to these four VLBA images at 15 GHz is plotted in Fig. 7. We detected fast relativistic speeds of β_app ~ 2–3c for six different SSPs in both the northern and southern limbs of the jet at z_obs ~ 6.5 mas and z_obs ~ 16 mas. In addition to this, eight SSPs exhibit apparent superluminal speeds of β_app ~ 1.5c, and subluminal speeds (β_app ~ 0.5c).

The displacements traced in the 15 GHz VLBA images do not provide enough material for tracing the velocity field in the jet. However, one can see that full-tracks VLBA observations at 15 GHz made at ~2 month intervals should be sufficient to obtain a reliable velocity field.

4. Discussion

4.1. Jet internal structure

The stacked cross correlation analysis discussed in Sect. 2.3 reveals a clear stratification of the flow, with a slow, subluminal component and a faster relativistic component detected in all three streams inside the jet. Both the slow and the fast velocity components exhibit similar speeds in each of the streams, which suggests that each of them reflects plasma motions occurring at similar radial separations (similar jet layers) in a transversally stratified flow. Since the fast velocity component is present in all three streams, it can be associated with the sheath, rather than the spine, of the jet. A faster spine is either not present in the flow (at least at the linear scales studied here), or it is deboosted and appears dimmer than the slower moving material of the sheath (in this case, VLBI observations at time intervals of smaller than 10 days are required to successfully disentangle this faster component from the two velocity components identified in this paper). In further analysis, we assume the latter possibility, which corresponds to the schematic view of the jet shown in Fig. 8.

The observed jet stratification, with a fast, accelerating sheath and a subluminal layer can be interpreted in several different ways. It may reflect a “two-fluid” jet (Sol et al. 1989) comprising a relativistic beam (or spine) originating either from the magnetosphere of the central engine, or from the inner part of the accretion disk, and a non-relativistic wind originating from outer parts of the accretion disk (Tsinganos & Bogovalov 2002). The viability of this scenario has been tested in RMHD simulations of AGN jets (e.g., Komissarov et al. 2007) which demonstrated that magnetic field lines originating from different parts of the accretion disk would result in substantial transverse stratification of the flow.

Alternatively, the fast sheath component may reflect the true speed of the flow (without substantial transverse stratification), while the slower one would then be interpreted as a pattern speed associated with current driven (CD) or Kelvin-Helmholtz (K-H) instability. In general, K-H instability pattern speed is expected to be considerably slower than the flow speed (Hardee 2000). It was also suggested that the subluminal apparent speeds of ~0.5c detected in the jet in M 87 at kiloparsec scale indicate a pattern speed from K-H instability (Lobanov et al. 2003).

4.2. Counter-jet and jet viewing angle

The presence of a counter-jet in M 87 jet has been suggested based on the structure detected in 15 GHz VLBA maps (Kovalev et al. 2007). A similar structural detail resembling a counter-jet is also observed in the 43 GHz images of M 87 (Ly et al. 2007; Walker et al. 2008). The WISE analysis performed on the 43 GHz VLBA data has revealed several SSPs located in the counter-jet structure and moving outward with speeds similar to the subluminal speed measured on the jet side at the same distances from the core. We will use this velocity measurement to constrain the viewing angle of the jet.

Relativistic Doppler beaming affects differently the flux density in the jet and the counter-jet. If the jet is seen at a viewing angle θ, the observed intensity ratio between the jet and counter-jet is given by: (2)with β_j, β_cj are the intrinsic speeds of the plasma in the jet and counter-jet respectively, γ_j, γ_cj are the corresponding Lorentz factors and α is the spectral index of the synchrotron radiation, I_ν ∝ ν^α. In a recent publication, α = −1 was estimated for the radio emission on sub-parsec scales in the M 87 jet (Hovatta et al. 2014).

We measure the intensity ratio ℛ by taking a longitudinal profile of emission brightness along the jet limbs. In order to reference properly the jet and counter-jet locations, special attention must be paid to the identification of the position of the core. Because of synchrotron self-absorption, the observed location of the jet core does not correspond exactly to the true base of the jet (Lobanov 1998). However, in the case of M 87, Hada et al. (2011) found that the central engine of M 87 and the radio core at 43 GHz are separated by a projected distance of only 41 ± 12 μas, which is small enough to consider, in our case, the core as the point of symmetry of the jet. The counter-jet component which yields a robust speed estimate is located in the northern limb of the counter jet. To preserve the symmetry, the jet to counter-jet intensity ratio is measured between the northern limb of the counter jet and the southern limb of the jet. To obtain the profile, flux density is measured at every pixel (0.035 mas; 1/8 of the beam size). The result is shown in Fig. 9. We note here again that both the intensity ratio and velocities are measured at the limbs of the jet, and therefore the effect of variable jet opening angle near the core would only be marginal and well below the statistical uncertainties of our measurements.

In order to avoid biasing by the bright core of the jet we exclude the inner 0.4 mas portion of the flow (corresponding to ≈1.4 FWHM of the restoring beam). At separations between 0.4 mas and 0.8 mas, we observe a sharp increase of the jet to counter-jet intensity ratio from ~7 to ~11 which may result from acceleration of the flow. In this region, robust displacement detections have been made in three components on the jet side (with speeds ranging from 0.16c to 0.24c) and in one feature on the counter-jet side with an apparent speed of 0.14c. These components are presented in Fig. 10.

The combination of the observed jet to counter-jet intensity ratio and the maximum and minimum measured apparent speeds provides conservative constraints on the viewing angle of the jet. This is illustrated in Fig. 11 which shows the ranges of the viewing angle and the intrinsic speed that can reconcile both the apparent speed and the measured jet to counter-jet intensity ratio:(3)This estimate assumes that both the jet and the counter-jet have the same intrinsic speed and brightness. In a stratified flow, it is however very likely that different layers are dominating the emission in the jet and the counter-jet sides. To take this into account, we consider the difference between the apparent speed found in the jet and the speed measured in the counter jet. The counter-jet component C1 is detected at distances of 0.5–0.6 mas from the core (see Fig. 10). In this region, we measure a jet to counter-jet intensity ratio ℛ = 9.5 ± 1.5. On the jet side, we can use a more robust velocity measurement obtained from the stacked cross correlation analysis discussed in Sect. 2.3. Between 0.5 and 1 mas from the core, we measure a jet apparent speed in the southern limb of the jet. Then, using Eq. (2), we obtain: (4)This estimate of the viewing angle is consistent with the constraint imposed by the maximum apparent speed β_app ~ 6c measured at HST-1 which requires a viewing angle θ ≲ 19°.

4.3. Jet rotation

4.3.1. Evidence and estimation of the jet rotation

The stacked cross correlation analysis has revealed substantial difference between the speeds measured in the northern and southern limbs of the jet (see Tables 2 and 3). At core separations of 0.5–1.0 mas, the significance of this difference reaches ~4σ, while it is above ~1.5σ for both the fast and slow velocity components detected between 1 and 4 mas.

This difference can be explained naturally by invoking jet rotation. The intrinsic speed, β, of the flow can be decomposed into poloidal (β_p = βcosα_rot) and toroidal (or azimuthal, β_φ = β sinα_rot) components, where α_rot is the angle between the intrinsic velocity and the jet axis. In this framework, for a clockwise rotation of the jet as viewed by the observer, the SSP components moving along the northern limb of the jet would be viewed at a larger viewing angle, increased by α_rot, while in the southern limb, the effective viewing angle would be reduced by α_rot. This effect would not apply to the central stream of the jet. However, since the velocities measured in the central stream can originate either from the near-side or the far-side of the jet, the effective viewing angle is increased by an angle ψ, which is zero when far-side and near-side of the jet contribute equally to the velocity detected in the central stream, and is bounded by the intrinsic instantaneous opening angle Θ of the jet. For the velocity vectors measured in the limbs, one need also to take into account the effect of the opening angle of the jet. The resulting geometrical setting is illustrated in Fig. 12.

Following this scheme, the effective viewing angles for the three streams in the jet can then be written as: (5)where θ_edge is the edge viewing angle corrected for the opening angle of the jet: (6)The intrinsic opening angle Θ can be estimated from the apparent opening angle measured in Fig. 6, using Θ = Θ_appsinθ.

Assuming that the flow has the same intrinsic velocity β in the three regions, the rotation angle α_rot can be determined by solving the system of two equations: (7)If the velocity in the central stream is also considered, an additional equation can be written: (8)A unique solution for Eqs. (7), (8) can be found with respect to the parameters α_rot, β, and ψ if the viewing angle of the jet, θ, is known. Based on the analysis of the counter jet feature (see Sect. 4.2), a viewing angle θ = 18° can be adopted for this purpose. In the subsequent analysis, we will also test independently this hypothesis considering an equal contribution of the far-side and near-side of the jet in the central stream velocity (ψ = 0).

We use a least-square fitting technique to solve Eqs. (7), (8). Robust estimates of the fit uncertainties are obtained from a Monte Carlo simulation with 10 000 trials using the input parameters drawn from Gaussian distributions with a mean and standard deviation provided by the fit value and the formal fit uncertainty for a given fit parameter. The resulting output distribution of the fit values also follows a Gaussian distribution, for each parameter of the fit. We therefore determine the uncertainty of the fit parameters by computing the standard deviation of this distribution.

The rotation parameters are calculated using the apparent velocities estimated from the SCC analysis in region A (0.5 mas ≤ z_obs ≤ 1 mas), and independently for the slow and fast component identified in region B (1 mas ≤ z_obs ≤ 3 mas). The results are summarized in Table 4, assuming a viewing angle θ = 18°. Evidence for a clockwise rotation of the jet is found in region A with a ~4σ significance, and for the fast component of region B. In this last case, we verify that ψ is compatible with the inequality of Eq. (5): ψ = −0.2 ± 5° while Θ = 2.7° at z_obs ~ 2.5 mas. Additionally we can obtain an independent estimate of the viewing angle of the jet θ assuming ψ = 0: (9)which supports very well our initial hypothesis. On the contrary, the slow component of region B is consistent with a counterclockwise rotation, which may indicate that this kinematic component reflects pattern motion rather than true physical rotation of the flow.

Plasma instability can produce pattern motions with a substantial azimuthal velocity component which can both co-rotate and counter-rotate with the flow rotation (Bodo et al. 1996; Porth 2013). This might be the case for the slow component detected in region B, which has already been associated with an instability pattern (see Sect. 4.1). A counterclockwise rotation of the instability pattern is also inferred for the jet in M 87 from the analysis of the internal structure of the flow on kiloparsec scales (Lobanov et al. 2003). We adopt therefore this explanation in our subsequent analysis of the physical conditions in the jet.

Jet rotation (irrespectively of its origin) can also explain the different broadening of the SCC peaks observed in Fig. 5 for the fast velocity component detected in Region B. The observed gradual increase in the peak width going from the northern limb to the southern limb agrees well with the hypothesis of a different effective viewing angle in the three parts of the jet. For the same intrinsic speed in all three streams of the jet, Fig. 13 demonstrates that a smaller effective viewing angle (southern limb), will result in a larger range of apparent velocities, while on the contrary, a larger effective viewing angle (northern limb) will result in a smaller range of apparent velocities.

The smaller effective viewing angle to the flow in the southern limb relative to the northern limb resulting from rotation would result in differential Doppler boosting. Instabilities and the potential slower outer layer make it difficult to accurately measure jet rotation from the brightness ratio between the two limbs. However, the mean brightness ratio between the northern and southern limbs is 0.56, and close to the 0.4 value predicted to accompany our estimated fast layer jet rotation.

4.3.2. Interpretation of the jet rotation

If the azimuthal components of the jet velocity determined above reflects physical rotation of the jet, they can be related to specific mechanisms of formation and acceleration of relativistic flows in AGN. In the context of an MHD jet launched by magneto-centrifugal forces, the initial toroidal velocity is associated either with the Keplerian speed at the launching location in the accretion disk (Blandford & Payne 1982, BP mechanism) or with the spin of the central engine (Blandford & Znajek 1977, BZ mechanism). Conservation of angular momentum will result in a slow decrease in the toroidal velocity beyond the light cylinder (Vlahakis 2015). Observational signature of jet rotation has been previously reported in jets from young stellar objects (YSO), based on spectroscopic radial velocity analysis (e.g., Anderson et al. 2003; Choi et al. 2011), and has been successfully associated with an MHD process, giving a constraint on the jet launching location. Jet rotation has so far not been directly detected in AGN jets, with the reported rotational speeds clearly related to the rotation of patterns inside the flow (e.g., Lobanov et al. 2003; Agudo et al. 2007, 2012; Cohen et al. 2015) resulting from plasma instability in the flow. If our measured rotation is associated with conserved angular momentum in the jet and if the jet is launched from the accretion disk via the BP mechanism, its direction should be the same as the rotational direction of the disk itself. Using spectroscopic HST measurements, Harms et al. (1994) found that the ionized gas associated with the accretion disk in M 87 rotates in a clockwise direction. This is the observed rotation direction found in region A and for the fast component in region B. It is also interesting to note that the toroidal velocity measured closer to the central engine is greater than that measured farther away, which is what the conservation of angular momentum predicts. We will now try to associate this observed rotation with the MHD properties of the flow, following Anderson et al. (2003).

In the approximation of a cold steady MHD jet, the equations of conservation of total (magnetic + matter) specific energy μ and angular momentum L can be combined, eliminating the toroidal magnetic field B_φ and the distribution of mass flux at the inlet of the flow η, and yielding the jet specific momentum J:(10)The jet specific momentum, J, is conserved along any particular field line and can be used to obtain an expression for the angular velocity, which is also a conserved quantity along a given field line: (11)At the origin of the outflow, we can assume that β_φ is well below the speed of light (Ωr₀ ≪ c, β_φ,in ≪ 1), which gives J ≃ γ_inc² ≃ c². Knowing the toroidal velocity and Lorentz factor of the flow, Eq. (11) can be used to estimate the angular velocity and the associated light cylinder radius (r_lc = c/ Ω). Assuming Keplerian rotation in the accretion disk, we can also obtain the launching location of a field line: (12)with M_BH denoting the mass of the central engine.

The resulting parameters obtained for region A and for the fast component in region B are summarized in Table 5 for a viewing angle of 18°. For completeness, the respective rotation parameters and angular velocity estimates for viewing angles between 15° and 21° are also provided in Fig. 14. Continuity of the flow between regions A and B, which can also be inferred from the VLBA images (see Fig. 2), indicates that both these regions can represent the same field line, and hence should yield similar underlying rotational properties of the flow. Indeed, the estimates listed in Table 5 for the two jet regions agree within the uncertainties, and a weighted average of the two measurements gives (13)Several potential difficulties in invoking jet rotation in the interpretation of some observed properties in YSO outflows have been discussed in the literature. Numerical simulations have shown that several disruptive effects, for instance shocks (Fendt 2011) or plasma instability (Staff et al. 2015) can indeed affect the rotational signature from YSO jets. These effects should be mitigated in our case, as our rotation measurements also include a fast component velocity which should be related to the inner layers of the jet and thus less affected by KH instability operating near the jet surface (however a body mode of the KH instability or current-driven instability could still have an impact on the observed jet rotation). We remark also that no re-collimation shocks are observed in the VLBA images of the M 87 jet and that the jet continues to be straight and steady up to Knot A at a de-projected distance of several kiloparsecs. We finally note that the location of the jet launching region that we derived from the angular velocity measurement is in the inner part of the accretion disk close to the innermost stable circular orbit (ISCO) of the black hole, which agrees well with the predication from numerical simulations of magnetized relativistic jets (Meier et al. 2001; Komissarov et al. 2007).

4.4. Jet collimation and acceleration

Analytical and numerical simulation MHD work has shown that the processes of jet collimation and acceleration are related and occur on the same spatial scales in relativistic flows. Most of the flow acceleration might result from the so-called “magnetic nozzle” effect (Li et al. 1992; Vlahakis 2015). This effect is at work whenever the flow collimates so that the distance between streamlines increases faster than the radius of the flow. This results in a small deviation of the local poloidal magnetic field from the mean magnetic field and this drives the acceleration. Since this deviation is small, the fluid is only slowly accelerated by a small residual force and the acceleration zone can extend over a large distance. Both numerical simulations and observational evidence suggest that jet acceleration continues up to distances of 10⁴–10⁵R_s from the central engine (Vlahakis & Königl 2004; Lobanov & Zensus 1996, 1999; Homan et al. 2015; Lee et al. 2008, 2016).

The ideal MHD equations describing plasma jets may be written as a set of two non-linear differential equations describing the equilibrium of forces perpendicular (Grad-Shafranov or transfield equation) and parallel (Bernoulli equation) to the magnetic surface (Camenzind 1986; Fendt & Ouyed 2004). A complete analytic solution of these equations can be obtained only in specific, restricted cases, (e.g., self-similar solutions). The force-free approximation is often used to describe jets in the near zone, inside the light cylinder (r ≲ r_lc; e.g., Blandford & Znajek 1977). At larger radial separations, different sets of approximations may need to be used.

In a recent paper, Lyubarsky (2009) obtained asymptotic solutions for the case of a cold, ideal, Poynting flux dominated (σ ≫ 1) MHD jet in the far-zone (r ≫ r_lc). In these solutions, different sets of approximations are used, depending on the decay of the external pressure confining the jet. In AGN, this confinement can be achieved through a variety of forces including gas pressure from the ambient medium, ram pressure from the outer, subrelativistic wind and the stress of a magnetic field anchored in the disk. In general, evolution of the confining pressure can be approximated by a power law: (14)Lyubarsky (2009) found two different collimation and acceleration regimes for the jet, depending on specific ranges of the power index κ.

If κ< 2, the jet is in the so-called equilibrium regime. The residual of the hoop stress and the electric force are counterbalanced by the pressure of the poloidal magnetic field so that, at any distance from the source, the structure of the flow is the same as the structure of an appropriate cylindrical equilibrium configuration. The jet shape is given by a power law r ∝ z^k, with k = κ/ 4 converging to a cylindrical shape at large distance. The acceleration is linear, meaning that it is proportional to the jet radius: γ ∝ r ∝ z^k.

If κ > 2, the jet is at first in a regime similar to the equilibrium regime described above. The jet shape is parabolic r ∝ z^k and the acceleration is linear. At some distance, the pressure of the poloidal magnetic field becomes small, and the flow could be viewed as a composition of coaxial magnetic loops. This regime is called non-equilibrium. The flow reaches ultimately a conical shape and the acceleration is determined by the curvature of the flow, γ ∝ z^(κ−2)/2.

We found in Sect. 2.4 that the evolution of the cylindrical cross-section of the jet in M 87 is described by a power law with an index k ~ 0.56–0.60. Asada & Nakamura (2012) found that the M 87 jet reaches a conical shape at kiloparsec scales, similarly to the κ > 2 case described above. To determine the acceleration profile, we compute the Lorentz factors for the velocity field obtained from WISE analysis at 7 mm (Sect. 2.1) and 2 cm (Sect. 3), using a viewing angle θ = 18°. We complement the radio measurements with the speeds determined from optical and radio images of the kiloparsec-scale jet in M 87 (Biretta et al. 1999; Cheung et al. 2007; Giroletti et al. 2012; Asada et al. 2014). The combined evolution of the Lorentz factor plotted in Fig. 15 indicates two different regimes for the acceleration. The acceleration is at first linear with γ ∝ r ∝ z^0.58 up to z ~ 8 mas (~10³R_s), and then it becomes weaker, with γ ∝ z^0.16, agreeing well with the measurements obtained at larger scales, up to the knot HST-1. The jet appears to enter the non-equilibrium regime (κ > 2) on these scales. One must however be cautious about this interpretation. The scaling relations obtained for the κ > 2 case are obtained for a conical jet. However the jet in M 87 transits to a conical shape only at z ~ 10³ mas (Asada & Nakamura 2012). An alternative explanation for this slower acceleration can be offered by an early saturation of Poynting flux conversion. In this case, the conversion should be nearly complete at z ~ 10³R_s, and acceleration would then continue slowly until full conversion is achieved.

In the following, the acceleration and collimation of the jet in M 87 is tested against the magnetically accelerated jet model. Two different methods are used to solve the MHD equations. The first method is based on the analytical approximations obtained by Lyubarsky (2009) which describes Poynting flux dominated flows and is mainly used to test the initial acceleration phase. A second method, that relaxes the assumption of a Poynting dominated jet, is used to investigate the second acceleration phase, the saturation of Poynting flux conversion and the transition to a kinetically dominated jet.

The velocity profile corresponding to the fast velocity component of the flow, which has been associated with the sheath of the jet in Sect. 4.1, is obtained by computing the mean of the fastest 10% of the speeds measured within individual bins of 0.1 mas in size. To improve the robustness of the velocity determination, resampling is done using bins of 0.2 mas for 0.4 ≤ z_obs< 2 mas, 0.3 mas for 2 ≤ z_obs< 4 mas and finally 0.6 mas for z_obs ≥ 4 mas. The result is plotted in Fig. 16. We ascribe these measured velocities to a single field line (located on a single magnetic flux surface) with a profile corresponding to the collimation profile of the sheath as determined in Sect. 2.4. Under this assumption, and assuming Ω = const along this field line, we can conveniently use variables scaled to the light cylinder radius, R = Ω r/c and Z = Ω z/c.

4.4.1. Asymptotic relations in the far-zone Poynting dominated approximation

In the ideal cold MHD approximation, one can further simplify the solutions by considering the additional approximation of a Poynting dominated jet (σ ≫ 1) in the far-zone (R ≫ 1,γ ≫ 1) (Lyubarsky 2009).

In the equilibrium case and for κ> 2 one can obtain the following relations for the jet shape: (15)with ϵ_p describing the ratio of the plasma pressure to the magnetic pressure: (16)where B₀ and p₀ are the characteristic magnetic field and the external pressure at the light surface, respectively. In the non-equilibrium case, the jet shape is asymptotically conical.

Evolution of the Lorentz factor can also be obtained, depending on the acceleration regime: (17)where is the local opening angle at the transition between the two regimes. From Fig. 16, we estimate that this transition occurs in the jet in M 87 around mas, which, in return, corresponds to .

The Lorentz factor and the collimation profile are fit simultaneously using the Levenberg-Marquardt least-squares minimization routine. Uncertainties in the fit parameters are determined from the standard error obtained from the estimated covariance matrix. The collimation profile is modeled using the relation from Eq. (15). The two phases of the acceleration, before and after , are modeled using the relations from Eq. (17) for the equilibrium and non-equilibrium case respectively.

Some modifications to the Lorentz factor are required for the fit to converge. An additional scaling, α_γ, was introduced modifying the Lorentz factor so that: γ^⋆ = a_γγ, and the Lorentz factor was phase shifted along the axial coordinate of the jet by , in order to be able to account for the fact that the maxima and minima found in the collimation profile do not correspond exactly to the maxima and minima found in the Lorentz factor profile. The Lorentz factor for the non-equilibrium case was obtained by Lyubarsky (2009) considering a non-oscillating conical jet shape. Figure 16 shows that, even after , the Lorentz factor is still oscillating. We account for this observation by extending, in our model, the oscillation to the non-equilibrium case, so that .

All parameters are allowed to vary freely, including the viewing angle θ and the angular velocity Ω. The best fit obtained under these conditions has a reduced χ² = 0.24. In detail, this corresponds to and for the collimation and Lorentz factor fit respectively. We recall that was found by modeling the collimation profile as a simple power law. This indicates an improvement, by a factor of ~2, of the quality of the fit provided by the Poynting flux approximation for the jet collimation profile. We note that the relatively small values obtained for the reduced chi-squared is most probably a consequence of our conservative choice for the errors in the jet radius and jet velocity.

The best fit is plotted in Fig. 17, and the derived parameters of the model are summarized in Table 6. It is remarkable that the viewing angle and the angular velocity agree, within the uncertainties, with the respective values found from the analysis of the counter jet (Sect. 4.2) and the jet rotation (Sect. 4.3). The oscillations in both the jet width and the Lorentz factor that we discussed in Sect. 2.4 are also particularly well reproduced by the model, and if we extend this solution, we find that the subsequent Lorentz factor located at z_obs = 6.8 mas and z_obs = 17.2 mas, coincide with the two regions of maximum apparent velocities identified in our WISE analysis of the MOJAVE observations (Sect. 3). Komissarov et al. (2015) have interpreted these oscillations as a standing magneto-sonic wave bouncing across the jet.

There are several reasons why a phase shift, φ_γ, and scaling factor , α_γ, are required in our modeling. A possible reason for the phase shift φ_γ is that the measured Lorentz factor may come from a streamline slightly inside the measured jet shape. In this case the maxima in the Lorentz factor will be offset from the maxima in the jet width, with a shift proportional to the jet radius (Norman et al. 1982). The need for the scaling factor, α_γ, may be found in the assumption made by Lyubarsky (2009) of a constant angular velocity for all field lines. This assumption is usually valid only for a jet launched from the magnetosphere. For the case of disk launching, we expect differential rotation and Tchekhovskoy et al. (2008) found γ<R under these conditions. We recall here also that relations were derived using assumptions, specifically R ≫ 1 and γ ≫ 1, which are only partially valid in our case.

While the Lorentz factor of the second acceleration regime is well modeled by the non-equilibrium case, it is not associated with transition to the conical jet shape that one should expect theoretically. The apparent opening angle at the location where the acceleration regime changes is , and therefore additional collimation is necessary to achieve the final apparent opening angle observed at kiloparsec scales. Early saturation of the Poynting flux could also result in a quenched acceleration, and would break the assumptions of a Poynting dominated flux used in this modeling. We will now investigate this possibility by using another method of solving the MHD equations.

4.4.2. Wind solutions

Solving simultaneously the full transfield and the Bernoulli equations is only possible in a restrictive number of cases. The transfield equation will determine the shape of the field lines while the Bernoulli equation governs the acceleration of the flow. It is however possible to solve only the Bernoulli equation by assuming a specific jet shape (e.g., Camenzind 1986; Fendt & Ouyed 2004; Toma & Takahara 2013). The resulting solution is called a wind solution. It is not strictly an exact MHD solution, and one should be careful not to over-interpret it. However it can be used in some specific cases, for example, for investigating the acceleration efficiency, or in our case, the acceleration profile.

In the context of a cold ideal MHD jet, the poloidal and toroidal velocities can be written as follows: (18)where B_p is the poloidal magnetic field, μ is the total (magnetic + matter) specific energy, γ_in and η are the Lorentz factor and the distribution of mass flux at the inlet of the flow, respectively. The quantity η is a field line constant and it is given by the following equation: (19)with ρ₀ describing the rest-mass density. We also define the ratio of the Poynting to the matter energy flux, also called magnetisation parameter, as: (20)The Bernoulli equation is obtained from: (21)Substituting β_p and β_φ using Eq. (18), one can reduce the Bernoulli equation to a quartic equation for γ. To solve it, one has to prescribe a model for the poloidal magnetic field B_p, which turns into a model for the flux function.

Toma & Takahara (2013) recently derived a set of models for B_p based on a realistic flux function Ψ(r,z) which describes parabolic field lines appropriate for M 87: (22)where ψ = Ψ/Ψ₀ and Ψ₀ is the total magnetic flux, and ζ(ψ), A(ψ), ϖ(ψ), and a(ψ) are arbitrary functions of ψ (Toma & Takahara 2013). The poloidal magnetic flux is then obtained recalling that B_p = (∇Ψ × φ) /r, and this yields (23)with primes denoting derivatives of the respective functions over ψ.

At r ≲ r_lc, the true shape of field lines will probably deviate from the shape derived from Eq. (22) where field lines are anchored into the accretion disk or the magnetosphere. However, the approximation should be valid in the region we are investigating.

In our case of interest, associating the jet sheath with the outer poloidal field line (ψ = 1) and using θ = 18° and Ω = 10^-6, we have ζ = ϖ = 0, A = 3.7 and a = 1 /k = 1.73, corresponding to the parabolic shape obtained by Asada & Nakamura (2012). We also adopt γ_in → 1.

With the assumptions specified above, the Bernoulli equation still contains 7 unknowns: the total specific energy μ, the distribution of mass flux at the inlet η, the total magnetic flux Ψ₀, and four parameters related to the geometry of the poloidal field lines (a′, A′, ζ′ and ϖ′). We can however isolate some parameters and treat different cases separately. We define three models corresponding to different configurations of the flux function:

• Model 1:

  a′ ≠ 0, A′ = 0, ζ′ ≠ 0, ϖ′ = 0,

• Model 2:

  a′ ≠ 0, A′ = 0, ζ′ = 0, ϖ′ ≠ 0,

• Model 3:

  a′ ≠ 0, A′ ≠ 0, ζ′ = 0, ϖ′ = 0.

The Bernoulli equation has two singularities at the Alfvénic point (AP) and at the fast magnetosonic point (FMP). The wind solution in which γ is growing toward infinity has to pass through both the AP and the FMP. The Bernoulli equation, being quartic for γ, admits four roots. To be valid, the wind solution needs to pass smoothly from one root to the other crossing the singularities. This provides an additional constraint on the solution, and for a given B_p(R) and η, a unique value of μ corresponds to a valid wind solution. Further examining the Bernoulli equation, we find that Ψ₀ and η can not be separated, and their ratio must therefore be treated as a single variable. Hence we end up with only three free variables for the three models described above.

A Nelder-Mead least-squares minimization routine is used to fit the Bernoulli equation to the measured Lorentz factor evolution. Within each iteration of the fitting algorithm, μ is determined so as to maintain a valid wind solution.

Satisfactory fits were obtained for Models 1 and 2, however we could not obtain a reasonable fit using Model 3. In order still to be able to analyze the impact of the parameter A′ on the solution, we assigned a small non-zero value to ϖ′ in this model (ϖ′ = 0.05). The results for the three models are summarized in Table 7 and corresponding wind solutions are plotted in Fig. 18.

Figure 15 suggests that the velocities obtained for knot HST-1 and from the 43 GHz VLBA images originate from the same field line. To test this hypothesis we also include to the fit the speed measured at HST-1 by Giroletti et al. (2012), assuming r_{HST - 1} = 26 mas (Asada & Nakamura 2012). This addition does not affect the reduced χ₂ of the fit, thus supporting our hypothesis.

Models 1 and 2 best reproduce the measured acceleration. Jet parameters derived with these models are also similar. In Model 3, a higher value of μ was required, and the resulting conversion from Poynting flux to kinetic energy is slower. In all three cases, the Alfvénic point is found to be located close to the light cylinder, in agreement with exact self-similar solutions (Vlahakis & Königl 2004). The fast magnetosonic point is located at r ~ 0.4 mas, corresponding to z ~ 1 mas, and in agreement with the expectation that γ_FMP ≃ μ^1/3.

A true solution of the Bernoulli equation for the jet in M 87 would most likely be represented by a combination of the three models. We also do not know the exact shape of the field line close to the central engine, and so we cannot reconstruct the initial acceleration which might have helped discriminate between the models, as suggested by the divergence of the poloidal and toroidal velocities observed in Fig. 18. Our result however strongly favors a moderate value of μ with an efficient transformation of Poynting flux to kinetic energy (Models 1 and 2).

4.4.3. Transition to a kinetically dominated jet

The wind solution obtained for the sheath suggests a total specific energy at the base of the jet μ ~ 6–10. In the first two models described in Sect. 4.4.2, conversion between magnetic and kinetic energy is efficient, and the equipartition condition σ ~ 1 is achieved at: (24)corresponding to (25)Results from MHD simulations of AGN jets (Komissarov et al. 2007) indicate r_eq ~ 30 r_lc. Considering the choice of total energy μ ~ 18 used for these simulations, our results for Model 1 and Model 2 are qualitatively similar to the results of Komissarov et al. (2007).

The conversion of Poynting flux to kinetic energy is less efficient for our Model 3, and equipartition is only achieved at r_eq ~ 40 mas ~ 300 r_lc. However this model is altogether less favorable as it does not reproduce the observed acceleration as well as the first two models.

4.5. Nature of the spine

In Sect. 4.1, we identified two primary velocity components in the jet. Further analysis has shown that the slow component reflects either pattern speed of a plasma instability or the bulk speed of an outer slower wind. The fast velocity component has been successfully associated with the bulk speed along a magnetic field line connected to the observed jet shape, and we have obtained a wind solution for it. An inner deboosted part of the jet could therefore be associated with an even faster speed along a field line located closer to the jet axis.

The limb brightened structure found in the M 87 jet essentially at all scales might be explained by the combination of an enhanced magnetic field at the jet edge, a faster speed close to the jet axis and a lower emissivity in the jet spine. At kiloparsec scales, kinematic analysis systematically found a higher apparent velocity in the optical domain (Biretta et al. 1999; Meyer et al. 2013) compared to the radio domain (Biretta et al. 1995; Cheung et al. 2007; Giroletti et al. 2012)², indicating that they trace different regions in the jet, and that the inner part (spine) of the jet is dominant at optical wavelength, as was also suggested earlier by Perlman et al. (1999).

In the stacked-epoch image of the jet in M 87, we measure an intensity ratio between the inner part (the spine) and the edge (the sheath) ℛ ~ 0.5 at z_obs ~ 3 mas. At this distance, we also observe γ_sheath ~ 2.7. Assuming that the lower intensity near the jet axis is mainly due to deboosting of the spine, we can calculate an estimate of the Lorentz factor of the spine necessary to reproduce the observed intensity ratio using: (26)Solving this equation, we obtain γ_spine ~ 6.5. If we combine this estimate with the speed observed at HST-1 in the optical and assume that they are connected to a single field line, we can find a wind solution that would satisfy these observations. To obtain an accurate estimate of the velocity at HST-1 in the optical, we take the weighted mean of the four apparent velocities measured by Biretta et al. (1999) at HST-1, which gives us γ_spine = 9.8 ± 1.9.

Starting with the best solution found in Sect. 4.4.2 using Model 1, we set Ω = 2.5 × 10^-6 s^-1 corresponding to a field line anchored to a launching location closer to the ISCO (r₀ ~ 2.8 R_s), μ = 12, and A = 6.6 so that r_spine(z) = 0.5r_sheath(z). We then look for a valid wind solution that would fit the two estimates of the Lorentz factor of the jet spine. The result, plotted in Fig. 19, indicates that while the sheath field line can be associated with the velocity measured in the radio at HST-1, the spine field line connects well the inferred velocity of the deboosted structure in 43 GHz VLBA maps and the velocity measured at HST-1 in the optical band.

While this model would explain the limb brightened structure at z_obs ≳ 2 mas, the decrease in spine speed close to the core would result in an accompanying increase in Doppler boosting of the spine. In the simplest interpretation a ridge brightened structure at z_obs ~ 1 mas with brightness ratio ℛ ~ 5 would occur. This is not observed in the 7 mm VLBA observations, or in the 3 mm Global VLBI Array (GMVA) maps (Hada et al. 2016). This lack of ridge brightening suggest a low spine emissivity as the reason for the observed limb brightening down to the smallest core separations.

4.6. Launching mechanism

The angular velocity of a sheath field line has been obtained independently from the observed flow rotation and from the modeling of jet acceleration and collimation. The two estimates agree remarkably well, within the uncertainty, with value Ω ~ 10^-6 s^-1. This parameter can be used to discriminate between the two major mechanisms that operate during jet launching. In the Blandford-Payne model, the jet is launched from the accretion disk and the field line angular velocity depends on the distance from the black hole at which this field line is anchored. If we assume a Keplerian rotation in the accretion disk, this angular velocity corresponds to r₀ ~ 5R_s. In the Blandford-Znajek mechanism, the angular velocity depends on the spin of the central black hole a_H. Following Tchekhovskoy (2015), the black hole angular speed is: (27)One can then obtain the angular speed of the respective field line,(28)Using 230 GHz VLBI observations, Doeleman et al. (2012) estimated the spin of the black hole in M 87 to be a_H ~ 0.6 based on a measure of the size of the smallest resolvable structure, identified as the ISCO. This value corresponds to Ω_F = 2.75 × 10^-6 s^-1 and is a factor three times larger than the angular velocity that we determined. The BP mechanism would therefore better account for our result, but a BZ launching cannot be firmly discarded owing to the uncertainty of this measurement. Nevertheless, the wind solutions obtained for the sheath and the spine imply μ_spine ≳ μ_sheath, which is a characteristic of differential rotation as found both analytically and from numerical simulations (Lyubarsky 2009; Komissarov et al. 2007). Thus we favor a disk launching at least for the sheath.

4.7. Mass-loss rate estimate

The mass-loss rate corresponds to the amount of matter extracted from the accretion disk and loaded into the jet. It can be obtained from: (29)The wind solution describing the acceleration of the sheath requires B₀/η ~ 2–10, with B₀ describing the characteristic magnetic field at the light cylinder. The mass loss rate can then be estimated recalling that : (30)Using 230 GHz VLBI observations, Kino et al. (2015) estimated the magnetic field at the base of the flow to be about 100 Gauss. Using this value for B₀, we get Ṁ ~ 10^-7–10^-8 M_⊙ yr^-1.

The mass accreted onto the black hole of M 87 is still undetermined. Using Chandra X-ray observations, a Bondi accretion rate of Ṁ_B = 0.12 M_⊙ yr^-1 was estimated by Di Matteo et al. (2003). However, the measured X-ray luminosity, L_X ~ 7 × 10⁴⁰ erg s^-1, suggests a mass-loss rate of Ṁ_acc ≪ Ṁ_B. Kuo et al. (2014) found an upper limit of Ṁ_acc< 1 × 10^-3 M_⊙ yr^-1 using Faraday rotation. This limit would yield Ṁ ~ 10^-4–10^-5Ṁ_acc.

Thus, the mass-loss rate is much lower than the total mass accreted onto the black hole, implying that only a small fraction of the accreted matter is transported into the jet.

5. Summary

In this paper, we have presented a detailed analysis of the two-dimensional kinematic evolution of the innermost part of the jet in M 87. This source has allowed us to investigate jet formation and propagation on scales of 10²–10³R_s. A WISE analysis (Mertens & Lobanov 2015, 2016) was performed on 43 GHz VLBA maps observed as part of the M 87 VLBA movie project (Walker et al. 2008; in prep.). We summarize here the main findings of this analysis:

• 1.

  We have obtained the first complete velocity field of an AGN jetat sub-parsec scale, revealing a structured and highly stratifiedjet, with a significant transverse velocity component.

• 2.

  The structure of the flow was also investigated using SWD decomposition. The shape of the streamline is parabolic, well described by a power law r ∝ z^0.6.

• 3.

  We observe an oscillation in flow expansion with a spatial period increasing with the distance as ~z_obs, which is correlated with a similar oscillating pattern in the apparent speed.

• 4.

  The viewing angle θ ≃ 18° was obtained using the speed of a moving component discovered in the counter jet, and the jet to counter jet intensity ratio. This viewing angle was also confirmed by the analysis of jet rotation and by the modeling of jet acceleration.

• 5.

  Flow stratification was analyzed using the SCC method (Mertens & Lobanov 2016), revealing a slow, mildly relativistic layer (β ~ 0.5c) associated either with an outer wind or with instability pattern speed, and a fast accelerating streamline (γ ~ 2.5 at z_obs ~ 3 mas)

• 6.

  Acceleration and collimation of the flow was modeled as an MHD jet, using two different methods. We found that the acceleration is at first linear γ ∝ r ∝ z^0.58, and then saturates, implying a total specific energy μ ~ 6, and that equipartition between Poynting flux and kinetic energy is reached at a distance z_eq ~ 3000 R_s, in excellent agreement with previous analytic and numerical simulation work. This suggests that the jet is already kinetically dominated at a distance of a few parsecs from the central engine.

• 7.

  We could not detect any features that could unequivocally be related to an inner spine. We estimated a spine Lorentz factor γ_spine ~ 6 that would deboost the spine emission sufficiently to explain the intensity ratio between sheath and spine. Subsequently we found a wind solution that would connect this Lorentz factor with the velocity observed at knot HST-1 in the optical regime. However, this solution breaks down close to the core where it would predict a ridge brightened instead of a limb brightened jet. This suggests a lower intrinsic emissivity for the jet spine.

• 8.

  Indication of jet rotation was found from the difference in velocities observed in the northern and southern limbs. This observation was used to obtain an angular velocity for the field line, Ω ~ 10^-6 s^-1 which, in turn, suggests that this streamline is most likely launched in the inner region of the accretion disk, at a distance r₀ ~ 5 R_s from the central engine. However, a BZ scenario in which the jet is launched from the magnetosphere of the central engine could not be excluded.

• 9.

  Finally, an estimate of the mass-loss Ṁ ~ 10^-7–10^-8 M_⊙ yr^-1 was obtained for the sheath, from the wind solution. This corresponds to Ṁ ~ 10^-4–10^-5Ṁ_acc.

Taken together, the WISE analysis of the jet in M 87 draws a coherent picture of a magnetically launched, accelerated, and collimated jet. Continuation and expansion of the work presented in this paper will enable testing the results obtained for the jet in M 87 and addressing other important aspects of the physics of relativistic flows. The methods that we have used to model the acceleration and collimation zone in the M 87 jet can be improved further to provide a more detailed physical description of the flow. One possibility would be to use a more realistic flux function obtained from numerical simulations. Another possibility would be to attempt to reproduce our results by performing full 2D or 3D MHD simulations using the physical parameters estimated for the jet in M 87. This would allow us to investigate the inner part of the jet, providing further clues about the reason for the edge brightened morphology.

The milliarcsecond-scale jet in M 87 has been observed at 15 GHz as part of the MOJAVE project. We demonstrated in this work that the fast stream line can be robustly detected at this frequency with VLBA observations made roughly every two months. Obtaining and analysing a well-sampled 15 GHz VLBA dataset on M 87 would enable us to extend the analysis of the acceleration profile up to z ~ 10⁴R_s. The high resolution provided by the GMVA observations at 86 GHz could also be critical in confirming the results we obtained. The structure of the flow could be probed with GMVA at distances as small as ~10R_s, providing an excellent probe of the innermost part of the jet.

---

Acknowledgments

This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2009). F.M. was supported for this research through a stipend from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. We sincerely thank Y. Kovalev for providing the VLBA images at 2 cm. F.M. thanks Y. E. Lyubarsky, C. Fendt, K. Toma and A. Tchekhovskoy for helpful discussions. We would also like to thanks the anonymous referee for the careful review that improved the quality of this paper. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

References

Appendix A: Intermediate-scales wavelet decomposition

Canonical wavelet transformations are performed on a set of scales ranging as powers of two of the initial, smallest scale used. The ability of the multiscale cross correlation (MCC; Mertens & Lobanov 2015) procedure to cross identify significant structural patterns (SSPs) in pairs of images relies on a smooth transition between scale j and scale j + 1 of the wavelet transform employed. This condition sets limits on detecting different speeds at different scales, and hence providing an accurate assessment of flow stratification traced by overlapping optically thin features with different intrinsic speeds (Mertens & Lobanov 2016). For any two adjacent scales of the wavelet transform, this limit is on the order of 2^j (Mertens & Lobanov 2015). As the velocity difference between two overlapping features becomes comparable or smaller than this limit, matching becomes more difficult. To improve cross identification of the individual features, one can construct a finer scaling function for which the wavelet decomposition will cover intermediate scales. This is performed by introducing an intermediate-scale wavelet decomposition (IWD). In a canonical wavelet transform, the coefficients w_j at a scale j of the wavelet decomposition contains information on spacial scales between 2^j−1 and 2^j. In the IWD, the coefficients at a given scale j contain information on spacial scales between 1.5 × 2^j−1 and 1.5 × 2^j. The MCC is then performed by inserting the scale j of the IWD between the scales j−1 and scale j of the SWD.

The IWD is realized in the WISE analysis by introducing an appropriate scaling function. The B-spline scaling function used in WISE has the following discrete low pass filter associated with it: (A.1)For the IWD, a modified B-spline 3 scaling function is used which has a discrete low pass filter in the following form: (A.2)The resulting sensitivity of the IWD to different structural scales is illustrated in Fig. A.1 which compares the responses of different wavelet scales to different spatial scales in the image, for the scaling functions B-spline and B-spline 3.

All Tables

All Figures