The Sub-Parsec Scale Jet in M 87

The Structure and Dynamics of the Sub-Parsec Scale Jet in M 87 Based on 50 Vlba Observations Over 17 Years at 43 Ghz

Abstract

The central radio source in M 87 provides the best opportunity to study jet formation because it has a large angular size for the gravitational radius of the black hole and has a bright jet that is well resolved by VLBI observations. We present intensive monitoring observations from 2007 and 2008, plus roughly annual observations that span 17 years, all made with the the Very Long Baseline Array at 43 GHz with a resolution of about 30 by 60 . Our high-dynamic-range images clearly show the wide-opening-angle structure and the counter-jet. The jet and counter-jet are nearly symmetric in the inner 1.5 milli-arcseconds (mas; 0.12 pc in projection) with both being edge brightened. Both show deviations from parabolic shape in the form of an initial rapid expansion and subsequent contraction followed by further rapid expansion and, beyond the visible counter-jet, subsequent collimation. Proper motions and counter-jet/jet intensity ratios both indicate acceleration from apparent speeds of to in the inner mas and suggest a helical flow. The jet displays a sideways shift with an approximately 8 to 10 year quasi-periodicity. The shift propagates outwards non-ballistically and significantly more slowly than the flow speed revealed by the fastest moving components. Polarization data show a systematic structure with magnetic field vectors that suggest a toroidal field close to the core.

galaxies: individual (M 87) — galaxies: jets — galaxies: active — radio continuum: galaxies — hydrodynamics — relativistic processes
\submitjournal

the Astrophysical Journal

\correspondingauthor

R. Craig Walker

0000-0002-6710-6411]R. Craig Walker \altaffiliationThe National Radio Astronomy Observatory (NRAO) is a facility of the National Science Foundation (NSF), operated under cooperative agreement by Associated Universities, Inc. (AUI). \affiliationNational Radio Astronomy Observatory, Socorro, NM 87801

\affiliation

Department of Physics & Astronomy, The University of Alabama, Tuscaloosa, AL 35487

\affiliation

Department of Physics, University of California, Santa Barbara, CA 93106-9530, USA

0000-0002-4245-2318]Chun Ly \affiliationSteward Observatory, 933 N Cherry Ave, Tucson, AZ 85721 \affiliationMMT Observatory, 933 N Cherry Ave, Tucson, AZ 85721

0000-0003-4740-8321] William Junor \affiliationISR-2, MS-B244, Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM 87545

1 Introduction

M 87 (Virgo A, NGC 4486, 3C 274) is a giant elliptical galaxy near the center of the Virgo Cluster that contains a very massive black hole and a prominent jet. The scale, in gravitational units per unit angle, is the best for any jet with observed structure making M 87 a prime target for studies of the jet launch region. At a distance to M 87 of  Mpc (Mei et al., 2007), pc and  mas yr. At this assumed distance, the black hole in the center of M 87 has a mass, determined from stellar dynamics, of M (Gebhardt, 2011, adjusted for a different assumed distance). This gives a scale of  R where the Schwarzschild radius  cm.1 With this scale, Very Long Baseline Array (VLBA; Napier et al., 1993) observations at 43 GHz (7 mm), with a resolution of 0.2 to 0.4 mas, probe structures with sizes on the order of 100 .

The well-known, kpc-scale jet in this galaxy is a prominent source of radio, optical, near-ultraviolet (NUV) and X-ray emission. There are remarkable similarities between the radio, optical and NUV emission on projected scales of kpc. The bright knots, the filamentary emission between the knots, and the bends and twists in the jet can easily be identified in the optical and NUV as well as the radio (Perlman et al., 2001). While X-ray images from Chandra (Marshall et al., 2002; Wilson & Yang, 2002) are not of comparable resolution to the radio or optical, the same bright knots and overall jet structure can easily be seen in the X-rays.

The M 87 jet has long been known to have a wide-opening-angle base on scales well below a mas (Junor, Biretta & Livio, 1999). On average the jet expands parabolically out to about the first bright knot (known as HST-1) at (projected distance pc) from the core (Asada & Nakamura, 2012; Hada et al., 2013), but with evidence for faster expansion at the smallest scales (see Hada et al., 2013, 2016). Older VLBI observations at 1.6 GHz indicate a full width half-maximum (FWHM) apparent opening angle of out to  mas (projected distance pc) from the radio core (Reid et al., 1989). Beyond HST-1 the jet expands conically (Asada & Nakamura, 2012), and Very Large Array (VLA) observations at 15 GHz indicate a FWHM apparent opening angle of from HST-1 to knot A located at (projected distance pc) from the radio core (Owen, Hardee & Cornwell, 1989). The radio core in M 87 is thought to be nearly coincident with the position of the black hole, unlike the situation in many blazars (e.g., Marscher et. al., 2008). Multi-frequency VLBA astrometric observations indicate an offset of only 41 as at 43 GHz (Hada et al., 2011), equivalent to about 6 .

Proper motions of jet components in the radio and optical range from subluminal to superluminal. A component located at  mas (projected distance pc), shows an angular motion of   mas yr that corresponds to a subluminal motion of (Reid et al., 1989). Kovalev et al. (2007) find motions of less than a few percent of the speed of light at 15 GHz for several components within 25 mas of the core with the VLBA. That result is based partly on the MOJAVE program (Lister & Homan, 2005), a recent reanalysis of which gives hints of acceleration and mild superluminal motion (Britzen et al, 2017). At higher frequencies, both subluminal and mildly superluminal (up to ) speeds, with evidence for acceleration, have been seen based on the 43 GHz VLBA data to be presented in this paper (Mertens et al., 2016, hereafter MLWH) and based on KaVA (KVN and VERA Array (Korean VLBI Network and VLBI Exploration of Radio Astronomy)) data at 22 GHz (Hada et al., 2017). Slow motions were seen by Hada et al. (2016) with 4 epochs, two at 86 GHz and two at 43 GHz. With one exception, those components were at core distances of about 1 mas or less. The fastest observed proper motion, found in the optical band at HST-1, is superluminal with (Biretta, Sparks & Macchetto, 1999). The fastest observed radio proper motions, also at HST-1, are superluminal with (Cheung, Harris & Stawartz, 2007), and and (both Giroletti et al., 2012). Giroletti et al. (2012) also find less reliable evidence for a component moving at , similar to the optical result. Asada et al. (2014) find evidence from their 1.6 GHz global VLBI data plus other data from the literature for acceleration from subluminal to superluminal speeds between about 160 mas and HST-1. In general, at and beyond HST-1 the radio proper motions are less than optical proper motions in the same region, and proper motions suggest jet deceleration (see Hardee & Eilek, 2011).

An important parameter for interpretation of the data is the jet angle to the line-of-sight. There is no one, simple way to determine that angle so the reasoning leading to the choice of the angle used in this paper is explained here at some length. The fastest observed proper motions and modeling can be used to constrain the viewing angle. The highest observed optical superluminal speed requires a jet viewing angle of . Recent model dependent values for the viewing angle are used by Nakamura & Meier (2014) to provide a self-consistent shock model for optical relativistic motions and particle acceleration at HST-1, and a value of found by Wang & Zhou (2009) using synchrotron spectrum model fitting to the kpc jet. MLWH found that a viewing angle of (see Table 6 in MLWH) provided a fit to the observed acceleration and collimation of the jet assuming a Poynting flux dominated approximation for the jet dynamics.

Less model dependent values for the viewing angle can be found in MLWH. A value of [see eq. 4 in MLWH] is obtained by MLWH between 0.5 and 1 mas from the core using apparent component speeds in the jet, , and the counter-jet, along with an intensity ratio, . An additional analysis in MLWH, based on the assumption that the significantly different speeds in northern and southern limbs of the jet are the result of jet rotation, leads to a viewing angle value of (see eq. 9 in MLWH). Thus, MLWH concluded that a viewing angle of represented a reasonable average based on these two techniques and Poynting flux dominated jet dynamics. We note here that the jet rotation and Poynting flux model viewing angle estimate lies above the maximum value allowed by the optically observed motion at HST-1, and the average viewing angle estimate lies uncomfortably close to the maximum viewing angle. Thus, we adopt the smaller viewing angle from MLWH of in order to relate angular scales to intrinsic scales, and point out that all the recent model dependent viewing angles lie within of this choice.

At a viewing angle of , intrinsic lengths are times longer than apparent (projected) lengths, and the intrinsic jet opening angle between HST-1 (intrinsic distance  pc) and knot A (intrinsic distance  pc) is . At that viewing angle, the observed optical and radio superluminal motions of and imply intrinsic Lorentz factors of () and () respectively.

The jet is observed to be edge-brightened from about 50 to 200 mas and also beyond . At angular distances out to  mas, the edge-brightened radio jet exhibits quasi-periodic deviation of the brightness centroid (Reid et al., 1989) that likely indicates brightness alternation from one side to the other. Beyond HST-1 the jet is similarly edge-brightened at 15 GHz and side-to-side alternation in the radio brightness profile is apparent from HST-1 through knot D located at (intrinsic distance  pc) from the 15 GHz radio core (Owen, Hardee & Cornwell, 1989). The jet also appears to contain twisted filaments from HST-1 to knot A (Lobanov, Hardee & Eilek, 2003) that are related to side-to-side variation in the brightness profile (Hardee & Eilek, 2011).

The jet structure seen in the optical and NUV is quite similar to that seen at 15 GHz, but there are important differences in detail. The optical/NUV emission is more concentrated in the knots and towards the jet axis, and the outer edges of the jet are less well defined than in the radio (e.g., Sparks, Biretta & Macchetto, 1996; Madrid et al., 2007). The major bright radio or optical knots can also be identified in the lower-resolution X-ray images, e.g., Perlman & Wilson (2005), but they can differ in position or structure from their radio or optical counterparts. At 15 GHz the projected magnetic field lies more or less along the edge of the jet flow (Owen, Hardee & Cornwell, 1989), which may be indicative of a shear layer (Laing, 1981). The optical/NUV knots also differ somewhat in polarization from the radio knots in that they tend to indicate magnetic field perpendicular to the jet flow, and suggest the presence of shocks (Perlman et al., 1999). In general, the broader radio profile, the more centrally concentrated optical/NUV structure, differences in the polarization structure, and the slower observed motions in the radio suggest velocity shear across the jet with the limb brightening observed in the radio being associated with the shear layer.

In this paper, we present results from monitoring observations of M 87 at 43 GHz using the VLBA. Images from 23 epochs in 2007 and early 2008 are used to study the structure and dynamics of the inner jet. An image showing the average structure was made by stacking the 23 images and is used for the structure analysis. In that image, the limb brightened jet is observed to 20 mas beyond the core. A faint structure, that extends slightly more than 1 mas beyond the core in the opposite direction, is speculated to be a counter-jet (Ly, Walker & Junor, 2007; Kovalev et al., 2007). We use transverse intensity slices to obtain the maximum intensity as a function of distance from the core and the width as a function of distance from the core. The individual images, and movies made with them, are used to study the rapid and complex changes near the core in an effort to extract velocities and possible acceleration. The 2008 observations showed a significant increase in flux density of the VLBA core that coincided with a flare seen at TeV energies as reported in Acciari et al. (2009). An early version of the 23-epoch average image, along with images from four individual epochs, was presented in that initial flare report. The flare led to an effort to catch additional cases of correlated radio and TeV flares, a result of which is that we have an ongoing series of roughly annual 43 GHz VLBA observations of M 87. We report results, based on those observations to 2016, along with legacy observations going back to 1999, of a study of changes to the jet envelope on time scales of years. Preliminary presentations of some of our data, including an initial movie, have appeared in Walker et al. (2008) and Walker et al. (2016) while an extensive analysis of the apparent velocity field in 2007 is given in MLWH.

The paper is organized as follows. In Section 2 we describe the observations. In Section 3, we present the observational results with subsections on the average structure (3.1), on the rapid structural evolution (3.2), on the polarized structure (3.3), on the long term evolution between 1999 and 2016 (3.4), and on the counter-jet structure based on an average image and on recent, high-sensitivity observations (3.5). In Section 4 we examine some of the possible implications of our observations for jet structure and collimation properties that will provide a framework for future theoretical and numerical modeling of jet acceleration and collimation regions. Finally, we summarize and further discuss our results in Section 5.

For a listing of the symbols used in this paper, see Table A1 in Appendix A.

2 The Observations

