© The Authors 2023

1. Introduction

M 87 is known as the best target for investigating the active galactic nucleus (AGN) jet formation mechanism due to its proximity (distance = 16.7 Mpc; Blakeslee et al. 2009) and its large central supermassive black hole (SMBH; Macchetto et al. 1997; Gebhardt & Thomas 2009; Gebhardt et al. 2011; Walsh et al. 2013). Recently, the Event Horizon Telescope (EHT) successfully imaged the first ever black hole shadow (Event Horizon Telescope Collaboration 2019, 2021) and firmly determined the mass of the central SMBH in M 87 (M_• = 6.5 ± 0.7 × 10⁹ M_⊙; Event Horizon Telescope Collaboration 2019). At the distance of 16.7 Mpc, 1 milli-arcsecond (mas) corresponds to ≈130 r_s ≈ 0.08 parsec, where r_s is the Schwarzschild radius of the SMBH. M 87 hosts a prominent radio jet that is detected from the radio to TeV γ-rays (e.g., Abramowski et al. 2012; EHT MWL Science Working Group 2021).

The current leading model for jet launching describes a jet that is driven by magnetic force and then accelerated by the relatively slow conversion of magnetic energy into kinetic energy (e.g., Blandford & Znajek 1977; McKinney 2006; Komissarov et al. 2007; Nakamura et al. 2018; Chatterjee et al. 2019). One way to test this scenario is to determine the radial profiles of the velocity field and the magnetic field along the jets. Over the last 15 years or so, the velocity field of the M 87 jet has been intensively explored via multi-epoch very long baseline interferometry (VLBI) observations that have allowed us to directly probe the kinematics of the jet components (e.g., Kovalev et al. 2007; Ly et al. 2007; Asada et al. 2014; Mertens et al. 2016). Furthermore, in order to avoid errors in the identification of knot motions due to low cadence, high-cadence observations were made and the global profile of the M 87 jet velocity field from 10² r_s to 10⁷ r_s was revealed in great detail by Park et al. (2019a). Interestingly, the obtained velocity field data showed an apparent discrepancy with those predicted by general relativistic magnetohydrodynamics simulations (e.g., McKinney 2006; Nakamura et al. 2018), posing a new challenge to the jet formation mechanism. Future VLBI observations with higher angular resolution and higher cadence may allow us to address this issue.

Meanwhile, there has been recent progress in constraining the magnetic field strength of the M 87 jet (for a review, see Hawley et al. 2015). In the bright knot HST-1, which is located in the ∼3.8 × 10⁵ r_s downstream of the SMBH, the magnetic field strength is constrained by the synchrotron cooling time in the X-ray energy band (Harris et al. 2003, 2009). In the region within ∼100 r_s, there has been major progress in millimeter and submillimeter VLBI observations, which successfully measured the size and the flux of the jet base (the radio core), where the optical depth for synchrotron self-absorption (SSA) becomes unity at the observed frequency (e.g., Doeleman et al. 2012; Hada et al. 2013, 2016). From this, the magnetic field strength of the radio core was estimated. The estimated strength suggests that the energy density of the magnetic field dominates the nonthermal electron energy density near the SMBH, which in turn suggests that magnetic fields are important for jet formation (Kino et al. 2014, 2015). The successful measurement of the frequency-dependent position of the core (the so-called core-shift) provides another way of constraining the magnetic field strength of the jet (Hada et al. 2011; Zamaninasab et al. 2014; Zdziarski et al. 2015). Most recently, from images of the central black hole’s photon ring and its polarization distribution obtained from EHT observations (Event Horizon Telescope Collaboration 2019, 2021), the magnetic field strength of the surrounding plasma was constrained.

However, our understanding of the magnetic field profile along the jet is currently limited. This is because only the magnetic field strengths of notable individual components have been estimated so far. The magnetic field strength between the radio core and HST-1 remains unknown. In particular, there has been no attempt to date to measure the magnetic field strength of the extended jet.

In this study we explore the magnetic field properties of the M 87 jet using spectral index maps. Although there have been many previous studies that have published spectral index maps of M 87 (Zavala & Taylor 2003; Dodson et al. 2006; Ly et al. 2007; Hada et al. 2011, 2012; Pushkarev & Kovalev 2012; Niinuma et al. 2014; Hovatta et al. 2014; Asada et al. 2016; Kravchenko et al. 2020), there has been no dedicated study trying to limit physical quantities using spectral index maps. In order to thoroughly constrain and characterize the radial profiles of physical quantities along the M 87 jet, we set up a dedicated observation program using the KVN and VERA Array (KaVA). After commissioning observations (Niinuma et al. 2014), M 87 has been observed as part of the KaVA AGN Large Program (hereafter, LP) since 2016.

This paper is organized as follows. In Sect. 2 we describe the multifrequency observations and data reduction of the KaVA LP and other archival data. In Sect. 3 we present the spectral index maps of M 87 between 22 and 43 GHz. In Sect. 4 we discuss a model of the spectral index distribution that can explain the observed spectral index profile of M 87 and constrain the magnetic field distribution of the jet. Furthermore, we compare our measurements to those from previous studies. Throughout the paper, we define the sign of the spectral index as S ∝ ν^+α.

2. Observations and data reduction

2.1. KaVA observations

M 87 observations using KaVA were performed nine times from February to June 2016. Each epoch consists of two sessions, at 22 and 43 GHz, which were observed within 1 to 2 days of each other. Each observation was allocated every 2−3 weeks. Among them, the data on June 1 were excluded due to antenna problems at Mizusawa and Ishigaki. As a result, eight out of the nine data sets are used for the spectral analysis. At each frequency, the total observation time at each epoch is approximately 7 h, and M 87’s on-source time is about 4 h and 30 min. More information about the observations is described in Park et al. (2019a). Data reduction was performed following the process described in Park et al. (2019a). The amplitude is calibrated via an a priori method using an opacity-corrected system temperature and the elevation-dependent gain curve for each telescope. The amplitudes were scaled up by a factor of 1.3 to correct for amplitude losses that occurred during correlation (Lee et al. 2015; Hada et al. 2017). However, since the same factor is applied at both frequencies, the spectral index maps do not change after the correction. Imaging with CLEAN and self-calibration was performed using the Difmap software package (Shepherd et al. 1994). When natural weighting is applied, the typical size of the full width at half maximum of the synthesized beam at 22 and 43 GHz is close to a circular shape with radii of 1.2 mas and 0.6 mas, respectively.

2.2. VLBA archival data

We also included archival quasi-simultaneous data from the Very Long Baseline Array (VLBA). A total of five observations, with two in 2010 and three in 2014 were obtained. Details of the observations for the data in 2010 and 2014 are described in Hada et al. (2011) and Hada et al. (2016), respectively. We re-performed the data reduction for these data. Initial data calibrations (a priori amplitude correction, atmospheric opacity correction, fringe-fitting, and bandpass) were performed using the Astronomical Image Processing System (AIPS; Greisen et al. 2003). For the 2014 observations using the RDBE digital backend system, the auto-correlation amplitude was further corrected after bandpass correction, as suggested in VLBA Memo #37. The subsequent image reconstruction was performed in Difmap based on the usual CLEAN/self-calibration procedure. The typical beam size is 0.8 × 0.4 mas at 22 GHz, and 0.45 × 0.22 mas at 43 GHz when natural weighing is applied.

3. Results

3.1. Spectral index maps

When observations are made at two frequencies in the same array, the low-frequency data are distributed over a shorter distance in the (u, v)−plane than the high-frequency data (e.g., see Fig. 1 in Park et al. (2019a) for the (u, v)-coverage of KaVA 22 and 43 GHz).

Creating spectral index maps using data with a nonidentical (u, v) range can cause the spectral index to become artificially steep in the extended region, where the sensitivity of the lower frequency is higher than that of higher frequency. Therefore, we excluded data on the long baseline of 43 GHz and the short baseline of 22 GHz, respectively. The resulting (u, v) range of the data is 33−170 Mλ for KaVA and 25−685 Mλ for VLBA. More details on the (u, v) range matching and its effects are described in Appendix A.

Then, we restored all the maps to the same circular beam of 1.2 mas × 1.2 mas (∼160 r_s), which is the comparable size of the synthesized beam of KaVA at 22 GHz. This beam size is 1.5 times larger than the original size of the VLBA’s synthesized beam. The selection of the size of the restored beam may affect the final spectral index map and further analysis. We explored this in Appendix B and found that the main conclusions are not changed.

