Faraday Rotation in the Jet of M87 inside the Bondi Radius: Indication of Winds from Hot Accretion Flows Confining the Relativistic Jet

Active galactic nuclei (AGNs) are powered by accretion of gas onto supermassive black holes at the centers of galaxies. It is now widely believed that there are two distinct modes of black hole accretion: cold and hot. A cold accretion flow forms an optically thick but geometrically thin disk, radiating thermal blackbody emission with the gas temperature in the range of 10⁴–10⁷ K (Shakura & Sunyaev 1973, see, e.g., Netzer 2013 for a review). On the other hand, hot accretion flows are thought to be optically thin but geometrically thick with a large portion of the gravitational binding energy of the accreted gas advected into the black hole (e.g., Ichimaru 1977; Narayan & Yi 1994, see, e.g., Yuan & Narayan 2014 for a review). The most critical factor in determining the accretion mode is the mass accretion rate ($\dot{M}$) relative to the Eddington rate (${\dot{M}}_{\mathrm{Edd}}$), or equivalently, the disk luminosity (${L}_{\mathrm{disk}}$) relative to the Eddington luminosity (L_Edd). Observationally, L_disk/L_Edd ≈ 1% is usually assumed to be a dividing line between the two accretion modes (e.g., Ghisellini et al. 2011; Heckman & Best 2014).

Most (≈98%) nearby AGNs spend their lives in a low accretion state, making them low-luminosity AGNs (LLAGNs, Ho 2008; Netzer 2013), which are thought to be powered by hot accretion flows. One of the representative models of hot accretion flows is the advection-dominated accretion flows (ADAFs, Ichimaru 1977; Narayan & Yi 1994, 1995a, 1995b), which is characterized by self-similar solutions with a density profile of ρ ∝ r^−1.5 and a constant mass accretion rate as a function of spherical radius (r). Two important properties found in ADAFs are that (i) the flows are convectively unstable and (ii) the Bernoulli parameter of the flow is positive, indicating that strong outflows are a natural outcome of hot accretion flows (e.g., Narayan & Yi 1994, 1995a). These properties led to two variants of ADAF, convection-dominated accretion flow (CDAF, e.g., Narayan et al. 2000; Quataert & Gruzinov 2000; Igumenshchev & Narayan 2002), and adiabatic inflow–outflow solution (ADIOS, e.g., Blandford & Begelman 1999, 2004; Begelman 2012), respectively.

A number of numerical simulations have been performed to better understand the dynamics of hot accretion flows (e.g., Stone et al. 1999; Igumenshchev & Abramowicz 2000; Machida et al. 2001; Igumenshchev et al. 2003; Pen et al. 2003, see Yuan et al. 2012b for a review). One of the most important findings consistently seen in those simulations is that the mass accretion rate decreases with decreasing radius, namely, ${\dot{M}}_{\mathrm{in}}(r)\propto {r}^{s}$ with s > 0, or equivalently, the density profile flatter than the one of ADAF self-similar solutions, i.e., ρ ∝ r^−q with q < 1.5. The CDAF model explains the inward decrease of ${\dot{M}}_{\mathrm{in}}$ with large fluxes of both inflowing and outflowing gas in turbulent convective eddies and predicts s = 1 and q = 0.5 (e.g., Narayan et al. 2000; Quataert & Gruzinov 2000; Igumenshchev & Narayan 2002). In the ADIOS model, the inward decrease of ${\dot{M}}_{\mathrm{in}}$ is due to a genuine mass loss via gas outflows; the model predicts 0 < s < 1 and $0.5\lt q=1.5-s\lt 1.5$ (Blandford & Begelman 1999, 2004). Values of s = 0.4–0.8 and q = 0.5–1 were preferentially found in simulations (see Yuan et al. 2012b and references therein), which is in general consistent with the ADIOS model. Indeed, both three-dimensional (3D) general relativistic magneto-hydrodynamic (GRMHD) simulations of hot accretion flows (Narayan et al. 2012) and two-dimensional (2D) simulations of hot accretion flows, including magnetic fields (Yuan et al. 2012a) showed that hot accretion flows are convectively stable, supporting that hot accretion flows can lose substantial mass via gas outflows (but see Bu et al. 2016a, 2016b for gas outflows on large spatial scales when the gravitational potential of the nuclear star cluster is included).

Nevertheless, knowledge of the properties of outflows from hot accretion flows has been limited due to the difficulty in tracing the actual outflows by discriminating them from turbulent motions. Yuan et al. (2015) used a "virtual particle trajectory" approach and overcame the difficulty in their 3D GRMHD simulations. They found that the outflows from hot accretion flows are dominant in the polar region, while inflows are filling in the equatorial regions, and the geometry of the outflows can be described as conical. Similar results were obtained in another GRMHD simulation in which the collimated and relativistic jet launched from a spinning black hole is surrounded by nonrelativistic gas outflows (Sadowski et al. 2013). We clarify the terminology of gas outflows with different physical properties: hereafter, jet refers to a highly magnetized, collimated, and relativistic gas outflow possibly launched from a spinning black hole (Blandford & Znajek 1977) or from the innermost region of an accretion disk (Blandford & Payne 1982), whereas wind refers to a moderately magnetized, un-collimated, and nonrelativistic gas outflow launched from the accretion flow.

Winds have been frequently observed in luminous AGNs for which cold accretion is thought to be operating (e.g., Crenshaw et al. 2003; Tombesi et al. 2010). However, it is challenging to confirm the presence of winds from hot accretion flows, i.e., in LLAGNs, because the winds are believed to be very hot and generally fully ionized (Yuan et al. 2018). Even though UV and X-ray absorption lines with high outflow velocities have been found in some LLAGNs (e.g., Tombesi et al. 2014), due to limited angular resolution it is unclear whether those outflows originate from the accretion flows or from outside regions (e.g., Crenshaw & Kraemer 2012). Accordingly, there have been attempts to directly determine the radial density profiles of hot accretion flows in a few nearby LLAGNs with X-ray observations. For example, Wong et al. (2011, 2014) presented a density profile of NGC 3115 broadly consistent with ρ ∝ r⁻¹ inside the Bondi radius, within which the gravitational potential energy of the central black hole is larger than the thermal energy of the gas, using Chandra X-ray observations. Russell et al. (2015) showed a similar density profile of ρ ∝ r⁻¹ for M87 inside the Bondi radius, and Russell et al. (2018) found a possible difference between the density profiles in the polar region, i.e., along the jet axis, with ρ ∝ r^−0.93, and in the equatorial region, with ρ ∝ r^−1.5 from Chandra observations. Although these results are consistent with the ADIOS model and possibly indicate the presence of winds in those LLAGNs, they were obtained near the Bondi radius; measurements of density profiles well inside the Bondi radius are needed for a firm conclusion. We note that there are some studies that favor the presence of winds in the supermassive black hole in our Galactic Center, Sagittarius A* (Sgr A*), using spectral energy distribution modeling (Yuan et al. 2003), modeling of the X-ray emission lines (Wang et al. 2013), numerical simulations reproducing the Fermi Bubbles possibly inflated by those winds (Mou et al. 2014), and modeling of the motion of the gas cloud G2 slowed down by a drag force (Gillessen et al. 2018).

Winds have important astrophysical implications. The actual rate of mass accreted onto the black hole could be substantially smaller than the accretion rate measured through X-ray observations at the Bondi radius (Bondi accretion rate, ${\dot{M}}_{\mathrm{Bondi}}$) due to the mass loss via winds. Therefore, a major factor in the faintness of LLAGNs might be the reduced mass accretion rate (Bower et al. 2003), not a very low radiative efficiency as usually assumed (Xie & Yuan 2012). Also, rotational energy of spinning black holes must be extracted efficiently to explain the observed high kinetic jet powers with small mass accretion rate (Nemmen & Tchekhovskoy 2015). Furthermore, winds have a large cross section and may regulate star formation in the host galaxies via momentum transfer (Yoon et al. 2018; Yuan et al. 2018).

