THE INNERMOST COLLIMATION STRUCTURE OF THE M87 JET DOWN TO ∼10 SCHWARZSCHILD RADII

In the present study, we aim to measure the width of the innermost jet. To this end, we fit a single elliptical Gaussian model to each image with the AIPS task JMFIT and derived the deconvolved parameters of the core region. Note that the derived Gaussian size from this simple modeling could yield a larger size than that of the true core (i.e., optical depth ∼1 surface) in the jet propagation direction because of blending of the optically thin jet part. Also, in the transverse direction to the jet axis, this method would give the total size of the true core plus any surrounding emission, if the core region has sub-structure that is not resolved by VLBA baselines in this direction (Dodson et al. 2006). However, here we are interested in measuring the entire width of the innermost jet.

The results are summarized in Table 1. Most of the derived values (particularly for the minor axes) are smaller than the beam size. Nonetheless, it is known that such sizes are measurable when the emitting region is sufficiently bright and robust self-calibration is possible using as many as 10 VLBA stations (eight at 86 GHz), such that the fringe visibilities can be accurately calibrated. For the M87 core, the derived sizes along the minor axes at 5, 8.4, 15, 23.8, 43, and 86 GHz correspond to amplitude decreases of 15%, 18%, 25%, 30%, 33%, and 30% at ∼60% of the longest baseline (80Mλ, 140Mλ, 240Mλ, 380Mλ, 700Mλ, and 1100Mλ, respectively, which yield effective angular resolutions at P.A. = 20°). These decreases are sufficiently larger than the typical VLBA amplitude calibration accuracy of ∼5% (a similar discussion is presented in Kellermann et al. 1998). At 86 GHz, previous Global Millimeter VLBI Array observations report similar sizes to the values listed in Table 1 (≲ 50 μas or 99 ± 21 μas; Krichbaum et al. 2006; Lee et al. 2008, respectively).

Table 1. Multi-frequency VLBA Observations of M87

Freq.	Date	Beam θbm	Irms	JMFIT Gaussian
Band		Size, P.A.			θmaj	θmin	P.A.	$\theta _{20^{\circ }} (\equiv W_{\rm c})$	$\frac{\theta _{20^{\circ }}}{\theta _{\rm bm, 20^{\circ }}}$
(GHz)		(mas)2, (deg.)	$\left(\frac{\rm mJy}{\rm beam}\right)$		(mas)	(mas)	(deg.)	(mas)
		(a)	(b)		(c)	(d)	(e)	(f)	(g)
2.3	2010 Apr 8	6.01 × 2.92, −3	0.79		3.26 ± 0.08	<1.06	286 ± 3	<1.06	<0.20
	2010 Apr 18	5.97 × 3.00, −4	0.66		3.47 ± 0.08	<1.06	287 ± 3	<1.06	<0.20
				Ave.	3.37	<1.06	287	<1.06	<0.20
5.0	2010 Apr 8	2.63 × 1.36, 1	0.43		1.36 ± 0.06	0.57 ± 0.16	285 ± 4	0.58 ± 0.17	0.25
	2010 Apr 18	2.70 × 1.31, 0	0.37		1.38 ± 0.06	0.56 ± 0.13	288 ± 3	0.56 ± 0.14	0.24
				Ave.	1.37	0.57	287	0.57	0.25
8.4	2009 May 23	1.48 × 0.55, −9	0.21		0.61 ± 0.07	0.38 ± 0.05	292 ± 5	0.38 ± 0.06	0.40
	2010 Apr 8	1.67 × 0.84, 2	0.37		0.84 ± 0.04	0.42 ± 0.06	295 ± 4	0.42 ± 0.07	0.29
	2010 Apr 18	1.63 × 0.85, −3	0.34		0.87 ± 0.05	0.40 ± 0.08	295 ± 4	0.40 ± 0.09	0.29
				Ave.	0.77	0.40	294	0.40	0.33