2.1 Data Aquisition

The primary goals of the 43 GHz VLBA M 87 project were to determine the apparent jet velocity and as much as possible about the dynamics of the jet close to the black hole. M 87 had been observed before, but with some combination of inadequate resolution to see very close to the core and time sampling intervals that were too long unless the motions are very slow (e.g., Reid et al., 1989; Junor & Biretta, 1995). For example, if the apparent motions are at the speed of light, a component would be moving at 4 mas per year which is about 20 beam widths for the VLBA at 43 GHz. Some other sources have apparent speeds of more than ten times that. If fast motions exist in M 87 near the core as they do on larger scales, they probably had not been detected because they were seriously undersampled. So a pilot project was conducted in 2006 during which time intervals of between three days and three months were sampled. The results were somewhat ambiguous, but it was concluded that the changes that were seen could be followed with a three week sampling interval. Such observations were made throughout 2007. During that period it became apparent that the sampling interval was too long — apparently related features moved two or more beam widths between observations and relating features was difficult. In early 2008, the observations were continued at about 1 week intervals. Despite the use of dynamic scheduling, the short interval forced the use of less optimal observing conditions, so the data quality in 2008 was not as high as in 2007. Many sessions had missing antennas which is not tolerated well when imaging a complex source with the sparse UV coverage of the VLBA. The primary data set for this paper is the first 11 epochs from 2007 and 12 epochs from 2008. One additional epoch, 2007 Nov. 02, is shown, but the remaining five 2007 epochs have technical problems so they are not used in this paper.

The correlation between a TeV flare and a radio flare noted in Section 1 led to an ongoing effort to catch other such flares in the hopes of further constraining the TeV emission location and mechanism. The TeV instruments could only observe at night which constrains the coordinated observations of M 87 to happen during the early part of each year. Observations were made with the VLBA once or twice each year to monitor the ambient condition of the radio jet. If a high energy flare was seen, additional observations were triggered. There were triggers in 2010 and 2016 but no clear correlated radio flare was found in our data. With our data plus an additional observation of their own that had high flux density and was almost coincident with the TeV flare, Hada et al. (2012) reported a weak radio flare associated with the 2010 TeV event. In 2012, there was a period of enhanced TeV emission that failed to reach the trigger level for our project, but for which observations on the VERA and the European VLBI Network detected a significant radio flare (Hada et al., 2014). A by-product of our 2010 trigger response is that we have a six-image stack for 2010. In 2016, most of the follow-up observations were too short to make good images, so only three full-track observations are available and one of those had poor observing conditions.

All of the observations reported here were made on the 10-station VLBA. Only the first two archival observations in 1999 and 2000 used additional antennas as noted in Table A2 in the Appendix. That table gives the dates, image noise, and peak and total flux densities from the individual observations that contributed to this paper (not for all of the observations that were made). The long-term analysis also used four average observations, listed in the Appendix in Table A3, one from each year with six or more individual observations.

The observations from 2004 and earlier are from the VLBA archive. They provide roughly annual information on M 87 at 43 GHz back to 1999. Several of those observations involved co-authors on this paper and several involved the use of M 87 as a phase reference source for astrometric studies or for observations of nearby weak Virgo Cluster sources. Literature references for the projects involved can be found in Table A2. Most of those observations involved less, sometimes significantly less, time on source and less bandwidth than the post-2006 observations so the image quality is generally not up to the same standard.

Starting with the pilot project in 2006, the observations involved a full track with data spanning most of the time that M 87 was visible from at least a few stations of the VLBA. The total time on M 87 at 43 GHz varied from about 5 to about 7 hours. The sources OJ 287, 3C 279, and occasionally 3C 84 or J1635+3808 were observed to use as polarization calibrators, fringe finders, and bandpass calibrators. Beginning in 2007, short phase-referencing scans on M 84 were observed in an attempt to measure the relative proper motions of M 84 and M 87 and to constrain any core wander. Such data are also available from 2001 October 12 and 2004 April 4 (Ly, Walker & Junor, 2007). To improve the quality of that astrometry, two to four roughly half-hour bursts of scans on sources around the sky were added to calibrate the atmosphere starting with 2008 March 12. Beginning in 2009, roughly 20% of the M 87 on-source time used 24 GHz to allow study of the spectral index and, possibly, Faraday rotation. The results from most of the polarization data, from the 24 GHz imaging, and the M 87/M 84 astrometry will be reported elsewhere.

Various changes have occurred to the VLBA hardware over the years spanned by this project. The most important is that the recording bandwidth has been increased both by increasing availability of recording media and, eventually, by a significant hardware upgrade. The bandwidths used are shown in Table A2. Between the pilot project, recorded at 128 Mbps (16 MHz per polarization) and the most recent observations, recorded at 2048 Mbps (256 MHz per polarization), the bandwidth improvements alone gave a sensitivity increase of a factor of four. The data were correlated at the Pete V. Domenici Science Operations Center in Socorro, New Mexico, on the original hardware correlator up through 2009. A new Distributed FX (DiFX) software correlator that provides significantly enhanced capabilities (Deller et al., 2011) began operation in time to be used for our 2010 and later epochs. A description of the status of the VLBA hardware at any given time, at least since 2009, can be found in the Observational Status Summaries such as Van Moorsel (2014) and others that can be found at the same site.

2.2 Data Reduction

The data were reduced with the Astronomical Image Processing System (AIPS; Greisen, 1995) following the usual procedures for VLBI data reduction including correction for instrumental offsets using the autocorrelations, bandpass calibration based on strong calibrator observations, and correction for atmospheric opacity based on the system temperature data. The a-priori amplitude calibration depended on the gains provided by VLBA operations, which are based on results from regular single-dish pointing observations of Jupiter, Saturn, and Venus averaged over many months. Additional calibrations were done to enhance the M 87/M 84 astrometry and are included here for completeness. The Earth Orientation Parameters were updated to values more accurate than those available at the time of correlation. For many of the epochs, corrections were made for the ionospheric delays based on global ionospheric maps produced by the Jet Propulsion Laboratory which are available at the NASA’s Crustal Dynamics Data Information System (CDDIS). Atmospheric delay corrections were made using AIPS task DELZN based on either delays measured on sources around the sky in the aforementioned 2 to 4 bursts of about a half hour duration or, for observations on or before 2008 March 06, on the residual fringe rates of the project sources and calibrators. These additional calibrations do not affect the images of M 87.

The images of M 87 were made using data that are both amplitude and phase self-calibrated. The flux scale for each epoch was set by normalizing the self-calibration gain adjustments to the a-priori gains for observations above elevation on only those antennas with good weather and instrumental conditions for that epoch. The flux densities are believed to be accurate to within about 10%, although the errors on the first few observations (especially 1999) may be higher. The phase self-calibration does not preserve source position so, for analysis of the dynamics, the images were aligned on the image peaks. During the radio flare of 2008, the astrometric observations of the M 87/M 84 relative positions suggest that there is a small core shift as new components appear ( Davies, F. et al., in preparation). For the 2008 event, the brightest we have seen by far, the effect was under 100 as while at other times it is less than a few tens of as, so the effect is small enough to ignore for purposes of this paper. For the imaging of this complex and fairly extended source, we found it important to use the robustness parameter (Briggs, 1995) to gain some of the advantages of both natural and uniform weighting together, and to be able to significantly vary the dirty beam, which seems to help convergence during the self-calibration and imaging iterations. We also found that significantly better images could be obtained using the multi-resolution capability of the imaging task IMAGR (Greisen, 2009). To aid comparison between epochs, nearly all images in this paper have had the point CLEAN components restored with a common beam of  mas, elongated in position angle . That resolution amounts to about with the adopted distance and black hole mass.

The observations were scheduled and correlated in a manner consistent with obtaining images of the polarized structure of M 87. The primary calibrators for this purpose were OJ 287 and 3C 279. These sources were used to determine the instrumental polarization using the standard procedures in AIPS meant for application to resolved calibrators of unknown structure. For the electric vector position angle (EVPA) calibration, we have utilized results from the Boston University Blazar project (Jorstad & Marscher, 2016, and references therein), which includes the two calibrators. That project has many more sources which aids the calibration significantly. We do not have simultaneous epochs, but past experience suggests that the absolute position angles will be the same when the images have the same total intensity structure, the same polarized intensity structure, and the same run of polarization angles along the source. When an epoch of the Blazar Project can be found with such a match for the total intensity and polarization structure for a calibrator, the EVPA calibration can be set to give the same polarization position angles for that calibrator. This method was used to calibrate the EVPA for both epochs for which we present polarization images later. Confidence in the method was gained from the fact that the two calibrators gave independent, but consistent calibrations for each epoch and the polarization angles for M 87 were consistent between the epochs. To our knowledge, this method of EVPA calibration has not been used before but should be useful for other projects.

For this project, obtaining a reliable polarization calibration proved to be possible, but very time consuming. It has only been completed to our satisfaction for two epochs, the results from which will be shown in Section 3.3. Those two epochs give similar results, so they likely reveal the nature of the polarization structure within about 0.5 mas of the core. To avoid significant additional delay in the publication of the total intensity results, we have chosen to defer a full polarization analysis. Once the rest of the polarization results are available, we should be able to study time dependencies and to detect the polarization farther down the jet. Also the 24 GHz data should allow study of Faraday rotation along the jet. Detection of polarization farther along the jet depends on higher sensitivity which can be obtained by stacking the 2007/2008 data and/or by using the recent, higher-sensitivity observations. Those recent observations are subject to frequency-dependent instrumental effects within our observing bands which, because of software limitations (since fixed), hindered our initial calibration attempts.

3 Observational Results

The results of our M 87 VLBA observational program are presented in this section. There are five sub-sections. Each sub-section begins by presenting images that address a different aspect of the information available from the observations. Then the images are analyzed to extract more quantitative information that should help understand the nature of the inner jet. Further discussion of what is learned from the data is deferred to Section 4.

Section 3.1 presents an average intensity image based on the 2007 and 2008 data, showing the persistent structure. An analysis based on slices is used to extract some of the detailed features of the jet shape which can shed light on the jet propagation near the core.

Section 3.2 shows the individual images from the 2007 and 2008 observations. These are the main movie observations, and the movies are available in the online version of the paper. An effort is made to measure component motions using traditional, component-tracking methods. This is not as sophisticated as what was done in MLWH, but provides added confidence in the detection of superluminal motion and of acceleration very near the core.

Section 3.3 presents the polarization results from the two epochs from 2007 for which we already have a successful polarization calibration. There is a systematic structure to the polarization near the core which suggests a toroidal field geometry.

Section 3.4 gathers the early observations, from before 2007, plus the results from ongoing efforts to catch more cases of activity related to TeV flares since 2009, to provide roughly annual images between 1999 and 2016. A movie based on these results is available in the online version of the paper. The jet is shown to move transversely over a period of several years. An analysis of that transverse motion, including derivation of a propagation speed of the pattern that is that sideways motion, is presented.

Section 3.5 uses the inner few mas of an average of nearly all our data over the full 17 years, along with images from two of the recent, high-sensitivity observations, processed specifically for high resolution, to delineate the structure of the counter-jet.

3.1 Average Intensity

A total of 23 images were used to make the movies showing the dynamics of the inner M 87 jet that are presented in Section 3.2. Those images have been averaged to improve the sensitivity and to show the persistent structure of the sub-parsec M 87 jet as seen at 43 GHz over a period of 1.2 years in 2007 and 2008. The detailed structure that varies from epoch to epoch is smeared out by the averaging process. The average image is shown in Figure 1.

Figure 1: A 23-epoch average radio image of the jet and counter-jet in M 87 based on data from 2007 and 2008. Angular to linear scales (in parsecs and in Schwarzschild radii ) are indicated for distances in the sky plane and for distances along the axis of the jet assuming it is oriented at to the line-of-sight. The beam with resolution  mas elongated in position angle is shown at the lower left. The off-source noise level is 62 Jy beam and the image peak is 0.83 Jy.

At a viewing angle of , 1 mas in the sky plane corresponds to an intrinsic (ie, de-projected) distance along the jet of 0.27 pc or about 480 . Structure on the counter-jet side, previously speculated to indicate a counter-jet (Kovalev et al., 2007; Ly, Walker & Junor, 2007), extends a little more than 1 mas to the other side of the radio core. There are three positions located at  mas,  mas, and  mas from the radio core where the jet appears to recollimate and then subsequently expand. The image, along with those of Figures 2 and 3, suggests an oblique filament crossing the jet associated with a brighter feature near the jet center at  mas from the radio core and possibly some fainter filaments crossing the jet at larger distances from the radio core. Overall the northern edge of the jet appears straighter than the southern edge. An average of all the epochs to be presented in Section 3.4 is not used for the analysis here because of smearing of the otherwise sharp edges caused by the transverse motions reported in that section.

