M87 black hole mass and spin estimate through the position of the jet boundary shape break

TL;DR

This paper introduces a novel method to estimate the mass and spin of the supermassive black hole in M87 by analyzing the transition in the jet boundary shape, supported by observational data and MHD modeling.

Contribution

The paper presents a new approach linking jet boundary shape transitions to black hole properties, combining observational data with magnetohydrodynamic modeling.

Findings

Black hole mass estimate is larger than previous estimates.

Jet boundary transition correlates with flow regime change.

Magnetic fields play a crucial role in jet dynamics.

Abstract

We propose a new method of estimating a mass of a super massive black hole residing in the center of an active galaxy. The active galaxy M87 offers a convenient test case for the method due to the existence of a large amount of observational data on the jet and ambient environment properties in the central area of the object. We suggest that the observed transition of a jet boundary shape from a parabolic to a conical form is associated with the flow transiting from the magnetically dominated regime to the energy equipartition between plasma bulk motion and magnetic field. By coupling the unique set of observations available for the jet kinematics, environment and boundary profile with our MHD modelling under assumption on the presence of a dynamically important magnetic field in the M87 jet, we estimate the central black hole mass and spin. The method leads us to believe that the M87 super massive black hole has a mass somewhat larger than typically accepted so far.

Tables

Table 1. Table 1: The non-dimensional parameters, which define the position of a “cabing” point, calculated for different magnetizations. The preferred values of σ M subscript 𝜎 M \sigma_{\rm M} , basing on M87 kinematics, are 5, 10, and 20.

${\sigma }_M$	$d_{\ast }$	${\tilde{p}}_{\ast }$
		(${10}^{-5}$)
(1)	(2)	(3)
5	33.6	1.39
10	52.4	1.02
20	79.8	0.75
30	82.0	0.60
40	115.9	0.59
50	134.2	0.51

Table 2. Table 2: Model and derived jet parameters.

${\sigma }_M$	$b$	$\Gamma {\theta }_j$	${\Gamma }_{\max }$	$M$	$a_{\ast }$	${\Psi }_0$	$W_j$	$R_L$	$R_L$
				$\left({10}^9{M}_{\odot }\right)$		(${10}^{33}G{\mathrm{cm}}^2$)	(${10}^{42}$ erg/s)	(pc)	($r_g$)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
5	2.045	$(0.103;\mathrm{\hspace{0.17em}0.024})$	3.4	$7.7\pm 2.7$	$0.17\pm 0.06$	2.9	1.0	0.018	47
10	2.050	$(0.127;\mathrm{\hspace{0.17em}0.033})$	5.1	$6.6\pm 2.1$	$0.22\pm 0.06$	2.8	2.3	0.012	35
20	2.070	$(0.179;\mathrm{\hspace{0.17em}0.057})$	8.6	$5.2\pm 1.5$	$0.26\pm 0.07$	2.9	5.7	0.008	32

Table 3. Table 3: Predicted jet and BH parameters for M = 6.5 × 10 9 ​ M ⊙ 𝑀 6.5 superscript 10 9 subscript M direct-product M=6.5\times 10^{9}\;{\rm M}_{\odot} .

${\sigma }_M$	$\Gamma {\theta }_j$	$\varphi$	$a_{\ast }$	$W_j$
				(${10}^{42}$ erg/s)
(1)	(2)	(3)	(4)	(5)
10	$0.059$	$3.1$	$0.21$	$2.3$
20	$0.060$	$3.1$	$0.32$	$5.7$

Equations

$$\frac{\mathrm{d}}{\mathrm{d}r}\left(\frac{B_{\phi }^2}{8\pi }+P\right)=0.$$

$$P_{\mathrm{ext}}=\frac{B_{\phi }^2}{8\pi }.$$

$$\tilde{p}=\frac{P_{\mathrm{ext}}}{{\left[{\Psi }_0/\left(2\pi R_{\mathrm{L}}^2{\sigma }_{\mathrm{M}}\right)\right]}^2}$$

$$\tilde{d}=\frac{d}{R_{\mathrm{L}}}$$

$$P_{\mathrm{ext}}(r)=P_0{\left(\frac{r}{r_0}\right)}^{-b}.$$

$$R_{\mathrm{L}}=\frac{d_{break}}{d_{\ast }({\sigma }_{\mathrm{M}})}.$$

$$r_{break}=r_0{\left[\frac{{\tilde{p}}_{\ast }({\sigma }_{\mathrm{M}})}{P_0}{\left(\frac{{\Psi }_0}{2\pi R_L^2{\sigma }_{\mathrm{M}}}\right)}^2\right]}^{-1/b}.$$

$$a_{\ast }=\frac{8(r_{\mathrm{g}}/R_{\mathrm{L}})}{1+16{\left(r_{\mathrm{g}}/R_{\mathrm{L}}\right)}^2}.$$

$${\Psi }_0=\varphi \sqrt{\dot{M}c}{r}_g,$$

$$\left(\frac{\dot{M}}{g/s}\right)=C_{\dot{M}}{\left(\frac{M}{10^9{M}_{\odot }}\right)}^2,$$

$$\begin{array}[]{l}\displaystyle\frac{M}{10^{9}M_{\odot}}=1.08\times 10^{2}\left(\frac{d_{break}/d_{*}(\sigma_{\rm M})}{\rm pc}\right)\left(\frac{r_{break}}{r_{0}}\right)^{-b/4}\times\\ \\ \displaystyle\sqrt{\frac{\sigma_{\rm M}}{\phi}}\left(\frac{P_{0}/\tilde{p}_{*}(\sigma_{\rm M})}{10^{-4}\,{\rm dyn/cm^{2}}}\right)^{1/4}\left(\frac{C_{\dot{M}}}{10^{24}\,{\rm g/s}}\right)^{-1/4}.\end{array}$$

$$W_{\mathrm{tot}}=\frac{c}{8}{\left(\frac{{\Psi }_0}{\pi R_{\mathrm{L}}}\right)}^2,$$