15	2009 Jan 7	1.04 × 0.47, −3	0.82		0.45 ± 0.02	0.28 ± 0.01	298 ± 5	0.28 ± 0.02	0.34
	2009 Jul 5	1.00 × 0.46, −8	0.63		0.39 ± 0.05	0.26 ± 0.03	305 ± 7	0.26 ± 0.06	0.35
	2010 Feb 11	1.00 × 0.43, −2	1.07		0.41 ± 0.02	0.21 ± 0.04	286 ± 2	0.21 ± 0.04	0.27
	2010 Apr 8	0.94 × 0.48, −5	0.59		0.47 ± 0.02	0.25 ± 0.04	300 ± 5	0.25 ± 0.04	0.33
	2010 Apr 18	0.92 × 0.49, −10	0.44		0.49 ± 0.02	0.24 ± 0.02	298 ± 3	0.24 ± 0.03	0.33
	2010 Sep 29	1.00 × 0.45, −2	0.72		0.42 ± 0.07	0.26 ± 0.03	307 ± 4	0.27 ± 0.08	0.35
	2011 May 21	0.93 × 0.43, −5	0.74		0.51 ± 0.05	0.28 ± 0.05	302 ± 5	0.29 ± 0.07	0.38
				Ave.	0.45	0.25	299	0.26	0.34
23.8	2010 Jan 18	0.60 × 0.29, −8	0.76		0.25 ± 0.04	0.20 ± 0.02	300 ± 3	0.20 ± 0.04	0.42
	2010 Apr 4	0.63 × 0.29, −4	0.68		0.27 ± 0.02	0.18 ± 0.03	313 ± 7	0.19 ± 0.04	0.38
	2010 Apr 8	0.54 × 0.29, −5	1.36		0.24 ± 0.03	0.18 ± 0.04	326 ± 6	0.20 ± 0.05	0.43
	2010 Apr 18	0.57 × 0.28, −13	1.12		0.25 ± 0.02	0.20 ± 0.03	327 ± 3	0.21 ± 0.04	0.54
	2010 May 1	0.62 × 0.28, −9	0.73		0.29 ± 0.02	0.19 ± 0.03	315 ± 6	0.20 ± 0.04	0.42
	2010 May 15	0.64 × 0.29, −12	0.80		0.31 ± 0.03	0.20 ± 0.04	317 ± 4	0.21 ± 0.05	0.48
	2010 May 30	0.62 × 0.28, −11	0.83		0.28 ± 0.02	0.16 ± 0.05	314 ± 5	0.17 ± 0.05	0.38
				Ave.	0.27	0.19	316	0.20	0.44
43	2009 Mar 13	0.29 × 0.13, −6	0.55		0.15 ± 0.01	0.13 ± 0.02	313 ± 9	0.13 ± 0.02	0.61
	2010 Jan 18	0.29 × 0.13, 2	1.01		0.11 ± 0.02	0.10 ± 0.02	262 ± 3	0.10 ± 0.03	0.39
	2010 Apr 8	0.26 × 0.13, −3	0.91		0.13 ± 0.02	0.11 ± 0.02	354 ± 4	0.13 ± 0.03	0.57
	2010 Apr 18	0.26 × 0.13, −8	1.12		0.13 ± 0.01	0.11 ± 0.01	358 ± 3	0.13 ± 0.01	0.63
	2010 May 1	0.29 × 0.14, −6	0.89		0.11 ± 0.02	0.11 ± 0.01	318 ± 19	0.11 ± 0.02	0.49
	2010 May 15	0.29 × 0.14, −9	0.95		0.12 ± 0.01	0.10 ± 0.02	254 ± 6	0.11 ± 0.02	0.48
	2010 May 30	0.30 × 0.15, −6	1.07		0.13 ± 0.01	0.09 ± 0.01	257 ± 8	0.10 ± 0.02	0.42
				Ave.	0.13	0.11	280	0.11	0.51
86	2007 Feb 18	0.25 × 0.08, −18	6.47		0.079 ± 0.021	0.065 ± 0.023	52 ± 23	0.074 ± 0.032	0.64

Notes. Columns: (a) uniformly weighted beam size; (b) image noise level; (c), (d), and (e) FWHM sizes of major/minor axes and position angle of the derived model; (f) projected FWHM size of the model at P.A. = 20°; (g) ratio of the Gaussian size divided by the beam size at P.A. = 20°.

Download table as:  ASCIITypeset image