During the self-calibration process, the absolute coordinate position of the source is lost and the brightest component in the image (the radio core) is shifted to the center of the map. Therefore, an additional step is required to align images of different frequencies. For the M 87 jet, the positions of the radio core at different frequencies have been measured using the phase-referencing technique (Hada et al. 2011). According to their asymptotic relation, r_RA = Aν^−α + B (α = 0.94 ± 0.09, A = 1.40 ± 0.16 and B = −0.041 ± 0.012), the position difference between the 22 GHz core and the 43 GHz core is $0.{036}_{-0.010}^{+0.013}$ mas in right ascension. Assuming the jet position angle is −72° (Walker et al. 2018), then the difference in declination is $0.{012}_{-0.003}^{+0.004}$ mas. We took into account this shift when aligning the maps. Recently, it has been shown that the core-shift of AGN jets could change with time (Plavin et al. 2019). However, the phase-referencing observations of the M 87 jet showed little change in 22−43 GHz core-shifts from 2010 to 2019 (≲10 micro-arcsec; Hada et al. 2012, 2014; Jiang et al. 2021), and thus we use a constant core-shift in our analysis.¹

After aligning the images of the two frequencies, the spectral index map is created by calculating the spectral index of each pixel using the following relation:

$$\begin{array}{c}\hfill \alpha =\frac{log(S_{{\nu }_2}/S_{{\nu }_1})}{log({\nu }_2/{\nu }_1)},\end{array}$$(1)

where ν₁ and ν₂ are the two observed frequencies, and S_ν1 and S_ν2 are the flux densities at each frequency. The error of the spectral index at each pixel was estimated according to the discussion in Kim & Trippe (2014); The error of the intensity of each pixel (i, j) is considered to be the sum of the systematic error and the thermal random noise (i.e., σ_Iν, ij = δ_νI_ν, ij + σ_rmsν). The factor of systematic amplitude error is empirically assumed as δ ∼ 10% (e.g., Hada et al. 2012; Niinuma et al. 2014; Hovatta et al. 2014; Cho et al. 2017). The thermal random noise (σ_rmsν) is obtained from the residual map after the CLEAN process for each individual epoch. Then, the error of the spectral index at each pixel is calculated as follows:

$$\begin{array}{c}\hfill {\sigma }_{\alpha ,ij}=\frac{1}{log({\nu }_2/{\nu }_1)}\times \sqrt{{\left(\frac{{\sigma }_{I_{{\nu }_1},ij}}{I_{{\nu }_1,ij}}\right)}^2+{\left(\frac{{\sigma }_{I_{{\nu }_2},ij}}{I_{{\nu }_2,ij}}\right)}^2}.\end{array}$$(2)

The spectral index error is typically σ_α ∼ 0.5 near the core and increases to σ_α ≳ 1 in the downstream jet. Because of the uncertainty in the core-shift, additional errors can occur during image alignment. We measured the amount of error due to image alignment for some epochs. As a result, σ_α ∼ 0.1 was found in the core region but it is almost negligible in the downstream jet (σ_α ≲ 0.02). Since the spectral index error is mainly resulting from the error of the intensity, the uncertainty of image alignment is not considered in this work.

The left and right panels in Fig. 1 are spectral index maps of the M 87 jet using 22 and 43 GHz images overlaid with 22 GHz contours observed by the VLBA and KaVA, respectively. All images have been rotated by −18° in order to align the jet central axis with the horizontal axis (e.g., Walker et al. 2018). The color indicates the spectral index between 22 and 43 GHz at each pixel. We blank pixels where the total intensity is less than 3σ_rmsν at one of the observed frequencies. In most maps, spectral indices can be obtained at distances up to 10 mas from the core, which is more than double than had been previously observed (Ly et al. 2007; Niinuma et al. 2014; Kravchenko et al. 2020). In particular, for some KaVA observations, spectral indices can also be obtained from structures located at ≈20 mas. However, we are not concentrating on this structure in this paper.

Fig. 1. Spectral index maps between 22 and 43 GHz obtained from VLBA and KaVA observations. All images have been rotated by −18°. The restoring beam size is 1.2 mas × 1.2 mas, drawn as a black circle in the bottom-left corner. Observing dates are shown to the left of each map. The contours represent the total intensity at 22 GHz. Contours start at 3σrms, increasing in steps of 2.

It is notable that in general, the spectral index morphology appears to be quite similar over the course of the observations. While there are differences from epoch to epoch, the general structures remain remarkably consistent.

In the spectral index maps, the M 87 jet shows a flatter spectrum at the core, while the extended jet has a steeper spectrum. These spectral distributions are commonly found in many AGN jets (e.g., O’Sullivan & Gabuzda 2009; Pushkarev & Kovalev 2012; Hovatta et al. 2014). We determined the weighted average spectral indices of the core from a 0.5  ×  0.5 mas box centered on the location of the peak flux density. We find α_core, KaVA = −0.25 ± 0.50 and α_core, VLBA = −0.20 ± 0.55, respectively. The spectral index values obtained from both facilities are identical within the error bars, giving an almost flat or slightly steep spectrum. The spectral index of the M 87 jet core at 22−43 GHz is more or less consistent with what has been previously reported (Hada et al. 2012; Niinuma et al. 2014; Kim et al. 2018a; Kravchenko et al. 2020), although it is not identical².

We note that there is a sudden spectral index increase at the edge of the jet for a few epochs (March 26, 2014, and May 7, 2014). However, these are possibly spurious features because they only appear intermittently in a few epochs and then disappear immediately. Such features are frequently found in many spectral index maps (e.g., O’Sullivan & Gabuzda 2009; Müller et al. 2011; Fromm et al. 2013; Hovatta et al. 2014; Boccardi 2015). One possible reason is the sparse sampling of the (u, v)-coverage, especially on short ranges, leading to strong variations in the outermost structure. In the following analysis, we exclude these high spectral index features at the edge of the jet.

3.2. Spectral index distribution along the M 87 jet

The spectral distribution along the jet can be studied in several ways. One method is to extract the spectral index values along the ridge line of the jet (e.g., Fromm et al. 2013; Hovatta et al. 2014). However, it is not straightforward to apply this method in our study because the M 87 jet is known to consist of double ridges (e.g., Hada et al. 2013, 2016). Instead, we constructed the spectral index distribution of the jet by taking the weighted average of the spectral indices in the direction perpendicular to the jet as it moves down the jet pixel by pixel:

$$\begin{array}{c}\hfill \overline{\alpha }=\frac{{\sum }_{i=1}^n({\alpha }_i{\sigma }_{\alpha ,i}^{-2})}{{\sum }_{i=1}^n{\sigma }_{\alpha ,i}^{-2}},\end{array}$$(3)

where α_i is the pixel value of the spectral index map, σ_α, i is the value obtained from the spectral index errors derived from Eq. (2), and n is the number of the pixels used for averaging. When calculating the weighted mean, we used pixel values located within one beam size from the jet axis (±1.2 mas). The reason is that it is wide enough to cover the spectral indices of the two ridges, while at the same time excluding edge structures with spurious features at some epochs (see Sect. 3.1). The error is obtained by multiplying the formal error of the weighted mean by $\sqrt{n}$ to compensate for the biases due to correlations between pixel values (Park et al. 2019b).

Figure 2 shows the radial distributions of the spectral index as a function of deprojected distance from the SMBH in units of r_s. Distributions obtained from KaVA and the VLBA are indicated by red and blue lines, respectively. The positional difference between the map center and the SMBH is corrected using the results of Hada et al. (2011). We assumed a viewing angle of 17° (Walker et al. 2018) to convert the observed projected jet distance to a deprojected distance, where 10 mas corresponds to ∼4500 r_s. The spectrum of the innermost jet is relatively flat as was previously described in Sect. 3.1. Downstream of the jet base the spectrum becomes steeper. At a distance of about 2 mas (∼900 r_s), the spectral index becomes α ∼ −0.7, which is similar to the average spectral index of the M 87 jet between 22−43 GHz reported in previous studies using the same VLBA data as used in this study (Hada et al. 2012, 2016). This value is also similar to the typical spectral index of jets observed in many AGNs (e.g., Pushkarev & Kovalev 2012; Hovatta et al. 2014).