Figure 2: A 3D slice profile rendering of the stacked 43 GHz VLBA image of M 87 shown in Figure 1. The amplitude scale is logarithmic to allow low level jet features to show without saturating the core region. The long axis (near vertical) is the right ascension offset. The short axis is the declination offset in mas. The image is rotated by to provide a good perspective on the edge-brightened jet structure.

A 3D slice profile rendering of the 23-epoch total intensity image shown in Figure 2 shows more quantitatively the edge-brightening seen in Figure 1. The southern edge is typically brighter than the northern edge by a few of tens of percent, but there is a wide scatter in the ratio with the north being brighter in a few locations.

In order to examine in detail the transverse structure and compare counter-jet to jet sides of the core, slices of the intensity image separated by 0.117 mas were made from  mas to 7.956 mas along the counter-jet and jet axis, i.e., separated by about half a beamwidth.

Figure 3: The locations of the sample intensity slices across the jet and counter-jet shown in Figures 4 and 5 are shown overlaid on a contour map version of the image in Figure 1. The contour levels, in mJy beam, are , 0.3, 0.6, 0.85, 1.2, 1.7, increasing from there by factors of . The convolving beam for the clean image is  mas elongated in position angle .

Figure 3 shows an intensity contour plot of the inner 10 mas of the image shown in Figure 1 with lines superimposed that mark the locations of the sample intensity slices shown in Figures 4 and 5. These sample slices are a subset of all of those calculated and used in the analysis. Here the slice axis is located at an angle of with respect to the R.A. axis. Figure 4 shows sample transverse intensity slices at 0.234 mas intervals in the inner  mas centered on the 43 GHz core.

Figure 4: Sample intensity slices at distances of 0.467, 0.702, 0.936 and 1.170 mas with respect to the radio core in the counter-jet (top row) and jet (bottom row) directions. The slice locations are shown in Figure 3. The northern edge of the jet and counter-jet is at positive values on the transverse (mas) scale.

Similar transverse edge-brightened structure is evident on both jet and counter-jet sides of the core, albeit on the counter-jet side inside  mas there is one position where the maximum intensity is interior to the bright ridgelines (at  mas; not shown in Figure 4) and one position where the intensity is centrally peaked (at  mas; shown in Figure 4). On the jet side all slices inside 1.2 mas are edge-brightened. The northern edge is brighter nearer the core along both jet and counter-jet in the two innermost slices. Along the jet the southern edge becomes brighter in the two slices farthest from the core, and there is weak indication of similar behavior along the counter-jet in the slice farthest from the core. In the 23-epoch intensity image we cannot follow the counter-jet farther from the core.

Figure 5 shows sample transverse intensity slices along the jet from  mas to  mas in 0.702 mas (approximately three beamwidth) intervals. The slices show the edge-brightening apparent in the 23-epoch image (Figure 1) and 3D slice profiles (Figure 2). Additionally, these slices show that the transverse intensity profile is more complicated than double edge peaked, in particular between and  mas (2.808 mas; shown in Figure 5) and also between and  mas (bottom row in Figure 5). The structures responsible for these more complicated intensity profiles might be associated with filaments crossing the jet diagonally. Diagonal filaments are suggested by faintly visible structure in Figure 1, and are somewhat easier to identify in the 3D slice profiles of Figure 2 and in the contour plot underlying Figure 3.

Figure 5: Sample intensity slices across the jet (top row) at positions 2.106, 2.808, 3.510 & 4.212 mas and (bottom row) at positions 4.914, 5.616, 6.318 & 7.020 mas with respect to the core as shown in Figure 3.

The intensities associated with the bright northern and southern jet edge ridgelines along with the maximum intensity at each transverse slice position are shown in the top panel in Figure 6.

Figure 6: Top panel: Maximum intensity (solid line) and maximum intensity associated with the northern (red triangles) and southern (green squares) ridgelines are plotted along with vertical lines indicating the first location not dominated by the unresolved core and subsequently the locations of change in the intensity decline along the jet. Vertical lines along the counter-jet side are placed at the same core distance as the first three lines along the jet. Bottom panel: Full width at half maximum (FWHM) and the ridgeline peak separation (RLPS) in mas along with vertical lines transferred from the top panel are shown. Simple linear fits to the apparent jet full opening angle, , for both FWHM and RLPS are indicated by the dotted lines. Angular scale is indicated along the bottom axis, and intrinsic scale for the jet at to the line-of-sight is indicated along the top axis.

Along the jet, the maximum intensity is always associated with one of the bright ridgelines. In the innermost  mas the intensity profile is dominated by the unresolved core although even at  mas the intensity on the jet side significantly exceeds the counter-jet side. Vertical lines in the intensity plot are placed across the jet at the first location where the unresolved radio core does not significantly affect the intensity ( mas), and at larger distances where the intensity decline with distance changes slope significantly. Three vertical lines across the counter-jet side are placed at the same radio core distance as the first three lines along the jet. Along the jet and counter-jet sides of the core the second vertical line at  mas is associated with a significant reduction in the intensity decline and the third vertical line at  mas is associated with a significant increase in the intensity decline.

The bottom panel in Figure 6 shows that change in the slope of the intensity decline can be identified with change in the expansion rate along the jet and the counter-jet side. The bottom panel shows the full-width half maximum (FWHM) width, the ridgeline peak separation (RLPS) between the intensity peaks associated with the edge-brightened ridgelines, and simple linear fits to the apparent jet full opening angle, , along the jet where the apparent full opening angle is defined as

where and are the change in jet width and in core distance in mas. The symmetry in the opening angle between the jet and counter-jet sides along with the similar edge-brightening confirms the identification of the faint structure as the counter-jet. In both jet and counter-jet we see an initial rapid opening, , followed by a region in which the width decreases slightly, . Subsequently both jet and counter-jet increase in width rapidly with , and the counter-jet becomes too faint to follow farther. The obvious symmetry on either side of the radio core in the jet and counter-jet structure at our 0.117 mas slice separation means that the radio core is located no more than half the slice interval, i.e., one quarter the beamwidth, and  mas from the central engine. Thus the radio core is located at   from the central engine, consistent with the results of Hada et al. (2011).

Following the second rapid increase in width the jet opening angle decreases first to and subsequently to . Up to about 3.9 mas from the radio core both FWHM and RLPS widths change their slopes at about the same radio core distance and to within the uncertainties is about the same for FWHM and RLPS determined widths. This spatially synchronized behavior is broken beginning at  mas where the FWHM width increases rapidly out to  mas without a similar increase in the RLPS width. The RLPS width shows a similar rapid increase beginning at  mas extending to  mas. It is unclear whether these differences between the two width measurements are significant or just reflect uncertainty in the measurement methods, especially given the weakness of the northern ridge relative to the center, to the large brightness difference between the sides, and changes occurring along the southern envelope in this region. Beyond 5.8 mas the RLPS width indicates , similar to the parsec- and kpc-scale jet opening angle. The implications of these structures will be discussed in Section 4.1.

3.2 Rapid Structure Evolution

In this section, we present individual epoch images and movies made from them that allow detailed studies of the structural evolution of the M 87 jet near the core. The structure is complex and rapidly evolving, and is not easily described as a series of moving features. The most convincing way to see that there is significant outward motion is to view the movies which give a strong visual impression of such motion. Despite the difficulties, we perform a traditional component motion analysis for the 2007 and 2008 data and show the results here. For a full velocity analysis of our 2007 data using more sophisticated techniques that can identify multiple, overlapping, velocity fields, see MLWH.

Intensity Images

Figure 7: Contour/gray scale plots of observations 1-6 of those that were made of M 87 with the VLBA at 43 GHz during 2007 at about three week intervals. The convolving beam size used was mas elongated along position angle . The contour levels (mJy beam) are , , , 1, 1.2, 1.4, 1.7, 2, 2.4, 2.8, 3.4, 4, 4.8, 5.7, 6.7, 8, 9.5, 11.3, 13.4, 16, 19, 23, 27, 32, 38, 45, 91, 181, 362, and 724. The image off-source noise level is between about 0.17 and 0.25 mJy beam. Components, identified by eye, are marked with crosses. Numbers are for components that appear to be related between epochs based on proximity and other cues in the structure.
Figure 8: Contour/gray scale plots of observations 7-11 of those that were made of M 87 with the VLBA at 43 GHz during 2007 at about three week intervals. The 11 images were used to make the movie shown in online animation Figure 9. The final image (bottom right panel) is from late in 2007. It is included for completeness but other data sets near it in time have problems and are not yet fully reduced so this image from late in 2007 was not used in the movie, the stack, and most of the analysis. The convolving beam size, contour levels, and marked components for all of the images shown here are the same as in Figure 7.
Figure 9: For the offline versions of this paper, this figure is a colorized version of one frame from the 2007 movie. Online at the journal, or in astro-ph ancillary file M87_VLBA_2007.mp4, the figure is the movie, showing the first 11 images in sequence, with pauses proportional in time to the interval between epochs (typically three weeks). The underlying images are the same as those in Figures 7 and 8. The series of ticks along the bottom indicates elapsed time and the moving box indicates which epoch is being displayed. The movie shows the actual observed epochs because confusing artifacts appear when interpolating undersampled data. The movie provides a clear visual impression of outward motion, but also of a rapidly evolving pattern that makes it hard to follow individual features over long periods.

Figures 7 and 8 show the total intensity contour plots from the 11 epochs at about three week intervals between 2007 January 27 (MJD 54127) and 2007 August 26 (MJD 54339). The line of squares next to the dates gives a visual representation for the time of each epoch, especially if the images are viewed as a movie. The filled square representing the particular epoch shown in each panel is enlarged. Over this period the standard deviation of the 43 GHz flux from the core (inner 1.2 mas) was , not significantly above the noise and uncertainties. There are significant changes to the intensity along the edges and down the centerline of the jet. Significant change is also seen in the intensity of the counter-jet. Some of the variation is the result of instrumental effects but other variation is the result of true intensity changes that may be related to motions along the jet. True structural changes are large enough that it is difficult to connect features from epoch to epoch when sampled at three week intervals.

Figure 10: Contour plots of the total intensity of M 87 at 43 GHz for first 6 epochs at about 1 week intervals observed in 2008. This set starts on 2008 January 21 and ends on 2008 February 20. The beam dimensions, contour levels, and meaning of the marks are the same as for Figure 7.
Figure 11: Contour plots of the total intensity of M 87 at 43 GHz for the last 6 epochs at about 1 week intervals observed in 2008. This set starts on 2008 February 25 and ends on 2008 April 5. The beam dimensions, contour levels, and meaning of the marks are the same as for Figure 7.
Figure 12: For the offline versions of this paper, this figure is a colorized version of one frame from the movie made from the 2008 data. Online at the journal, or in astro-ph ancillary file M87_VLBA_2008.mp4, this figure is a movie showing the 12 images in sequence, with pauses proportional in time to the interval between epochs (typically a week). The underlying images are the same as those in Figures 10 and 11. The series of ticks along the bottom indicate elapsed time and the moving box indicates which epoch is being displayed. The data quality was not as good as in 2007 because tight observing date constraints precluded waiting for better weather. Three epochs were so degraded that the images were not included, as evidenced by gaps in the time sequence. Despite these degradations, the movie still gives a clear visual impression of outward motion.

A movie of the color image versions of the 2007 images is shown in online animation Figure 9. Offline, that figure is a sample frame from the movie. Interpolation is not used for the movie because the typical motion in 2007 between epochs is about 2.5 times the beam width and interpolation leads to distracting artifacts in the structure or the noise depending on how the interpolating is done. Viewing the movie is the best way to see the jet motions.

Observations continued between MJD 54363 and 54464 at three week intervals, but most have issues that complicate processing. Those issues have been partially overcome for the observation on 2007 October 26 (MJD 54406), which is shown as the last panel in Figure 8. This epoch is sufficiently isolated in time from the other epochs that it is not used in the time series analysis or the 23-image stack. It is included in the 49-image stack in Section 3.5 on the counter-jet. Technically this image is of interest because the observations were made at four frequencies (39.4, 41.3, 43.1, and 44.8 GHz) in an attempt to improve the UV coverage. The benefit of four frequencies was not fully realized because the two stations giving the longest baselines were degraded — one was out for maintenance, the other had poor weather.