Another important role played by winds is their effect on the collimation of active galactic nucleus (AGN) jets. It has been a long-standing problem how jets in AGNs can be highly collimated and accelerated to nearly the speed of light. It is widely accepted that the acceleration and collimation zone in AGN jets are co-spatial and located within about 10⁵ Schwarzschild radii (r_s, Marscher et al. 2008). MHD models predict that magnetic fields can accelerate AGN jets to relativistic speeds if the jets are systematically collimated (e.g., Vlahakis 2015). It is difficult for the jets to be confined by themselves (e.g., Eichler 1993; Begelman & Li 1994; Komissarov et al. 2007) and an external confining medium is necessary to produce the observed highly collimated jets. Previous theoretical studies suggest that winds are the primary candidates for this medium (Tsinganos & Bogovalov 2002; McKinney & Gammie 2004; Bogovalov & Tsinganos 2005; De Villiers et al. 2005; Gracia et al. 2005; Globus & Levinson 2016; Nakamura et al. 2018).

M87 serves as a unique laboratory for studying AGN jets and their formation, collimation, and acceleration thanks to its proximity with a distance of 16.7 Mpc (Mei et al. 2007) and its extremely massive black hole with a mass of M_BH = (3.5–6.6) × 10⁹ M_⊙ (Gebhardt & Thomas 2009; Gebhardt et al. 2011; Walsh et al. 2013). Accordingly, this source has been studied extensively especially on scales corresponding to the jet acceleration and collimation zone. One of the most notable results is the discovery of an edge-brightened jet structure with a systematic collimation of the jet on scales ≳100 r_s (Junor et al. 1999). The large-scale collimation profile shows a transition from a semi-parabolic jet with z ∝ R^1.7, where z and R denote the jet distance and the jet radius, respectively, to a conical jet at a transition location near the Bondi radius (Asada & Nakamura 2012). The precise constraint on the location of the black hole by core-shift analysis (Hada et al. 2011) together with the source size measured with the Event Horizon Telescope (EHT) at 1.3 mm (Doeleman et al. 2012) allowed the constraint of the innermost collimation profile. The profile is consistent with a parabolic geometry (Nakamura & Asada 2013) but shows indication of a slight deviation from the larger scale profile (Hada et al. 2013, see also Hada et al. 2016; Mertens et al. 2016; Kim et al. 2018a; Walker et al. 2018).

There has been growing evidence for gradual acceleration of the jet inside the Bondi radius as well, though the scale on which bulk jet acceleration occurs is a matter of debate. Observations of HST-1, a peculiar feature that consists of a quasi-stationary component from which superluminal components are emerging and is the location of the multiwavelength flare observed around 2005 (Cheung et al. 2007), show superluminal motions with velocities larger than 6c (with c being the speed of light) at optical wavelengths (Biretta et al. 1999), and with velocities of $\approx 4c$ at radio wavelengths (Cheung et al. 2007; Giroletti et al. 2012). Asada et al. (2014) found a systematic acceleration of the jet at a distance of ≈10⁵ r_s, supported by the slow velocities obtained on smaller scales with the Very Long Baseline Array (VLBA) at 15 GHz (Kovalev et al. 2007). However, as already noted in Kovalev et al. (2007), the observed one-sideness of the jet at a distance of only ≈3 mas from the radio core is difficult to explain with sub-luminal motions at the same distance. Other studies suggest that the jet acceleration occurs on a much smaller scale (Mertens et al. 2016; Hada et al. 2017; Walker et al. 2018) and constraining the acceleration profile at various jet distances is still ongoing (J. Park et al. 2018, in preparation).

The observation of jet collimation and acceleration on the same spatial scales is consistent with the scenario that the jet is collimated by an external medium with a relatively shallow pressure profile, which results in gradual acceleration of the jets in an MHD process (Komissarov et al. 2009; Lyubarsky 2009). However, it has not been possible to either probe the external medium with observations or to verify the general picture of jet collimation and acceleration. In this study, we investigate Faraday rotation, the rotation of the plane of linear polarization by intervening magnetic fields, in the jet of M87. When linearly polarized emission passes through a magnetized medium, Faraday rotation occurs. The amount of rotation of the electric vector position angle (EVPA), Δχ, is related to the Faraday rotation measure (RM) via Δχ = RMλ², where λ is the wavelength. RM is proportional to the integral of the product of free electron density (n_e) and line of sight component of the magnetic field (B) along the path from emitter to observer (l), meaning $\mathrm{RM}\propto \int {n}_{e}(l)B(l){dl}$ (e.g., Gardner & Whiteoak 1966). Thus, observations of the Faraday rotation of polarized jets can probe the magnetized medium between the jet and the observer, i.e., the external medium. Unfortunately, the jets in nearby LLAGNs are usually very weakly polarized (see, e.g., Bower et al. 2017 for further discussion), and the Faraday rotation observations have been limited to specific emitting regions in some sources such as Sgr A* (e.g., Bower et al. 2003, 2018; Marrone et al. 2006, 2007; Liu et al. 2016), 3C 84 (e.g., Taylor et al. 2006; Plambeck et al. 2014; Nagai et al. 2017; Kim et al. 2018b), and M87 (e.g., Zavala & Taylor 2002; Kuo et al. 2014). In this work, we obtain RM values at various locations in the M87 jet by exploiting multifrequency VLBA data from multiple epochs and present the radial RM profile of the jet between 5000 and 200,000 r_s. Then, we test the conjecture that winds are launched from hot accretion flows and serve as the external confining medium of the jet using the RM data.

We searched the VLBA archive for data suitable for a study of linear polarization and Faraday rotation in the M87 jet. We selected those data in which (i) different sub-bands are sufficiently separated in wavelength, (ii) both parallel and cross-hand visibilities are available, and (iii) M87 is observed as a primary target in full-track observing mode. Using these criteria, we are left with one data set at 8 GHz, four data sets at 5 GHz, and many data sets at 2 GHz. We note that there are multifrequency VLBA data obtained quasi-simultaneously in seven different sub-bands from 8.1–15.2 GHz in the literature (Zavala & Taylor 2002), which we could not find in the VLBA archive. Therefore, these data are not included in our analysis but we show that our results are consistent with the results of their work in Section 3.4. We found that the distribution of RM in the jet in different data sets at 2 GHz are more or less the same and chose three data sets among them for which all 10 VLBA antennas are available and the weather was good. We show the list of the eight VLBA archive data sets we analyzed and the basic information for each observation in Table 1. In total, we analyzed eight different polarization data sets of M87 taken by the VLBA (one at 8 GHz, four at 5 GHz, and three at 2 GHz).

A standard data post-correlation process was performed with the NRAO's Astronomical Image Processing System (AIPS, Greisen 2003). We corrected ionospheric dispersive delays using the ionospheric model provided by JPL, antenna parallactic angles, and instrumental delays using scans on bright calibrators. Amplitude calibration was performed by using the antenna gain curves and system temperatures with an opacity correction. We performed global fringe fitting with a solution interval between 10 and 30 s, assuming a point-source model. Bandpass calibration was performed by using scans on bright calibrators. The cross-hand R-L phase and delay offsets were calibrated by using scans on bright calibrators. We used the Caltech Difmap package (Shepherd 1997) for imaging and phase and amplitude self-calibration. We determined the feed polarization leakage (D-terms) for each antenna and for each sub-band by using the task LPCAL (Leppänen et al. 1995) in AIPS with a total intensity model of the D-term calibrators. We used OQ 208 or M87 for the D-term correction (Table 1) because of their very low degree of linear polarization (usually ≲1%).