$${\dot{M}}_{\mathrm{eject}}\approx \frac{W_{\mathrm{tot}}}{{\sigma }_{\mathrm{M}}{c}^2},$$

Full text

††pagerange: M87 black hole mass and spin estimate through the position of the jet boundary shape break–References

1 Introduction

The object Messier 87 (also known as NGC 4486 and Virgo A; hereafter – M87 in short) is a super-giant elliptical galaxy. At the redshift of $z=0.0043$111NASA/IPAC Extragalactic Database, http://ned.ipac.caltech.edu, accessed 2019.03.14, M87 is one the closest galaxies with active galactic nuclei (AGN). Long before identification as an AGN, the object attracted attention as the first jet, discovered in optical observations a century ago (Curtis, 1918). This jet, later detected in radio emission, has become a test bench for major models of AGN phenomena. Together with the Crab Nebula, M87 was one of the first celestial objects which facilitated the role of synchrotron emission in astrophysics (Shklovsky, 1958).

The bright radio jet in M87 is one of extragalactic structures with best studied morphological properties on the angular scales from arcminutes down to sub-milliarcseconds. The external medium in the inner area of M87 is also best studied among AGN of various classes. M87 is the only galaxy with measurements of particle number density and a temperature of ambient medium at the distance to the central source $\sim 10^{5}$ gravitational radii, which is very close to a Bondi radius (Di Matteo et al., 2003; Russell et al., 2015). There is an extensive information on kinematics and jet transversal structure (see, e.g., Mertens et al. 2016, Asada & Nakamura (2012), Nakamura et al. (2018), Hada et al. (2018), Lister et al. (2019)). All these observing data make the jet in M87 an ideal object for application of the theoretical models which connect the physical properties of the jet and its ambient medium.

A change in M87 jet shape along its extension has been first reported by Asada & Nakamura (2012). It was shown that the power index $k$ in the dependence of jet width $d\propto r^{k}$ on the deprojected distance $r$ from the “central engine” along the jet changes at a $r\sim 100$ pc from $k\approx 0.6$ at small distances to $k\approx 0.9$ at large ones. Later the same “cabing” jet boundary shape behaviour was discovered for 1H0323+342 by Hada et al. (2018), and the position of a break in this source suggested that the mass of a central black hole in 1H0323+342 might be underestimated (Hada et al., 2018). The jet geometry transition was also reported for NGC 6251 (Tseng et al., 2016), for NGC 4261 (Nakahara et al., 2018), and for Cyg-A (Nakahara et al., 2019). As demontstrated recently by Kovalev et al. (2019), a similar morphological pattern in jet shape (“cabing point”) is observed in ten nearby AGN.

The modern AGN paradigm associates many of their manifestations with the presence of a super-massive black hole (SMBH) as the major galactic gravitator. The SMBH mass defines the appearance of AGN and their major observable characteristics. Currently available estimates of the SMBH mass in M87 are based on a variety of measurements and corresponding interpretations. Over the past two decades these values were reported within the range from $M=(3.2\pm 0.9)\times 10^{9}M_{\odot}$ to $M=9.5^{+0.22}_{-0.23}\times 10^{9}M_{\odot}$ (Macchetto et al., 1997; Gebhardt & Thomas, 2009; Gebhardt et al., 2011; Walsh et al., 2013; Oldham & Auger, 2016) based on the dynamical behaviour of various constituencies of galaxy population in the SMBH’s gravitational field.

In this paper we propose a method of estimating BH mass for core-jet AGN that involves another SMBH manifestation – a powerful relativistic jet launched from the circumnuclear area of the source. The method is based on MHD modelling (Beskin et al., 2017) of a jet boundary shape and matching the model’s “cabing” point in the jet shape with its observed position. We estimate the central BH mass and spin independently following the theoretical model by Beskin et al. (2017), and using the measurements of jet parameters in M87: the ambient pressure, the plasma flow kinematics, jet opening angle, and the position of a jet shape break – the “cabing” point.

The paper is organized as follows. We describe the multifrequency observational data used to recover the M87 jet shape on the scales from $10^{-1}$ pc to $10^{4}$ pc and the “cabing” point position. In Section 3 we discuss briefly the MHD model that allows us to reconstruct the observed break in a jet shape for the smooth ambient pressure profile as well as all the model assumptions and the values needed to obtain the black hole mass and spin rate. In Section 4 we define the method of fitting the jet profile by two power-laws, and in Section 5 we discuss the errors. We present the results in Section 6.

Throughout the paper, we use the $\Lambda$CDM cosmological model with $H_{0}=71$ km s*-1* Mpc*-1*, $\Omega_{m}=0.27$, and $\Omega_{\Lambda}=0.73$ (Komatsu et al., 2009).

2 Observational data

We use the multi-frequency radio interferometry data, reported by Asada & Nakamura (2012), Hada et al. (2013), and Hada et al. (2016), and collected in the paper by Nakamura et al. (2018). For each data set we use distance along the jet taking into account the error for the core data, and a radius of a detected feature with the error in radius determination (see Figure 1). The data sources and their thorough description are as follows. The 1.8 GHz data are obtained with MERLIN (Asada & Nakamura, 2012). The 2.3, 5.0, 8.4, and 22 GHz data come from the Very Long Base Array (VLBA) as reported by Hada et al. (2013). The 15 and 43 GHz VLBA data have been reported by Asada & Nakamura (2012) and Hada et al. (2013). The 86 GHz data set is provided by the VLBA–High Sensitivity Array (HSA) observations Hada et al. (2016). Nakamura et al. (2018) use the luminosity distance $D_{\rm L}=16.7$ Mpc (Blakeslee et al., 2009).

The VLBA core data at frequencies 5.0, 8.4, 15.4, 23.8, 43.2, 86.3 GHz are dedscribed by Hada et al. (2013) and at frequencies 43 and 86 GHz – by Nakamura & Asada (2013); Hada et al. (2016). The Event Horizon Telescope (EHT) core data at 230 GHz obtained by Doeleman et al. (2012); Akiyama et al. (2015). However, we do not use these core data to fit the jet boundary form in the parabolic domain due to large errors in the determination of the core position along the jet due to the core shift estimates (Hada et al., 2011). We note however, that, described below, the major contribution in the BH mass estimate is provided by the jet boundary data on the scales larger than those of the core.