We estimated the standard errors of the derived Gaussian parameters for each datum as follows. Formally, statistical size uncertainties purely based on signal-to-noise ratios (size of the fitted model divided by its peak-to-noise ratio; e.g., Fomalont 1999) result in very small values for M87 (a level of a few μas or smaller) because the core is bright at each frequency (peak-to-noise >70 ∼ 1000). However, in practice, the model parameters are more strongly affected by imperfect CLEAN/deconvolution processes under limited uv coverage. Then, we divided each individual data set into three subsets with equal integration times and repeated the deconvolution and JMFIT processes for each subset individually. Through this procedure, we obtained three sets of quasi-independent fitting results for each epoch; the rms scatter can be calculated for each parameter (i.e., major/minor axes and P.A.). These values in the scatter are adopted as realistic errors for each model parameter in Table 1.

As an additional check, we conducted a fake visibility analysis to examine how precisely sizes smaller than beam sizes can be recovered (a similar analysis is described in Taylor et al. 2004); using the AIPS task UVCON, we created fake visibilities that are equivalent to the derived Gaussian parameters at each frequency with similar uv coverage as the actual observations. We produced 10 independent visibility data sets at each frequency by adding random noise at a level seen in the actual observations, repeated the CLEAN/deconvolution and JMFIT procedures on each data set, and calculated the rms scatter of the recovered model parameters at each frequency. We confirmed that these scatter values were smaller than one-third of the quoted errors in Table 1 at frequencies smaller than 43 GHz and less than one-half of the quoted errors at 86 GHz.

We note that JMFIT derived quite a small size (less than 15% of the beam size) for the minor axis of the 2 GHz core, which indicates a marginal constraint. Therefore, in Table 1 we instead set one-fifth of the beam size as an upper limit, which corresponds to a ∼5% amplitude decrease at the longest baseline in this direction.

In Figure 2, we show Gaussian size projected along P.A. = 20° as a function of frequency (hereafter, denoted W_c).⁷ This projection is perpendicular to the overall jet axis, corresponding to the direction of the jet width. We found W_c to be clearly frequency dependent, becoming smaller as frequency increases. To determine the averaged frequency dependence of W_c, we fit a power-law function to this plot using the data at 5, 8.4, 15, 23.8, 43, and 86 GHz (the 2 GHz data were excluded because of their upper limits). We found the best-fit solution to be W_c(ν)∝ν^ξ where ξ = −0.71 ± 0.05. Interestingly, when this frequency dependence is extended toward higher frequency, its extrapolated size appears to result in a similar value to the size measured by a recent 230 GHz VLBI experiment (Doeleman et al. 2012), which used a circular Gaussian fit. For reference, we also show a Gaussian size distribution projected along the jet propagation direction, but this distribution does not seem to fit the 230 GHz core size.

Zoom In Zoom Out Reset image size

In Figure 3, we show the radial profile of the M87 jet width as a function of the (de-projected) distance along jet (assuming a jet inclination angle i = 15°; Biretta et al. 1999; Perlman et al. 2011). Here, we investigated the jet width profile using the following two procedures: (1) measurements of W_j(r); for the region where the jet is clearly resolved into a two-humped shape at each frequency, we made transverse slices of the jet every 0.01 ∼ 0.5 mas along the jet and each of the two-humped structures was fit with a double-Gaussian function. We then defined the separation between the outer sides of the half-maximum points of the two Gaussians as the jet width. When the jet at a particular frequency becomes single-humped toward the upstream region, we measured the jet width using higher frequency images because the jet is clearly resolved again into a two-humped shape at higher resolution. By using the images between 2 and 43 GHz, such measurements were repeated over a distance of ∼10⁵ R_s down to ∼100 R_s along the jet. This process is basically the same as that used in AN12, but we exclude measurements for the single-humped region. We aligned the radial positions of the W_j(r) profiles between different frequencies by adding the amount of core shift measured in H11 (described below). This amount and the associated positional error at each frequency contribute only tiny fractions of the distance where W_j(r) was measured at each frequency (a level of 10^{−1 ∼ −2} at 43 GHz to 10⁻³ at 2 GHz), so the horizontal error bars for W_j(r) are not shown in Figure 3. At 86 GHz, we could not perform reliable measurements of the jet width because the edge-brightened jet was only marginally imaged at a level of 2σ–3σ. (2) measurements of W_c(r); closer to the central engine, we furthermore constructed a radial distribution of W_c. Because our previous astrometric study (H11) measured locations of the cores as 41, 69, 109, 188, 295, and 674 μas at 43, 23.8, 15, 8.4, 5, and 2 GHz from the convergence point of the core shift (in R.A. direction or P.A = 270°), we can set their de-projected distances along the jet (P.A. = 290°) as 24, 40, 63, 108, 170, and 388 R_s for i = 15°, respectively. For the 86 and 230 GHz cores, we can also set their de-projected positions as 12 and 5 R_s by assuming the same asymptotic relation (r_c∝ν^−0.94; H11) upstream of the 43 GHz core. Here, we assume that the central engine is located at the convergence point of the core shift specified in H11.