The EVPA calibration was performed by comparing the integrated EVPAs of the VLBI maps of the calibrators with the EVPAs obtained in contemporaneous single dish polarization observations of the University of Michigan Radio Astronomy Observatory (UMRAO), or by using 3C 286 for which a stable integrated EVPA of ≈33° is known at the frequencies of our interest (Perley & Butler 2013), if available. However, we note that EVPA calibration is not critical for our purpose because the expected amount of EVPA correction for different sub-bands is almost the same. For example, we present the RM map and EVPA as a function of λ² at the map center of one of the calibrators in BC210B session, 0716+714, in Figure 1. Even though we could not perform EVPA correction for this epoch (see Table 1), the difference in EVPAs in different sub-bands is much smaller than the error bars and the obtained RM value is consistent with the previous measurements with the VLBA (Hovatta et al. 2012). We check the RM of the calibrators in all the data we analyzed to ensure that the detected RM for M87 is not due to potential errors in EVPA calibration but is intrinsic to the source itself.

Zoom In Zoom Out Reset image size

We obtained RM values at various positions in the M87 jet from measuring EVPAs in different sub-bands (intermediate frequencies) in each data set (see Table 1). We considered four error sources in EVPA: random error, systematic error induced by imperfect CLEAN procedures, by imperfect D-term calibration, and by imperfect EVPA calibration. We present the details of error estimation in Appendix A. For obtaining RM maps, we first convolved the maps in different sub-bands with the synthesized beam of the sub-band at the lowest frequency. Then, we fitted a linear function to the EVPAs from different sub-bands versus λ² for each pixel where the linear polarization intensity exceeds 1.5σ in all sub-bands, with σ being the full uncertainty, including D-term errors and CLEAN errors (Hovatta et al. 2012). We discuss the significance levels of the observed rotation measures (RMs) in Appendix B. We fitted the EVPA data several times, including potential nπ rotations, and used the fit that provided us with the lowest χ² value.

3.1. RM Maps

In Figure 2, we present example RM maps overlaid on the total intensity distribution of the jet for one observation at each frequency. The EVPA as function of λ² at three different locations of the jet is shown with good λ² fits. We obtained good λ² fits for the other RMs measured at different locations as well (some of them are shown in Figure 3) and also in the other five data sets not presented in Figure 2. We omitted the jet and RMs at ≈900 mas from the core at 2 GHz for better visualization; those data are presented in Figure 3. At lower frequencies, Faraday rotation is observable in more outward regions of the jet due to longer cooling times of the jet plasma. At higher frequencies, Faraday rotation is observable closer to the compact upstream emission thanks to better angular resolution and less depolarization. We note that the RM distributions are patchy at all frequencies because significant linear polarization is detected only in some parts of the jet, possibly due to substantial depolarization in the other parts. We also note that it is unlikely that those patchy RMs are artifacts because we found that the RMs in different epochs at the same observing frequency are detected in similar locations of the jet (see Appendix C).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Radial RM Profile

To obtain a radial RM profile along the jet, we calculated spatially binned RM by taking the weighted mean of all values in each separated region of the map with similar RM values. A priori, taking a weighted mean over a part of a map assumes that all individual pixels are independent from each other, which is not the case here. Pixels values are correlated across the extension of a resolution element (here, the synthesized beam). Thus, we first calculated a mean value, then its formal error (which assumes all pixels to be uncorrelated), and then multiplied this formal error by $\sqrt{n{{\rm{\Sigma }}}_{\mathrm{FWHM}}/{{\rm{\Sigma }}}_{\mathrm{RM}}}$, where n is the number of the pixels used for taking the mean, Σ_RM the size of the map region with RM values, and Σ_FWHM the area within the FWHM of the synthesized beam. We present the mean distance from the black hole, the binned RM values, and corresponding RM errors in Table 2. This data is used for our further analysis.

Table 2.  Binned RM Values

Session	Proj. Dist. (mas)	RM $(\mathrm{rad}\,{{\rm{m}}}^{-2})$	${\sigma }_{\mathrm{RM}}(\mathrm{rad}\,{{\rm{m}}}^{-2})$
(1)	(2)	(3)	(4)
BJ020A	10.87	−10163.47	3965.69
	15.55	−6381.18	2278.09
	18.35	−3421.25	541.23
	20.97	−3391.10	849.69
	24.69	−6054.02	615.75
	30.23	−5688.09	804.28
	33.37	−4489.37	2327.39
	51.38	−3908.49	1239.84
	66.79	−3374.41	1852.78
BJ020B	20.07	−3628.40	1049.11
	25.39	−4495.79	472.80
	30.05	−5698.63	3064.56
	43.83	−1932.51	1169.71
	59.31	−2022.79	1056.72
	156.90	−627.96	603.85
	163.58	−110.87	312.19
	171.72	−77.88	194.93
BC210B	30.73	−5414.65	1117.57
	41.89	−3101.13	1317.59
	50.59	−2229.21	1535.19
	92.36	−742.15	478.03
	160.65	−259.05	232.85
	170.49	161.17	201.12
BC210C	23.18	−2736.92	779.55
	92.35	−1345.83	919.99
	156.32	−45.10	414.45
	170.03	−32.60	136.35
	183.02	−501.42	505.80
BC210D	23.47	−3135.18	1653.81
	49.63	−2119.48	693.08
	168.70	−185.79	221.31
BH135F	149.44	−247.30	288.74
	170.00	81.37	198.19
	315.35	449.35	285.42
	347.12	434.49	196.92
	370.33	61.64	137.53
	873.80	1242.69	36.47
	899.69	1340.05	143.18
BC167C	149.35	−240.64	592.94
	169.17	12.77	233.97
	187.30	−88.17	460.41
	316.27	439.96	547.68
	326.85	314.65	358.40
	345.39	362.52	262.42
	368.23	−39.88	508.98
	865.47	1270.15	102.53
	883.03	1269.38	44.13
BC167E	170.44	6.70	217.93
	317.98	−59.62	388.00
	327.55	41.67	292.40
	346.71	−14.73	394.56
	368.51	167.98	446.63
	384.56	209.87	460.24
	866.06	1276.49	117.15
	883.04	1221.10	48.23

Note. (1) Project code of VLBA observations. (2) Mean projected distance from the black hole of the region where the RMs are measured, in units of milliarcseconds. (3) Binned RM values in units of rad m⁻². (4) 1σ errors of the binned RMs. All RM values are those before subtracting 130 rad m⁻² and the RM errors are before adding 300 rad m⁻² in quadrature to the uncertainties (see Section 3.3).

Download table as:  ASCIITypeset image

In Figure 4, we present the absolute values of RM as a function of de-projected distance from the black hole in units of r_s. We assumed a black hole mass of M_BH = 6 × 10⁹ M_⊙ (Gebhardt et al. 2011), a viewing angle of 17 deg⁸ (Mertens et al. 2016), and a distance between black hole and radio core as estimated by core-shift analysis of the M87 jet (Hada et al. 2011) to convert the observed projected jet distance from the radio cores into the de-projected distance from the black hole. Remarkably, the RM decreases systematically along the jet over nearly two orders of magnitude in distance (from 5000–200,000 r_s) inside the Bondi radius (3.6 × 10⁵ r_s; Russell et al. 2015). Our results substantially improved the radial RM profile of the M87 jet that was previously limited to a specific jet location at ≈20 mas from the core obtained in the pioneering RM study of the M87 jet (Zavala & Taylor 2002). The sign of the RM is preferentially negative inside the Bondi radius, except in the outer jet region (at distance of ≈2 × 10⁵ r_s) where RM errors are comparable to the RM values, which makes the RM sign ambiguous. However, at the location of HST-1 (at ≈4 × 10⁵ r_s), the observed RMs suddenly increase by a factor of ≈10 compared with those at ≈2 × 10⁵ r_s and their signs are always positive, which is opposite to the signature observed in the inner jet region (Figure 4). This result is in good agreement with the previous measurements by the Very Large Array (VLA) observations (Chen et al. 2011). Since we focus on the behavior of RMs inside the Bondi radius in this paper, we briefly discuss the results of RMs at HST-1 in Section 4.6 and more detailed results will be presented in a forthcoming paper (J. Park et al. 2018, in preparation).