Fig. 2. Evolution of the spectral index between 22 and 43 GHz along the deprojected distance from the SMBH. Weighted mean values across the jet are used in this distribution (see Sect. 3.2). The projected distance is displayed on the upper axis in mas, and the corresponding deprojected distance is displayed on the lower axis in rs. Red lines represent the spectral index distribution from 2016 KaVA observations and blue lines the spectral index distribution from VLBA observations from 2010 and 2014. The shaded area is the 1σ error of the spectral index.

Two interesting features are found in Fig. 2 in the region farther than ∼2 mas. First, the spectral indices obtained from both VLBA and KaVA decrease rapidly with distance, reaching a steep spectrum of α ∼ −2.5. This high rate of decline in the spectral index with distance and the resulting very steep spectrum have not been reported in previous works on the M 87 jet. Second, there are hints that the spectral index stops declining at ∼6 mas in the KaVA data. Interestingly, a constant spectral index distribution after a certain distance is often found in many other AGN jets (e.g., Kovalev et al. 2008; Hovatta et al. 2014; Haga et al. 2015; Boccardi 2015; Lisakov et al. 2017; Pushkarev et al. 2019; Baczko et al. 2019; Park et al. 2021). However, it is not clear that this is also seen in the VLBA data. We cannot completely rule out the possibility that the spectral index may continue to decline in this region. There is a time difference of about ∼2 years between the KaVA and the VLBA observations. Therefore, the spectral distribution of the jet in this region could vary with the observation time.

4. Discussion

4.1. Comparison with previous works

In Sect. 3 we were able to obtain spectral index distributions of the M 87 jet up to 10 mas (∼4500 r_s) from the core for the first time thanks to the good image quality of KaVA and the VLBA at 22 and 43 GHz. An important characteristic of the distribution is that the spectral index decreases between 2−6 mas (∼900 r_s – ∼2500 r_s) from α ∼ −0.7 to α ∼ −2.5. In the KaVA data, it does not decrease after ∼6 mas.

The change of the spectral index can be obtained from the change in the energy distribution of nonthermal electrons N(γ), since the slope of the electron distribution (p) and the spectral index is related by α = (p + 1)/2 (Rybicki & Lightman 1979). The evolution of N(γ) as a function of time can be investigated by solving the transfer equation (or continuity equation) including nonthermal electron injection and energy losses.

One previous study discussing the synchrotron spectrum of parsec-scale AGN jets has proposed two possibilities for spectral steepening using analytical solutions of the transfer equations for two nonthermal electron injection scenarios (Hovatta et al. 2014). (1) In the case of a constant and continuous injection of nonthermal electrons, the synchrotron spectrum has a break at γ_break, and the spectral index steepens by 0.5 before and after the break. (2) If there is no injection after the initial injection, the cutoff energy γ_max decreases with time. At the cutoff energy, the synchrotron spectrum decreases exponentially, and thus the spectral index decreases to an arbitrarily low value.

However, both scenarios cannot explain the properties observed in the spectral index distribution of the M 87 jet (Fig. 2). Firstly, the observed decrease in the spectral index is much greater than 0.5. Also, the decrease in the spectral index may not continue and may stop at a certain value, remaining at that value after that.

Therefore, in this study, we attempt to explain the properties of the observed spectral index distribution by numerically solving the transfer equation. Section 4.2 describes the transfer equation, and assumptions for physical parameters to solve the equation numerically, including magnetic field, jet radius, bulk jet velocity, and electron injection. In Sect. 4.3, the importance of electron injection is discussed in order to explain key features of the observed spectral index distribution. In Sect. 4.4 we try to constrain the magnetic field strength and the electron injection of the M 87 jet at the distance of 2−10 mas by comparing the observed spectral index distribution to the model. Then, in Sect. 4.5, we compare our estimate for the magnetic field strength with previous estimations.

4.2. Spectral index model

In the frame of a cross section of the jet co-moving with the jet’s flow, the transfer equation is written as (e.g., Ginzburg & Syrovatskii 1964; Longair 2011; Blasi 2013)

$$\begin{array}{c}\hfill \frac{N(\gamma ,\tau )}{\partial \tau }+(\mathrm{\nabla }·v)N(\gamma ,\tau )+\frac{\partial }{\partial \gamma }[b(\gamma ,\tau )N(\gamma ,\tau )]=Q(\gamma ,\tau ).\end{array}$$(4)

Here, N(γ, τ) is the number density of the nonthermal electrons as a function of energy γ and time τ, v is the velocity of the system, Q is the nonthermal electron injection rate, and b is the energy loss rate. In the co-moving frame, the bulk velocity along the jet direction is zero and only the expansion (or contraction) rate along the radial direction remains. In our discussion, we assume that adiabatic losses and synchrotron losses are the dominant energy loss processes. Then, b can be written as (e.g., Rybicki & Lightman 1979; Longair 2011)

$$\begin{array}{c}\hfill b(\gamma ,\tau )=\frac{\mathrm{d}\gamma }{\mathrm{d}\tau }=-b_{\mathrm{adi}}\gamma -b_{\mathrm{sync}}{\gamma }^2,\end{array}$$(5)

where $b_{\mathrm{adi}}=\frac{1}{R}\frac{\mathrm{d}R}{\mathrm{d}\tau }$ and $b_{\mathrm{sync}}=\frac{4}{3}\frac{{\sigma }_T}{m_{\mathrm{e}}c}\frac{B^2}{8\pi }$ is the coefficient of the adiabatic losses and the synchrotron losses, respectively. Here, R is the radius of the cross section, and B is the magnetic field strength in the co-moving frame. We assumed that the pitch angle distribution is rapidly isotropized during synchrotron cooling³.

In order to numerically solve Eqs. (4) and (5) along the jet, several physical quantities should be known as a function of distance: the bulk jet velocity profile (Γ(z), where Γ is the Lorentz factor of the bulk jet velocity), the jet radius profile (R(z)), the magnetic field strength profile (B(z)), and the nonthermal electron injection function (Q(z)). Among them, we take the bulk jet velocity and the jet radius as a function of jet distance from previous observations; Γ(z)∝z^0.16 (Park et al. 2019a), R(z)∝z^0.56 (Asada & Nakamura 2012; Hada et al. 2013, 2016; Nakamura et al. 2018). The magnetic field at this distance is assumed to be dominated by the toroidal component, in which case the magnetic field strength is inversely proportional to the radius of the jet as the magnetic flux is conserved (Marscher et al. 2010). Then, the magnetic field strength in the co-moving frame is (Lyutikov et al. 2005; Zamaninasab et al. 2014)

$$\begin{array}{c}\hfill B(z)\propto \frac{1}{R\mathrm{\Gamma }}=B_i{\left(\frac{z}{z_i}\right)}^{-0.72},\end{array}$$(6)

where B_i is the strength of the magnetic field at z_i.

The energy spectrum of the nonthermal electron injection function, Q(γ), is generally assumed to be a power law with a slope p_inj limited to the minimum and maximum energies, γ_inj, min and γ_inj, max. In addition to this, we assume that the number of nonthermal electron injections varies with distance with a power of −q:

$$\begin{array}{c}\hfill Q(\gamma ,z)=Q_0{\gamma }^{p_{\mathrm{inj}}}{\left(\frac{z}{z_i}\right)}^{-q},\\ \hfill {\gamma }_{\mathrm{inj},\min }<\gamma <{\gamma }_{\mathrm{inj},\max },\end{array}$$(7)

where Q₀ is the nonthermal electron injection rate at z_i for electrons of γ = 1. For example, if q = 0, the same amount of nonthermal electrons as Q(z_i) is injected along the jet continuously. For q = ∞, there is no additional injection after the initial injection at z_i. The above two cases resemble the two scenarios of the spectral steepening suggested by Hovatta et al. (2014). If 0 < q < ∞, nonthermal electrons are injected continuously, but the injection volume decreases with distance. There are indeed theoretical models that suggest that electron acceleration along the jet is possible; one is the boundary shear acceleration between the jet and the surrounding matter (Ostrowski 1998), while the other is Kelvin-Helmholtz instability driven magnetic reconnection (e.g., Sironi et al. 2021). In Sect. 4.3, we discuss the effect of the electron injections on the spectral index distribution.