Zoom In Zoom Out Reset image size

We confirmed that the edge-brightened region can be well expressed as a parabolic structure: W_j(r)∝r^{0.56 ± 0.03} (solid line in Figure 3). This result is in good agreement with the finding of r^{0.58 ± 0.02} in AN12. For the region around ∼100 R_s, where the independent measurements of W_j and W_c overlap with each other, W_c at 5 and 8.4 GHz smoothly connect with W_j at 43 GHz. The combination of the core size (W_c∝ν^ξ, where ξ = −0.71 ± 0.05) and core position (r∝ν^α, where α = − 0.94 ± 0.09; H11) yields the following radial dependence of W_c: $W_{\rm c}(r) \propto r^{\frac{\xi }{\alpha }} = r^{0.76\pm 0.13}$ (dotted line in Figure 3), which is slightly steeper than that of the outer jet W_j(r), although the uncertainty is still large. In the present work it is still difficult to distinguish whether the values of W_c at 5, 8.4, 15, and 23.8 GHz are on the solid line or the dashed line due to their positional and size uncertainties. On the other hand, W_c at 43 and 86 GHz tends to be below the solid line. At 230 GHz, the exact profile cannot be discriminated because the data are totally dominated by positional uncertainties.

We note that the two methods used for W_j(r) (a double-Gaussian fit to each slice image) and W_c(r) (based on a single two-dimensional elliptical Gaussian model for each two-dimensional image) are different from each other. Nevertheless, the observed consistencies of the values between the two methods were confirmed in the overlapping region, indicating that W_c is actually a good tracer for the width of the innermost jet region.

Theoretically, the shape of a magnetized jet is determined by the detailed force balance across the poloidal field lines. This balance is described by the trans-field equation, which was originally derived by Okamoto (1975) for a steady, axisymmetric, rotating magnetized flow in a gravitational field without external pressure. Later, magnetic hoop stresses associated with toroidal field lines were invoked and played a major role in the realization of global collimation of magnetized jets (Blandford & Payne 1982; Heyvaerts & Norman 1989; Sakurai 1985; Chiueh et al. 1991).

Here, we observationally confirmed that the M87 jet is well characterized with a parabolic collimation profile between ∼100 and ∼10⁵ R_s. This result is consistent with the prior work of AN12, where these authors also found a transition into a conical shape above ∼10⁶ R_s. Regarding the formation of a parabolic shape, recent theoretical studies indicated the importance of external pressure support at the jet boundary. Okamoto (1999) analytically showed that hoop stresses alone would not produce a global collimation and numerical studies also suggest a similar conclusion (e.g., Nakamura et al. 2006; Komissarov et al. 2007, 2009; Toma & Takahara 2013). Komissarov et al. (2009) show that when the external gas pressure follows p_ext∝r^−a where a ≲ 2, the jet maintains a parabolic shape as W_j∝r^a/4. On the other hand, for a > 2, the jet eventually becomes conical due to insufficient external support. If the observed radio emission of M87 traces the exterior edge of a magnetized jet, the measured width profile suggests a ∼ 2. As a source of such a confinement medium, AN12 propose an interstellar medium bounded by the gravitational influence of the central black hole, such as a giant advection-dominated accretion flow (Narayan & Fabian 2011). We add that a purely hydrodynamic jet can also produce a gradual parabolic collimation of M87 (e.g., Bromberg & Levinson 2007).

Interestingly, at the end of the collimation region, both AN12 and Bromberg & Levinson (2007) suggest HST-1 (a peculiar knot at a deprojected distance ≳120 pc or 2 × 10⁵ R_s) to be a stationary shocked feature resulting from the overcollimation of the inner jet, as originally proposed by Stawarz et al. (2006) and Cheung et al. (2007) to explain the broadband flaring activity. While our recent kinematic observations of HST-1 show clear evidence of superluminal motions of the overall structure, a weak, slowly moving feature was also found upstream of HST-1 (Giroletti et al. 2012), which indeed could be related to a recollimation shock.