Zoom In Zoom Out Reset image size

3.3. Contribution of RM Sources outside the Bondi Radius

3.4. Variability

Our data are obtained in different periods from 1996 to 2014, so RM variability might affect the results. We also included the results of a previous study of RM of the M87 jet (Zavala & Taylor 2002) for investigating potential RM variability. One can divide our data into four time groups, obtained in 1995–1996, 2000.48, 2006–2008, and 2013–2015. We show the absolute values of RM from different groups with different colors as a function of distance from the black hole in the left panel of Figure 5. The data obtained in different periods do not show significant deviation from each other. We also present the RM values as a function of time obtained in four different jet distance ranges, 15–40, 40–70, 100–200, and 200–400 mas, with different colors (the right panel of Figure 5). The mean values from different groups in the same jet distance range are consistent with each other within 1σ in almost all cases, suggesting that there is no significant temporal variability in RM. However, one exception is the case of positive RMs detected at ≈25 mas from the radio core in 2000.48 presented in the literature (Zavala & Taylor 2002). This value is larger than others obtained at similar jet distance by a factor of ≈2, and its sign is opposite. It is reasonable to consider that the positive RMs might be locally transient and not related to a global behavior of RM of the M87 jet because (i) the region of positive RMs is much smaller than that of negative RMs at a similar jet distance by a factor of ≈20 (Zavala & Taylor 2002) and (ii) positive RMs are not detected in other epochs and at other jet distances except in the outer jet region where RM errors are comparable to the RM values, which makes the RM sign ambiguous.

Zoom In Zoom Out Reset image size

3.5. The Faraday Screen

3.5.1. Internal Faraday Rotation and Depolarization

If the Faraday rotating electrons are intermixed with the synchrotron-emitting jet plasma, internal Faraday rotation can occur. Burn (1966) showed that the complex polarization (${ \mathcal P }$) of a synchrotron-emitting uniform slab with a purely regular magnetic field (see Sokoloff et al. 1998 for the case of a nonuniform or an asymmetric slab) is given by

$$\begin{eqnarray}&&{ \mathcal P }\equiv Q+{iU}={p}_{0}I\displaystyle \frac{\sin \phi {\lambda }^{2}}{\phi {\lambda }^{2}}{e}^{2i\left({\chi }_{0}+\tfrac{1}{2}\phi {\lambda }^{2}\right)},\end{eqnarray} \tag{ 1 }$$

where Q, U, and I are intensity in Stokes Q, U, and I maps, respectively, p₀ is the intrinsic fractional polarization, χ₀ the intrinsic EVPA, and ϕ the Faraday depth. However, internal Faraday rotation in sources with more realistic geometries and magnetic field structures usually results in deviation from a λ² law after total rotations ≳45° (Burn 1966; Sokoloff et al. 1998; Homan 2012).

We tested whether the observed degree of polarization and Faraday rotation can be explained with Equation (1) or not. We compared the degree of linear polarization expected in this model, ${p}_{{\rm{L}},\mathrm{internal}}={p}_{0}| \mathrm{sinc}(2\mathrm{RM}{\lambda }^{2})|$, with the observed one, p_L,obs. We assumed p₀ ≈ 0.75 because this is the maximum allowed degree of linear polarization for optically thin synchrotron radiation (Pacholczyk 1970). In Figure 6, we present p_L,obs/p_L,internal as a function of de-projected distance from the black hole. Most of the data points are much larger than unity, indicating that internal Faraday rotation in a uniform slab permeated by a regular magnetic field is not responsible for the observed jet RM. In addition, we frequently measure EVPA rotations larger than 45° with good λ² fits at various locations in the jet at all frequencies as shown in Figure 3. The fact that we could not find any statistically significant difference in RMs obtained at different frequencies at a given distance also supports an external origin (Figure 4). Thus, the systematic decrease of RM shown in Figure 4 must originate from the magnetized plasma outside the jet (external Faraday rotation).

Zoom In Zoom Out Reset image size