Given B_i and q, Eqs. (4) and (5) are numerically solved and N(γ, τ) is obtained, and the spectral index distribution α(z) between the two observed frequencies is obtained from N(γ, τ). Throughout the paper, the initial and final distance z_i and z_f are set as ∼900 r_s and ∼4500 r_s, which is the deprojected distance of 2 mas and 10 mas of the M 87 jet, respectively. The minimum and maximum energy of the injection function, γ_inj, min and γ_inj, max, is 10 and 10⁵, and the energy slope of the injection function is p_inj = −2.4 (α_inj = −0.7), which is similar to the average value of the observed spectral index at 2 mas (Fig. 2). As Q₀ does not affect the spectral index we arbitrarily selected Q₀ as 1.00 × 10²⁰. Details of the model are described in Appendix D.

4.3. The effect of the nonthermal electron injections on the spectral index distribution

If we assume that the flattening of the spectral index seen in Fig. 2 is real, we can test the effect of the nonthermal electron injections on the spectral index distribution. If this assumption is wrong, we discuss the implications of this in Sect. 4.4. We made a number of spectral index distribution models by adopting different values of q. Here, we set the four different cases of the nonthermal electron injection function in the jet to be q = 0, 5, 10, and ∞, while the magnetic field distribution is assumed to be the same B_i with a value of 0.3 G.

Figure 3 summarizes the results of the test models. The left panel shows the evolution of electron energy distribution N(γ) as a function of the jet distance drawn as a vertical dashed line in the right panel. The color of the distribution represents the value of q used in the model. A blue and a red vertical line represent the electron energies emitting synchrotron radiation at 22 and 43 GHz in the observer’s frame. The slope of the nonthermal electron injection is p_inj = −2.4. The right panel of Fig. 3 shows the spectral index distribution between 22−43 GHz versus jet distance, α_{22 − 43 GHz}(z), obtained from N(γ, z). The horizontal dashed line represents the spectral index of the injection function (α_inj = −0.7).

Fig. 3. Results of the modeled electron energy distributions (N(γ)) and the spectral index distributions as a function of distance along the jet. This is done to test the effect of the nonthermal electron injection rate in the jet. Left: Modeled electron number density distributions with different values of q = 0, 5, 10, and ∞. The different colors represent different values of q in the model: q = 0 represents continuous new electron injections in the jet, q = ∞ no new electron injections in the jet, and q = 5 and 10 new electrons injected into the jet but the volume decreases with distance. In all cases, we set the same magnetic field distribution (Bi = 0.3 G). The slope of the nonthermal electron injection is pinj = −2.4. The vertical red and blue lines represent the energies of the nonthermal electrons, which emit at 22 and 43 GHz. The location down the jet where the calculations were made is written in the top right corner. Right: Modeled spectral index distributions as a function of distance, which correspond to the slope of the electron number density distributions between the vertical lines in the left panel. The vertical dashed black line is the location down the jet where the calculations were made. The horizontal dashed line is the slope of the injection function (αinj = −0.7). An animation of this figure is available.

In the absence of electron injection (q = ∞), cutoff energy is shown in N(γ, z) where no higher energy electrons can be found. The cutoff energy moves to lower energy over distance. At z ∼ 2000 r_s (middle left panel), the cutoff energy is located in between γ(ν_obs = 43 GHz) and γ(ν_obs = 22 GHz). Therefore, α_{22 − 43 GHz}(z) diverges to −∞ around this distance. On the other hand, in the constant injection case (q = 0), the electrons still exist beyond the cutoff energy due to the newly injected one. Therefore, the spectral index changes only slightly with distance⁴.

In the case of q = 5 and 10, the amount of newly injected electrons decreases with distance, thus the slope of N(γ, z) at energies higher than the cutoff energy is steeper than in the case of q = 0. As a result, the energy distribution becomes a broken power law, where the position of the break (γ_break) is the same as the cutoff energy in the q = ∞ case. Interestingly, α_{22 − 43 GHz}(z) for q = 5 and 10 show similar trends to the observed distributions, that is, α decreases until ∼5 mas and no longer decreasing or slowly decreases after that. This is because γ_break first passes through γ(ν_obs = 43  GHz) and then γ(ν_obs = 22  GHz), and accordingly, the slope between these two frequencies shifts from the first slope of the broken power-law to the second. The larger q, the steeper the second slope of N(γ), so the change in α will be greater. For example, if q = 10 then |Δα| ∼ 2, which is far greater than the expectation by Hovatta et al. (2014; |Δα| = 0.5), and similar to our observations.

From this experiment, we found that our model can reproduce the characteristic trend of the spectral index distribution seen in the KaVA observations for cases where the amount of nonthermal electron injections have a distance dependence. The combination of nonthermal electron injections and energy losses creates a broken power-law energy distribution in which the slope beyond γ_break is determined by the change in the amount of injection with distance (i.e., the parameter q of the electron injection function). The larger the q, the greater the slope beyond γ_break and the change in α. The continuous injection and non-injection cases proposed by Hovatta et al. (2014) correspond to two extreme cases in our model (q = 0 and ∞, respectively). In the former case, the change in α is the smallest, and in the latter, the break energy becomes the cutoff energy with no electrons above it, making the spectral index diverge.

4.4. Constraining the magnetic field strength and nonthermal electron injection along the M 87 jet

Now we constrain the physical quantities of the M 87 jet using our numerical model. In order to compare the modeled spectral index distributions with the observations quantitatively, we define an “allowed region” based on the observational spectral index distributions with three parameters: (1) The spectral index range at z_i (α_i), (2) the location of the end of the steepening region (z_s, f), and (3) the spectral index range after the steepening region (α_f). The gray shaded region in Fig. 4 is the allowed region of the spectral index of the M 87 jet drawn with parameters of −1.2 < α_i < −0.3, z_s, f = 5.5 mas, −3.0 < α_f < −1.5. More than 90% of the observed profiles lie within the allowed region.

Fig. 4. Constraining the magnetic field strength and nonthermal electron injection of the M 87 jet. Left: Modeled spectral index distributions constrained by the observations. The distributions must lie in the gray shaded allowed region (see Sect. 4.4 for details). Colored lines are the modeled spectral index distributions. Different colors represent different initial magnetic field strengths, Bi. The dashed and dotted red line is a modeled spectrum index distribution with q = 10 and q = 11, respectively, with Bi = 0.3 G. Right: Parameter space (Bi, q). The black region is the allowed values of Bi and q constrained by the observations.

We then create many modeled spectral index distributions for different values of the parameters B_i and q. The parameters B_i were set from 0.1 G to 1.2 G in 0.1 G increments, and the parameter q was set from 5 to 13. We limited the parameters by excluding models that exist outside the allowed region. The colored lines in the left panel of Fig. 4 are the modeled spectral index distributions that lie completely inside the allowed region. The color indicates the initial magnetic field strength used in the model. In a strong magnetic field, the steepening region is located closer to the core, and the spectral index decreases less. As the magnetic field strength becomes weaker, the spectral index decreases over a longer distance, and the change in the spectral index are greater. For the same magnetic field strength, q essentially determines the final spectral index as discussed in Sect. 4.3. For example, for the same B_i = 0.3 G, the model with q = 11 (red dotted line) have steeper spectral index distribution than the q = 10 case (red dashed line). By constraining the modeled distributions to the allowed region, B_i is constrained to 0.3 G < B_i < 1.0 G and q is constrained to 7 < q < 11. The right panel of Fig. 4 shows the allowed parameter space of B_i and q.

We note that our method is highly sensitive to the choice of the allowed region. In Sect. 3.2 we found that a constant spectral index distribution from the VLBA data is less clear, so it is possible that the spectral index distribution could continue to decay farther down the jet (blue lines in Fig. 2). This would suggest two possibilities: (1) Spectral index flattening still occurs but farther downstream. In this case, the model indicates a weaker B_i. (2) The spectral index declines indefinitely. In this case, this would suggest that no electrons are injected beyond ∼6 mas. Considering that there are several years between the VLBA and KaVA observations, there is a possibility that the magnetic field strength and the electron injection rate could change with time. Since 2017, we have been monitoring the M 87 jet with the East Asia VLBI Network (EAVN; Cui et al. 2021) at 22 and 43 GHz quasi-simultaneously with better sensitivity and (u, v)-coverage than KaVA. These observations could help resolve these possibilities.

4.5. Radial magnetic field distribution of the M 87 jet

In Sect. 4.4 we successfully inferred the magnetic field strength from the modeled spectral index distributions in the M 87 jet from 2−10 mas (∼900 r_s − ∼4500 r_s in the deprojected distance). Therefore, the strength of the magnetic field evolves a function of distance as