Detailed procedures of estimating the jet width are described in Asada & Nakamura (2012) and Hada et al. (2013). In short, we made transverse slices of the jet at various distances from the core. For each slice, we fitted a double-Gaussian function (if the slice is clearly resolved into a two-humped shape, which applies to most of the slices) or a single Gaussian (if the slice is single-peaked). We then defined the separation between the outer sides of the half-maximum points of the two Gaussians as the width of the jet at each distance (for the single Gaussian case, its deconvolved FWHM was taken as the jet width). Finally, the jet radius ($d$) at each distance was defined as a half of the jet width.

3 Black hole mass and spin determination

As was already stressed, the uniqueness of the M87 jet is in the availability of direct information not only on the jet boundary shape, but also on the ambient pressure $P_{\rm ext}$ in close vicinity of the jet “cabing” region (Young et al., 2002; Di Matteo et al., 2003; Russell et al., 2015). Below we show that this additional information gives us the possibility to determine such key parameters of the “central engine” as the total magnetic flux $\Psi_{0}$ in the jet and the radius of the light cylinder $R_{\rm L}=c/\Omega$. In turn, magnetically arrested disk (MAD) assumption allows us to decouple the mass $M$ of the supermassive black hole and its spin parameter $a_{*}$. We designate the distance along the jet and the jet radius as $r$ and $d$, respectively. The function $d(r)$ determines the jet boundary shape. The position and radius of the “cabing” point at which the jet shape changes from parabolic to conical is designated as $r_{break}$ and $d_{break}$, respectively.

Below we use a model of the transversal structure of a jet based on the now generally accepted MHD theory within the framework of the approach of the Grad-Shafranov (GS) equation (Heyvaerts & Norman, 1989; Pelletier & Pudritz, 1992; Heyvaerts, 1996). More precisely, on its one-dimensional cylindrical version, when a second-order partial differential GS equation can be reduced to two first-order ordinary differential equations (Beskin, 1997; Lery et al., 1998, 1999; Beskin & Malyshkin, 2000; Beskin & Nokhrina, 2009). This approach has well proven itself for both non-relativistic and relativistic flows. In particular, just within this approach, it has been predicted theoretically that in a parabolic magnetic field, effective particle acceleration becomes possible (Beskin & Nokhrina, 2006). Earlier, on the basis of solutions for quasi-spherical outflow, it was believed that effective acceleration in a magnetically dominated wind is impossible (Michel, 1969; Kennel et al., 1983; Bogovalov, 1997). Later this conclusion has been repeatedly reproduced by numerical simulations (see, e.g., McKinney 2006; Narayan et al. 2007). Among other things, demonstrating the full consistency of a semi-analytical modelling with the numerical simulations, the existence of a denser core along a jet axis was obtained by Beskin & Nokhrina (2009). It was corroborated by independent numerical models (Komissarov et al., 2007; Tchekhovskoy et al., 2009; Porth et al., 2011). Finally, it was also shown by  Beskin & Zheltoukhov (2013) how asymptotic relations obtained in the framework of the one-dimensional approach used in our present work make it possible to reproduce convincingly the results of numerical simulations for the black hole magnetosphere obtained by McKinney et al. (2012).

In what follows we use the most developed version, in which we assume that an electric current $J$ locked inside the jet (Beskin et al., 2017). In this model, the flow velocity and magnetic and electric fields vanish at the jet edge $d(r)$. In this case, the current sheet at the edge is absent. In numerical modelling, such a structure has been known for non-relativistic trans-sonic flows (Romanova et al., 2009). Recently this structure was reproduced for relativistic outflows as well (Bromberg & Tchekhovskoy, 2016). The fall of a flow bulk motion Lorentz factor down to unity at the jet boundary is clearly seen in the numerical simulations by Nakamura et al. (2018), in accordance with the assumption used here.

Assuming that the flow remains supersonic up to the very boundary of a jet, one can write down the force balance at the jet boundary as

[TABLE]

Here $B_{\varphi}$ is a toroidal magnetic field, which dominates the poloidal field $B_{\rm p}$ outside the light cylinder, and $P$ is a jet plasma pressure that transits smoothly into the pressure of the ambient environment. Integrating this equation through the thin boundary layer where the gradient of the gas pressure balances the magnetic stress, we obtain

[TABLE]

Indeed, as was shown by Kovalev et al. (2019), even for finite temperature the magnetic pressure dominates the force balance inside a jet up to the very thin boundary layer. The importance of a gas pressure at the jet boundary is supported by numerical simulations in (Nakamura et al., 2018).

Here we must emphasise the key difference of our model comparing to other ones (see, e.g., Lyubarsky 2009). We note that in the framework of the approach considered here it is necessary to specify five integrals conserved on magnetic surfaces (energy density flux, angular momentum density flux, angular velocity of field lines, entropy, and mass-to-magnetic flux ratio). For the major part of the jet, we use standard values prescribed by the condition of a smooth crossing of the singular surfaces (Alfvénic and fast magnetosonic). However, near the outer boundary of the jet, the integrals were chosen in such a way that the condition of the total zero longitudinal current within the jet $J(d)=0$ was satisfied. As was already stressed, such a structure of the integrals of motion corresponds to the results of numerical simulation (Romanova et al., 2009; Bromberg & Tchekhovskoy, 2016).

Thus, the solution obtained by Beskin et al. (2017) provides that the major part of an electric current is locked inside the bulk jet volume, and only a residual electric current $J_{\rm res}$, that defines $B_{\varphi}$ in Equation 2, is left in the outer thin jet layer. This implies that the characteristic toroidal magnetic field $B_{\varphi}=2J_{\rm res}/cd$, which constitutes the major pressure at the jet edge, is much lower than in the models without an electric current closure.