However, internal Faraday rotation might still be responsible for depolarization. As already noted in Section 3.1, at many locations of the jet linearly polarized emission is not detected in our data, making the RM distributions patchy. In general, depolarization originates from either internal Faraday rotation or spatial variations in the RMs of the external Faraday screen on scales smaller than the resolution of the observations (e.g., Burn 1966; Tribble 1991; Sokoloff et al. 1998; Homan 2012). The depolarization mechanism of AGN jet emission has been extensively investigated recently, thanks to observations with large bandwidths (e.g., O'Sullivan et al. 2012, 2017; Hovatta et al. 2018; Pasetto et al. 2018), or VLBI observations at many different observing frequencies (e.g., Kravchenko et al. 2017). Investigating the depolarization mechanism of the M87 jet is difficult for us because we have a limited number of observing frequencies with relatively short λ² spacings available. However, we found that the data collected in the BJ020A and BJ020B sessions could be combined because their observing dates and frequencies are relatively close to each other (Table 1).

We obtained the RM map as described in Section 2 after considering a core-shift effect (e.g., Lobanov 1998) between 5 and 8 GHz by employing 2D cross correlation of the optically thin emission regions in the image plane (Croke & Gabuzda 2008) and present the map in the left panel of Figure 7. We note that the results are not significantly affected by the core-shift. Significant RMs were detected in small parts of the jet because linear polarization at 5 GHz has not been detected in most parts of the jet in the inner jet region (at distances less than ≈60 mas), where the jet emission could be detected at 8 GHz. Nevertheless, an RM patch was detected at ≈25 mas from the core over a region with a size comparable to the beam size. In the right panel of Figure 7, we present Q/I, U/I, $p\equiv \sqrt{{Q}^{2}+{U}^{2}}/I$, and χ as a function of λ² in this region.

Zoom In Zoom Out Reset image size

In order to investigate the depolarization mechanism, we tried to model the Stokes I, Q, and U intensity simultaneously at different wavelengths, known as the qu-fitting technique (e.g., Farnsworth et al. 2011; O'Sullivan et al. 2012). We used a model for the complex polarization, which includes the effect of depolarization due to random magnetic fields (σ_RM) and ordered magnetic fields (ΔRM), given by

$$\begin{eqnarray}&&{ \mathcal P }={p}_{0}{{Ie}}^{2i({\chi }_{0}+\mathrm{RM}{\lambda }^{2})}{e}^{-2{\sigma }_{\mathrm{RM}}^{2}{\lambda }^{4}}\mathrm{sinc}{\rm{\Delta }}\mathrm{RM}{\lambda }^{2}\end{eqnarray} \tag{ 2 }$$

(Sokoloff et al. 1998). We followed a recent study that detected a very high RM of (3.6 ± 0.3) × 10⁵ rad m⁻² in the quasar 3C 273 with Atacama Large Millimeter Array (ALMA) observations at 1 mm (Hovatta et al. 2018) and fitted Equation (2) with σ_RM = 0 (model 1, the black dotted lines in the right panel of Figure 7) and with ΔRM = 0 (model 2, the red dashed lines) to the data points. The best-fit parameters are (p₀ = 0.10 ± 0.01, ΔRM = 532 ± 62 rad m⁻², χ₀ = 184° ± 6°, RM = −5195 ± 43 rad m⁻²) and (p₀ = 0.10 ± 0.01, σ_RM = 171 ± 25 rad m⁻², χ₀ = 184° ± 6°, RM = −5194 ± 43 rad m⁻²) for models 1 and 2, respectively. Both models can explain the data well with the reduced chi-square ${\chi }_{r}^{2}\equiv {\chi }^{2}/\mathrm{dof}$, where dof is the degree of freedom, of 0.66 and 0.64 for models 1 and 2, respectively. This is due to the sparse sampling of the data in the λ² space, which prevented us from solving the degeneracy. Nonetheless, the observed depolarization at ≈25 mas from the core is likely due to a gradient in RM by ≈532 rad m⁻² either in the jet or in the external Faraday screen across the beam or due to random magnetic fields with σ_RM ≈ 171 rad m⁻² in the external screen (Sokoloff et al. 1998; O'Sullivan et al. 2017; Hovatta et al. 2018; Pasetto et al. 2018). We also obtained good λ² fits for the EVPA rotation larger than 4π, supporting an external origin of the observed RM. The observed RM of ≈−5194 ± 43 rad m⁻² for model 2 is consistent with that obtained in the same location by using only BJ020A (8 GHz) data,−5535 ± 1226 rad m⁻², within 1σ and BJ020B (5 GHz) data, −4469 ± 431 rad m⁻², within less than 2σ. The deviation larger than 1σ in the latter case might be due to a non-negligible time gap of ≈18 days between the two data sets.

3.5.2. A Jet Sheath

If the Faraday screen is placed in the immediate vicinity of the jet, e.g., like a sheath surrounding the jet as claimed for other distant AGNs (e.g., Zavala & Taylor 2004; Jorstad et al. 2007; O'Sullivan & Gabuzda 2009; Hovatta et al. 2012; Park et al. 2018), then one expects significant RM gradients across the jet with a possible change of the sign of the RM; this is seen in numerical simulations (Broderick & McKinney 2010). This signature has indeed been frequently observed in the jets of many blazars (e.g., Asada et al. 2002, 2008; Gabuzda et al. 2004, 2015, 2018; Hovatta et al. 2012). The transverse RM gradients are related to toroidal magnetic fields in the jet and/or in the sheath, which can be naturally produced in the inner part of the accretion disk and/or in the black hole's ergosphere. These magnetic fields play a crucial role in launching and powering of relativistic jets (Meier 2012). MHD theories predict that poloidal magnetic fields that are dominant near the jet base become rapidly weak at larger distance and toroidal fields become dominant relatively far from the black hole (e.g., Vlahakis & Königl 2004; Komissarov et al. 2009).

However, for M87 the observed sign of RM is negative almost everywhere inside the Bondi radius (Figure 4). Furthermore, we found that there is no significant difference between the RMs on the northern and southern jet edges at a given distance, and the RMs appear to vary only as function of radial distance (see Appendix D). Recently, linear polarization structure of the core of the M87 jet at 43 GHz has been revealed, showing the inferred magnetic field vectors wrapped around the core⁹ (Walker et al. 2018). This suggests that toroidal fields might be dominant already on scales of ≈100 r_s, which makes it difficult to explain the observed single (negative) RM sign and no significant difference in RMs between the north and south edges with the Faraday screen consisting of a jet sheath.

We checked whether the observed RMs can be explained by the sheath model or not if poloidal magnetic fields are somehow dominant in the sheath at distance ≳5000 r_s, as indicated by a recent study of time variable RM in the radio core of a nearby BL Lac object Mrk 421 (Lico et al. 2017). We assumed (i) the same parabolic geometry of the sheath as that observed for the jet, i.e., z ∝ R^1.73 (Asada & Nakamura 2012; Nakamura & Asada 2013), with the radius of the outer boundary of the sheath being twice the radius of the jet (see the left panel of Figure 8), (ii) a constant velocity of the sheath at different distances, (iii) no reversal in the magnetic field direction along the line of sight, and (iv) the sheath consisting of nonrelativistic cold plasma. These assumptions led us to the scaling relations of n_e(z) ∝ R⁻² ∝ z^−1.16 and ${B}_{p}(z)\propto {R}^{-2}\propto {z}^{-1.16}$, with R being the radius of the sheath and B_p the poloidal magnetic field strength. We integrated $\mathrm{RM}\propto \int {n}_{e}(l)B(l){dl}$ numerically along each line of sight for each RM data point between the jet boundary and the sheath boundary (see the left panel of Figure 8) and fitted this function to the data points at different distances with a coefficient left as a free parameter. The best-fit model is indicated by the dashed line in Figure 4.

Zoom In Zoom Out Reset image size

3.5.3. Hot Accretion Flows

We use the scaling relations ${n}_{e}{(r)={n}_{e,\mathrm{out}}(r/{r}_{\mathrm{out}})}^{-q}$ with 0.5 ≤ q ≤ 1.5 and B(r) = B_out(r/r_out)⁻¹, where r is the radial distance from the black hole and n_e,out and B_out are the electron number density and the magnetic field strength at r_out, respectively. The former is based on self-similar solutions of hot accretion flows (Blandford & Begelman 1999; Yuan & Narayan 2014). The latter is based on the assumption that toroidal magnetic fields are dominant in the accretion flows (e.g., Hirose et al. 2004). We note that we are restricted to 1D scaling relations due to the limitations of the 2D accretion flow models, including non-negligible magnetic fields currently available, especially at small polar angles which is of our interest because of the small jet viewing angle (e.g., Mosallanezhad et al. 2016; Bu & Mosallanezhad 2018). In other words, we assume here that the quantities of the flows would be spherically symmetric for regions with a polar angle smaller than the jet viewing angle of 17°.

We employed $\mathrm{RM}=8.1\times {10}^{5}\int {n}_{e}(l)B(l){dl}$ (RM in units of rad m⁻², n_e in units of cm⁻³, B in units of Gauss, and l in units of parsec; Gardner & Whiteoak 1966) for "cold" nonrelativistic plasma, which applies to the relatively large spatial scales probed in this study (Yuan & Narayan 2014). We also performed numerical integration along each line of sight between the jet boundary and the Bondi radius (see the right panel of Figure 8, see also Section 3.3 for a discussion of the potential contribution by gas outside the Bondi radius). The result of fitting this function to the observed RM values measured inside the Bondi radius is indicated by the solid line in Figure 4 with the best-fit parameter of q = 1.00 ± 0.11, which indicates ρ ∝ r⁻¹, with ρ being the mass density. We could also obtain ${n}_{e,\mathrm{out}}{B}_{\mathrm{out}}$ from the fitting and when using ${n}_{e,\mathrm{out}}\approx 0.3\,{\mathrm{cm}}^{-3}$ at the Bondi radius measured by the X-ray observations (Russell et al. 2015), we obtain ${B}_{\mathrm{out}}\,=(2.8\pm 0.8)\times {10}^{-6}\,{\rm{G}}$.

4.1. Jet Sheath versus Hot Accretion Flows

In Section 3.5, we considered three different sources of Faraday rotation, (i) the jet itself, (ii) a sheath surrounding the jet, and (iii) hot accretion flows. Given that the observed EVPA rotations are larger than 45° at various locations in the jet with good λ² scalings, and the observed degree of linear polarization is usually much higher than that expected in the internal Faraday rotation model, we excluded the scenario (i) in Section 3.5.1. Although the hot accretion flows model (the solid line in Figure 4) apparently fits the data better than the sheath model (the dashed line in Figure 4), a statistical analysis is necessary to properly determine the better model. In Table 3, we present the values of reduced chi-square (${\chi }_{r}^{2}$) and Bayesian information criterion (BIC) obtained in the best fit for each model. The BIC is defined as $\mathrm{BIC}\equiv -2\mathrm{ln}{{ \mathcal L }}_{\max }+k\mathrm{ln}N$, where ${{ \mathcal L }}_{\max }$ is the maximum likelihood and $-2\mathrm{ln}{{ \mathcal L }}_{\max }$ is equivalent to the χ² value for the best-fit model in case of Gaussian errors (when neglecting a constant term), k the number of free parameters in the model, and N the number of data points used in the fit. The BIC allows one to compare the goodness of fit of different models with different numbers of free parameters (Schwarz 1978). The difference between the BIC values (ΔBIC) for two models quantifies how strongly one model is preferred over the other one, where a model with a lower BIC value is more favored by the data. Conventionally, 0 < ΔBIC < 2 represents weak evidence, 2 < ΔBIC < 6 positive evidence, 6 < ΔBIC < 10 strong evidence, and 10 < ΔBIC very strong evidence (e.g., Jeffreys 1961; Kass & Raftery 1995; Mukherjee et al. 1998; Liddle 2004). The value of BIC for the hot accretion flows model is smaller than that for the jet sheath model by ≈24 (Table 3), indicating that the former is strongly favored by the data.

Table 3.  Comparison of the Models

Model	ne Profile	B Profile	${\chi }_{r}^{2}$	BIC
(1)	(2)	(3)	(4)	(5)
Jet sheath	${n}_{e}(z)\propto {z}^{-1.16}$ (fixed)	$B(z)\propto {z}^{-1.16}$ (fixed)	1.73	86.8
Hot accretion flows	${n}_{e}(r)\propto {r}^{-1.00\pm 0.11}$ (fit)	$B(r)\propto {r}^{-1}$ (fixed)	1.16	62.4

Note. (1) Model applied to the RM data. (2) Density profile. The definition of z and r is explained in Figure 8. (fixed) means that the fixed profile is used in the model, whereas (fit) means that the index in the power law is left as a free parameter in the fitting. (3) Magnetic field strength profile. (4) Reduced chi-square for the best fit. (5) BIC for the best fit. The number of data points used in the fitting is 49.

Download table as:  ASCIITypeset image

We note that the above conclusion is based on the results obtained by using several assumptions on the jet sheath. For example, we assumed that the sheath geometry is the same as the jet, which may not be true. When we relax this assumption and leave the power-law index in the width profile of the sheath as a free parameter, i.e., z_sheath ∝ R_sheath^η, and fit the sheath model to the data points, then we obtain the best fit with η = 2.49 ± 0.17. This indicates that the sheath is more strongly collimated than the jet, which is unlikely because the inner part (closer to the axis) of streamlines is thought to be more collimated than the outer part for collimated outflows (e.g., Komissarov et al. 2007, 2009; Tchekhovskoy et al. 2008; Nakamura et al. 2018). Or, if we assume that toroidal fields are dominant in the sheath and fix the sheath geometry, we obtain a relatively good fit with a BIC value comparable to that of the hot accretion flows model. However, as noted in Section 3.5.2, it is difficult to explain the absence of a systematic difference between the RMs on the south and north edges in this case.

An alternative scenario is that the Faraday screen consists of dense clouds with ordered magnetic fields that are entrained by the jet (suggested by Zavala & Taylor 2002). The volume filling factor of these clouds, if they exist, is expected to be very small and this might explain why the RMs are detected in only small parts of the jet. Although we could not exclude this possibility, the observed depolarization at longer wavelengths does not seem to support this scenario. We present the distribution of the degree of linear polarization overlaid on the contours of total intensity emission for one observation at each frequency in Figure 9. At higher observing frequencies, the distribution of significant linear polarization becomes more continuous and the degree of linear polarization becomes higher at a given distance, notably at ≈20 and ≈170 mas from the core. This suggests that the Faraday screen consists of a continuous and extended medium such as winds but significant depolarization in large parts of the jet makes the observed patchy RM distributions especially at lower frequencies. We will investigate the depolarization mechanism at various locations of the jet with dedicated multifrequency polarimetric observations in the near future, which will allow us to identify the Faraday screen more rigorously.

Zoom In Zoom Out Reset image size

Taken as a whole, we conclude that attributing the Faraday screen to hot accretion flows is most consistent with the data presented in this paper and we discuss the results obtained by applying the hot accretion flows model hereafter.

4.2. Winds and the Faraday Screen

4.3. Jet Collimation by Winds

4.4. Misalignment

The dominance of a single RM sign for M87 implies that the background light source, i.e., the jet, exposes only one side of the toroidal magnetic loops in the Faraday screen. This situation can be realized when there is a misalignment between the jet axis and the symmetry axis of the toroidal field loops (Figure 10). This is another indication for winds or inflows as the dominant source of Faraday rotation because the jet sheath is tightly attached to the jet and cannot be tilted relative to the jet axis. Since the jet is highly collimated and narrow (Junor et al. 1999; Asada & Nakamura 2012; Doeleman et al. 2012), only a slight misalignment by ≈5° can result in observations of a fixed RM sign over a large distance range. Such small misalignments seem to be quite common in hot accretion flows even when the magneto-spin alignment effect, an alignment of the accretion disk and jets with the black hole spin by strong magnetic fields near the black hole, operates (McKinney et al. 2013).

Zoom In Zoom Out Reset image size

We note that it is unlikely that poloidal magnetic fields are responsible for the observed RMs of M87 because in that case one expects ρ ∝ r⁰ from B ∝ r⁻², which is impossible to explain with the accretion models currently available (Yuan & Narayan 2014). However, there is indication of non-negligible poloidal fields as well as toroidal fields—resulting in helical magnetic fields—in the jet environment of other distant AGNs, which results in transverse RM gradients with no sign changes (e.g., Asada et al. 2002; Zamaninasab et al. 2013; Gómez et al. 2016; Gabuzda et al. 2018, see also Section 3.5.2). The existence of non-negligible poloidal fields was indicated even for the M87 jet at HST-1 from the observed moving knots with both fast and slow velocities, which could be explained by quad relativistic MHD shocks in a helical magnetic field permeating the jet (Nakamura et al. 2010; Nakamura & Meier 2014). In Sections 3.5.2 and 4.1, we explained that poloidal magnetic fields might be very weak at distances ≳5000 r_s probed in this study and we concluded that hot accretion flows and winds are more probable to be the Faraday screen than the jet sheath. However, if the jet experiences recollimation, which may lead to formation of standing shocks (e.g., Daly & Marscher 1988; Gómez et al. 1995; Agudo et al. 2001; Mizuno et al. 2015; Martí et al. 2016; Fuentes et al. 2018), then the strength of poloidal fields could be substantially enhanced. Indeed, the width of HST-1 is significantly smaller than expected from the parabolic (conical) width profile inside (outside) the Bondi radius (Asada & Nakamura 2012), which has been explained with a hydrodynamic recollimation shock (e.g., Stawarz et al. 2006; Bromberg & Levinson 2009; Asada & Nakamura 2012). Also, the core of blazars is often identified with a recollimation shock (e.g., Daly & Marscher 1988; Marscher 2008; Cawthorne et al. 2013). This may explain the presence of non-negligible poloidal fields in the sheath of blazar jets and in HST-1, but not in the M87 jet inside the Bondi radius.

4.5. Mass Accretion Rate

4.6. RM at HST-1

4.7. EHT Observations

We studied Faraday rotation in the jet of M87 with eight VLBA data sets. We found that the magnitude of RM systematically decreases with increasing distance from 5000 to 200,000 r_s. Our work leads us to the following principal conclusions:

• 1.  
  We found that the degree of linear polarization in the jet is usually much higher than that expected in the case of internal Faraday rotation in a uniform slab with regular magnetic fields. In addition, we found that EVPA rotations are larger than 45° at various locations in the jet and always follow λ² scalings, which is difficult to reproduce with internal Faraday rotation in a synchrotron-emitting region with a realistic geometry and magnetic field structure. We conclude that the systematic decrease of RM must originate from the magnetized plasma outside the jet, supporting an external Faraday rotation scenario.
• 2.  
  We found that the observed sign of RM is predominantly negative inside the Bondi radius, without indication of significant difference in RMs detected on the north and south edges. The observed radial RM profile is difficult to explain with a sheath surrounding the jet permeated by poloidal magnetic fields being the Faraday screen. This implies that the Faraday screen consists of hot accretion flows, not of the jet sheath.
• 3.  
  We applied hot accretion flows model to the RM data points and obtained a best-fit function consistent with $\rho \propto {r}^{-1}$. This result is in good agreement with the ADIOS model in which substantial winds, nonrelativistic un-collimated gas outflows, are launched from hot accretion flows. The winds are likely surrounding the highly collimated relativistic jet and probably a dominant source of the observed RMs (Figure 10). However, we see indication for non-negligible contribution of inflows to the observed RMs as well.
• 4.  
  The density profile we obtained leads to the pressure profile of the winds, an external medium surrounding the jet, which is ${p}_{\mathrm{gas}}\propto {r}^{-5/3}$. This profile is consistent with a scenario in which the jet is substantially collimated by the winds, resulting in gradual acceleration of the jet in an MHD process. This is in agreement with the observed gradual collimation and acceleration of the jet inside the Bondi radius.
• 5.  
  The negative RM sign preferentially found inside the Bondi radius indicates that the jet exposes only one side of the toroidal magnetic loops in the Faraday screen (Figure 10). We conclude that the jet axis and the wind axis are misaligned with respect to each other. Since the jet is narrow, a slight misalignment by only ≈5° can lead to a fixed RM sign at distances ≳5000 r_s. According to recent GRMHD simulations (McKinney et al. 2013), such a (small) misalignment seems to be common in hot accretion flows, depending on the history of gas accretion, even when the magneto-spin alignment effect operates.
• 6.  
  The mass accretion rate can be substantially lower than the Bondi accretion rate due to the winds; we obtained ${\dot{M}}_{\mathrm{BH}}=1.6\times {10}^{-4}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, assuming a radial self-similarity of the density profile. This leads to a radiative efficiency of 3.8% at ${\dot{M}}_{\mathrm{BH}}/{\dot{M}}_{\mathrm{Edd}}=1.2\times {10}^{-6}$, which indicates that the radiative efficiency is not as small as usually assumed and the faintness of the nucleus of M87 is mainly due to the reduced mass accretion rate. Also, we obtained a jet production efficiency of ≳110%, implying that extraction of rotational energy of a spinning black hole might be at work in a MAD state.
• 7.  
  The RM at HST-1, located outside the Bondi radius, is larger by an order of magnitude and shows the opposite sign compared to the RM profile inside the Bondi radius. We conclude that this might be related with a recollimation shock that possibly forms in HST-1.

We conclude with several caveats that need to be addressed in future studies. We used simple 1D self-similar solutions for the density and magnetic field strength in the hot accretion flows model, while it is unclear whether this is valid or not. Studies of 2D solutions of hot accretion flows showed a breakdown of spherical symmetry (e.g., Mosallanezhad et al. 2016; Bu & Mosallanezhad 2018), though the behavior of physical parameters measured close to the jet axis is poorly constrained yet. We assumed $\approx 130\,\mathrm{rad}\,{{\rm{m}}}^{-2}$ for the contribution of the diffuse gas in M87 outside the Bondi radius based on the results of RM studies of the large-scale jet but this could be uncertain. We assumed that the same radial density profile holds for the polar region and for the equatorial region to estimate the mass accretion rate. This may not be true as seen in a recent study by Russell et al. (2018), though their results are obtained relatively close to the Bondi radius. We conclude that a sheath surrounding the jet is unlikely to be the Faraday screen based on the fact that the RMs detected on the southern and northern sides of the jet at a given distance are similar to each other. However, we could not test whether there are significant transverse RM gradients in the jet due to limited sensitivity and/or substantial depolarization. We plan to perform polarimetric observations with high sensitivity and having both short and long λ² spacings to constrain the origin of Faraday rotation more robustly and to investigate the depolarization mechanism in the near future.

We thank the referee for constructive comments, which helped improve the paper. The VLBA is an instrument of the Long Baseline Observatory. The Long Baseline Observatory is a facility of the NSF operated by Associated Universities, Inc. This research has made use of data from the University of Michigan Radio Astronomy Observatory, which has been supported by the University of Michigan and by a series of grants from the National Science Foundation, most recently AST-0607523. This work made use of the Swinburne University of Technology software correlator, developed as part of the Australian Major National Research Facilities Programme and operated under licence (Deller et al. 2011). We acknowledge financial support from the Korean National Research Foundation via Global PhD Fellowship Grant No. 2014H1A2A1018695 (JP) and Basic Research Grant No. NRF-2015R1D1A1A01056807 (S.T.). M.K. acknowledges the financial support of the JSPS KAKENHI program under Grant Nos. JP18K03656 and JP18H03721.

In this appendix, we present the details of error estimation for linear polarization quantities. We used the following equations to estimate the errors in linear polarization quantities (Roberts et al. 1994; Hovatta et al. 2012):

$$\begin{eqnarray}&&{\sigma }_{P}=\displaystyle \frac{{\sigma }_{Q}+{\sigma }_{U}}{2}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\sigma }_{\mathrm{EVPA}}=\displaystyle \frac{{\sigma }_{P}}{2P},\end{eqnarray} \tag{ 4 }$$

where σ_P, σ_EVPA, σ_Q, and σ_U are uncertainties in the polarized intensity, EVPA, Stokes Q, and U data, respectively. σ_Q and σ_U are estimated by adding different noise terms in quadrature, i.e.,

$$\begin{eqnarray}&&{\sigma }^{2}={\sigma }_{\mathrm{rms}}^{2}+{\sigma }_{\mathrm{Dterm}}^{2}+{\sigma }_{\mathrm{CLEAN}}^{2}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\sigma }_{\mathrm{Dterm}}=\displaystyle \frac{{\sigma }_{{\rm{\Delta }}}}{{\left({N}_{\mathrm{ant}}{N}_{\mathrm{IF}}{N}_{\mathrm{scan}}\right)}^{1/2}}{\left({I}^{2}+{\left(0.3{I}_{\mathrm{peak}}\right)}^{2}\right)}^{1/2}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\sigma }_{\mathrm{CLEAN}}=1.5{\sigma }_{\mathrm{rms}},\end{eqnarray} \tag{ 7 }$$

where σ_rms, σ_Dterm, and σ_CLEAN denote rms noise, D-term errors, and CLEAN errors, respectively. We estimated the rms noise in the residual maps after the CLEAN procedure in Difmap by shifting the maps by about a hundred times the beam size, corresponding to an off-center rms noise. We note that the rms noise in Stokes Q and U maps and in different sub-bands are similar and provide an average of them for each data set in (3) in Table 4. σ_Δ is the D-term scatter (provided in (4) in Table 4, see below), N_ant the number of antennas, N_IF the number of IFs, N_scan the number of scans with independent parallactic angles, and I_peak the peak of the total intensity map. We present the stations participating in the observations in Table 4. In our case, N_IF = 1 because we analyzed each sub-band data separately. We assumed N_scan = 8 because all our data sets observed M87 as a primary target in a full-track observing mode. We assumed the error from imperfect EVPA calibration of 3^◦ because relatively small errors are expected, as can be seen in Figure 1 (see also Section 2), and added this error to Equation (4) in quadrature.

Table 4.  Information About Data and Errors

Project Code	Stations	rms Error (mJy/beam)	D-term Scatter (%)	Pfalse	Nfalse	Nobs
(1)	(2)	(3)	(4)	(5)	(6)	(7)
BJ020A	VLBA	0.139	0.26	8.72 × 10−5	38	3081
BJ020B	VLBA	0.128	0.25	1.34 × 10−6	0	2824
BC210B	VLBA, −MK, −OV	0.090	0.41*	<3.64 × 10−7	0	3427
BC210C	VLBA, −MK, −KP	0.092	0.41*	<3.96 × 10−7	0	4947
BC210D	VLBA, −KP, $-\tfrac{1}{2}$ PT	0.074	0.41	2.04 × 10−5	7	1345
BH135F	VLBA	0.174	0.41*	2.50 × 10−7	0	1125
BC167C	VLBA	0.173	0.41*	<2.62 × 10−7	0	863
BC167E	VLBA	0.175	0.41*	1.33 × 10−5	2	775

Note. (1) Project code of VLBA observations. (2) VLBA stations participating in the observations (MK: Mauna Kea; OV: Owens Valley; KP: Kitt Peak; PT: Pie Town). (3) Averages of off-center rms errors in Stokes Q and U maps in units of mJy/beam. (4) Scatter in the D-terms obtained with different sources in units of %. We could not derive reliable D-term scatters in some data sets (marked with *) and thus assumed that the errors for these data sets are similar to that of session BC210D, 0.41% (see Appendix A for more details). (5) Probability of detecting false RMs with the RM values and ${\chi }_{r}^{2}$ similar to the observed ones (Appendix B). P_false is obtained from integrating the false-alarm probability distribution function (FPDF) between the minimum and maximum observed RMs for each data set presented in Table 2. The symbol < in front of the values for some sessions means that we could not find any pixel of false RM and we provide an upper limit. (6) Number of pixels of false RMs expected to be seen in the jet. (7) Number of pixels of observed RMs in the jet.

Download table as:  ASCIITypeset image

We estimated the D-term scatters by comparing the D-terms obtained from different sources. However, this was not always possible because some data sets did not have more than one source suitable for D-term calibration or because of a small number of scans (less than three) on other D-term calibrators. Specifically, we obtained reliable D-terms for BJ020A and BJ020B using three sources, OQ 208, OJ 287, and M87, because all of these are suitable for D-term calibration, i.e., either weakly polarized or moderately polarized but having compact geometries, and they are observed in multiple scans over large parallactic angle ranges (see e.g., Roberts et al. 1994; Aaron 1997; Park et al. 2018 for details of D-term calibration). We present the D-terms of different antennas obtained from different sources in Figure 11 (top for BJ020A and middle for BJ020B). The scatter in the D-terms obtained by using different sources is about ≈0.25% for both left-handed circularly polarized (LCP) and right-handed circularly polarized (RCP) data.

Zoom In Zoom Out Reset image size

For the BC210 data sets, we could obtain a reasonably small scatter of ≈0.4% only for the BC210D data (the bottom panel of Figure 11), using M87 and 0716+714, because the number of scans on 0716+714 is at most two or three and two antennas were missing in the other two data sets. However, the D-terms measured by using M87 for BC210B and BC210C data sets are likely quite reliable because we obtained clear linear polarization in all different sub-bands that are consistent with the results of BJ020B and BC210D (Figure 13). Thus, we assumed that the D-term scatters for these data sets are similar to that of BC210D and used 0.41%.

For the L band (2 GHz) data sets, we could not obtain the D-term scatters because only three sources, M87, 3C 273, and 3C 286, were observed. M87 can serve as a good D-term calibrator thanks to its very low degree of linear polarization. However, the other two sources show quite strong (≳10%) linear polarization over large extended jet regions and thus they are not suitable for D-term calibration. We assumed that the D-term scatters of these data sets are the same as those of BC210D, i.e., 0.41%. This is because the D-term scatters of the VLBA tend to be larger at higher observing frequencies (e.g., Gómez et al. 2002).

We obtained RM for each pixel where the linear polarization intensity exceeds 1.5σ in all sub-bands, with σ being the full uncertainty (Equations (3) and (5)). As there are at least four independently processed sub-bands per data set, the total (Gaussian) probability of false detection of RM is <3.2 × 10⁻⁴. However, linear polarization intensity does not follow a Gaussian probability distribution for a small signal-to-noise ratio (S/N; Wardle & Kronberg 1974; Trippe 2014). Thus, we need to carefully check the potential chance of detection of artifacts in the observed RMs. This kind of test has been done by performing extensive simulations in previous studies (e.g., Roberts et al. 1994; Hovatta et al. 2012; Algaba 2013; Mahmud et al. 2013). They generate simulated data sets with the known polarized intensity distributions, e.g., a uniform fractional polarization and EVPA across the source's total intensity structure, and add errors introduced by various effects discussed in Appendix A. The significance level can be inferred from the number of simulated data sets where the input polarized model is distorted.

However, this approach might not apply to our study because the observed linear polarization is very patchy in all data sets possibly due to substantial depolarization (Section 3.5.1). We present an alternative approach to infer the significance levels of the observed RMs. If the criterion of 1.5σ cutoff is not strict enough and this introduces many false RMs in the jet, then one should expect to see many similar RMs outside the jet region (where the total intensity emission of the jet is not significant) as well. This is because all the error sources, i.e., random errors, CLEAN errors, and D-term errors can be distributed across the entire map, not specific to the jet region. Although the D-term errors depend on the total intensity (Equation (6)), this intensity is usually smaller than ≈30 mJy/beam where significant RMs are observed in the jet. Thus, the second term in Equation (6) is always dominant, and it would be fair to compare the observed RMs in the jet with the false RMs outside the jet region generated by errors.

For the regions outside the jet, we computed the number of pixels that satisfy the following two conditions: (i) polarized intensity above 1.5σ is detected in all sub-bands and (ii) ${\chi }_{r}^{2}\lesssim 1.1\mbox{--}1.5$ are obtained for the λ² fit to the EVPAs, similarly to the observed RMs (We note, however, that the results are not significantly changed when we did not consider ${\chi }_{r}^{2}$.) This calculation was done by using the maps having similar fields-of-view to those of the jet to properly compare with the observed jet RMs and to avoid the bandwidth-smearing and the time-average smearing effects. We obtained histograms of false RMs and divided them by the total number of pixels outside the jet region in the maps, which can serve as the FPDFs of detecting RM.

In Figure 12, we present the FPDF for the 8 GHz data as an example. The probability of detecting false RMs with −10,163 ≲ RM ≲ −3374 rad m⁻² (the range of observed RMs at 8 GHz, see Table 2) with good λ² fits, obtained from integrating the hatched region, is 8.72 × 10⁻⁵ (P_false). Accordingly, one can expect to detect false RMs in the jet region in approximately 38 pixels (N_false), while significant RMs are detected in more than 3000 pixels (N_obs). We present the values of P_false, N_false, and N_obs for all data sets in Table 4. The values of P_false for the other data sets, obtained by integrating the FPDF between the minimum and maximum observed RMs (Table 2) for each data set, are even smaller, resulting in very small or zero N_false.

Zoom In Zoom Out Reset image size

In Table 5, we present P_false obtained by using five different S/N cutoffs. When 1σ cutoff is used, P_false values are non-negligible, up to ≈2 × 10⁻³ for BJ020A data set. However, P_false decreases rapidly as the cutoff level increases for all data sets, becoming smaller than ≈9 × 10⁻⁵ with 1.5σ cutoff. Therefore, we conclude that almost all of the observed RMs obtained by using the 1.5σ cutoff is intrinsic to the source.

We present the RM maps for the whole 2 and 5 GHz data sets in Figure 13. The RMs in different epochs at the same observing frequency are detected in similar locations of the jet, notably at ≈170 and ≈320–370 mas from the core at 2 GHz and at ≈20–30 and ≈150–200 mas from the core at 5 GHz. We note that the difference in the locations of some RMs might be due to relatively large time gaps (3 months to 17 yr) between different data sets, given that the jet is known to move relativistically already at distances less than ≈10 mas from the core (Mertens et al. 2016; Hada et al. 2017; Walker et al. 2018).

Zoom In Zoom Out Reset image size

In Figure 14, we present the absolute values of RM as a function of de-projected distance for the RMs detected on the northern and southern jet edges with different colors. We determined whether the observed RMs are located in the north or south edges by comparing the position of the RMs with the brightness centroid of the transverse intensity profile at the given distances. We found that out of 49 regions where significant RMs are detected, 10 are located in the southern jet edges.

Zoom In Zoom Out Reset image size