In Figures 7 and 8, crosses indicate the visually-determined positions of local maxima in the total intensity. An attempt to relate these features from epoch to epoch was made by blinking rapidly back and forth between the epochs. The attempt is hampered, as noted above, by the large motions between epochs. But a variety of visual cues including the peak location, and also aspects such as major dips in the ridge line, made it possible to identify a significant number of components. Each feature that appeared related in three or more successive epochs was given a number. The motions of these components are analyzed in Section 3.2.2.

Figures 10 and 11 show total intensity contour plots of M 87 for 12 epochs at intervals of about 1 week from 2008 January 21 to 2008 April 05. These plots are similar to the images shown in Figures 7 and 8 for the 2007 data, although the observing interval is significantly shorter. Feature identification was done in the same manner as for the 2007 data and the results are included in our motion analysis. In principle, with the shorter interval, it should have been easier to relate features between epochs. But, as noted in Section 2, many of the images were of lower quality than those in 2007 so feature identification was as or more difficult. A movie of the color image versions of the 2008 images is shown in online animation Figure 12. Offline, that figure is a sample frame from the movie. Due to the lower quality of many of the images, the impression of motion is not as strong as with Figure 9, but the impression of motion is still present.

Serendipitously, TeV flaring occurred during the first two weeks of February 2008 and it is around this time that a significant long term increase in the core intensity at 43 GHz began. By April 2008 the core intensity had increased by 60% relative to the pre-TeV flaring average. Over the same time period the jet edges first brighten near to the core with brightening extending outwards over time. The increased brightening along the jet edges inside 1 mas can be seen by comparing the panel for January 21 in Figure 10, from before the TeV flaring, to the panel for April 5 in Figure 11. It is even more clear in the difference images shown in Figure 3 of Acciari et al. (2009).

Component Motion Analysis

Figures 7-8 and 10-11 show a jet structure that is complex and rapidly evolving.

Figure 13: Panel A: Positions of components that lie along the southern edge of the M 87 jet as indicated in Figures 7, 8, 10, and 11 as a function of time (For reference, 2007 January 1 is MJD 54101) are marked by circles whose area is proportional to the component peak flux density. Where three components have the same number in successive epochs, a fit was done for speed and position. Each such fit is represented by a line segment in the figure. Panel B: A histogram of the number of line segments in Panel A as a function of apparent speed . There are concentrations of speeds near zero, mostly from positions near the core, and concentrations at c, although intermediate and some faster speeds are also seen. Panel C: A plot of apparent speed with core distance (mid-point of line segment). The dotted line is the right boundary of the exclusion region resulting from the requirement that the first of the three points in a line segment be resolved from the core. The boundary shown is for the three-week cadence of the 2007 observations. As is apparent from looking at Panel A, higher velocities are at larger core distances. For distances less than about 2 mas, an acceleration region, as identified by MLWH, is seen.

The details of the structure, especially in the regions farther from the core, are somewhat uncertain because of the high dynamic range required to image the jet in the presence of the bright core, the complexity of the overall structure being imaged with limited UV coverage, and the unfortunate loss of one or more antennas in some epochs, especially during 2008. The fact that many of the marked features in the figures are not numbered is an indication that the structure is not easily described by long-lived discrete “components”. However, we can present a component motion analysis based on the marked and numbered features.

Figure 14: A figure identical to Figure 13, but using data from components along the northern edge of the jet. See the caption of Figure 13 for details. The speeds seen in both edges are discussed in the text along with comparison with results from MLWH. Compared to the southern edge, the northern edge shows more of a preponderance of intermediate speeds with a broad peak near . As along the southern edge, an acceleration region is seen at distances less than about 2 mas.

Figures 13 and 14 show the results from such an analysis of the features seen in Figures 7 through 11. The two figures show results for the southern and northern edges of the jet, respectively. There were insufficient features for a similar analysis along the jet centerline or in the counter-jet. In each figure, Panel A shows a time series of the positions of all the marked features. When three features have the same number in successive epochs the feature is assumed to correspond to the same moving component. The feature is then connected by a line segment where the line segment is the result of a linear fit to the position and apparent speed of the component. The other two panels are based on those line segment results. Panel B shows a histogram of fitted speeds for the line segments. Panel C shows the speed plotted as a function of core distance where the position of the mid point of a line segment is used for the core distance. There is an exclusion region for fast components near the core because the first feature closest to the core must be resolved from the core. The dotted line indicates the boundary of the exclusion region assuming that the first point must be 0.4 mas from the core and that the observing cadence is once per three weeks, which is appropriate for 2007. For 2008 with an observing cadence of once per week, the boundary of the exclusion region is steeper.