For the cylindrical geometry, the Grad-Shafranov and Bernoulli equations, describing a full MHD flow, become a set of ordinary differential equations easily solvable numerically. The cylindrical flow solution reproduces accurately the axisymmetric flow solution if the derivatives along the jet are negligible. It was shown by Nokhrina et al. (2015) that the solution, obtained within the cylindrical geometry, is applicable for formation of the jet structure for the highly collimated flows. This allows us to use the cylindrical approach to the problem of modelling a well collimated jet. In non-dimensional variables the solution of these equations depends only on the Michel’s magnetization parameter $\sigma_{\rm M}$, which is defined as the ratio of Poynting flux to the plasma rest mass energy flux at the base of a flow. Integrating the system of two ordinary differential equations describing internal structure of a jet (see Beskin et al. (2017) for more details), we obtain a non-dimensional external pressure

[TABLE]

as a function of non-dimensional jet radius

[TABLE]

for different initial magnetizations $\sigma_{\rm M}$. Here we use the natural inner scale for both poloidal and toroidal magnetic field $B_{scale}=\Psi_{0}/(2\pi R_{\rm L}^{2}\sigma_{\rm M})$, written through the total magnetic flux in a jet $\Psi_{0}$. The corresponding scale for pressure is $B_{scale}^{2}$. The integration of MHD equations for the given integrals (see Beskin et al. 2017) provides the numerical factor, that relates this pressure scale with the corresponding jet inner pressure, needed to balance the ambient pressure.

Below we assume a power-law dependence of the ambient pressure on the distance from the central source:

[TABLE]

Here $P_{0}$ is the ambient pressure amplitude at the distance $r_{0}$ from the BH. The exponent $b$ attains values between 1 and 2.5. The largest value 2.5 corresponds to the supersonic regime of a Bondi accretion of a gas described by the adiabatic equation of state $P\propto n^{\gamma}$ with $\gamma=5/3$. However, the recent theoretical studies of a gas accretion onto SMBH provide smaller values $b\in(1.0,\,2.1)$ (Quataert & Narayan, 2000; Narayan & Fabian, 2011), and the recent observations by Park et al. (2019) favour $b\lessapprox 2.0$. Thus, within our model we are able to determine the jet boundary shape $d(r)$ for the given ambient pressure profile $P_{\rm ext}(r)$.

As was shown by Beskin et al. (2017), the obtained jet boundary dependence $d(r)$ has a pronounced break in the domain, where the flow transits from magnetically-dominated regime to the quasi-equipartition of plasma bulk motion kinetic energy density and the energy density of magnetic field. For the pressure profile predicted by the Bondi accretion model with $b\approx 2$, we obtain a clear transition from a parabolic to conical shape consistent with the results by Asada & Nakamura (2012) and Nakamura et al. (2018).

In our semi-analytical solution, the non-dimensional jet radius $d_{*}$ and ambient pressure $\tilde{p}_{*}$ are defined as functions of the Michel magnetization parameter $\sigma_{\rm M}$. The results of the simulations are presented in Table 1. These simulations provide the position of the “cabing” point. The essence of our method is in comparing the jet’s geometry at the “cabing" point as obtained in the simulations with the observed shape of the jet.

Using now Equation 4 we obtain for the light cylinder radius

[TABLE]

On the other hand, Equation 5, rewritten for the “cabing” point, together with Equation 3, allows us to find the total magnetic flux $\Psi_{0}$ in a jet with the measured $r_{break}$:

[TABLE]

Here the pressure amplitude $P_{0}$ at the distance $r_{0}$ is known from the observations, while $\sigma_{\rm M}$ and $\tilde{p}_{*}$ — from the modelling.

The results presented above are direct outcomes of MHD modelling of the jet structure (Beskin et al., 2017). For MHD models, the intrinsic length scale is the light cylinder radius $R_{\rm L}$, not the gravitational radius $r_{\rm g}=GM/c^{2}$. Thus, the position and radius of the cabing point depends on both the BH mass and its spin. Indeed, $R_{\rm L}$ and $r_{g}$ can be related for the maximum BH energy extraction rate condition $\Omega_{\rm F}=\Omega_{\rm H}/2$ (Blandford & Znajek, 1977). Here $\Omega_{\rm F}$ is a field lines rotational velocity, and $\Omega_{\rm H}$ is a BH angular velocity. Introducing the BH spin $a_{*}\in[0;\;1]$, we obtain the relation between $a_{*}$, $r_{\rm g}$ and $R_{\rm L}$:

[TABLE]

Gravitational radius may be recovered if we assume that the total magnetic flux $\Psi_{0}$ is locked with the mass accretion rate $\dot{M}$ (Narayan et al., 2003). The numerical simulations by Tchekhovskoy et al. (2011) provide the following dependence:

[TABLE]

with $\phi\sim 50$ in Gaussian units (Tchekhovskoy et al., 2011) for a disk being in a magnetically arrested state (MAD). The same relation holds for a standart and normal evolution disk (SANE), with lower values of $\phi$ (Narayan et al., 2012). An analysis of a sample of 76 radio-loud sources (Zamaninasab et al., 2014) gave the same result with $\phi=(52\pm 5)\Gamma\theta_{j}$, where $\Gamma$ is a Lorentz factor of a bulk flow, and $\theta_{j}$ is a jet half-opening angle. For $\Gamma\theta_{j}\ll 1$ the disk state is SANE, not MAD (see discussion in section 6). For the Bondi accretion, we use the expression for an accretion rate that depends on the mass of a central BH. It is defined by the relation

[TABLE]

where $C_{\dot{M}}$ depends on the ambient gas particle number density and temperature (see Di Matteo et al. (2003) for more detail). Substituting now Equation 6, Equation 9, and Equation 10 into Equation 7, we finally obtain the following expression for the BH mass in M87:

[TABLE]

We stress that we do not use in this formula any results that were obtained under a priory assumption on the BH mass. Similarly, we use the expression Equation 10 and do not use a direct estimate of the accretion rate.