For the first time, we have revealed a detailed collimation profile down to r ∼ 10 R_s by investigating the multi-frequency properties of the radio core. An intriguing point here is that the measured collimation profile suggests a possible tendency of a wider jet opening angle around ∼10 to ∼100 R_s from the central engine.

Since the two methods used to make our width measurements switch around r ∼ 100 R_s, one could speculate that the profile change is related to some systematic effects due to the different methods. To check the two-humped jet shape even closer to the core (r ≲ 100 R_s) more clearly, we created a 43 GHz image convolved with a slightly higher resolution (a circular beam 0.14 mas in diameter). The image is shown in Figure 4. More than ∼0.4 mas downstream of the core (de-projected distance ∼220 R_s for i = 15°), where the jet is parabolic on a logarithmic distance scale, the two ridges are already oriented in a similar direction (dark blue region in Figure 4). On the other hand, the opening angle made by the two ridges appears to broaden more rapidly within ∼0.3 mas of the core (the region with contours in Figure 4), resulting in a more radially oriented structure near the base. Such a tendency is consistent with the observed possible transition from the solid line to the (steeper) dashed line around r ∼ 100 R_s in Figure 3.

Zoom In Zoom Out Reset image size

A transition of the jet collimation profile near a black hole is actually suggested in some theoretical aspects. In the framework of relativistic MHD jet models, most of the energy conversion from magnetic to kinetic occurs after the flow passes through the fast-magnetosonic point (the "magnetic nozzle" effect; Li et al. 1992; Vlahakis & Königl 2003). Beyond this point, the magnetized jet starts to be collimated asymptotically into a parabolic shape because the increasing plasma inertia winds the field lines into toroidal directions and thus amplifies the hoop stresses (e.g., Tomimatsu & Takahashi 2003). The radius of the fast point is typically a few times the light cylinder radius R_lc (Li et al. 1992), where R_lc is of the order of (1 ∼ 5) R_s (Komissarov et al. 2007). Thus, if the M87 jet is magnetically launched, the observed possible transition of the jet shape around 10 ∼ 100 R_s could be explained by this process. Moreover, the jet on this scale is likely to have complicated interactions with its surrounding medium, including accretion flows, coronae, and disk winds (e.g., McKinney 2006). The geometries of these processes and the local pressure balance at the jet boundary would affect the initial jet shape. Alternatively, a change in the jet shape could occur as an apparent effect due to projection, if the jet inclination angle is not constant down to the black hole (McKinney et al. 2013).

It is interesting to note that the time-averaged dependence ν^−0.71 of W_c appears to connect with the 230 GHz core size. However, this apparent connection should be treated with caution; the uv coverage used in Doeleman et al. (2012) yields the highest (∼3000Mλ) angular resolution along the jet direction for M87, with ∼5 times shorter projected baselines transverse to the jet. Thus, the derived size of ∼40 μas using their circular Gaussian fit could be more weighted by the structure along the jet, unless the brightness pattern of the jet base (when projected toward us) is actually close to circular. To clarify the exact relationship of the core size at the higher frequency side, the addition of north–south baselines in the future Event Horizon Telescope array is crucial, which can be realized by including Chilean stations such as ASTE, APEX, and ALMA (Broderick et al. 2009; Doeleman et al. 2009).

The results presented here have shed new light on crucial issues that should be addressed more rigorously in future observations: where is the exact location of the profile change between ∼10 and ∼100 R_s and how does it change (e.g., a sharp break or a gradual transition)? In addition, simultaneous observations at multiple frequencies are important; the dynamical timescale of our target region is of the order t_dyn ∼ 10 R_s/c ∼ 7 days, so the jet structure (size and position of the core) can be variable on this timescale (e.g., Acciari et al. 2009). To address these issues, we are currently conducting new high-sensitivity VLBA observations at 86 GHz in combination with quasi-simultaneous sessions at lower frequencies, which allows more robust investigations of the jet structure within 100 R_s and thus tests some specific models more quantitatively. Finally, we also stress the significance of future imaging opportunities with RadioAstron at 22 GHz or global-VLBI at 43 and 86 GHz including ALMA baselines, because these facilities will provide images with drastically improved resolutions, especially in the direction transverse to the M87 jet.