The data show a wide scatter in apparent speeds. The speeds range from slightly negative to near . Near the core, the speeds are slower with faster speeds at larger core distances. The histograms show concentrations near zero and near for both edges with an additional concentration near along the northern edge. The feature closest to the core is moving very slowly if at all, but that feature is generally just the inner part of the jet that is not well resolved from the core. Toward the end of the data set, in 2008, the innermost northern feature (#35) does move away more decisively. That feature appears to be related to the radio flare in the core region that occurred around the time of the TeV flares seen between MJD 54497 and 54509, and is shown in difference images in Acciari et al. (2009). Here this feature is seen to have accelerated to at  mas from the core. For features farther from the core, there is a general trend for speeds to increase out to  mas, reaching speeds albeit with large scatter. Simply averaging all the points more than 1.8 mas from the core for both years gives results for the north and south ridges of and with standard deviations (SD) of and , respectively. These averages correspond to 5.44 and 8.30 mas yr. But averages may be too simplistic if, for example, additional acceleration is involved.

Our results can be compared to those presented for the 2007 data by MLWH. That study used the Wavelet-based Image Segmentation and Evaluation (WISE) velocity field analysis combined with a stacked cross correlation analysis (Mertens & Lobanov, 2015) to analyze the velocity field in a manner that was able to pick out overlapping fast and slow motions. Velocities from the WISE analysis, as shown in Figure 3 of MLWH, show a wide scatter with predominantly slow speeds near the core and speeds over a wide range up to nearly beyond about 2 mas. Visual inspection doesn’t show clear evidence for multiple components. The velocity plots shown in our Figures 13 and 14 look similar in character. It was the stacked cross correlation analysis in MLWH that helped delineate the velocity structure of the jet. The analysis found that, at core distances between 0.5 and 1 mas, the motions are in the range of with the higher speeds in the north. Between 1 and 4 mas, there are two velocity components. The slower component has (north) and (south). The faster component has (north) and (south). Speeds were also determined along the center of the jet with values found between those of the northern and southern edges.

Figure 15: Magnetic field vectors overlaid on total intensity contours and a polarized intensity grey scale image for the 43 GHz VLBA observations made on 2007 January 27 (MJD 54127) and 2007 May 10 (MJD 54230). The contour levels are 1, 2, 4, 8, 16, 32, 64, 128, 256, and 512 mJy beam. The convolving beam, with a resolution of  mas, elongated along position angle , is shown at the lower right. The peak total flux densities are 666 and 660 mJy beam, respectively, while the peak polarized flux densities, at slightly different positions from the peak total flux densities, are 6.15 and 4.91 mJy beam, respectively. The percent polarizations at the positions of the peak polarized intensity are 1.5% and 1.1% at the two epochs and rises to about 4% at about 0.4 mas along the jet ridges. The cutoff levels for plotting of the polarization vectors are 0.5 mJy beam for the total intensity (2.5 and 3.1 times the RMS noise for the two epochs) and 1.0 mJy beam for polarized intensity (4.3 and 5.0 times the RMS noise of the Q and U images for the two epochs). Both epochs show a consistent polarization structure near the core. On the jet side, the magnetic vectors wrap around the core. There is a minimum in the polarized intensity near the peak of the total intensity where the polarzation direction jumps by about .

The stacked cross correlation analysis found a clear division at about between the fast and slow components. For comparison, using only 2007 data (the only data used by MLWH), our visually-determined components at core distances greater than 1.8 mas, have an average with SD of for the northern edge and with SD of for the southern edge. Similar averages for the WISE data in MLWH Figure 3 for all points beyond 1.8 mas are and . The WISE analysis clearly picks out more slow components than our visual analysis. Restricting the average to components that have individual speeds above to try to focus on the fast component seen by MLWH, the averages for our data are and while they are and for the MLWH Figure 3 data. Thus, the data sets give essentially the same result for the faster components. In fact, they are remarkably close given the limitations of the data and of our visual analysis method. For further discussion of the velocity results, the reader is referred to Section 4.2 and to MLWH.

3.3 Polarization Images

Images of the polarized emission of a radio jet can provide much information about the magnetic fields in the jet. For that reason, all of our observations included the necessary cross-hand correlations and calibration observations to allow imaging of the polarization. Difficulties were encountered with the calibration and a full polarization analysis has been deferred. But two epochs, 2007 January 27 and 2007 May 10, were successfully calibrated and imaged. The polarization structure within 3 mas of the core is shown in Figure 15 for those two epochs. Magnetic field vectors are shown under the assumption that they are rotated from the measured electric field vectors and that there is no significant relativistic aberration or rotation measure. These assumptions will need to be checked because rotation measures high enough to matter (greater than about 5000 radians m) have been observed at other positions and frequencies in M 87 (Zavala & Taylor, 2002; Kuo et al., 2014). It may be possible to extract rotation measures from our data but that analysis has not yet been completed. The magnetic field vectors are superimposed on total intensity contours and a grey scale image indicating the polarized intensity. The polarization percentage varies from near zero close to the core where the field direction flips, to about 4% about 0.4 mas from the core along the southern edge of the jet, and along the northern edge of the jet in the first epoch. Farther along the jet, and closer to the core along the northern edge in the second epoch, the polarized emission is too weak for reliable detection in these two data sets. The peak polarized flux densities (noted in the figure caption) are seen about 0.15 mas southwest of the core. At those locations, the polarization percentages are 1.5% and 1.1% at the first and second epochs. On the counter-jet (east) side of the core, within the wings of the beam, the polarization is 1% to 2%.

The magnetic field orientation shows a consistent pattern at the two epochs. This might be expected as there were no large changes to the intensity structure during this period. On the jet side of the core the magnetic field vectors are approximately perpendicular to the local jet axis. Rotation of the vectors between north and south edges reflects the wide jet opening angle. On the counter-jet side of the core, the magnetic field vectors are approximately parallel to the counter-jet axis, again rotation between north and south edges reflects the wide opening angle. While the polarization angle change near the core may reflect the magnetic field structure, it is also possible to get such an effect in VLBI sources as a result of the transition from optically thick to optically thin emission. A related opacity effect is expected to displace the radio core from the actual black hole location. Such a shift for 43 GHz in M 87 was measured by Hada et al. (2011) to be as. Given that the polarized emission on the counter-jet side of the core is within the wings of the convolving beam, it is possible that it originates in the portion of the jet between the black hole and the radio core, i.e., not from the counter-jet. As noted above, the polarized intensity drops rapidly away from the core but the percentage polarization rises over the short distance for which it is measured reliably. There is some suggestion that the magnetic field vectors are more flow aligned along the northern and southern edges of the jet at  mas from the core, although this result needs to be confirmed with lower noise polarization data.

Magnetic field vectors that are wrapped around the core, reflecting the initial large jet opening angle, suggest a toroidal jet magnetic field geometry. However, modeling is needed to take into account the relatively small viewing angle to the line of sight, along with possible optical depth, Faraday rotation, and relativistic aberration effects. The polarization data from our observations that have not yet been reduced should provide much higher sensitivity, allowing detection of polarization farther from the core. The data taken since 2009 will provide 24 GHz data to allow measurement of Faraday rotation. These improvements should eventually provide more information on the magnetic field structure.

3.4 Long-Term Evolution

This project began as an effort to measure possible superluminal motions near the core of M 87 by sampling significantly faster than had been done in previous projects. But slower, longer-term sampling has the potential to provide information on aspects of the jet behavior other than component speeds. For example, relatively slow changes in the jet envelope or direction could be observed, providing information about jet propagation and the external medium. As described in Section 2, we have, somewhat fortuitously, roughly annual observations covering 17 years. Prior to our pilot project, a few archival observations done for other reasons are available. Between 2009 and 2016, annual observations were made as part of an effort to observe radio counterparts to TeV flares. In this section, we first present the images from these roughly annual observations and subsequently determine displacement motions from these images. A clear finding is that the jet moves transversely on time scales of several years.

Intensity Images

The total intensity evolution of the sub-parsec scale radio jet in M 87, spanning a 17 year period from 1999 to 2016, is shown in Figures 16, 17, and 18. Approximately annual observations at 43 GHz involving the VLBA are shown.

Figure 16: First five previously published epochs between 1999 to 2004 (Junor, Biretta & Livio, 1999; Ly, Walker & Junor, 2007) of the roughly annual observations of M 87 made between 1999 and 2016. All images have been rotated by . The white lines indicating the jet edges on the 2007 image are there to aid visual detection of changes at different epochs. The images in Figures 16, 17, and 18 are available as a movie in online animated Figure 19.
Figure 17: Middle six epochs of the roughly annual observations of M 87 made between 1999 and 2016. The 2006.3, 2007.4, 2008.1, and 2010.3 images are actually stacks (noise-weighted means) of 6, 11, 12, and 6 images, respectively, taken within a few months of the indicated mean time. As in Figure 16, the white lines indicate the jet edges on the 2007 intensity image. Recall that the stacks will smear out variable structure and emphasize persistent structures. The images in Figures 16, 17, and 18 are available as a movie in online animated Figure 19.

Where multiple observations were made in less than a year, average images are shown. Image quality improves significantly with time as the VLBA hardware improved and as more average images became available. See Table A2 for details of the observations, including the bit rate and the off-source image noise.

Figure 18: Final five epochs of the roughly annual observations of M 87 made between 1999 and 2016. This group consists of single-epoch images (not stacks). The final 4 are based on observations made after the VLBA bandwidth upgrade. With typically 8 times the bandwidth or nearly three times the sensitivity of the older observations, the individual epochs have noise levels comparable to the stacks in Figure 17. As single-epoch images, transient details are not smeared. Again, the white lines indicate the jet edges on the 2007 intensity image. The images in Figures 16, 17, and 18 are available as a movie in online animated Figure 19.

The images in Figure 16 are early observations and have limited sensitivity. In some of the very early observations, the actual observing time was short because M 87 was a phase calibrator, not the primary target of the archival observations shown here. Average images are available from the pilot project in 2006, the main monitoring observations in 2007 and 2008, and from the trigger response in 2010.

Significant transverse displacement of the jet, especially between about 2 and 8 mas from the core, is seen in Figures 16, 17, and 18, particularly in the epochs after 2011.1 relative to earlier epochs. Lines have been drawn on the figures that follow the bright jet edges in the 2007.4 average image. These lines help make the northern shift of the jet clear, especially in the later epochs. The transverse displacements are discussed more extensively below.

Figures 16, 17, and 18 also show that the relative brightness of southern and northern edges experiences significant changes over the 17 year period in the innermost few mas. There are epochs when the northern side is brighter, and epochs when the southern side is brighter. Which is brighter can be different as a function of position along the jet. It does appear that, when there is a transverse shift, the side toward which the shift happens brightens in the innermost mas. The best case is for the large shift to the north which started in 2011.1. For that epoch, the northern side is brighter by a factor of 2 to 3 near 2 mas. The northern side remains brighter as the shift propagates along the jet, but, as the jet starts to return toward its original position, the southern side becomes brighter.

Displacement Motion

The results of an effort to measure the propagation speed of the jet offset are shown in Figures 20 and 21. The transverse offset of the measured center of the jet relative to a nominal jet center line, that extends at position angle from the core, is shown in Figure 20 for core distances at 2 to 8 mas at 1 mas increments (the offset is actually in the y direction in the rotated images). The offsets show a general drift north superimposed on what appears as a quasi-sinusoidal variation. The patterns for each distance from the core show a lag, indicating a propagation speed that is addressed in this section. For core distances at and beyond 3 mas, the offsets were measured from images that had been blurred along the jet direction (convolved with a beam of mas elongated along ) to reduce effects of fine scale structure. At 2 mas, the full resolution images were used. Slices perpendicular to the jet were made at the core distances noted.

Figure 19: In the online version of this paper at the journal, or in astro-ph ancillary file M87_VLBA_Long_Term.mp4, this figure is a movie made from the roughly annual images shown in Figures 16, 17, and 18. It gives a visual sense of the changes that occurred in the source between 1999 and 2016, including the side-to-side displacements. Those shifts are described in detail in Section 3.4.2 in terms of a linear shift and a sine wave, both functions of time, both of which propagate down the jet at a pattern speed. A fit to the offsets of the jet center from a nominal center line gives the parameters of these patterns and their pattern speed, which is fast, but not as fast as the faster component motions. A pattern speed is also derived using a correlation analysis giving a similar result. For the offline versions of this paper, the image from 2013.0 is shown.

The offsets of the jet center from the PA line were determined from the slices by the following algorithm. First the peak intensity south of the PA line was determined. Then, the locations of the most southern interpolated points at one quarter and one half of the that peak intensity were found. The average of those two points was treated as the position of the southern side and one half of the difference was treated as an “error bar”. The process was repeated for the northern side using the most northern points at a quarter and a half of the peak value north of the PA line. The separate determinations of the north and south peaks prevents spurious results when one ridge is significantly brighter than the other — a frequent case. The center was taken to be the average position of the northern and southern side points and the number used for weighting in the analysis was the quadratic sum of the “error bars”. In practice, these “error bars” did not vary much between points and, to avoid clutter, only those for the 2 mas and 8 mas points are shown in Figure 20.

Figure 20: Measured offsets of the M 87 jet center from a line extending from the core along position angle of . The data used are the images shown in Figures 16, 17, and 18 with one extra observation (2016 April 10 in Table A2) included at the end. Offset measurements are based on finding the outermost half and quarter power points on each side, relative to the peak on that side. “Error” bars are based on the separation of the two points on each side and do not vary much so only those for the points at 2 and 8 mas are shown. Dashed lines are the results, for different core distances as indicated by color, of a least squares fit for a model, described in the text, whose main components are a linear relation and a sine function, both with time, both propagating down the jet. The “errors” were used as the weights in the fit. The apparent propagation speed, with its formal error from the fit, is . The bottom panel shows the residual data which are the measured data points with the fitted model subtracted.

Two methods were used to determine a speed of propagation of the transverse displacement pattern. The first method is a least squares fit for the parameters of a simple equation that roughly describes the data in Figure 20. The purely empirical equation is not based on any particular physical model and is given by

(1)

where is the offset from the nominal centerline in mas, is the date in decimal years, is the number of years since 2000.0 (to reduce the chances of numerical problems), and is the time a feature seen at time at core distance was at the core assuming linear motion from the core at speed . The form of used to parameterize the equation is not meant to imply that linear pattern motion occurs or exists very close to the core. The fitted parameter results for , , , , , , , and are given in Table 1. The offset, , for each core distance given by Equation 1 is shown as dashed lines in Figure 20 with the colors matching the data colors. The bottom panel shows the residuals after subtracting the model from the data.

Parameter Description Fit value Std. Error Units
Offset at , 0.024 mas
Slope with distance 0.005 0.004 mas offset per mas core distance
Linear change with time 0.022 0.002 mas yr
Angular pattern speed 3.362 0.68 mas yr
Sine wave amplitude at 0.061 0.017 mas
Index of sine wave amp dependence on 0.69 0.16
Sine wave frequency. 0.0973 0.003 Cycles mas
Phase at 31 15 Degrees
Table 1: Fitted Parameters of the Jet Displacement Model

A linear change in offset with core distance at a fixed time can be created in at least two ways. First, an error in the choice of the position angle assumed for the jet center relative to the core would produce such an effect because the resulting linear offset would grow with core distance. Second, if the offset changes with time at a constant rate, and that change propagates down the jet at finite speed, the same linear change in offset with distance will be produced. Both effects are included as terms in Equation 1. This would not normally work for the fit because the terms are completely correlated. However, we assume that the pattern propagation speeds for the linear change in the jet center offset and for the sine wave are the same. This causes the pattern propagation speed to be controlled by the sine wave, leaving the error in the choice of jet position angle to absorb any additional linear change in offset with distance. This assumption decorrelates the terms. The fitted result that the linear change with distance () is very small suggests that the assumed position angle was properly chosen for the reference epoch (2000.0).

This first method gives a period result of years within the errors. However visual inspection of Figure 20, concentrating on the 2 and 3 mas data points, suggests a shorter period. These two distances are the only ones measured from images before 2006 and have an extra seven-year time baseline. Indeed, a fit to just the 2 and 3 mas data gives a period result of years. Nevertheless, adding the data from farther down the jet strongly discourages the shorter period. The differences seem to result from the 2006 and 2016 data — the extremes of the time range for the more distant data. As a check of the 2016 March 14 data, and an indication of the reliability of the rest of the data, a 2016 April 10 epoch (the result of a high energy trigger and included at the end of Table A2), not shown in Figure 18, was included for the fit. The results using this 2016 April 10 epoch show good agreement with those ending with the 2016 March 14 epoch. A few more years of good data, covering more than a full period, should provide some clarification. In any event, this method gives an apparent pattern speed of mas yr for . At a viewing angle the apparent pattern speed would imply an intrinsic pattern speed of .

Figure 21: The correlation functions used to determine the propagation speed of the side-to-side variations in the M 87 jet center position. These provide an alternate determination to that from the least squares fit shown in Figure 20. The data used are the subset of those shown in Figure 20 which were taken between 2006 and 2016. Data for each core distance, measured at 1 mas intervals in core distance, were correlated with that of each other distance as a function of time lag. The time axis is scaled by dividing by the separation in core distance of the pair so that the axis becomes the propagation time per mas (years mas). The correlation function for each pair of core distances is identified by a unique combination of color and line type as shown in the legend. For most pairs, the line type indicates the initial core distance of each pair, while the color indicates the separation in core distance of the pair. The  mas pair does not conform to the pattern because the plotting software used provided too few line types. The average over all the correlation curves is shown by the bold black line. The weighted average of the peaks of the individual curves is  year mas while the peak of the average curve is at year mas.

In the second method used to measure the displacement pattern speed, for each core distance, the measured data points at the different epochs were linearly interpolated to a dense time series and each pair of time series (different distances) was cross correlated. This was done by calculating the “Pearson’s r” statistic for time lags between and 5 years. The normalization is done separately for each lag to try to minimize the impact of end effects since the maximum lags are a noticeable fraction of the overall range. The correlations should peak at a time lag related to the pattern propagation speed. This method of measuring the propagation speed does not assume any particular functional form for the offsets so it is insensitive to the time scale of any oscillations. The correlation functions are plotted in Figure 21. The correlations are plotted against time lag divided by the core-distance separation of the two time series being correlated. This turns the time lag axis into units of propagation time per mas, which is an inverse speed. The actual speed is not a natural unit for plotting the correlation functions because there is an infinity and sign flip in speed at zero lag. The correlations were then averaged as a function of the inverse speed with the result shown by the bold black line in Figure 21. The peak of the average curve is at 0.350 year mas which is 2.86 mas yr or . An average speed was also calculated by measuring the peaks of each correlation function, again scaled by the distance, and doing a weighted mean. For the weights, the width to the 90% amplitude points were used. That gives reasonable relative weights and does not significantly affect the calculated standard error which takes into account the scatter of the points. The result, determined in this manner, is  year mas which is mas yr or . At a viewing angle the apparent pattern speed (weighted average) would imply an intrinsic pattern speed of .

The speed measured from the correlation analysis () is slower than, but within the errors of, the speed measured with the least squares fitting (). This match between the speeds from the two methods lends credence to the results. In either case, the speed is significantly less than the component speeds measured in the component motion analysis ( for the fast components). Thus the transverse shift of the jet is not propagating down the jet ballistically. Another indication that the offset is not ballistic is that it does not grow linearly with distance. The amplitude deviation of the fitted sine wave from the nominal axis grows with distance by only a power-law index of .

3.5 The Counter-Jet

Essentially all of the images presented in this paper show evidence for a counter-jet. A counter-jet can provide important clues to the nature of the jet. Firstly, its mere existence shows that the M 87 jet is two-sided as it is launched, which is commonly assumed, but rarely proven in superluminal sources. This provides evidence against flip-flop models (e.g., Rudnick & Edgar, 1984) where the jet is one-sided at any given time. Secondly, the ratio of brightness of the jet and counter-jet, if it is assumed to be the result of Doppler boosting, gives a measure of the jet speed. Thirdly, similarities or differences in the structure of the jet and counter-jet help distinguish features related to fundamental dynamical structures in the jet and surrounding medium from more ephemeral, time-dependent structures based on details of the stability of the jet generation mechanism.

While all of our images show the counter-jet, those based on data that predates the recent upgrades to the VLBA did not reliably delineate the structure of the counter-jet. The counter-jet feature is weak and located in the region near the bright core where artifacts from residual calibration errors are most pronounced. Stacked images showed that it has a structure roughly symmetric with the main jet and some of the recent, high-sensitivity, observations show the structure more clearly. Here we first push the sensitivity at the normal resolution as far as possible with our data by stacking all of the images, with weighting based on the off-source noise. Subsequently, we explore the structure close to the core by adjusting imaging parameters to emphasize resolution in the creation of images from two of the best high-sensitivity epochs.

An image based on stacking 49 images from the epochs listed in Table A2 in the appendix (excluding the last epoch), is shown in Figure 22.

Figure 22: A noise-weighted average image made from the first 49 epochs of the 50 epochs listed in Table A2. The image shows the average counter-jet structure at 43 GHz as far from the core as is allowed by our data. The convolving beam was the usual  mas, elongated along position angle . The contour levels, in mJy beam, are , 0.2, 0.4, 5.6, 8, 11.3, 16, increasing from there by factors of . The peak is 726 mJy beam and the off-source RMS noise level is 47 Jy beam. Blue lines (jet side) and red lines (counterjet side) show the locations where intensity measurements were made to give the sidedness data shown in Figure 25. The lines on the jet side follow the ridges while those on the counter-jet side are the jet-side lines reflected through the core.

Only the inner few mas are shown to emphasize the counter-jet and because the side-to-side variations discussed in Section 3.4 blur the otherwise sharp-edged structure on larger scales. On these small scales, the side-to-side variation-induced blurring is relatively small compared to the resolution. This image allows measurement of the sidedness intensity ratio to  mas from the core, provides cross jet resolution and shows the edge-brightening similar to that seen on the jet side. The blue (jet side) and red (counterjet side) lines crossing the core in Figure 22 show the locations of the pixels whose intensities are used to measure the sidedness that will be discussed in Section 4.2.2. These pixels are chosen to follow the north and south ridges on the jet side. On the counter-jet side, each pixel chosen is simply at the symmetrically opposite position from the corresponding pixel on the jet side. This builds in an assumption that the jet/counter-jet structure is symmetric which the image shows is close to, but not precisely, true. Recall that, as with any stack spanning times longer than the time scales for changes in the source, the detailed changing structure is smoothed out, leaving the persistent structure.

Our recent higher-sensitivity, wide-bandwidth observations have provided an opportunity to delineate the structure of the counter-jet at individual epochs in more detail than was previously possible.

Figure 23: Mild super-resolution 2013 and 2016 images using the wide bandwidth system on the VLBA. The convolving beam of mas, elongated north-south, is shown in the lower right corner. It is % the size of the fitted beam in the north-south direction and % in the east-west direction. The contour levels, in mJy beam, are , , , 1, 2, 2.8 and increase from there by factors of . Jagged blue (jet side) and red (counterjet side) lines connect pixel center points where intensity measurements were made for the sidedness data used in Figure 25. Green lines on the counter-jet side in the 2013 image connect counter-jet points plotted at 90% of the core distance of corresponding jet side points. Those green lines are dotted where they overlap the red lines. Note the symmetry between the jet and counter-jet, especially in 2013. At very low levels there are residual calibration artifacts near the core, particularly in the 2016 image.

To examine the counter-jet structure close to the core with two of the best data sets, 2013 January 12 and 2016 March 16, we have used imaging parameters that lead to higher resolution than our other images at the cost of degraded image quality on larger scales. This was accomplished by using a “robustness” (Briggs, 1995) approximating uniform weighting. Uniform weighting emphasizes the long baselines even though they are a relatively small fraction of the data. Then the point source CLEAN components were convolved with a beam smaller than that fitted to the dirty beam. This “super-resolution” increases the noise level and can be risky if overdone. Here the fitted beams were mas elongated in PA in 2013 and mas elongated in PA in 2016. The final convolving beam used was mas in PA . Thus the super-resolution reduced the beam area by about 50%, mostly in the north-south direction. The peak flux densities are 0.62 and 0.48 Jy beam in 2013 and 2016, respectively.

The resulting fairly mild super-resolution images of the inner jet from the two data sets are shown in Figure 23. The images show a structure in which the counter-jet is a reflection of the main jet, but at lower intensity. The bright edges of the main jet that exit the core at an apparent opening angle of about have corresponding bright edges on the counter-jet side. As with the stacked image, the blue lines (jet side) and red lines (counderjet side) show the locations where measurements of intensity along the jet were made for the sidedness data discussed in Section 4.2.2. The blue and red lines follow the jet ridge line on the jet side and are reflected across the core for the counter-jet side. They are jagged because measurements were made at pixel centers.

The northwest jet edge and the southeast counter-jet edge, especially on 2013 January 12, show a distinct, sharp bend in both jet and counter-jet at about 0.44 mas from the core. Such a feature is not common in our other images, so its appearance in both jet and counter-jet suggests a common origin. If the bend were moving outwards at a significant fraction of the speed of light, one would expect to see the bend closer to the core on the counter-jet side because of differential light travel time. That allows limits to be placed on the propagation speed of the bend. The green lines in the top panel in Figure 23 (2013 January 12 epoch) show a % core-distance reduction compared to the red lines (which are symmetric with the blue lines on the jet side). This small reduction does appear to be a somewhat better match to the counter-jet ridge line. A 10% difference requires that the speed associated with bend motion be no more than about . This inner jet and counter-jet region is where the slowest moving components are seen (see Figures 13 and 14). On the other hand, there is a significant sidedness intensity ratio in this region (see Section 4.2.2) which suggests higher material speeds. Thus, we conclude that material is moving through the bend.

The results of the sidedness measurements and their implications are discussed in Section 4.2.2.

4 Flow Dynamics and Emission

In this section, we explore what can be learned about the flow within the jet and about the structure of the jet based on the observational results of Section 3. In particular:

Section 4.1 examines the regions of expansion and recollimation of the jet envelope in the inner 7 mas of the jet, shown previously in Section 3.1 based on the 23-epoch average. Here we provide intrinsic opening angle values based those data, on the 49-epoch average image data of Section 3.1 and on the images presented to study the counter-jet in Section 3.5. This information can be used to provide an indication of the interaction between the jet and the external medium. It is distinct from Section 4.3 which is concerned with the internal flow structure of the jet.

Section 4.2 addresses the implications of the proper motion measurements and of the ratio of the intensities of the jet and counter-jet for jet flow speed and acceleration. Additionally, the different proper motions and intensities along north and south sides of the jet and counter-jet are used to explore poloidal versus helical flow models.

Section 4.3 addresses implications for the internal jet spine/sheath structure from the observed slice profiles, and discusses potential differential Doppler boosting effects resulting from helical flow. These issues are distinct from Section 4.1 which addresses the jet envelope.

Section 4.4 explores the implications of propagation down the jet of the long-term changes in jet transverse position. Here we explore the constraints placed on the jet to external medium interface that would lead to the observed propagation assuming that non-ballistic propagation can be described by the Kelvin-Helmholtz helical mode.

4.1 Expansion & Recollimation

Here the shape of the jet, especially regions of expansion and recollimation, is explored. That shape can provide constraints on the interaction of the jet and the external medium. In Figure 24 we redisplay the FWHM and RLPS widths from Figure 6, which are based on the 23-image stack of 2007 and 2008 data, and average the two to give a more reliable indicator of the behavior of the jet width as a function of the intrinsic distance in units of . As implied in the earlier discussion of Figure 6, there are sufficient uncertainties in the separate width measurements, especially in the region of greatest differences between about 5 and 6 mas, that we will not attempt to interpret those differences in terms of jet physics. Linear fits to the average width, indicated by the dotted line, are made between the vertical lines which are shifted up to  mas relative to the locations indicated in Figure 6 in order to better accommodate the linear fitting. The resulting intrinsic full opening angle associated with the linear fits is defined as where the viewing angle and and are the change in apparent width and apparent distance. Beyond 7.4 mas the dotted line in Figure 24 indicates the intrinsic opening angle of the jet beyond HST-1 from Owen, Hardee & Cornwell (1989). Our results suggest three collimation regimes: (1) an initial rapid expansion out to () followed by jet contraction over a distance (), (2) a second equally rapid expansion beginning at () with recollimation occurring over a distance (), and (3) a third slower expansion beginning at () with recollimation still occurring at () where intrinsic distances in assume a viewing angle of .

Figure 24: Average (AVGW) of Full width at half maximum (FWHM) and the ridgeline peak separation (RLPS) determined widths (cyan dash–dotted line) along with linear fits to the intrinsic opening angle, (blue dotted line), between the vertical lines from the 23-epoch image. First and third vertical lines are shifted by  mas and  mas relative to the locations indicated in Figure 6 in order to better accomodate the linear fitting. Numerical labels show the intrinsic full opening angle associated with the linear fits between the vertical lines. Beyond 7.4 mas the dotted line and label indicate the intrinsic opening angle of the kpc-scale jet (see Section 1). The width along the right axis and the intrinsic core distance along the upper axis for the jet oriented at to the line of sight are given in units of .

On average our expansion is consistent with the parabolic structure out to about HST-1 found by Asada & Nakamura (2012). Our results confirm the initial rapid parabolic expansion found by Hada et al. (2013) inside . Our results also show that the structure consists of a subsequent contraction followed by additional expansions and recollimations beyond that is more complicated than the slower single parabolic expansion beyond a few hundred found in Hada et al. (2013).

The 49-epoch image (Figure 22), that shows the jet and counter-jet and is based on 17 years of data, also reveals different opening angle regimes for the jet. Motions blur the structure in the 49-epoch image but the image still reveals a high initial opening angle near the core, followed by a narrower opening angle indicating slower expansion with distance. If we use the blue lines to determine the opening angle, we find inside 0.620 mas and beyond 0.620 mas and out to  mas. The corresponding intrinsic opening angles are and , respectively.

The mild super-resolution single epoch 2013 and 2016 images (Figure 23) that show the inner jet and counter-jet structures provide a more detailed view of the inner structure not blurred by motions. Now it is possible to see a rapid initial expansion, followed by a slower expansion in both jet and counter-jet. In the 2013 image the blue lines indicate inside 0.5 mas and between and  mas. In fact the initial opening angle inside 0.5 mas must be even larger: we find an initial opening angle of if we connect the location of the sharp bend in the lines along the south counter-jet and north jet to the core position. The corresponding intrinsic opening angles are inside 0.5 mas and between and  mas. In the 2016 image the blue lines indicate inside 0.4 mas and between and  mas. Again if we connect the core position to the sharp bend in the south counter-jet and north jet we find an initial opening angle of . The corresponding intrinsic opening angles are inside 0.4 mas and between and  mas.

At even smaller scales than we resolve, Hada et al. (2016), using data at 86 GHz, find a region of jet expansion with apparent opening angle of inside 0.2 mas. They also find a subsequent contraction at a core distance of about mas. At a viewing angle of , the corresponding intrinsic opening angle is inside 0.2 mas.

The multiple expansion and collimation regions indicate a jet launching region in which equilibrium parabolic expansion is achieved only after 5 to 7 mas, i.e., several thousand . The non-equilibrium behavior found here and at smaller scales by Hada et al. (2016) should provide important information for theoretical/numerical modeling of the jet launching region. For example, we may be seeing the interaction between a black hole launched jet spine, velocity shearing sheath region and disk wind cocooning medium (see also Hada et al., 2016) as well as the extent of the jet launching and acceleration region. We note that the initial intrinsic opening angle near the jet base is such that a portion of the jet flow comes directly into the line-of-sight at a viewing angle of to the jet axis. That some jet flow can be directly into the line-of-sight near to the jet base has consequences for models of the location and mechanism responsible for the TeV and associated radio flaring observed from M 87.

4.2 Jet Flow

We now turn to the flow dynamics of the jet as revealed by the component motion measurements and by the observed sidedness ratios described in Section 3.2. The component motions represent the material flow speed if their origin is some feature of the underlying flow, such as a density enhancement, that moves with the jet material. They can also be slower than the underlying material if they are a pattern speed, such as might occur with instabilities, shocks, or interactions with external influences of some sort. The components are unlikely to be faster than the material speed. The sidedness ratio is most likely the result of Doppler boosting and is a direct consequence of the material flow speed. A sidedness ratio higher than implied by the true material speed would require some asymmetry in the physical parameters of the jet/counter-jet system. A sidedness ratio lower than expected for the faster material motions also requires some complexity to the jet, such as the presence of some slower material that is contributing some of the emission. In the presence of Doppler boosting, which is assumed when there are superluminal apparent motions or high sidedness ratios, there can also be beaming-induced brightness differences between similar regions with different flow angles. In this section, the speeds implied by the component motions and sidedness are explored. Differences in speed and brightness between the sides of the jet are interpreted in terms of the possible flow patterns such as a helical flow. This discussion complements the extensive proper motion study along the jet using the WISE analysis technique that can be found in MLWH.

Proper Motion Implications

Apparent jet motions obtained from the wavelet analysis given in MLWH Tables 2 and 3 indicate averages of (north) and (south) between 0.5 and 1 mas. Assuming pure poloidal flow along the jet edges, this would indicate intrinsic speeds of (north) and (south) at a viewing angle of . Between 1 and 4 mas there are two velocity components. The slower component has (north) and (south). The faster component has (north) and (south). Following the same viewing angle and flow assumptions, the slower velocity component would indicate (north) and (south), and the faster velocity component would indicate (north) and (south).

It is hard to understand the velocity differences between opposite jet edges found assuming a purely poloidal flow. On the other hand, different values along opposite jet edges can naturally result from flow rotation about the jet axis. With rotation, the flow on each edge will be at different angles to the line-of-sight. Faster apparent motions will be seen along the edge where the flow is at an angle to the line-of-sight that comes closer to the angle at which apparent motions are maximized.

The critical angle at which apparent motions are maximized, , for is . This angle is greater than the poloidal flow angle to the line-of-sight along the jet edges, which is slightly larger than the jet axis angle to the line-of-sight (see Equation 4). Thus, the higher apparent speeds along the northern jet edge between 0.5 and 1 mas and the higher apparent speeds along the northern jet edge of the faster component between 1 and 4 mas would indicate clockwise flow rotation about the jet axis. Here the “clockwise” terminology for the sense of rotation is for an observer looking towards M 87 along the jet (not counter-jet) axis and thus remains the same for a jet and counter-jet that share a common source of rotation. Such rotation places the northern edge flow closer to the critical angle. The higher apparent speeds along the southern edge of the jet for the slower component between 1 and 4 mas would indicate counter-clockwise flow rotation for that component.

If the faster component motion along the northern edge is interpreted as resulting from clockwise rotation, then the speed of the implied clockwise helical flow can be given approximately by

(2)

where and are the apparent north (N) and south (S) edge motions, and provides an estimate of the difference in flow angle relative to pure poloidal flow along the jet edges at line-of-sight angle . Here indicates clockwise flow rotation, and we have made the simplifying assumption that the flow angle to the line-of-sight can be expressed as . A little manipulation allows us to solve for which is found from

(3)

In Equations (2) & (3), is the viewing angle to the jet edges and is related to the viewing angle to the jet axis by

(4)

where here is the viewing angle to the jet axis and is the intrinsic half opening angle.

Between 0.5 and 1 mas the intrinsic half opening angle declines from about to about (see Figure 24). Here we use an average of so the viewing angle to the jet edges is . Now using the apparent edge motions reported above from MLWH so that and , we find that

provides an indication of intrinsic clockwise helical flow less than 1 mas from the core. The likely range of values was deterimined by offsetting and , separately, by plus and minus one sigma, then taking the geometric sum of the positive changes and of the negative changes. The results were sufficiently symmetric to use an average absolute offset as the quoted error.

Between 1 and 4 mas, declines from to (see Figure 24) but is at the smaller angle over two thirds of this range. Here we use a distance-weighted average of so that the viewing angle to the jet edges is . Now using the faster apparent component motions from MLWH so that and , we find that

provides an indication of intrinsic clockwise helical flow a few mas from the core. Thus, in a flow rotation scenario we find evidence for a clockwise helical flow accelerating from inside 1 mas to at larger distances with the flow becoming more poloidal at larger distances.

For the slower-velocity components from MLWH, if the edge motions, and , are interpreted as resulting from counter-clockwise rotation, then the flow angle and speed of the implied counter-clockwise helical flow beyond 1 mas is and . Here we note that angles closer to the poloidal case () result in a smaller value for the speed . Our results for the fast and slow component flow angles allow for poloidal flow of both components within less than 1.5 times the errors. Nevertheless, our results suggest the possibility of a complex velocity shearing region with opposite sense of rotation for the fast and slow components found by MLWH.

Jet/Counter-jet Intensity Ratio Implications

The ratio of the intensities of the jet and counter-jet provides information about the jet flow, under an assumption of an inherently symmetric system, that is independent of the component motions. The average image of M 87 shown in Figure 22 and in the super-resolution images of M 87 shown in Figure 23 have been used to obtain intensity ratios between counter-jet and jet edges that are symmetrically opposite the core. The resulting intensity ratios are shown in the upper panels in Figure 25 where the data points are at image pixel positions along the blue and red lines shown in the images. The ratios, and are between intensities at positions on the counter-jet (cj) side and at positions on the main jet (j) side that are symmetrically opposite the core. The positions were chosen to be along the northern (N) and southern (S) ridges on the jet side. The structures on the counter-jet side did not determine the positions, although the images show that the symmetry assumption is good but not exact. In Figure 25 the ratio , which is the intensity ratio choosing the maximum counter-jet and jet edge intensities at each core distance independent of symmetry considerations, is included for comparison purposes. Recall that the intensity ratios are dominated by the beam until the jet and counter-jet features are resolved from the core.

Figure 25: Upper panels: Counter-jet/jet intensity ratios at image pixel positions along the blue and red lines in the M 87 images of Figures 22 and 23. Each set is the ratio between intensities at positions on the counter-jet (cj) side and the main jet (j) side that are symmetrically opposite the core. For comparison a black line indicates the intensity ratio choosing the maximum counter-jet and jet edge intensities at each position. The positions were chosen to be along the northern (N) and southern (S) ridges on the jet side. The shaded region indicates where south counter-jet and north jet ratios are dominated by the image resolution and core. The intrinsic scale in units of () along the jet is shown at the top. Middle panels: Intrinsic jet speed implied by the ratios shown in the upper panel at the viewing angle appropriate to the jet edges, . An angle to the line-of-sight of the jet axis and a spectral index of are assumed. The vertical lines indicate the location where the jet opening angle changes significantly, seen mainly as a kink in the lines along the northern jet edge. Bottom panels: Intrinsic flow speed, , and flow angle, , relative to pure poloidal motion along the jet edges.

Since the beam is elongated more along the inner south ridge on the counter-jet side and north ridge on the jet side, these intensity ratios will remain high farther from the core, as is evident in all the intensity ratio panels. For example, the 49-epoch ratio stops its steep core-dominated decline at  mas whereas the ratio stops its steep core dominated decline at  mas, and that the 2013 ratio stops its steep core dominated decline at  mas whereas the ratio stops its steep core dominated decline at  mas. Also the 2016 and ratios do not appear beyond core dominated influence until  mas. The region of core influence is shaded in the figure panels. Beyond the region of core influence, the fact that the regions sampled are symmetric about the core, as is the beam, should ensure that the ratios are not significantly sensitive to the beam shape.

The persistence of larger ratios between the south counter-jet and north jet edges, , than between the north counter-jet and south jet edges, , can be interpreted in terms of Doppler boosting and suggests that the northern side of the counter-jet and jet is deboosted relative to the south side of the counter-jet and jet. In Figure 25 the middle panels show the computed values of a purely poloidal speed along the north counter-jet edge and south jet edge, , and along the south counter-jet edge and north jet edge, , required by the intensity ratios. Speeds are computed from

(5)

where, as previously stated, Values are computed using the viewing angle to the jet axis , the intrinsic half opening angle , and the spectral index (; appropriate to the kpc-scale jet in the radio band; Frazer Owen, private communication). In the panels the vertical line indicates the location where the opening angle changes. For all computations, the intrinsic opening angles for the average 49-epoch image are assumed to be inside 0.620 mas and beyond 0.620 mas. The intrinsic opening angles for the 2013 and 2016 images are assumed to be: (2013) inside 0.5 mas and beyond 0.5 mas and (2016) inside 0.4 mas and beyond 0.4 mas.

For an assumed purely poloidal flow, the intensity ratios associated with the 49-epoch intensity image require a higher speed along the north counter-jet and south jet edges and lower speed along the south counter-jet and north jet edges. These relative speeds are not consistent with the faster motions along the north side of the jet and slower motions along the south side of the jet found from component motions by MLWH. A resolution to this inconsistency for poloidal flow will be suggested in the discussion of helical flows below. The trend of intensity ratios with core distance provides evidence for speed changes, and, in particular, acceleration in resolved regions near the core. On average the flow implied by the intensity ratios increases from to between and  mas ( to  ), then slowly increases to at 1.2 mas ( ) and finally accelerates to at  mas ( ). At larger distance the counter-jet intensity approaches the level of the noise and imaging artifacts so the intensity ratios cannot be determined from our data. The counter-jet/jet intensity ratio implies that the flow speed increases most rapidly in the regions of about 0.4 to 0.7 mas and about 1.2 to 1.4 mas where the two rapid expansions of counter-jet and jet occur (see Figure 24). Recall that the counter-jet width behaves similarly to the jet width within the mas of the core where it can be seen.

The average purely poloidal flow behavior inferred from the intensity ratios for the 2013 and 2016 super-resolution images are very different from each other. In part the differences may result from small separation from the core along with noise and imaging artifacts particularly in the 2016 image. Nevertheless, the 2016 image intensity ratios imply an intrinsic flow that on average rapidly accelerates from at 0.3 mas to at  mas that exceeds the speed found for the 49-epoch average image at  mas. The better quality 2013 image intensity ratios imply a much more complex behavior with an intrinsic flow that on average first increases from at 0.3 mas to at  mas, then decreases to at 0.5 mas where the opening angle changes from rapid to a slower expansion. At larger distance the flow accelerates back to at 0.6 mas. Notice that typical speeds from the 2013, 2014 and 49-epoch images at 0.5 mas are in the range 0.4 to .

We can compare typical poloidal flow speed values inferred from the counter-jet/jet intensity ratios to the proper motions associated with bright features within the jet shown in Figures 13 and 14 and the apparent motions from the wavelet analysis reported by MLWH. The proper motions shown in Figures 13 and 14 indicate at core distances  mas with motions up to at core distances  mas (somewhat less on the north side). Figure 16 in MLWH also shows clearly that the velocities of the fastest 10% of features increase from less than at  mas to over beyond  mas. Thus, tracking individual components and the wavelet analysis indicate acceleration in the same region indicated by the counter-jet/jet intensity ratios. In particular, at a viewing angle of the wavelet analysis indicates acceleration along the jet edges from at  mas to at  mas from the core. Thus, for an assumed purely poloidal flow the intensity ratios provide qualitative but not quantitative agreement with the speeds and acceleration found by the component tracking and by the wavelet analysis.

In Section 4.2.1 we found that faster northern edge motions along the jet can be associated with helical flow. We examine the intensity ratio implications for helical flow in the lower panels in Figure 25. Here we compute the difference in flow angle, , relative to pure poloidal motion along the jet edges, and the flow speed, , required to produce the observed counter-jet/jet intensity ratios, where again we have made the simplifying assumption that the flow angle to the line-of-sight can be expressed as . In these lower panels the flow speed is computed from

(6)

where

A little manipulation allows us to solve for which is found from

(7)

and as previously, , indicates clockwise helical flow. Our helical flow analysis here depends on both north and south edge counter-jet and jet intensity ratios and thus can only be applied if both edge ratios are beyond the region of core influence as shown by the shaded region in Figure 25. These equations look very similar to Equations 2 and 3 but are in terms of intensity ratios rather than component speeds.

The assumption of helical flow solves the edge flow issues associated with the intensity ratios raised by the assumption of purely poloidal flow. Now flow will appear faster on the northern jet edge and slower on the southern jet edge as observed. The helical flow speeds and acceleration associated with the 49-epoch image intensity ratios are comparable to the values indicated by the component tracking and by the wavelet analysis inside 1 mas, i.e., acceleration from to between 0.6 and 0.8 mas, albeit with values of , considerably larger than suggested by our estimate of based on average edge motions inside 1 mas. Inside 1 mas the 49-epoch image intensity ratios probe the counter-jet and jet region where the first rapid expansion ceases and the jet recollimates.

Speeds and acceleration from component tracking and the wavelet analysis exceed the values suggested by the 49-epoch image counter-jet/jet intensity ratios at core distances beyond 1 mas, although by 1.4 mas has decreased to a value of about that is comparable to our estimate based on average edge motion between 1 and 4 mas. A possible explanation for slower motion based on intensity ratios beyond 1 mas is the slow components detected in the wavelet analysis. Those components would not be subject to strong beaming so adding the emission from them and from the fast components would mimic an intermediate speed, which is what is seen. Recall that the intensity ratio is set by the material flow speed, in contrast to component speeds which can be set by patterns moving at different speeds from the material. While it is possible for a pattern to move faster than the material speed, e.g., flow plus sound or Alfvén speed, it is unlikely that pattern effects can explain why the counter-jet is brighter than would be expected for the observed component speeds.

At very small angular scales the super-resolution image intensity ratios probe the first rapid counter-jet and jet expansion region. The better quality 2013 intensity ratios would indicate a helical flow that accelerates from to between 0.4 and 0.6 mas in the initial rapid counter-jet and jet expansion region that extends to  mas (see Figure 24). This provides our only evidence for rapid acceleration associated with rapid counter-jet and jet expansion and towards flow values consistent with the faster component motions found by the WISE analysis between 1 and 4 mas. The poorer quality 2016 intensity ratios could be taken to indicate the creation of a slower moving counter-clockwise helical flow, i.e. , beyond 0.55 mas, consistent with the slower moving component motions found by the WISE analysis between 1 and 4 mas. In any event, it is clear that at these small angular scales the situation is complex and any quantitative analysis is subject to noise and imaging artifacts in addition to any real temporal changes in the flow structure.

4.3 Jet Structure

We now turn from jet shape and flow speeds to examine what can be learned about the internal structure of the jet from its intensity distribution. The pronounced edge brightness is a clear sign that the jet is not uniform and that the main emission comes from the edge. This raises questions such as how thick is the emitting region and what might be in the center. Also Doppler boosting implied by the relativistic motions and helical jet flow have implications for the intensity distribution of the jet that are compared with the imaging results. Interpreting the brightness distribution has to begin with an understanding of what part of the jet is actually being observed along a given line of sight in a jet oriented close to the line-of-sight.

Lines of Sight

The line-of-sight through the jet axis passes through the jet nearside, center and farside at significantly different core distances. This complicates understanding the implications of the transverse intensity profiles for sheath and spine emission structure. It also complicates the use of mid-jet as opposed to jet edge proper motions for understanding the jet flow and Doppler beaming effects. For example, at a viewing angle of 17°, a sightline that passes through the jet axis at 1 mas from the core, passes through the farside and nearside of the jet at  mas from the core, respectively. A typical line-of-sight through the jet axis crosses the nearside of the jet from 2 times (close to the core) to 1.5 times farther from the core than it crosses the farside of the jet. As a consequence the average farside edge intensity is 3 (close to the core) to 2 times the nearside edge intensity where the sightline leaves the jet (see Figure 6 to compare northern or southern ridgeline intensities at different distances from the core). If emissivity is relatively uniform, the intensity where a sightline crosses the jet axis might be the average of the nearside and farside “edge” intensities along the sightline. This intensity difference along a sightline through the jet axis is problematic for determination of the nearside, center and farside location of observed centerline structures. We will address this issue in what follows.

Sheath & Spine

Transverse slice profiles (see Figures 4 and 5) indicate that the intensity at the jet centerline can be less than half that of the bright edges. Production of such edge-brightened profiles requires a high sheath and low spine effective emissivity. If we assume a cylindrical jet and that emission comes only from a sheath of thickness , where and are the radius of the jet and spine respectively, then the sheath path length through the jet axis for viewing angle can be written as . The path length through the sheath along a line-of-sight through the jet edge at (i.e., the midpoint of the sheath), can be written as

Thus, a cylindrical jet with uniform sheath emissivity and a jet edge intensity twice the center intensity, as suggested by the transverse slice profiles, means

The actual jet is expanding, so the intrinsic emissivity and sheath thickness will change with the jet’s radius and core distance. Also sightlines cross the jet edges and jet axis at very different core distances. However, given that we found the average edge intensity at the distance where a line-of-sight crosses the jet centerline to be on the order of the nearside and farside average, we can still conclude that the jet consists of a low apparent emissivity jet spine and high apparent emissivity jet sheath with observed emission largely coming from a sheath with thickness on the order of . Note that these arguments apply once the jet ridges are resolved from each other and the jet expansion is moderate — say at core distances beyond about  mas.

For intensity profiles to be as edge-brightened as they appear despite the path length along sightlines through the jet spine being on the order of three times the path length through the nearside and farside jet sheath, the “effective” spine emissivity must be at least an order of magnitude less than the “effective” sheath emissivity. Intrinsic spine emission can be Doppler deboosted relative to intrinsic sheath emission by a factor

To achieve an order of magnitude deboosting, , for with the intrinsic sheath speed and Lorentz factor of and appropriate between 1 and 4 mas, would require a spine Lorentz factor and speed of and . This required minimum gives an apparent motion at of . This lower limit is consistent with the seen at HST-1 in the optical by Biretta, Sparks & Macchetto (1999) and the of the fastest component seen in the radio by Giroletti et al. (2012). Thus, it is possible that a high speed Doppler deboosted spine lies within the slower moving sheath.

Rotation, Expansion & Doppler Boosting

The presence of clockwise helical flow in the sheath would lead to differential Doppler boosting across the jet. If we consider the values found above for the region between 0.5 and 1 mas (1 and 4 mas region in parentheses) of (2.9°) and (0.924) with (17.3°), then we would predict an enhanced Doppler boost of the southern jet edge relative to the northern jet edge in the range

where the numerical value range corresponds to the results for 0.5 to 1 mas and 1 to 4 mas in that order.

This predicted differential Doppler boosting is lower closer to the core. A comparison of ridgeline intensities from Figure 6 shows a southern ridgeline intensity equal to or greater than the northern ridgeline intensity at core distances  mas. However, the northern ridgeline intensity is higher at core distances  mas in the 23-epoch image. While the ratio is significantly less close to the core, it would still predict a brighter southern edge. It seems more likely that the maximum edge intensity is a function of change in jet direction associated with long-term epoch evolution (see Figures 16 to 18), and jet edge intensity differences are not likely the result of differential Doppler boosting resulting from rotation.

There will also be differential Doppler boosting between the nearside (ns) and farside (fs) of the jet along lines-of-sight through the jet axis because, thanks to the jet opening angle, the sides are at different angles to the line-of-sight. If we consider with opening angle between 0.5 and 1 mas, and between 1 and 4 mas with then we would predict an enhanced Doppler boost of the nearside relative to the farside in the range

where the smaller value is closer to the core, and and . This amount of differential Doppler boost between nearside and farside can partially compensate for higher farside intrinsic intensities along lines-of-sight through the jet axis.

Clockwise helical flow around the jet will lead to apparent transverse motion across the jet centerline with northwards apparent transverse motions across the nearside of the jet and southwards apparent transverse motions across the farside of the jet. MLWH find a tendency towards northwards apparent transverse motions and this indicates a bias towards seeing motions on the nearside of the jet (see Figures 2 and 3 in MLWH).

4.4 Long-Term Pattern Motion

In Section 3.4, we showed that the transverse position of the jet shifts on a quasi-periodic about 8 to 10 yr timescale. This period provides a characterization of the variation timescale but does not have an adequate time span to claim true periodic motion such as would be expected from precession. In this section, we discuss the implications of these shifts for the properties of the interface between the jet and external medium if their propagation is described by the Kelvin-Helmholtz helical mode.

The structural variations observed are consistent with a jet base that changes direction on the timescale noted above, leading to a variation in flow direction that appears to propagate down the jet at apparent speeds where the range is set by the two different methods to determine the speed used in Section 3.4.2. At a viewing angle of , the intrinsic pattern motion range is . The intrinsic pattern motion is significantly greater than the slower intrinsic component motions from MLWH of but is significantly less than the faster intrinsic component motions from MLWH of at distances beyond 1 mas. MLWH suggested that the fast and slow component motions associated with the bright jet edges indicated velocity shear associated with the jet-external medium interface. Thus, we suggest that the pattern can be thought of as moving through an external (e) medium which itself can be moving with flow speed but that the pattern moves more slowly than the internal jet (j) flow speed which is .

In what follows let us assume that the jet with flow speed (), determined by the MLWH faster moving jet components beyond 1 mas, is separated by a sharp velocity discontinuity from an external medium with flow speed (), determined by the MLWH slower moving components beyond 1 mas. We assume a direction change at the jet base on a  yr timescale and that this direction change propagates with an intrinsic pattern speed of .

A non-ballistically moving bulk displacement of the jet can be described by the properties of the Kelvin-Helmholtz helical mode (see Hardee, 2007) provided the displacement frequency is less than the resonant frequency (see Equation C1). See Appendix C for a brief review of the important equations. Numerical study shows that the helical mode above resonance does not lead to bulk displacement of the jet as a whole (Hardee et al., 2001). Below resonance the helical mode propagation properties are relatively insensitive to the presence of a velocity shear layer (see Birkinshaw, 1991). The important propagation properties are revealed by analytical approximations for helical wave propagation along a cylindrical jet with purely poloidal flow and magnetic field inside and outside a sharp velocity shear edge (see Hardee, 2007)). A simple linear displacement at the jet base can be described by two counter-rotating equal amplitude helical displacements. Provided flow and magnetic helicity are relatively small beyond 1 mas, approximations to pattern motions assuming no helicity should prove adequate for estimating the jet to external medium interface properties required to produce the observed pattern motion.