The described above method is based on an assumption that the real jet boundary, determined by the condition $\Psi=\Psi_{0}$, corresponds to the visible jet boundary. The latter is obtained as a cut at half maximum of intensity. Although in general case their coincidence might not be exact, it holds for our jet transversal structure model. We assume that the synchrotron self-absorbed emission is produced by highly relativistic plasma with an energy distribution $dn=k_{e}\Gamma^{-p}d\Gamma$. The emitting particle number density amplitude $k_{e}$ is either equal or proportional to the total local particle number density in a jet (Lobanov, 1998; Nokhrina et al., 2015). The intensity depends on the emission $\rho$ and absorption $\kappa$ coefficients for a synchrotron emission (Ginzburg & Syrovatskii, 1965). They, in turn, are defined by the plasma conditions: the particle number density of emitting plasma and magnetic field amplitude. In case of an optically thick part of a flow, the intensity $I$ is depends on the magnetic field roughly as $\propto B^{-1/2}$, while in the optically thin region $I\propto nB^{(p+1)/2}$. In both cases, the profiles of $n$ and $B$ are such (Chernoglazov et al., 2019), that the intensity grows towards the jet boundary, falling rapidly only in a very thin layer in its vicinity. As the flow is relativistic, the Doppler factor also affects the received intensity. For a high bulk Lorentz factor, the observer may be out of a cone of emission and receive the suppressed intensity, as can be seen in the Doppler maps by Chernoglazov et al. (2019). Thus, we expect that the observed jet boundary corresponds indeed to the model jet boundary $d(r)$. The effect of a jet slowing down at the boundary is expected in real jets and supported by the numerical simulations by McKinney (2006); Dexter et al. (2012); Nakamura et al. (2018).

4 Break in the M87 jet shape

As was reported by Asada & Nakamura (2012), the M87 jet boundary shape changes from approximately parabolic ($d\propto r^{0.5}$) to approximately conical ($d\propto r$). Our modelling predicts such the transition (the cabing point) as the flow accelerates from initially magnetically dominated regime to the energy equipartition. The change in a jet boundary shape occurs without a change in an ambient pressure profile. Thus, in order to compare theoretical predictions with observational data, we need to pin the observed position of the cabing point by approximating the jet boundary shape by two power laws to determine the SMBH mass.

The procedure of fitting the two power laws is as follows. We use the MERLIN and VLBA imaging data as described in section 2 at frequencies 1.8, 2.3, 5.0, 8.4, 15.0, 22.0, 43.0, and 86.0 GHz (Asada & Nakamura, 2012; Hada et al., 2013, 2016; Nakamura et al., 2018). For each frequency, we have a set of measured jet radii $d$, de-projected (for the viewing angle of $14^{\circ}$) distances along the jet $r$, and an error in determination of $d$. Figure 1 represents the $d(r)$ dependence for observational data their fit by two power laws. The first guess is that the change the power law index (cabing) occurs at the distance corresponding to the data obtained at 2.3 GHz. After this rough guess we divide the full sample into two sets (“parabolic” and “conical”) choosing a point from the 2.3 GHz sample as a boundary between them. For each such choice, we fit the power law parameters for two data sets of the full sample. We choose from the resulting set of possible approximations one that minimises the standard error in the expected conical domain. However, we also find the position of cabing point for every cut inside the 2.3 GHz sample, and use them to estimate the error in BH mass and spin values due to possible uncertainty in the cabing point determination.

The obtained jet shape break for the “best” choice of a sample cut is at $r_{break}=43.41$ pc with the corresponding jet width radius $d_{break}=0.60$ pc. The power laws are: $d_{\rm pc}=0.07\,r_{\rm pc}^{0.57}$ for the parabolic domain, and $d_{\rm pc}=0.02\,r_{\rm pc}^{0.90}$ for the conical domain. The result of this fit is presented in Figure 1. Here $d_{\rm pc}$ is a jet radius measured in pc, and $r_{\rm pc}$ is a distance along a jet in pc.

We use the full data at 1.8 GHz in contrast with Nakamura et al. (2018), who excluded the farthest 4 points as a suspected jet wiggle that drives the conical domain fit to be more shallow. We have checked how deleting these points alter the results. We note that, indeed, the conical domain fit becomes steeper: $d\propto r^{0.92}$, with the cabing position moving to $r_{break}=45.16$ pc and $d_{break}=0.62$ pc. However, this changes the final results for the mass and spin at the level $\sim 0.1\%$. In fact, this demonstrates the robustness of our result. The fit in the parabolic domain is very well defined. The final expression for the mass Equation 11 has a term $d_{break}/r_{break}^{b/4}$, which varies very slowly as the jet boundary shape break follows the nearly parabolic trend (parabola holds, conical domain changes its slope). Because of this, we do not need to exclude the points that possibly reflect the local jet wiggle.

5 Error budget

There are four major sources of errors in the mass determination by the method presented here: (i) errors from determining the cabing position; (ii) errors in the jet half-width; (iii) errors in $C_{\dot{M}}$ determination; (iv) errors due to uncertainty in $\Gamma\theta_{j}$ estimate.

In the available data set, the cabing position depends on our choice of attributing the data points from 2.3 GHz sample to parabolic or conical domains. We calculated the cabing position for every point from 2.3-sample being an expected boundary between parabolic and conical domains. We observe that, starting from our “best” choice (77 of 394), the cabing position stays more or less constant until about 250 of 394 points of 2.3 GHz sample are attributed to the parabolic domain. We calculate the mean and standard deviation for this “plato” sample of cabing positions and jet half-widths. The result is: $d_{break}=0.62\pm 0.02$ pc and $r_{break}=44.7\pm 1.9$ pc.

We have performed bootstrapping to model how errors (see Figure 1) in a jet half-width determination affect the cabing position for the “best” choice of dividing the sample into parabolic / conical domains. The mean values with standard deviations are $d_{break}=0.61\pm 0.02$ pc and $r_{break}=44.5\pm 1.9$ pc.