$$\begin{array}{c}\hfill B(z)=(0.3-1.0\mathrm{G}){\left(\frac{z}{2\mathrm{mas}}\right)}^{-0.72},\\ \hfill 2\mathrm{mas}<z<10\mathrm{mas},\end{array}$$(8)

where exponent −0.72 comes from Eq. (5).

Figure 5 shows the magnetic field strength of the M 87 jet as a function of the distance from the SMBH. The magnetic field distribution obtained by our models and the observations is drawn as a black area. The figure also summarizes the magnetic field strengths at various locations of the M 87 jet using observations over a wide range of wavelengths. The cross represents the theoretically expected magnetic field strength at the black hole surface based on Blandford-Znajek power estimation (Blandford et al. 2019)⁵. A recent estimation by Kino et al. (2022) also indicates a similar magnetic field strength in this region. The circles represent magnetic field strengths in radio cores at various VLBI frequencies, obtained based on SSA theory (Kino et al. 2014, 2015; Hada et al. 2012, 2016; Reynolds et al. 1996; EHT MWL Science Working Group 2021)⁶⁷. Here, the distance from the core to the SMBH and its error are determined using the result of core-shift (Hada et al. 2011). The error of the magnetic field strength is calculated by applying 25% error to the radio core size and 5% error to the flux of the core (Hada et al. 2013). For the data shown as a ranged value (Kino et al. 2014, 2015), the error bars of the field strength is not displayed. The diamonds indicate field strengths obtained with VLBI observations using the core-shift method (Zamaninasab et al. 2014; Jiang et al. 2021)⁸, light curve modeling (Acciari et al. 2009), brightness temperature (Kim et al. 2018b), and EHT observations (Event Horizon Telescope Collaboration 2019, 2021). The field strengths at HST-1 estimated by applying the synchrotron cooling time to the variability of the X-ray band are drawn as squares (Harris et al. 2003, 2009).

Fig. 5. Magnetic field distribution of the M 87 jet. The black area is the distribution of the magnetic field strength estimated in this study. The solid lines with the gray shaded region are the extrapolation of the upper and lower limits of the result in the upstream direction (B(z)∝z−0.72). The dashed lines are the magnetic field distribution which assumes a steeper slope of the jet radius profile in the region closer to the SMBH (B(z)∝z−0.92). The solid lines without a shaded region are the upper and lower bounds of the magnetic field distribution, assuming fast acceleration in the region closer within 2 mas (B(z)∝z−1.12). The dotted line is the extrapolation of the maximum values of the result to the downstream region. Markers are estimated magnetic field strength from previous studies. The blue cross is the theoretically expected magnetic field strength in a Blandford-Znajek jet (Blandford et al. 2019). The circles are magnetic field strength values obtained in VLBI cores derived using SSA theory (Kino et al. 2014, 2015; Hada et al. 2012, 2016; Reynolds et al. 1996; EHT MWL Science Working Group 2021). The diamonds are values obtained from VLBI observations using various methods other than SSA theory (Acciari et al. 2009; Kim et al. 2018b; Zamaninasab et al. 2014; Event Horizon Telescope Collaboration 2019, 2021; Jiang et al. 2021). The squares are the magnetic field strengths at HST-1 (Harris et al. 2003, 2009). The horizontal dotted-dashed lines are the magnetic field strengths estimated by fitting the broadband SED of the jet using a single-zone SSC emission model (Abdo et al. 2009; MAGIC Collaboration 2020; EHT MWL Science Working Group 2021).

The solid lines with gray shaded region are the extrapolation of the upper and lower bounds of the magnetic field distributions estimated in this study (B(z)∝z^−0.72; Eq. (8)). The dashed lines represent the distribution that reflects a steeper slope of the jet radius profile in the region close to the SMBH as R(z)∝z^−0.76 (accordingly, B(z)∝z^−0.92; Hada et al. 2013). Interestingly, our estimation of B(z) on 900 − 4500 r_s scale smoothly connects to values based on synchrotron emission down to the EHT observing scale, keeping B(z)∝z^−0.72. The consistency of the slope of the magnetic field distributions with distance may imply that the majority of the magnetic flux of the jet near the SMBH is preserved without serious dissipation out to ∼4500 r_s.

Mertens et al. (2016) reported that the M 87 jet shows linear acceleration in the region closer than 2 mas, meaning that the bulk jet velocity profile in this region is proportional to the jet radius (Γ(z)∝R(z)∝z^0.56). The solid lines without shaded regions are the upper and the lower bounds of the magnetic field distribution estimated by applying it. In this case, due to the fast acceleration, the magnetic field strength at 2 mas should increase to 0.8 − 2.4 G in order to achieve a similar amount of spectral index steepening. Also, the slope of the magnetic field distribution becomes steeper (B(z)∝z^−1.12) in the region closer than 2 mas. As a result, a very strong magnetic field distribution is inferred that deviates from the magnetic field strengths estimated by observations and theory.

Extending our scaling relation downstream (dotted line) appears to be considerably higher than the magnetic field strength of HST-1. HST-1 is a complex structure in which re-collimation occurs, and the magnetic field strength of the upstream jet before HST-1 is expected to be smaller than that of HST-1. Therefore, the magnetic field strength measured in HST-1 should be considered as an upper limit. Given that our extrapolated maximum value to the location of HST-1 is roughly an order of magnitude higher than the upper limits set by X-ray variability studies in HST-1, this suggests that assumptions such as the conservation of the magnetic flux along the jet may be wrong. For example, a recent large-scale jet simulation showed that pinch instability can lead to magnetic flux dissipation along the jet (Chatterjee et al. 2019).

The horizontal dotted-dashed lines are the magnetic field strengths estimated by fitting a single-zone synchrotron self-Compton (SSC) emission model to the broadband spectral energy distribution (SED) of the jet, specifically focusing on the very high energy (VHE) γ-rays (Abdo et al. 2009; MAGIC Collaboration 2020; EHT MWL Science Working Group 2021). Comparing these values directly with the magnetic field distributions in our study, the estimated location of the VHE emission is in the downstream jet between z ∼ 4500 r_s and ∼10⁵ r_s, which differs significantly from the most favored location of VHE flare activities (a few tens of r_s from the SMBH; see, e.g., Acciari et al. 2009; Abramowski et al. 2012; Hada et al. 2014; Akiyama et al. 2015). One possible explanation is that VHE γ-rays are emitted from the outer layer of the radio-emitting region (e.g., Tavecchio & Ghisellini 2008). Indeed, a recent simultaneous broadband SED study found that the VHE emitting region is not overlapping with the radio-emitting region, and requires a larger area (EHT MWL Science Working Group 2021). However, there are various physical processes that suggest the possibility of the VHE emission not only originating in the innermost jet base but also in the downstream jet (e.g., the HST-1 knot; see Abramowski et al. 2012; Rieger & Levinson 2018, and references therein). The origin and location of the VHE emission in the M 87 jet are not yet conclusive.

5. Summary and conclusions

In this work we have investigated the spectral index distribution of the M 87 jet using KaVA LP observations from 2016 and archival VLBA observations from 2010 and 2014. This is the first systematic study of the spectral properties in the innermost regions of the M 87 jet at 22 and 43 GHz, and it allows us to investigate the spectral index out to ∼10 mas from the jet base. Our study can be summarized as follows:

1. The overall spectral index morphology of the M 87 jet does not appear to have changed significantly over the course of these observations.

2. The observed spectral index distributions of the M 87 jet from the KaVA observations show a rapid steepening from ∼2 mas until ∼6 mas from α ∼ −0.7 to −2.5. Beyond ∼6 mas, the spectral steepening stops and remains at the same value until ∼10 mas. The trend of the spectral index distribution obtained from the VLBA data showed similar rapid spectral index decay, but the constant spectral index beyond ∼6 mas is less clear. This suggests that the spectral index distribution of the M 87 jet may change with time.

3. To understand the details of the spectral evolution in the M 87 jet, we modeled the synchrotron spectral index by solving the transfer equation and calculating the evolution of the energy distribution of the nonthermal electron densities in the jet. We find that if the amount of the nonthermal electron injections decreases with distance, the spectral index distribution reproduces the characteristics of the observed distributions.