The  yr timescale for the change in direction at the jet base implies an angular frequency of  rad s. An estimate for the conditions required for the lower limit to the resonant frequency to exceed this angular frequency can be obtained by assuming that, in the “sonic” limit, the sound speeds in the jet and external medium are the same or that, in the “Alfvénic” limit, the Alfvèn speeds in the jet and external medium are the same. Here the “sonic” limit means that the sound speed greatly exceeds the Alfvén speed (weakly magnetized flow scenario) and the “Alfvénic” limit means that the Alfvén speed greatly exceeds the sound speed (strongly magnetized flow scenario). In these scenarios the resonant frequency (see Equation C2)

at say 3 mas where 0.9 mas  cm provided . We have used the FWHM value for the jet width at 3 mas from Figure 24 to estimate , and is the sound speed in the sonic limit or the Alfvén speed in the Alfvénic limit. This required lower limit to the sound speed or Alfvén speed in the jet and external medium at about 7% of lightspeed lies well below the maximum possible sound speed or Alfvén speed of or , respectively. Thus, it is clearly possible for a quasi-periodic  yr timescale jet displacement to be described by the low to resonant frequency helical mode propagation properties.

Pattern motion, , is a function of the jet speed , the external medium speed , and the enthalpy ratio, between the jet and external medium. Here the enthalpy where , and are the rest mass density, pressure and adiabatic index, respectively.

If pattern motion is described by helical mode propagation below the resonant frequency, the required enthalpy ratio is given by Equation C5. If magnetic fields are small, then in the sonic limit where and for (), and (). If the jet is similar in density to the surrounding medium unless the faster moving jet components are not truly representative of the underlying jet flow. A much faster jet flow of () such as required to Doppler deboost the spine and produce the apparent superluminal motions in the optical at HST-1 implies a jet over an order of magnitude less dense than the surrounding external medium. If magnetic fields are large, then in the Alfvénic limit where and . The Alfvénic limit suggests somewhat higher jet density relative to the unmagnetized case as in general .

If pattern motion is more accurately described by helical mode propagation near the resonant frequency, the required enthalphy ratio is given by Equations C6. The principal difference with the low frequency form lies in the square in the velocity term so that where