We conclude that the typical errors arising from VLBI and MERLIN data are $\pm 0.02$ pc for $d_{break}$ and $\pm 1.9$ pc for $r_{break}$.

The errors provided by Di Matteo et al. (2003) give the error in the total expression for mass around 1% due to errors in particle number density $n=0.170\pm 0.003\;{\rm cm^{-3}}$ and temperature $kT=0.80\pm 0.01\;{\rm keV}$ measurements.

However, the scatter in numerical values for $\Gamma\theta_{j}$ make the major contribution into the error budget, being finally of about an order higher than the errors due do $r_{break}$ and $d_{break}$ position modelling and the errors in measurements of $C_{\dot{M}}$. Thus, we present the result for BH mass as an interval of values corresponding to obtained by our modelling interval for $\Gamma\theta_{j}$ with the errors (i)–(iii). Also we directly give the mean value for BH mass, with the total error including the uncertainty due to $\Gamma\theta_{j}$ values. See the next Section.

6 M87 black hole mass

We use the closest to the central BH pressure and density measurements by Russell et al. (2015): $kT=0.91$ keV and $n_{e}=0.31\;{\rm cm^{-3}}$ at approximately $r_{0}=0.22$ kpc. This gives the pressure amplitude $P_{0}=0.45\times 10^{-9}\;{\rm dyn/cm^{2}}$ at $r_{0}=0.22$ kpc from the BH.

The Bondi mass accretion rate was obtained by Di Matteo et al. (2003) basing on measurements of density and temperature of a hot interstellar medium (ISM) using the observed X-ray emission at distances $\lesssim 100$ pc from the black hole. We use the Equation 10 with the obtained numerical value $C_{\dot{M}}=7\times 10^{23}\;{\rm g/s}$.

We choose the value of an exponent $b$, defining the pressure profile, so as to fit the observed jet shape. For example, for $b=2.07$ and $\sigma_{\rm M}=20$ the model predicts $d\propto r^{0.57}$ for a parabolic domain and $d\propto r^{0.82}$ for a conical domain. In fact, we do not fit precisely both power–laws, describing the observational data in section 4. We set the exponent $b$ so as to fit the parabolic domain (see Figure 1). In this case, the conical domain of our model still fits the data within the error bars.

We also need to choose the initial magnetization parameter $\sigma_{\rm M}$. The magnetization $\sigma_{\rm M}$ defines the maximum bulk flow Lorentz factor, that can be achieved by the flow if all the Poynting flux energy is converted to the bulk plasma kinetic energy. It was shown (see, e.g., Beskin & Nokhrina 2006; Tchekhovskoy et al. 2009; Komissarov et al. 2009; Lyubarsky 2009), that the plasma in highly collimated outflows accelerates effectively only up to approximately $\Gamma\sim\sigma_{\rm M}/2$. Further downstream, the acceleration continues very slowly. Thus, the observed Lorentz factors in M87 can provide us with the estimate for $\sigma_{\rm M}$. Mertens et al. (2016) obtained the Lorentz factors of the order of $\Gamma\sim 3$ at $r\sim$ a few parsecs. The detected by Biretta et al. (1999) Lorentz factors at few hundred parsec is $\sim 10$. We present here results for three values of the magnetization parameter $\sigma_{\rm M}$: 5,10, and 20, which are consistent with the discussed above observed bulk flow Lorentz factors. For these three models we calculate the predicted jet shape profile and find the cabing point position.

There are theoretical as well as observational constraints on $\Gamma\theta_{j}$. It was discovered by Tchekhovskoy et al. (2009) and Komissarov et al. (2009) that the condition $\Gamma\theta_{j}<1$ corresponds to the casual connectivity across a jet, ensuring the effective plasma acceleration up to equipartition. Komissarov et al. (2009) showed that the approximate equality $\Gamma\theta_{j}\approx 1$ should hold for the power-law acceleration regime in a jet. On the other hand, the observations provide the median value $\Gamma\theta_{j}=0.17$ (Pushkarev et al., 2017). The high resolution data obtained by Mertens et al. (2016) also allows to estimate this value for M87 jet specifically. Measurements and analysis by Mertens et al. (2016) provide the apparent opening angle $\theta_{\rm app}$ at the distances $\sim 0.3-4.0$ pc varying from about $18^{\circ}$ closer to the BH to $\approx 7^{\circ}$ further downstream. The intrinsic opening angle depends on the apparent opening angle as $\theta_{j}=\theta_{\rm app}\sin\varphi/2$, where we use the same viewing angle $\varphi=14^{\circ}$ (Nakamura et al., 2018) as was used for the de-projection for the result self-consistency. It gives $\theta_{j}\approx 0.038$ at $r=0.3$ pc and $0.015$ at $r=4.0$ pc. The Lorentz factor at the same scales varies (Mertens et al., 2016) from roughly $1.2$ to $\approx 3$. This provides $\Gamma\theta_{j}\sim 0.046$ at $r=0.3$ pc and $\Gamma\theta_{j}\sim 0.044$ at $4$ pc, the resultant value being much smaller than theoretical upper boundary for this value. In this paper we use the results of our modelling of a jet structure to bound the possible values of $\Gamma\theta_{j}$. We calculate the maximum Lorentz factor across a jet and the jet shape boundary derivative $d/r=\tan\theta_{j}$ for each $r$. We observe that the parameter $\Gamma\theta_{j}$ does not stay constant along the jet (as was first observed by Komissarov et al. (2009)). It starts at the value $\approx 0.1$ in the parabolic region and runs down up to approximately the cabing point, where it assumes a constant value, corresponding to the maximal Lorentz factor, attained by the jet for the given magnetization, multiplied by the roughly constant opening angle of conical domain.We use the interval of values for this parameter in our Equation 11 for the BH mass determination. The scatter in this parameter provides the major contribution into errors in the final result for the BH mass. We should note that $\Gamma\theta_{j}$, obtained within our modelling is consistent with the result by (Mertens et al., 2016), but differs strongly from the assumption $\Gamma\theta_{j}=1$ by Zamaninasab et al. (2014).