4. We compare the model of the spectral index distributions with the allowed region, which is defined based on the observations. As a result, we find that the initial magnetic field strength at 2 mas of 0.3 G < B_i < 1.0 G and electron injection function with 7 < q < 11 reproduce the observed spectral index distribution appropriately. We conclude that the distribution of the magnetic field strengths as inferred from the model at 2−10 mas (∼900 r_s − ∼ 4500 r_s in the deprojected distance) in the M 87 jet is B = (0.3−1.0 G)(z/2 mas)^−0.72.

5. We compared the magnetic field strengths of the M 87 jet with those estimated in previous observational studies. The extrapolation of the magnetic field distribution is in good agreement with the estimates by VLBI in the inner jet region, indicating that the majority of the magnetic flux of the jet near the SMBH is preserved without serious dissipation up to ≲4500r_s. However, when extrapolating to the HST-1 region, our maximum magnetic field strength is roughly an order of magnitude higher than the upper limits found from X-ray variability studies.

In this work, the magnetic field strength of the M 87 jet at a distance of several thousands r_s has been estimated for the first time. This study showed that, if the bulk jet velocity profile and radius profile are known, the magnetic field strength and nonthermal electron injection of the extended jet can be inferred from the spectral index profile. Consequently, it provides the magnetic field profile of the jet, which is an important physical quantity to validate the launching model of the jet.

Movie

Movie 1 associated with Fig. 3 Access here

Acknowledgments

We thank the anonymous referee for constructive comments to improve the quality of this paper. This work is based on observations made with the KaVA, which is operated by the Korea Astronomy and Space Science Institute and the National Astronomical Observatory of Japan. This research was supported by the Korea Astronomy and Space Science Institute under the R&D program supervised by the Ministry of Science and ICT. H.R., B.W.S., and A.C. acknowledge support from the KASI-Yonsei DRC program of the Korea Research Council of Fundamental Science and Technology (DRC-12-2-KASI). B.W.S. is grateful for the support by the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (MSIT) of Korea (NRF-2020K1A3A1A78114060). A.C. acknowledges support from the National Research Foundation of Korea (NRF), grant Nos. 2022R1A2C100298211, and 2022R1A6A1A03053472. J.P. and I.C. acknowledge financial support from the Korean National Research Foundation (NRF) via Global Ph.D. Fellowship Grant 2014H1A2A1018695 and 2015H1A2A1033752, respectively. J.P. is supported by an EACOA Fellowship awarded by the East Asia Core Observatories Association, which consists of the Academia Sinica Institute of Astronomy and Astrophysics, the National Astronomical Observatory of Japan, the Center for Astronomical Mega-Science, the Chinese Academy of Sciences, and the Korea Astronomy and Space Science Institute. I.C. acknowledges financial support by the Spanish Ministerio de Economía y Competitividad (grants PID2019-108995GB-C21). K.H. acknowledges support from the Mitsubishi Foundation grant No. 201911019. This work is partially supported by JSPS/MEXT KAKENHI (grants 18K03656, 18H03721, 19H01943, 18KK0090, 21H01137, 21H04488). S.S.-S. is supported by JSPS KAKENHI grant No. 21K03628. K.L. and S.T. acknowledge financial support from the National Research Foundation (NRF) of Korea through grant 2022R1F1A1075115. J.O. was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2021R1A6A3A01086420). J.-S.K. has been supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2016R1A5A1013277 and 2020R1A2C1007219). J.-Y.K. acknowledges support from the National Research Foundation (NRF) of Korea (grant no. 2022R1C1C1005255). S.-S.L. was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MIST), (2020R1A2C2009003). Y.C. is funded by the China Postdoctoral Science Foundation (2022M712084). R.-S.L. is supported by the Key Program of the National Natural Science Foundation of China (grant No. 11933007), the Key Research Program of Frontier Sciences, CAS (grant No. ZDBS-LY-SLH011), the Shanghai Pilot Program for Basic Research – Chinese Academy of Sciences, Shanghai Branch (JCYJ-SHFY-2022-013), and the Max Planck Partner Group of the MPG and the CAS. Z.S. acknowledges support from the Major Program of the National Natural Science Foundation of China (Grant No. 11590780, 11590784), and Key Research Program of Frontier Sciences, CAS (grant No. QYZDJ-SSWSLH057).

References

Appendix A: Effect of the (u, v)-range matching on the spectral index distribution

As described in Sect. 3.1, it is important to match the (u, v) ranges of the two frequencies when producing spectral index maps. Considering the typical (u, v) ranges for KaVA and VLBA at 22 and 43 GHz given in Table A.1, we used data within 33−170 Mλ for KaVA and 25−685 Mλ for VLBA. This cuts out the data on average ∼ 21 % (∼ 46 %) for KaVA 22 GHz (43 GHz), and ∼ 5 % (∼ 8 %) for VLBA 22 GHz (43 GHz), respectively. As a result, the image noise increases on average by ∼ 7 % (∼ 33 %) for KaVA 22 GHz (43 GHz), and ∼ 0 % (3 %) for VLBA 22 and 43 GHz, respectively.

The upper panel of Fig. A.1 compares the spectral index distributions before (red lines) and after (blue lines) matching the (u, v) ranges. This shows the distributions using the full (u, v) range have steeper spectra at longer distances. The difference between these two distributions – the artificial steepening along the jet distance (Δα_arti) – is shown in the lower panel. The gray lines are from individual epochs and the thick black line is their average. We also calculated the expected Δα_arti under the assumption that the loss of sensitivity for large-scale structures at the shortest baseline is ∼sin(πθ/θ_max)/(πθ/θ_max), (e.g., Thompson et al. 2017), where θ is the size of the structure and θ_max is the angular size corresponding to the shortest baseline. Here, we assume θ as jet widths and adopt it from Hada et al. (2013), and θ_max of KaVA and VLBA at 22 and 43 GHz are summarized in Table A.1. After calculating the expected Δα_arti for KaVA and VLBA, they were weighted and averaged by the number of observations to calculate the final value. This value is displayed in a red line in the lower panel.

Fig. A.1. Effect of the (u, v)-range matching on the spectral index distribution. The upper panel shows spectral index distributions between 22 and 43 GHz before (red lines) and after (blue lines) matching the (u, v) range. Distributions using the full (u, v) range have steeper spectra at longer distances. Differences between the two distributions (i.e., artificial steeping due to inconsistent (u, v) ranges) are indicated by gray lines in the lower panel. The thick black line is the average difference, and the red line is the expected artificial steepening along the jet. See text for details.

The averaged Δα_arti increases with distance, close to the expected value, and is ∼ − 0.2 at ∼ 10 mas. Individual cases show Δα_arti ∼ −0.6 at a distance of ∼ 7 mas in severe cases. Even taking this into account, the observed spectral index decreases more significantly (Δα ∼ −2). Therefore, we conclude that the decrease in the spectral index with the distance of the M 87 jet is intrinsic.

Appendix B: Effect of beam size on the spectral index distribution

Typical beam sizes of VLBA 22 and 43 GHz are 0.8 × 0.4 mas and 0.45 × 0.22 mas, respectively, which are much smaller than the typical beam sizes of KaVA. In Fig. B.1 the spectral index distribution between the VLBA 22 and 43 GHz maps was convolved into two different beams: 0.8 × 0.8 mas and 1.2 × 1.2 mas. The radius of each beam is equal to the major axis of the VLBA and KaVA beams at 22 GHz, respectively. Compared to the spectral index distribution using the KaVA beam, the distribution using the VLBA beam has a greater scatter with some local fluctuations. As the beam size increases, the local structures are averaged and the overall distribution is smoother. However, it can be seen that the global trend of the spectrum does not change with the size of the beam.

Fig. B.1. Spectral index distributions between VLBA 22 and 43 GHz maps convolved with beams of different sizes. The blue lines are convolved with a circular beam of 1.2 × 1.2 mas, where the radius is the same as the major axis of the KaVA 22 GHz beam. The red lines are convolved with 0.8 × 0.8 mas, which is equal to the major axis of the VLBA 22 GHz beam. The distribution with the VLBA beam has a larger scattering than the distribution with the KaVA beam.

Appendix C: Alternative explanation of the spectrum steepening

In our model, we assumed a relationship of α = (p + 1)/2 between the slope of the electron distribution and the synchrotron spectrum (Rybicki & Lightman 1979). This assumption is maintained when the electron energy corresponding to the observation frequency is within a sufficiently wide range of γ in the optically thin region. The observation frequencies 22 and 43 GHz correspond to γ ∼ 10^{1.6 − 2.0} and are located within the electron energy range of the injection function (γ_inj, min = 10¹, γ_inj, max = 10⁵). However, without electron injection, γ_max decreases over time and thus the above assumption may not hole. In this case, the synchrotron spectrum has a different shape depending on the synchrotron emission model we use.