The model parameters and results are presented in Table 2. The model parameters that we set are in columns (1)–(2): the initial jet magnetization and the exponent $b$ set to fit exactly the parabolic domain jet boundary shape. The calculated parameters of the central BH and jet are in columns (3)–(10). We calculate within our model non-dimensional jet shape break parameters $d_{*}$ and $\tilde{p}_{*}$ (see Table 1), the interval for $\Gamma\theta_{j}$, and the maximum Lorentz factor of bulk motion, attained by the flow. We calculate the BH mass using the Equation 11. We use the measured values for $C_{\dot{M}}$, $P_{0}$, and $r_{0}$. We use the results of our fitting the observational data for $d_{break}$ and $r_{break}$ (see section 4), which are consistent with the results by Asada & Nakamura (2012). We also put the model parameters $\sigma_{\rm M}$, $d_{*}$, $\tilde{p}_{*}$, and our estimates for $\Gamma\theta_{j}$. For the result for the SMBH mass in Table 2 we present the median value, obtained for each magnetization for the interval $\Gamma\theta_{j}$ with an error due to uncertainty in this parameter. The same is for the BH spin $a_{*}$, which we find using Equation 8. The result for the total magnetic flux obtained using Equation 7 depends on the model parameters and pressure measurements only. To calculate the total jet power we use the expression (Nokhrina, 2018)

[TABLE]

that relates jet power $W_{j}$ with the Poynting flux power at the jet base. This formula neglects the initial power in plasma kinetic energy, which is justified for sufficient magnetizations.

We obtain the different values of BH mass for different magnetizations:

$M\in(5.0\pm 0.3;\;10.4\pm 0.6)\times 10^{9}\,M_{\odot}$ for $\sigma_{\rm M}=5$;

$M\in(4.4\pm 0.3;\;8.7\pm 0.5)\times 10^{9}\,M_{\odot}$ for $\sigma_{\rm M}=10$;

$M\in(3.8\pm 0.2;\;6.7\pm 0.4)\times 10^{9}\,M_{\odot}$ for $\sigma_{\rm M}=20$.

Here the mass interval is due to calculated interval for $\Gamma\theta_{j}$, and errors are due to errors in jet half-width determination, errors in determination of the cabing point, and errors in mass accretion rate estimate. The corresponding intervals for the BH spin:

$a_{*}\in(0.11\pm 0.01;\;0.22\pm 0.01)$ for $\sigma_{\rm M}=5$;

$a_{*}\in(0.15\pm 0.01;\;0.28\pm 0.02)$ for $\sigma_{\rm M}=10$;

$a_{*}\in(0.19\pm 0.01;\;0.33\pm 0.02)$ for $\sigma_{\rm M}=20$.

7 Results and discussion

Within the jet model with an electric current locked inside a jet (Beskin et al., 2017) we obtained a clear transition from parabolic to conical jet boundary shape for the ambient pressure given by Equation 5. The break in a jet form occurs as the flow transits from magnetically dominated regime to the rough equipartition between plasma bulk motion kinetic and magnetic field energy. We propose to associate the positions of a model break with the observed one to obtain the jet and BH parameters. Together with an assumption for the dynamically important magnetic field near the BH and measurements of a pressure amplitude, mass accretion rate, and kinematics, we are able to estimate the BH mass, spin, and total magnetic flux in a jet.

In order to obtain the jet shape break position we use the VLBA and MERLIN data collected by Nakamura et al. (2018). We should note that not all data from that paper were used here to obtain the result. First of all, our aim is to pin the cabing point position, so the data that is as uniform as possible and that represents the jet boundary shape is of the most importance. We exclude the core data due to errors in core position (Hada et al., 2011). We also do not use the European Very Long Base Interferometry Network (EVN) data, because these data cover mainly the HST-1 complex around $r\sim 100$ pc, and may not reflect the jet shape behaviour, but rather the special features of HST-1 itself.

The resultant mass depends on the initial magnetization $\sigma_{\rm M}$. For all the values of $\sigma_{\rm M}$, the obtained mass is much bigger than that obtained by Walsh et al. (2013) using a spectral analysis of gas velocity dispersion. For the high magnetization $\sigma_{\rm M}$, it is somewhat consistent with the result by Gebhardt et al. (2011), based on both gas velocity dispersion measurements and stellar dynamics. The median value for $\sigma_{\rm M}=10$ is close to the result of Oldham & Auger (2016), obtained basing on the analysis of stellar and cluster dynamics, although the scatter in mass in our paper is bigger due to the strong dependence of the result on the parameter $\Gamma\theta_{j}$. But the magnetization $=5$, favoured by the kinematics detected in (Mertens et al., 2016), points to even bigger BH mass value.

The total magnetic flux depends very weakly on $\sigma_{\rm M}$, and its value of the order of $10^{33}\;{\rm G\,cm^{2}}$ is consistent with the results by Nokhrina (2018). In contrast, the BH spin depends on $\sigma_{\rm M}$, as the light cylinder radius $R_{\rm L}$ depends on $d_{*}$ only. This gives the scatter in $a_{*}$ estimate from $0.11$ to $0.33$ for different magnetizations. Thus, we have obtained the moderate spin parameter of the order of $0.1-0.3$ for M87 SMBH. The numerical simulations (Tchekhovskoy et al., 2011; McKinney et al., 2012) favour the spin $>0.5$ in order to obtain the jet power of the order of $\dot{M}c^{2}$. On the other hand, semi-analytical and numerical modelling of a BH spin evolution (King et al., 2008; Barausse, 2012; Volonteri et al., 2013; Sesana et al., 2014) predict moderate spins $a_{*}\in(0.1,\,0.7)$ for low redshift $z<2$ AGN, with BH residing in elliptical galaxies tending to have smaller spins, which is consistent with our result.