The most widely used synchrotron emission models were proposed in Kardashev (1962) and Pacholczyk (1970; KP model) and Jaffe & Perola (1973; JP model). In the JP model, which is used in our spectral index model, the scattering due to turbulence in the magnetic field randomizes the electron pitch angles on shorter timescales relative to the radiative lifetime (e.g., Hardcastle & Krause 2013). Consequently, the cutoff energy is determined only by magnetic field strength, not by the pitch angle.

In the KP model, however, the pitch angle of individual electrons remains constant over the lifetime, and electrons with a pitch angle parallel to the magnetic field do not lose energy. Hence, the cutoff energy depends not only on the magnetic field strength but also on the pitch angle. Assuming an isotropic pitch angle, the synchrotron spectrum does not drop exponentially but has a break with the spectral index ranging from α_inj at low frequencies to 4/3α_inj − 1 at high frequencies (Kardashev 1962)⁹.

Figure C.1 shows a schematic picture of the spectral index distribution using the KP and JP models in the absence of electron injection. Applying an initial slope of α_inj = −0.7, the KP type energy loss predicts the slope at energies higher than the break energy as 4/3α_inj − 1 ≈ −1.933, which is within the allowed region defined from observations. Therefore, this scenario may also explain the spectral index steepening.

Fig. C.1. Schematic diagram of the spectral index distributions without electron injection, assuming different synchrotron emission models. The blue line assumes the JP model, which is used in this study. The red line shows the spectral index distribution with the KP model. Assuming that the injection spectral index is αinj = −0.7, the expected index after steepening is 4/3αinj − 1 ≈ −1.933, which is within the allowable range obtained from observation.

In summary, we suggest two scenarios to explain the observed spectral index distribution. (1) JP-type energy loss and particle injection with distance dependence. (2) KP-type energy loss and no electron injection. Both scenarios make an electron energy distribution a shape of a broken power law, and hence it is possible to explain the observed trend of the spectral index distributions. Further study is needed to discuss which scenario is more plausible on the M 87 jet. One of the expectations is that the brightness distribution between the JP model scenario and the KP model scenario is different since the JP model scenario requires continuous electron injection and the KP model scenario does not. Another distinguishing point between the JP and the KP scenarios is the fractional polarization. The JP model needs a fully random field plus an isotropic pitch angle. This will lead to very low fractional polarization. The KP model predicts linear polarization (depending on the magnetic field geometry). Therefore, fractional polarization could be a good indicator of the synchrotron loss model.

Appendix D: Details of the spectral index model

The general transfer equation of the electron number density distribution is written as follows (e.g., Ginzburg & Syrovatskii 1964; Longair 2011; Blasi 2013):

$$\begin{array}{c}\hfill \frac{\partial N}{\partial \tau }-\mathrm{\nabla }·\left(D\mathrm{\nabla }N\right)+\mathrm{\nabla }·\left(vN\right)+\frac{\partial }{\partial \gamma }\left(bN\right)=Q,\end{array}$$(D.1)

where D is the diffusion coefficient, v is the velocity of the system, b is the energy loss rate, and Q is the nonthermal electron injection rate. In our discussion, the frame of reference is a cross section of the jet co-moving with the jet’s flow, then v ⋅ ∇N = 0. We further neglect the diffusion of the electron. Then, Eq. (D.1) becomes Eq. (4).

To solve the general transfer equation numerically, we change the partial derivative equation into a series of ordinary differential equations (e.g., Arfken & Weber 2005). First, we differentiate b(γ), (Eq. 5) with respect to γ,

$$\begin{array}{c}\hfill \frac{\partial b}{\partial \gamma }=-b_{\mathrm{adi}}-2b_{\mathrm{sync}}\gamma .\end{array}$$(D.2)

Putting Eq. (D.2) into Eq. (D.1) and using the relation $\mathrm{\nabla }·v=\frac{3}{R}\frac{\mathit{dR}}{d\tau }=3b_{\mathrm{adi}}$, we get

$$\begin{array}{c}\hfill \frac{\partial N}{\partial \tau }+\frac{d\gamma }{d\tau }\frac{\partial N}{\partial \gamma }=Q+2(b_{\mathrm{sync}}\gamma -b_{\mathrm{adi}})N.\end{array}$$(D.3)

The left-hand side of Eq. (D.3) is the same as a total derivative of N(γ, τ) with respect to τ. Consequently, the transfer equation changes into a series of ordinary differential equations:

$$\begin{array}{cc}& \frac{\mathit{dN}}{d\tau }=2(b_{\mathrm{sync}}\gamma -b_{\mathrm{adi}})N+Q\hfill \end{array}$$(C.4.1)

$$\begin{array}{cc}& \frac{d\gamma }{d\tau }=-b_{\mathrm{adi}}\gamma -b_{\mathrm{sync}}{\gamma }^2,\hfill \end{array}$$(C.4.2)

which are solvable numerically. The number density of electrons N(γ) changes with time with respect to Eq. (C.4.1), and at the same time, each electron in N(γ) loses its energy by following Eq. (C.4.2), (same as Eq. 5).

To solve the problem, we additionally need to know the relationship between the distance of the reference frame from the SMBH and the time of the frame (i.e., the proper time). This is because the physical quantities in the equation were obtained as a function of distance, while the transfer equation deals with the change with time. When τ is the proper time and z is the distance, τ and z are in the following relation:

$$\begin{array}{c}\hfill \frac{d\tau }{\mathit{dz}}=\frac{d{t}_{\mathrm{obs}}}{\mathit{dz}}\frac{d\tau }{d{t}_{\mathrm{obs}}}=\frac{1}{\mathrm{\Gamma }\beta c},\end{array}$$(C.5)

where t_obs is the time in the observer’s frame, $\beta =\sqrt{1-{\mathrm{\Gamma }}^{-2}}$ is the bulk jet speed in units of the speed of the light and Γ is the bulk Lorentz factor. Here, we use the equation of apparent jet speed ${\beta }_{\mathrm{app}}c\equiv \frac{\mathit{dz}}{d{t}_{\mathrm{obs}}}sin\theta =\frac{\beta csin\theta }{1-\beta cos\theta }$ and Doppler factor $\delta \equiv \frac{d\tau }{d{t}_{\mathrm{obs}}}=\frac{1}{\mathrm{\Gamma }(1-\beta cos\theta )}$ (e.g., Ghisellini 2013). The relation between τ and z is then calculated by integrating Eq. (C.5):

$$\begin{array}{c}\hfill {\int }_{{\tau }_i}^{{\tau }_f}d\tau ={\int }_{z_i}^{z_f}\frac{\mathit{dz}}{c\sqrt{{\mathrm{\Gamma }}^2-1}}.\end{array}$$(C.6)

The distance range of our calculation is from z_i = 2 mas to z_f = 10 mas, which correspond to z_i ≈ 889 r_s and z_f ≈ 4446 r_s in deprojected units assuming a viewing angle of 17° (Walker et al. 2018). Applying the observed bulk jet velocity field of the M 87 jet (Park et al. 2019a) and setting τ_i = 0s, the travel time is τ_f ∼ 2.2 × 10⁸s (∼ 6.9 years). If we apply the velocity field with faster acceleration in the region closer than 2 mas (Mertens et al. 2016), the travel time between 2 − 10 mas is ∼ 2.4 years.

Now we solve the Eqs. (C.4.1) and (C.4.2) to obtain the electron number density distribution as a function of time, N(γ, τ) ¹⁰, which is further converted to N(γ, z) using Eq. (C.6). From N(γ, z), the slope of the number density distribution between γ(ν_obs = 22 GHz) and γ(ν_obs = 43 GHz) along the jet distance, p_{22 − 43GHz}(z), is obtained by using the relation of the observed frequency of synchrotron radiation and the energy of nonthermal electrons (Rybicki & Lightman 1979):

$$\begin{array}{c}\hfill {\nu }_{\mathrm{obs}}=\frac{3eB}{4\pi m_{\mathrm{e}}c}{\gamma }^2\delta .\end{array}$$(C.7)

Finally, the spectral index distribution α_{22 − 43GHz}(z) is obtained by using the relation α = (p + 1)/2 (Rybicki & Lightman 1979).

Appendix E: Notes on the injection function Q(γ, z)

The electron injection function in our calculation is assumed to be a power-law distribution with a sharp cutoff at γ_inj, min = 10 and γ_inj, max = 10⁵. The minimum and maximum range in the electron energy distributions are constrained by the observations. The M 87 jet has been detected in VLBI observations at frequencies as low as 1.6 GHz (e.g., Cheung et al. 2007), and as low as ∼20 MHz using Low-Frequency Array (LOFAR; e.g., de Gasperin et al. 2012). From Eq. (C.7), the electron energy corresponding to synchrotron radiation at 20 MHz is γ(ν_obs = 20MHz)∼10^{0.1 − 0.4} when B_i = 0.3 − 1.0G and Γ = 2 is assumed. Likewise, the detection of X-ray emission in the M 87 jet between the core and HST-1 could constrain the maximum energy of the electron distribution (e.g., Perlman & Wilson 2005; EHT MWL Science Working Group 2021). From Eq. (C.7), the X-ray synchrotron emission requires γ(ν_obs = 300PHz)∼10^{5.2 − 5.5}. Perlman & Wilson (2005) suggested that even higher energies of nonthermal electrons of γ ∼ 10^6 − 8 is required in order to emit X-ray synchrotron emission. We also found that electrons with energies greater than 10⁵ have little effect on the spectral index distribution at 22−43 GHz. Figure E.1 shows the spectral index distribution model using different γ_inj, max in the injection function. For all models, q = 10 and B_i = 0.3G were applied. As γ_inj, max increases, the spectral index distributions gradually converge, and when γ_inj, max is greater than 10^4.5, the spectral index distribution does not change. This is because the synchrotron cooling time of the electrons with energy γ ≳ 10⁵ is ${\tau }_{\mathrm{sync}}=\gamma {\left(\frac{d\gamma }{d\tau }\right)}^{-1}\lesssim 8.6\times {10}^4\mathrm{s}$, which is much shorter than the travel time (∼2.2 × 10⁸s). From this experiment, we found that, for γ_inj, max > 10⁵, the spectral index distribution remains unchanged regardless of the shape of the high-energy tail of the injection function (e.g., a sharp cutoff or an exponentially decreasing). Therefore, in this study, we assume the injection function as a power-law energy distribution and the sharp cutoff at γ_inj, max = 10⁵.

Fig. E.1. Comparison of model spectral index distributions using nonthermal electron injection functions with different γinj, max. Distributions with different colors represent different values of γinj, max. For all distributions, the magnetic field strength was assumed to be Bi = 0.3G. As γinj, max increases, the spectral index distributions converge and there is little change between the models with γinj, max = 104.5 and γinj, max = 106.

All Tables

All Figures

	Fig. 1. Spectral index maps between 22 and 43 GHz obtained from VLBA and KaVA observations. All images have been rotated by −18°. The restoring beam size is 1.2 mas × 1.2 mas, drawn as a black circle in the bottom-left corner. Observing dates are shown to the left of each map. The contours represent the total intensity at 22 GHz. Contours start at 3σrms, increasing in steps of 2.
In the text

Fig. 2. Evolution of the spectral index between 22 and 43 GHz along the deprojected distance from the SMBH. Weighted mean values across the jet are used in this distribution (see Sect. 3.2). The projected distance is displayed on the upper axis in mas, and the corresponding deprojected distance is displayed on the lower axis in rs. Red lines represent the spectral index distribution from 2016 KaVA observations and blue lines the spectral index distribution from VLBA observations from 2010 and 2014. The shaded area is the 1σ error of the spectral index.

In the text

Fig. 3. Results of the modeled electron energy distributions (N(γ)) and the spectral index distributions as a function of distance along the jet. This is done to test the effect of the nonthermal electron injection rate in the jet. Left: Modeled electron number density distributions with different values of q = 0, 5, 10, and ∞. The different colors represent different values of q in the model: q = 0 represents continuous new electron injections in the jet, q = ∞ no new electron injections in the jet, and q = 5 and 10 new electrons injected into the jet but the volume decreases with distance. In all cases, we set the same magnetic field distribution (Bi = 0.3 G). The slope of the nonthermal electron injection is pinj = −2.4. The vertical red and blue lines represent the energies of the nonthermal electrons, which emit at 22 and 43 GHz. The location down the jet where the calculations were made is written in the top right corner. Right: Modeled spectral index distributions as a function of distance, which correspond to the slope of the electron number density distributions between the vertical lines in the left panel. The vertical dashed black line is the location down the jet where the calculations were made. The horizontal dashed line is the slope of the injection function (αinj = −0.7). An animation of this figure is available.

In the text

Fig. 4. Constraining the magnetic field strength and nonthermal electron injection of the M 87 jet. Left: Modeled spectral index distributions constrained by the observations. The distributions must lie in the gray shaded allowed region (see Sect. 4.4 for details). Colored lines are the modeled spectral index distributions. Different colors represent different initial magnetic field strengths, Bi. The dashed and dotted red line is a modeled spectrum index distribution with q = 10 and q = 11, respectively, with Bi = 0.3 G. Right: Parameter space (Bi, q). The black region is the allowed values of Bi and q constrained by the observations.

In the text

Fig. 5. Magnetic field distribution of the M 87 jet. The black area is the distribution of the magnetic field strength estimated in this study. The solid lines with the gray shaded region are the extrapolation of the upper and lower limits of the result in the upstream direction (B(z)∝z−0.72). The dashed lines are the magnetic field distribution which assumes a steeper slope of the jet radius profile in the region closer to the SMBH (B(z)∝z−0.92). The solid lines without a shaded region are the upper and lower bounds of the magnetic field distribution, assuming fast acceleration in the region closer within 2 mas (B(z)∝z−1.12). The dotted line is the extrapolation of the maximum values of the result to the downstream region. Markers are estimated magnetic field strength from previous studies. The blue cross is the theoretically expected magnetic field strength in a Blandford-Znajek jet (Blandford et al. 2019). The circles are magnetic field strength values obtained in VLBI cores derived using SSA theory (Kino et al. 2014, 2015; Hada et al. 2012, 2016; Reynolds et al. 1996; EHT MWL Science Working Group 2021). The diamonds are values obtained from VLBI observations using various methods other than SSA theory (Acciari et al. 2009; Kim et al. 2018b; Zamaninasab et al. 2014; Event Horizon Telescope Collaboration 2019, 2021; Jiang et al. 2021). The squares are the magnetic field strengths at HST-1 (Harris et al. 2003, 2009). The horizontal dotted-dashed lines are the magnetic field strengths estimated by fitting the broadband SED of the jet using a single-zone SSC emission model (Abdo et al. 2009; MAGIC Collaboration 2020; EHT MWL Science Working Group 2021).

In the text

Fig. A.1. Effect of the (u, v)-range matching on the spectral index distribution. The upper panel shows spectral index distributions between 22 and 43 GHz before (red lines) and after (blue lines) matching the (u, v) range. Distributions using the full (u, v) range have steeper spectra at longer distances. Differences between the two distributions (i.e., artificial steeping due to inconsistent (u, v) ranges) are indicated by gray lines in the lower panel. The thick black line is the average difference, and the red line is the expected artificial steepening along the jet. See text for details.

In the text

Fig. B.1. Spectral index distributions between VLBA 22 and 43 GHz maps convolved with beams of different sizes. The blue lines are convolved with a circular beam of 1.2 × 1.2 mas, where the radius is the same as the major axis of the KaVA 22 GHz beam. The red lines are convolved with 0.8 × 0.8 mas, which is equal to the major axis of the VLBA 22 GHz beam. The distribution with the VLBA beam has a larger scattering than the distribution with the KaVA beam.

In the text

Fig. C.1. Schematic diagram of the spectral index distributions without electron injection, assuming different synchrotron emission models. The blue line assumes the JP model, which is used in this study. The red line shows the spectral index distribution with the KP model. Assuming that the injection spectral index is αinj = −0.7, the expected index after steepening is 4/3αinj − 1 ≈ −1.933, which is within the allowable range obtained from observation.

In the text

Fig. E.1. Comparison of model spectral index distributions using nonthermal electron injection functions with different γinj, max. Distributions with different colors represent different values of γinj, max. For all distributions, the magnetic field strength was assumed to be Bi = 0.3G. As γinj, max increases, the spectral index distributions converge and there is little change between the models with γinj, max = 104.5 and γinj, max = 106.

In the text