The result for the total jet power, obtained with Equation 12, is consistent with the estimates of an average jet power $W_{j,av}\sim 3\times 10^{42}\;{\rm erg/s}$, needed to evacuate the inner cavities (Young et al., 2002). This result is also marginally consistent with the jet power obtained by Levinson & Globus (2017) within a model of the recollimation shock in HST-1 due to a jet interaction with a disk outflow. However, the other theoretical modelling by Stawarz et al. (2006) predict higher jet power $\sim 10^{43}-10^{44}\;{\rm erg/s}$ needed to feed the radio lobes (Owen et al., 2000). This may be an indication that the rough estimate of jet power by Beskin (2010) without a numerical factor is more robust, providing for M87 total jet power a few of $10^{43}\;{\rm erg/s}$. Indeed, the factor $1/8$ does provide the correlation of a magnetic flux with the averaged over large period of time power (Nokhrina, 2018), but it also depends on the particular choice of MHD integrals.

We are able to fit the the observed jet boundary shape with the theoretical curve in the parabolic domain with $b\approx 2.05-2.07$, which is consistent with Bondi accretion flow models Quataert & Narayan (2000); Narayan & Fabian (2011). However, the direct measurements of a particle number density in ISM either by X-ray observations (Russell et al., 2015), or by modelling the Faraday rotation measure on the ambient medium (Park et al., 2019), provide $n\propto r^{-1}$, which corresponds to smaller $b$ for adiabatic flow. This caveat may be solved if the temperature rises closer to the central source, as was predicted by Quataert & Narayan (2000) and discussed in Russell et al. (2015).

Our model means that the more or less effective plasma acceleration takes place up to the cabing point, or, approximately, up to HST-1 (Asada & Nakamura, 2012). This is different from modelling by Mertens et al. (2016), in which the acceleration saturation is set at the distance approximately $4$ pc from the jet base. However, the longer acceleration domain obtained within our model is consistent with the observed by Biretta et al. (1999) Lorentz factors of the order of $10$ at the HST-1 (the cabing region (Asada & Nakamura, 2012)). The kinematics obtained by radio interferometric measurements (Mertens et al., 2016; Lister et al., 2019) with low detected Lorentz factors favour the smaller magnetizations $\sigma_{\rm M}=5\div 10$, with the predicted bigger central BH mass. On the other hand, the optical observations of velocities at $\sim 100$ pc scales by Biretta et al. (1999) favour $\sigma_{\rm M}=20$.

As was stressed, the method for determining BH mass proposed above is based on the existence of a statistical dependence Equation 9, which relates the total magnetic flux $\Psi_{0}$ to the accretion rate ${\dot{M}}$. In cases where the accretion rate ${\dot{M}}$ can be found independently, the procedure for determining the mass may be changed. In particular, note that the relation Equation 10 provides the accretion rate ${\dot{M}}\approx 0.2-0.4\,M_{\odot}$ yr*-1* for the masses we estimated for M87. We plan to address the question of mass determination for the other sources with the detected jet boundary shape break, in particular 1H0323+342 (Hada et al., 2018), in the future work.

Using the definition of a magnetization parameter $\sigma_{\rm M}$, one can rewrite the mass ejection rate ${\dot{M}}_{\rm eject}$ in a jet in the form

[TABLE]

where $W_{\rm tot}$ Equation 12 is the total energy losses in a jet. This value does not depend on the assumed jet composition. For $\sigma_{\rm M}\sim 10$ we obtain the reasonable value of the mass ejection rate ${\dot{M}}_{\rm eject}\sim 10^{-4}M_{\odot}$ yr*-1*. This mass loss rate in a jet corresponds to the mass density $\rho$ through the given jet cross section $\dot{M}_{\rm eject}=\rho cS$. At the distance 1 pc from the “central engine” (typical distance where the particle number density is calculated through the core shift effect) the jet radius $\sim 0.1$ pc. If the mass density is defined by the electrons, than the particle number density at 1 pc is of the order of 100 ${\rm cm}^{-3}$, which is in agreement with another independent evaluation by core-shift data Nokhrina et al. (2015). For the protons the appropriate particle number density is 1800 times less. Therefore, this result points at the mainly electron-positron composition of the M87 jet.

7.1 Corroboration by the EHT results

The brand new EHT results (EHT Colaboration et al., 2019) provide the BH mass in M87 as $M=(6.5\pm 0.7)\times 10^{9}\;{\rm M}_{\odot}$. This result is in agreement with the choice of $\sigma_{\rm M}=10,\;20$, in accordance with the observed kinematics (Biretta et al., 1999; Mertens et al., 2016).

The BH mass value $M=6.5\times 10^{9}\;{\rm M}_{\odot}$ corresponds to the BH and jet properties, listed in Table 3.

We predict that the jet is highly casually connected $\Gamma\theta_{j}\ll 1$. The disk state is far from the MAD ($\phi\sim 50$), and obtained value $\phi\approx 3$ (in Gaussian units) suggests the SANE disk state. We also predict the moderate spin of the order of $0.2-0.3$. This value has not been probed by the (EHT Collaboration & et al., 2019) modelling. The total jet power corresponds to the highest obtained by the EHT collaboration results, being closer to the estimates obtained in the previous works (Stawarz et al., 2006; Owen et al., 2000; Young et al., 2002; Levinson & Globus, 2017). Again, this power may be higher by the factor of about four, but this needs further investigation.

The proposed method of estimating the BH mass and spin, total magnetic flux in a jet, and total jet power, may prove to be a powerful instrument in probing the BH physics. It is in full accordance with the EHT results and multitude of previous studies of BH environment, as well as jet morphology and kinematics. At the same time, this instrument presented here requires resolution of jet boundary shapes on the scale of tens of parsecs, or $10^{5}$ gravitational radii, which is an attainable goal for the modern VLBI systems. We also note that the comparison of the “traditional” cm–dm–wavelength VLBI results discussed in this work and EHT results on M87 will offer a powerful calibration method for future interpretation of high-resolution studies in many AGN.

Acknowledgements

We thank the anonymous referee for suggestions which helped to improve the paper. This research was supported in part by the 5-100 Russian Academic Excellence Project (Agreement number 05.Y09.21.0018) and by the Russian Foundation for Basic Researches (grant 17-02-00788).