Free Access
Issue
A&A
Volume 632, December 2019
Article Number A2
Number of page(s) 16
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201936150
Published online 21 November 2019

© ESO 2019

1. Introduction

More than a century ago, a bright jet in the Virgo cluster was discovered by Curtis (1918). The jet is connected to a nucleus that resides in the center of M 87, which is an elliptical galaxy. Since its discovery, the jet of M 87 has been subject to extensive radio observations (Bolton et al. 1949; Mills 1952; Baade & Minkowski 1954; Turland 1975; Owen et al. 2000; de Gasperin et al. 2012), and various millimeter observations; 7 mm (43 GHz) (Junor et al. 1999; Ly et al. 2004; Walker et al. 2008, 2018; Hada et al. 2011, 2013, 2016), 3 mm (86 GHz) (Krichbaum et al. 2006; Hada et al. 2013; Kim et al. 2018), and 1.3 mm (228 GHz) (Doeleman et al. 2012). At millimeter-wavelengths, the radio emission shows a source morphology that is consistent with a jet launched from the putative supermassive black hole in the center of the radio core with a mass of MBH  =  6.2 × 109M (Gebhardt et al. 2011) and at a distance of d  =  16.7 Mpc (Mei et al. 2007). This black hole is one of the primary targets of global millimeter very long baseline interferometry (mm-VLBI) observations by the Event Horizon Telescope Collaboration, which has the aim to spatially resolve black-hole shadows (Falcke et al. 2000; Goddi et al. 2017), and succeeded in the case of M 87 (EHT Collaboration 2019a,b,c,d,e,f). The shadow of a black hole is a depression of flux in the radiation field surrounding the black-holes event horizon, for a non-rotating black hole its size on the sky is given by 2 27 G M / ( c 2 D ) $ 2\sqrt{27}G M/(c^2 D) $, with G the Gravitational constant, M the black-hole mass, c the speed of light, and D the distance to the black hole. Due to the large set of observations available across the electromagnetic spectrum (see e.g., Prieto et al. 2016) and the event horizon scale mm-VLBI observations, it is possible to use the M 87 jet as a laboratory to study jet launching and particle acceleration.

Since the discovery of M 87, relativistic jets have been studied in great detail in theory. The analytical model by Blandford & Königl (1979) describes an isothermal jet model that can explain the observed flat radio spectra of jets. They recover the observed relation between source size (r) and frequency (ν) to be r ∝ ν−1. An addition to this model was made by Falcke & Biermann (1995), who connected the accretion rate to the jet.

Broderick & Loeb (2009) modeled M 87 with an analytic, force-free jet model. Their best-fit model is consistent with 43 GHz observations. The model parameters include a black-hole spin of a* = Jc/GM2 = 0.998, a viewing angle of i = 25°, and a jet foot-point at r = 10 rg, where the gravitational radius rg is defined as rg = GM/c2. The disk consists of both thermal and accelerated electrons, but the fraction of accelerated electrons is low (around one percent). Inside the jet, only an accelerated electron population is present. Their model uses a black-hole mass of M = 3.4 × 109M (Walsh et al. 2013).

General-relativistic magnetohydrodynamical (GRMHD) simulations are often used to study the dynamics of accretion flows. Next, we review some of the earlier GRMHD based models of the M 87 jet. The first model of M 87 based on GRMHD simulations was presented by Dexter et al. (2012), who computed synthetic synchrotron maps based on a high-spin GRMHD simulation. Their models included a thermal electron population in the disk and a power-law based electron population in the jet. Their best-fit model, at an inclination of 25°, showed counter-jet dominated emission, meaning that most of the radiation detected by the observer originates in the jet that points away from the observer. Dexter et al. (2012) obtained a mass-accretion rate of ≈ 10−3 M yr−1, and a power-law index of the non-thermal electron distribution function of p = 3.25 − 3.5, where they used a constant electron-to-proton temperature ratio of Tp/Te = 10.

Mościbrodzka et al. (2016) used GRMHD simulations and a Monte Carlo-based radiative-transfer code to model the full spectral energy distribution (SED) of an accreting supermassive black hole from radio to X-ray, as well as ray-traced images of the accretion flow at 43, 86 GHz, and 230 GHz. A thermal distribution function of the electrons was assumed across the simulation domain, and the electron physics was modeled by coupling the ion-to-electron temperature as a function of plasma β = P/Pmag, where P is the gas pressure and Pmag the magnetic pressure. The electrons were thermally distributed both in the disk and the jet. The authors obtained a mass accretion rate of ≈ 9 × 10−3 M, a favored inclination angle of 20° or 160° and an ion-to-electron temperature ratio in the disk of Ti/Te = 100. Smaller values of the ion-to-electron temperature ratio resulted in an excess of X-ray emission. The 230 GHz images showed counter-jet dominated emission. Subsequently, Mościbrodzka et al. (2017) performed polarized radiative transfer calculations of the jet launching foot point of M 87 to obtain Faraday rotation measurements. It is shown that the best-fit jet-dominated model from Mościbrodzka et al. (2016) recovers the observed 1% polarization fraction and rotation measure of the jet base in M 87.

Recently, Ryan et al. (2018) performed 2D-axisymmetric two-temperature GRMHD simulations that include radiative cooling. The authors conclude that radiative cooling is important in the inner region (r <  10 rg) of the accretion flow, and that the black-hole mass of M = 6.2 × 109M and spin a* = 0.9375 simulation recovers the observed radio and X-ray emission and image size at 230 GHz. The jet-opening angle in their model at lower frequencies is too narrow compared to the millimeter-observations of the jet base in M 87 and the model assumes a thermal electron population in the entire simulation domain. Chael et al. (2019) also performed a two-temperature radiative GRMHD model of a magnetically arrested disc (MAD; Narayan et al. 2003; Tchekhovskoy et al. 2011). The model recovers observables such as jet-opening angle, image size, core shift, and radio SED. This model also considers a thermal electron population and, therefore, does not fit the NIR and optical emission.

In 2019, the Event Horizon Telescope published its first set of results, showing an asymmetric ring-like structure in the radio core of M 87 at 228 GHz. This ring-like structure is evidence for the existence of a black hole shadow and consistent with predictions from GRMHD models (EHT Collaboration 2019a,b,c,d,e,f). A detailed comparison of GRMHD models with the data can be found in EHT Collaboration (2019e). The main assumption in these models that we want to address in this work is that the electron distribution function is taken to be thermal in the entire simulation domain.

All of the models have in common that they are based on GRMHD simulations that use spherical polar coordinates with a radial grid that is logarithmically spaced. Such a grid has the advantage of high resolution close to the event horizon but introduces a polar axis that needs careful treatment, potentially resulting in numerical issues that affect the jet outflow. GRMHD codes often track only the dynamically important ion fluid, with no direct knowledge of the electrons available. One of the open questions in modeling the electromagnetic radiation emerging from accreting black holes is, therefore, the shape of the distribution function of the radiatively important electrons. The often-made assumption that the electrons in the full simulation domain are in a thermal-relativistic Maxwell–Jüttner distribution potentially breaks down in regions where non-ideal effects are important.

These non-ideal effects are expected to be strongest in the highly magnetized regions of the jet, where they can be associated with magnetic reconnection accelerating electrons to very large energies. In the case of M 87, features of electron acceleration are observed in the NIR and optical wavebands (see e.g., Prieto et al. 2016, and references therein). We, therefore, need a distribution function that describes the electrons that are not in thermal equilibrium. Particle ensembles that are not in thermal equilibrium can be described in the framework of Tsallis statistical mechanics (Tsallis 1988). In this framework, the κ-distribution function plays a key role. In Fig. 1 we show that the κ-distribution function is a combination of a thermal core at low values of the Lorentz factor γ, which asymptotically turns into to a power-law with power-law index p = κ − 1 for large γ values. In the limit of κ → ∞, the κ-distribution becomes the Maxwell–Jüttner distribution function (Rezzolla & Zanotti 2013). The κ-distribution function is observed at a variety of astrophysical systems such as the solar wind, solar magnetosphere, Jovian magnetospheres, planetary nebula, and many others (see for a review Pierrard & Lazar 2010 and references therein).

thumbnail Fig. 1.

Example of the κ-distribution (orange) for a κ value of 3.5, and for comparison a Maxwell–Jüttner distribution (black), and a power-law distribution (yellow) with p = 2.5. The κ-distribution function is a combination of a powerlaw with a thermal core.

In Davelaar et al. (2018a), we introduced a κ-jet model for the accreting black hole in the center of the Milky Way, Sagittarius A* (Sgr A*). This model is a combination of a thermal and a κ-distributed electron population. In the accretion disk, we set the eDF to a thermal distribution function, while in the jet we use a mix of thermal and κ-distributed electrons. The ratio between the two species is a free parameter of the model. In the case of Sgr A*, we found that ∼5 − 10% of the electrons is κ-distributed in the event of flares, and they are negligible in the quiescent state. The non-thermal prescription used in this model was an electron distribution with a constant power-law in the outflow of the simulation domain with a fixed power-law index.

To improve the model we here connect the electron-acceleration parameters to information from local kinetic plasma simulations. Kinetic plasma, or particle-in-cell (PIC) simulations, are capable of resolving the micro-physics scales that GRMHD simulations cannot reach. Although local, these type of simulations can provide first-principle parametrizations of particle-acceleration processes. For our model, we consider a parametrization of the power-law index for trans-relativistic reconnection as found by Ball et al. (2018). Reconnection is known to be an efficient particle accelerator in magnetized environments (see e.g., Sironi & Spitkovsky 2014; Guo et al. 2014; Sironi et al. 2015; Werner et al. 2016, 2018; Petropoulou et al. 2016; Werner & Uzdensky 2017). Besides this parametrization, we also extended our model with an injection radius, which corresponds to the footpoint of the jet where electron acceleration can become important.

In this work, we apply thermal and κ-jet models to the accreting black hole in M 87. The dynamics of the accretion flow are drawn from GRMHD simulations performed in Cartesian–Kerr–Schild (CKS) coordinates. This prevents numerical artefacts and directional biases of the jet caused by the presence of a polar axis, this will be studied in detail a future work. In addition, the use of adaptive mesh refinement (AMR) allows us to capture the instabilities in the jet sheath, and, at the same time, to resolve the magneto-rotational instability (MRI) in the disk. We use the results of this simulation to generate SEDs, synthetic synchrotron maps (images), and optical-depth maps of the jet-launching region in M 87. We extend the general-relativistic-ray-tracing (GRRT) code RAPTOR, rendering it compatible with AMR data structures. We fit synthetic SEDs obtained from our GRRT simulations to observational data.

The plan of the paper is as follows: in Sect. 2 we describe our GRMHD simulation setup, as well as the electron model that we use in our radiative-transfer calculations. In Sect. 3 we compute SEDs, synchrotron and opacity maps, source sizes, and core shifts. In Sect. 4, we compare our results to previous works and observations. We summarize our results in Sect. 5.

2. Methods

In this section, we describe the GRMHD simulation setup, the coordinates used to simulate the accretion flow and radiation transport, and introduce our electron-physics model.

2.1. GRMHD simulations

The dynamics of the accretion flow onto the black hole are simulated using the Black Hole Accretion Code (BHAC, Porth et al. 2017), which solves the GRMHD equations

μ ( ρ u μ ) = 0 , ( 1 a ) $$ \begin{aligned} \nabla _{\mu } (\rho u^{\mu }) = 0, \qquad \qquad \mathrm{(1a)} \end{aligned} $$

μ T μ ν = 0 , ( 1 b ) $$ \begin{aligned} \nabla _{\mu } T^{{\mu \nu }} = 0, \qquad \qquad \mathrm{(1b)} \end{aligned} $$

μ F μ ν = 0 , ( 1 c ) $$ \begin{aligned} \nabla _{\mu }\, ^{*}\!F^{\mu \nu } = 0, \qquad \qquad \mathrm{(1c)} \end{aligned} $$

where ∇μ denotes the covariant derivative, ρ the rest-mass density, uμ the fluid 4-velocity, Tμν the energy-momentum tensor of the combined perfect fluid and electromagnetic fields, and *Fμν the dual of the Faraday tensor (Fαβ).

The system is closed by the ideal-MHD approximation corresponding to a plasma with infinite conductivity Fμνuν = 0, and by the equation of state of an ideal fluid (see, e.g., Rezzolla & Zanotti 2013) h ( ρ , P ) = 1 + γ ̂ γ ̂ 1 P ρ $ h(\rho,P)= 1 +\frac{\hat{\gamma}}{\hat{\gamma}-1}\frac{P}{\rho} $, where h and P are the specific enthalpy and gas pressure in the fluid frame, and the adiabatic index γ ̂ = 4 / 3 $ \hat{\gamma}=4/3 $. The simulation is initialized with a Fishbone–Moncrief torus (Fishbone & Moncrief 1976) with its inner radius at 6 rg and its pressure maximum at 12 rg. The initial configuration of the magnetic field is a single poloidal loop described by the vector potential Aϕ ∝ max(ρ/ρmax − 0.2, 0). The initial density and pressure are normalized so that ρmax = 1. The initial magnetic field is also normalized such that the ratio between maximum gas pressure Pmax and maximum magnetic pressure Pmag, max is Pmax/Pmag, max = 100. The disk is, therefore, weakly magnetized. In order to break the initial equilibrium state and accelerate the development of the MRI, we add 5% “white noise” random perturbations to the pressure. This triggers the MRI, which transports angular momentum and allows accretion onto the black hole (Balbus & Hawley 1991).

The black-hole’s dimensionless spin parameter was set to be a* = 0.9375, where J is the angular momentum. For this value of a* the inner horizon is at r ≈ 1.34799 rg.

2.2. AMR grid in Cartesian–Kerr–Schild coordinates

The simulation is performed on a Cartesian (rectangular) grid. The covariant metric gμν of a rotating black hole in Cartesian–Kerr–Schild (CKS) coordinates is given by (Kerr 1963; Rezzolla & Zanotti 2013)

g μ ν = η μ ν + f l μ l ν , $$ \begin{aligned} g_{\mu \nu } = \eta _{\mu \nu } + f l_\mu l_\nu , \end{aligned} $$(2)

where ημν = ( − 1, 1, 1, 1) is the Minkowski metric, and

f = 2 r 3 r 4 + a 2 z 2 , ( 3 a ) $$ \begin{aligned} f = \frac{2r^3}{r^4 + a_*^2 z^2},\qquad \qquad \mathrm{(3a)} \end{aligned} $$

l ν = ( 1 , r x + a y r 2 + a 2 , r y a x r 2 + a 2 , z r ) , ( 3 b ) $$ \begin{aligned} l_\nu = \left(1, \frac{rx+a_* y}{r^2 + a^2}, \frac{ry-a_* x}{r^2 + a_*^2}, \frac{z}{r}\right),\qquad \qquad \mathrm{(3b)} \end{aligned} $$

where r is given by

r 2 = R 2 a 2 + ( R 2 a 2 ) 2 + 4 a 2 z 2 2 , $$ \begin{aligned} r^2 = {\frac{R^2 - a_*^2 + \sqrt{(R^2 - a^2)^2 + 4a_*^2z^2}}{2} }, \end{aligned} $$(4)

and

R 2 = x 2 + y 2 + z 2 . $$ \begin{aligned} R^2 = x^2+y^2+z^2. \end{aligned} $$(5)

All units of length are scaled by the gravitational radius rg which is given by rg = GM/c2. In the limit of R ≫ a*, the radius r → R. The contravariant metric is given by

g μ ν = η μ ν f l μ l ν , $$ \begin{aligned} g^{\mu \nu } = \eta ^{\mu \nu } - f l^\mu l^\nu , \end{aligned} $$(6)

where lν is given by

l ν = ( 1 , r x + a y r 2 + a 2 , r y a x r 2 + a 2 , z r ) · $$ \begin{aligned} l^\nu = \left(-1, \frac{rx+ay}{r^2 + a_*^2}, \frac{ry-a_*x}{r^2 + a_*^2}, \frac{z}{r}\right)\cdot \end{aligned} $$(7)

The use of AMR allows us to increase the resolution in regions of interest during runtime. The decision to refine is made based on the Löhner scheme (Löhner 1987), which quantifies variations of the density and the plasma magnetization σ, defined as σ = b2/ρ, where b = b μ b μ $ b=\sqrt{b^{\mu}b_{\mu}} $ is the magnetic-field strength in the fluid frame. The code is allowed to refine up to a maximum level of refinement that depends on the location in the computational domain; greater levels of refinement are allowed in the regions where the jet is expected to form and the disk is expected to reside. This distinction is made based on a radius r and polar angle θ, and for the jet this region is between θ <  15° or θ >  165°. The maximum allowed refinement level as a function of radius and polar angle are shown in Table 1. The base resolution of the grid is 96 × 96 × 192 cells in x, y, and z-directions, respectively. The simulation domain is x ∈ ( − 500 rg, 500 rg), y ∈ ( − 500 rg, 500 rg) and z ∈ ( − 1000 rg, 1000 rg). We simulate up to tf = 104rg/c, which corresponds to 37.5 orbital periods of the accretion disk at the pressure maximum. At the end of the simulation, the domain contains around 70 million cells.

Table 1.

Maximum AMR refinement radii in rg for the different AMR levels.

2.3. Ray tracing in AMR CKS grid

In order to perform general-relativistic ray-tracing calculations in Cartesian coordinates within the block-based AMR data structure of BHAC, it has been necessary to extend our general-relativistic ray-tracing code RAPTOR (Bronzwaer et al. 2018). In particular, the initial conditions for the rays, also called the “virtual camera”, employ a tetrad basis in which the initial wave-vectors are described (Noble et al. 2007), a description of the implementation of this in RAPTOR can be found in Davelaar et al. (2018b). The tetrad camera uses a set of trial vectors to generate a tetrad basis by using a Gramm-Schmidt orthogonalization procedure. In spherical coordinate systems, the trial vectors are unit vectors pointing along the t, r, θ, ϕ-directions. In our case, we have to transform this into Cartesian coordinates. The virtual camera is constructed at a position (xc, yc, zc) in space which is computed based on the following parameters: (i) the radial distance between the camera and the black hole rc; (ii) the inclination with respect to the black hole spin axis θc; (iii) the azimuthal angle around the black hole spin axis ϕc. The tetrad trial vectors can then be defined as

t 0 μ = ( 1 , 0 , 0 , 0 ) , ( 8 a ) $$ \begin{aligned} t^\mu _{0} = (1, 0, 0, 0), \qquad \qquad \mathrm{(8a)} \end{aligned} $$

t 1 μ = ( 0 , sin ( θ c ) cos ( ϕ c ) , sin ( θ c ) sin ( ϕ c ) , cos ( θ c ) ) , ( 8 b ) $$ \begin{aligned} t^\mu _{1} = (0, -\sin (\theta _{\rm c})\cos (\phi _{\rm c}), -\sin (\theta _{\rm c})\sin (\phi _{\rm c}),-\cos (\theta _{\rm c})),\qquad \qquad \mathrm{(8b)} \end{aligned} $$

t 2 μ = ( 0 , sin ( ϕ c ) , cos ( ϕ c ) ) , ( 8 c ) $$ \begin{aligned} t^\mu _{2} = (0, -\sin (\phi _{\rm c}),\cos (\phi _{\rm c})),\qquad \qquad \mathrm{(8c)} \end{aligned} $$

t 3 μ = ( 0 , cos ( θ c ) cos ( ϕ c ) , cos ( θ c ) sin ( ϕ c ) , sin ( θ c ) ) . ( 8 d ) $$ \begin{aligned} t^\mu _{3} = (0, -\cos (\theta _{\rm c})\cos (\phi _{\rm c}), -\cos (\theta _{\rm c})\sin (\phi _{\rm c}),-\sin (\theta _{\rm c})).\qquad \qquad \mathrm{(8d)} \end{aligned} $$

The choice of trial vectors results in a right-handed basis where the observer is facing the black hole.

The integration of the geodesic equations is done by solving the second-order differential equation

d 2 x α d λ 2 = Γ μ ν α d x μ d λ d x ν d λ · $$ \begin{aligned} \frac{\mathrm{d}^2 x^{\alpha }}{\mathrm{d}\lambda ^2} = -\Gamma ^{\alpha }_{\ \mu \nu } \frac{\mathrm{d}x^{\mu }}{\mathrm{d}\lambda } \frac{\mathrm{d}x^{\nu }}{\mathrm{d}\lambda }\cdot \end{aligned} $$(9)

where Γ μ ν α $ \Gamma^{\alpha}_{\ \mu \nu} $ are the connection coefficients, xα is the geodesic position, and λ is the affine parameter. We use a fourth-order Runge–Kutta algorithm, where the connection coefficients are evaluated using a finite-difference derivative of the metric.

The step-sizing for the geodesic integration in RAPTOR was adopted since it relies on spherical logarithmic coordinates. First we compute a required step-size based on the geodesic wave-vector

d λ x = Δ / ( | k x | + δ ) , ( 10 a ) $$ \begin{aligned} \mathrm{d} \lambda _{x} = \Delta \ / \left( \left| k^x \right| + \delta \right), \qquad \qquad \mathrm{(10a)} \end{aligned} $$

d λ y = Δ / ( | k y | + δ ) , ( 10 b ) $$ \begin{aligned} \mathrm{d} \lambda _{y} = \Delta \ / \left( \left| k^y \right| + \delta \right), \qquad \qquad \mathrm{(10b)} \end{aligned} $$

d λ z = Δ / ( | k z | + δ ) , ( 10 c ) $$ \begin{aligned} \mathrm{d} \lambda _{z} = \Delta \ / \left( \left| k^z \right| + \delta \right), \qquad \qquad \mathrm{(10c)} \end{aligned} $$

d λ geod = R | d λ x | 1 + | d λ y | 1 + | d λ z | 1 , ( 10 d ) $$ \begin{aligned} \mathrm{d} \lambda _{\rm geod}=\frac{R}{{\left| \mathrm{d}\lambda _{x} \right|}^{-1} + {\left| \mathrm{d}\lambda _{y}\right|}^{-1} + {\left| \mathrm{d}\lambda _{z} \right|}^{-1}},\qquad \qquad \mathrm{(10d)} \end{aligned} $$

where kx, y, z are the wave-vector components in the x, y, z directions, δ is a small real number to prevent divisions by zero, and Δ is a scale factor for the step-size (typically Δ ≈ 0.01). Then we compute a required step-size based on the AMR cell size dx

k max = max ( k x , max ( k y , k z ) ) , ( 11 a ) $$ \begin{aligned} k_{\rm max} = \max (k^x,\max (k^y,k^z)),\qquad \qquad \mathrm{(11a)} \end{aligned} $$

d λ grid = d x n k max , ( 11 b ) $$ \begin{aligned} \mathrm{d} \lambda _{\rm grid} = \frac{\mathrm{d}x}{n k_{\rm max}},\qquad \qquad \mathrm{(11b)} \end{aligned} $$

where n sets the amount of steps per cell. We typically use at least two steps per cell. We then compare both the geodesic and AMR based step-sizes and use the smallest of the two to ensure convergence; dλ = min(dλgeod, dλgrid).

For the radiative-transfer part of the ray-tracing calculation, we need the plasma variables at the location of the geodesics. We interfaced RAPTOR with the AMR data structure of BHAC, and reconstruct the full AMR grid. The BHAC AMR block-based data structure is parsed by the code. When we integrate the geodesics we use a nearest-neighbor approach to interpolate the grid-based plasma variables to the geodesics.

2.4. Electron model and radiative-transfer model parameters

Since GRMHD simulations are scale-free, we have to re-scale the plasma variables from code units to c.g.s. units. Units of length are scaled with ℒ = rg, while units of time are scaled with 𝒯 = rg/c, the mass unit is set by ℳ. Estimates of the mass of M 87 are used to constrain the length and time units, we use a mass of M = 6.2 × 109M (Gebhardt et al. 2011), the mass used in this work is slightly smaller than the mass of M = (6.5 ± 0.7)×109M reported in EHT Collaboration (2019a), but the used value for the black hole mass is within the ermror margins. The mass unit ℳ, which sets the accretion rate, however, is unknown. It is, therefore, a fit parameter. The mass unit is directly proportional to the accretion rate via cgs = simℳ𝒯−1, where sim is the accretion rate in simulation units. In order to scale the relevant plasma quantities to c.g.s units, the following scaling operations are performed: ρ0 = ℳ/ℒ3, u0 = ρ0c2, and B 0 = c 4 π ρ 0 $ B_0=c\sqrt{4\pi \rho_0} $.

As mentioned before, our GRMHD simulation only simulates the dynamically important protons. Therefore, we need to parametrize the electron properties, such as their distribution functions, densities, and temperatures, in post-processing. The plasma is assumed to be charge-neutral, so that ne = np throughout the domain. For the electron temperature we employ the parametrization of Mościbrodzka et al. (2016):

T ratio = T p / T e = R low 1 1 + β 2 + R high β 2 1 + β 2 , ( 12 a ) $$ \begin{aligned} T_{\rm ratio}= T_{\rm p}/T_{\rm e} = R_{\rm low}\frac{1}{1+\beta ^2} + R_{\rm high} \frac{\beta ^2}{1+\beta ^2}, \qquad \qquad \mathrm{(12a)} \end{aligned} $$

Θ e = U ( γ ̂ 1 ) m p / m e ρ T ratio , ( 12 b ) $$ \begin{aligned} \Theta _{\rm e} = \frac{U(\hat{\gamma } - 1) m_{\rm p}/m_{\rm e}}{\rho T_{\rm ratio}},\qquad \qquad \mathrm{(12b)} \end{aligned} $$

where mp is the proton mass, me is the electron mass, U is the internal energy, which is related to h and ρ via h ( ρ , U ) = 1 + γ ̂ U ρ $ h(\rho,U)= 1 +{\hat{\gamma}}\frac{U}{\rho} $, therefore, P = U(γ − 1), where P is the gass pressure. The dimensionless electron temperature Θe can be re-scaled to c.g.s units via T = Θemec2/kb, where kb is the Boltzmann constant. The parameters Rlow and Rhigh are free parameters of the model; Rlow sets the temperature ratio in the jet, where β ≪ 1, and Rhigh sets the temperature ratio in the disk where β ≫ 1 (Fig. 2).

thumbnail Fig. 2.

Parameterization of the κ parameter. κ as function of β and σ as found by Ball et al. (2018). A high value of κ corresponds to steep particle spectra with power-law index p = κ − 1. We overplotted contours of constant κ in black.

For the electrons’ energy-distribution function, we follow a similar recipe as described in Davelaar et al. (2018a). we use the relativistic isotropic κ-distribution function for the electrons, which is given by (Xiao 2006)

d n e d γ = N γ γ 2 1 ( 1 + γ 1 κ w ) ( κ + 1 ) , $$ \begin{aligned} \frac{dn_{\rm e}}{d\gamma } = N \gamma \sqrt{\gamma ^2 -1} \left( 1 + \frac{\gamma -1}{\kappa w}\right)^{-(\kappa + 1)}, \end{aligned} $$(13)

where γ is the Lorentz factor of the electrons, κ is the parameter that sets the power-law index p via p = κ − 1, w sets the width of the distribution function, and N is a normalization factor such that the electron distribution function contains ne electrons.

The width w of the κ distribution sets the amount of energy in the distribution. In the case that κw ≫ 1 the total energy in the κ distribution is given by

E κ = 3 κ κ 3 n e w . $$ \begin{aligned} E_{\kappa } = \frac{3 \kappa }{\kappa -3} n_{\rm e} w. \end{aligned} $$(14)

We couple this energy to the energy present in a thermal distribution (Ethermal = 3neΘe) and add a source term based on the magnetic energy

E κ = 3 κ κ 3 n e w = 3 n e Θ e + ϵ B 2 8 π , $$ \begin{aligned} E_{\kappa } = \frac{3 \kappa }{\kappa -3} n_{\rm e} w = 3 n_{\rm e} \Theta _{\rm e} + \tilde{\epsilon } \frac{B^2}{8\pi }, \end{aligned} $$(15)

here ϵ $ \tilde{\epsilon} $ is used to join smoothly between between the κ-distribution and the magnetic energy. After a bit of algebra, we can rewrite the width as

w = κ 3 κ Θ e + ϵ κ 3 6 κ m p m e σ . $$ \begin{aligned} w = \frac{ \kappa -3 }{\kappa } \Theta _{\rm e} + \tilde{\epsilon } \frac{ \kappa -3 }{6 \kappa } \frac{m_{\rm p}}{m_{\rm e}} \sigma . \end{aligned} $$(16)

In the limit of σ ≪ 1, the κ-distribution energy is set by the thermal energy, while in the magnetized regime the energy is set by the magnetic energy. The ϵ $ \tilde{\epsilon} $ parameter is set by

ϵ = ϵ 1 2 ( 1 + tanh ( r r inj ) ) , $$ \begin{aligned} \tilde{\epsilon } = \epsilon \frac{1}{2}\left(1 + \tanh ( r - r_{\rm inj} )\right), \end{aligned} $$(17)

where rinj is the injection radius from which the magnetic energy contributes to the w parameter, and ϵ is the base value for radii larger than rinj; hereafter, we will consider two cases: where ϵ is zero or non-zero. Here, ϵ is a parameter expressing the equipartition between the energy in the distribution function and that available in magnetic energy (see, e.g., Falcke & Biermann 1995; Markoff et al. 2008; Prieto et al. 2016).

The power-law index of the electrons distribution functions (eDFs) is based on sub-grid particle-in-cell (PIC) simulations of trans-relativistic reconnection by Ball et al. (2018), who simulated two-dimensional reconnection layer (Harris sheath) for an electron-ion plasma for multiple values of the plasma β and of the magnetization σ. One of the benefits of this type of plasma simulation is that one obtains eDFs from first principles. In Ball et al. (2018) these eDFs are then used to fit the power-law index p as a function of β and σ as

p = A p + B p tanh ( C p β ) ( 18 a ) $$ \begin{aligned} p = A_p + B_p \tanh \left(C_p \beta \right) \qquad \qquad \mathrm{(18a)} \end{aligned} $$

A p = 1.8 + 0.7 / σ ( 18 b ) $$ \begin{aligned} A_p = 1.8 + 0.7/ \sqrt{\sigma }\qquad \qquad \mathrm{(18b)} \end{aligned} $$

B p = 3.7 σ 0.19 ( 18 c ) $$ \begin{aligned} B_p = 3.7\, \sigma ^{-0.19}\qquad \qquad \mathrm{(18c)} \end{aligned} $$

C p = 23.4 σ 0.26 . ( 18 d ) $$ \begin{aligned} C_p = 23.4\, \sigma ^{0.26}.\qquad \qquad \mathrm{(18d)} \end{aligned} $$

These fits are obtained for 10−4 <  β <  1.5 and 0.1 <  σ <  7.2, which corresponds to the typical values that we find in the jet sheath, which is the main source of synchrotron emission in our jet-models.

Finally, we exclude all emission from regions where σ >  5.0, this is what we call the jet spine. These regions are unreliable for modeling since the thermal energy in highly magnetized regions is unreliable in GRMHD simulations. We also exclude all emission from regions where floor values are applied, these regions typically resides inside the magnetized jet. This results in three regions inside our simulation domain; the disk, the jet sheath, and jet spine. The disk resides where σ is much smaller than one and plasma β is large than one, the jet sheath resides where σ is of the order unity and β is smaller than one. In the case of our κ-jet model we set the electron distribution function to a relativistic κ-distribution function into the disk and jet sheath and no electrons are present in the jet spine.

The emission and absorption coefficients for the thermal electron distributions are taken from Leung et al. (2011), and in the case of the κ-distribution, the fit formula taken are from Pandya et al. (2016).

2.5. SED cut-off

The SED of M 87 shows a clear cut-off in flux around ν = 1015 Hz (Prieto et al. 2016). We will consider three potential sources for this cut-off.

First, we assume that the cut-off is caused by synchrotron cooling in the jet, which becomes important when the synchrotron-cooling time of the electron is comparable with the typical dynamical time. Under these conditions, the cooling (cut-off) frequency is given by

ν cool = 18 π σ T 2 m e c 2 e B 3 z jet , $$ \begin{aligned} \nu _{\rm cool} = \frac{18 \pi }{\sigma _T^2} \frac{m_e c^2 e }{B^3 z_{\rm jet}}, \end{aligned} $$(19)

where σT is the Thomson cross-section, and zjet the position along the jet.

Second, we assume that the break takes place at the synchrotron burn-off limit, that is, at the maximum energy that a particle can gain while emitting synchrotron radiation. The maximum Lorentz factor in this case is

γ max = 3 m e 2 c 4 E 4 π e 3 B 2 , $$ \begin{aligned} \gamma _{\rm max} = \sqrt{\frac{3m_{\rm e}^2 c^4 E}{4\pi e^3 B^2}}, \end{aligned} $$(20)

where E is the electric field, and the cut-off frequency is then given by

ν cut-off = 3 2 γ max 2 ν c , $$ \begin{aligned} \nu _{\text{ cut-off}} = \frac{3}{2} \gamma _{\rm max}^2 \nu _c, \end{aligned} $$(21)

with νc = eB/(2πmec).

Finally, we assume that break is given by the Hillas criterion (Hillas 1984), stating that the maximum Lorentz factor achievable can be estimated by equating the gyration radius of the electron and the size of the acceleration region L. This results in a maximum Lorentz factor of

γ max = eBL m e c 2 , $$ \begin{aligned} \gamma _{\rm max} = \frac{eBL}{m_{\rm e} c^2}, \end{aligned} $$(22)

which results in a cut-off frequency of ν ≈ 1015 Hz after using Eq. (22) in (21). In this way, we can also we can estimate the typical size L of the acceleration-region

L = 4 π ν cut off m e 3 c 5 e 3 B 3 4.5 × 10 7 cm ( ν cut off / 10 15 Hz ) ( B / 1 G ) 3 · $$ \begin{aligned} L = \sqrt{\frac{4\pi \nu _{\rm cut-off} m_{\rm e}^3 c^5}{e^3 B^3}} \approx 4.5 \times 10^{7} \mathrm{\,cm} \sqrt{\frac{(\nu _{\rm cut-off}/10^{15}\,\mathrm{Hz})}{{(B/1~\mathrm{G})^3}}}\cdot \end{aligned} $$(23)

Interestingly, the maximum size L can be interpreted as the size of plasmoids as was done by Petropoulou et al. (2016) and Christie et al. (2019) for blazars.

3. Results

In this section, we present the results of our GRMHD simulations and how the SEDs they produce can be compared with the observational data. at three observational relevant frequencies at two inclinations. We also show how we can compute from the synthetic images the source size and core shifts, and how they compare with the observations.

3.1. Structure of the accretion disk and jet in the AMR simulation

A representative snapshot of the simulation is shown in Fig. 3. The simulation produces a well-resolved relativistic jet up to the edge of the simulation domain at 1000 rg in the z-direction. At z = 40 rg the jet diameter is resolved with 160 cells, and with 32 cells at z = 1000 rg. The accretion rate through the event horizon is shown in the left panel of Fig. 4 as a function of time; the accretion rate initially increases sharply, and then settles around 0.2 (in code units) at later times. The jet and wind power are instead shown in the right panel of Fig. 4; both of the quantities are calculated by performing the integral over the constant r = 100 rg surface

E ˙ = 0 2 π 0 π ( T t r ρ u r ) χ ( · ) g d θ d ϕ $$ \begin{aligned} \dot{E} = \int _0^{2\pi }\int _{0}^{\pi } (- T^r_t - \rho u^r)\, \chi _{(\cdot )}\, \sqrt{-g}\,\mathrm{d}\theta \,\mathrm{d}\phi \end{aligned} $$(24)

thumbnail Fig. 3.

GRMHD simulation snapshot. Left panels: slice in the (x, z) and (x, y) planes of the density in code units. Middle panels: slice along the (x, z) and (x, y) planes of the magnetization parameter b2/ρ, over-plotted with the grid block sizes. Right panels: slice along the (x, z) and (x, y) planes of the dimensionless ion temperature. Shown with a black circle is the location of the event horizon.

thumbnail Fig. 4.

Mass accretion rate and jet and wind power. We computed global quantities of the Cartesian GRMHD simulation, accretion rate in code units as a function of time (left panel), and jet and wind power in code units as a function of time (right panel).

where the function χ( ⋅ ) selects only material in the jet, wind, or disk following the setting

χ jet = ( b 2 / ρ > 1 or μ > 2 ) ( 25 a ) $$ \begin{aligned} \chi _{\rm jet} = (b^2/\rho >1~\mathrm{or}~ \mu > 2) \qquad \qquad \mathrm{(25a)} \end{aligned} $$

χ wind = ( not χ jet and h u t > 1 ) ( 25 b ) $$ \begin{aligned} \chi _{\rm wind} = (\mathrm{not}~\chi _{\rm jet}~\mathrm{and}~ -h u_t > 1)\qquad \qquad \mathrm{(25b)} \end{aligned} $$

χ disk = ( not χ jet and not χ wind ) , ( 25 c ) $$ \begin{aligned} \chi _{\rm disk} = (\mathrm{not}~\chi _{\rm jet}~\mathrm{and}~\mathrm{not}~ \chi _{\rm wind}),\qquad \qquad \mathrm{(25c)} \end{aligned} $$

and μ denotes the energy flux normalized to the rest-mass energy in the radial direction μ=( T t r ρ u r )/(ρ u r ) $ \mu=(-T^r_t-\rho u^r)/(\rho u^r) $. Hence, the jet is defined as the region which reaches asymptotic Lorentz factors of at least 2. The optional condition b2/ρ >  1 also selects the flow in the inner axial region, where the Poynting flux necessarily vanishes. The disk wind is then the remaining unbound material and the disk itself is composed of the bound material.

Analytic work on radial profiles of relativistic jets was performed by Blandford & Königl (1979) and subsequently by Falcke & Biermann (1995). In these Blandford–Königl jet-models, the electron density decreases as a function of radius as ρe ∝ r−2, the magnetic field strength as B ∝ r−1, and the equipartition electron temperature in the jet is constant. The temperature in the disk is set by the virial theorem, and follows Te ∝ r−1. To compare our simulations with these analytical formulae, we compute averages on spherical shells at different fixed radii of the electron density ρe, magnetic-field strength B, and electron temperature Te. This is done by performing the following integral

q ( r ) = 1 Δ t ( 0 2 π q ( t , r , θ , ϕ ) γ ( r , θ ) d θ d ϕ 0 2 π g ( r , θ ) d θ d ϕ ) d t . $$ \begin{aligned} q(r) = \frac{1}{\Delta t} \int \left(\frac{\int \int _0^{2\pi } q(t,r,\theta ,\phi ) \sqrt{-\gamma (r,\theta )}\,\mathrm{d}\theta \,\mathrm{d}\phi }{\int \int _0^{2\pi } \sqrt{-g(r,\theta )}\,\mathrm{d}\theta \,\mathrm{d}\phi } \right)\,\mathrm{d}t. \end{aligned} $$(26)

where γ(r, θ) is the determinant of the three metric The integral in the θ-direction depends on the local plasma criteria. We consider two regions of interest; a jet sheath, for which 0.1 <  σ <  5.0, and the accretion disk, σ <  0.1. The time average, on the other hand, is performed using snapshots of the simulation between t = 5000 rg/c and t = 104rg/c, with a total of hundred snapshots. The computed radial profiles are shown in Fig. 5 and are over-plotted with the analytic predictions (Blandford & Königl 1979; Falcke & Biermann 1995). The equipartition electron temperature in the jet (right panel) shows a flat profile up to 200 rg, followed by is an increase of temperature that correlates with the break in the profile of the electron density. The break is caused by de-collimation of the jet, whose origin could be due to the limited initial size of the torus. The wind emitted by the disk effectively acts as a collimation agent; however, because of its limited size, the collimation stalls at radii r >  200 rg.

thumbnail Fig. 5.

Lineprofiles. To compare our simulation with analytical prediction for plasma quantities by Blandford & Königl (1979) and Falcke & Biermann (1995) in the jet we computed radial profiles. We show radial profiles of the dimensionless electron density (left), magnetic-field strength (middle) and electron temperature (right). Black lines correspond to jet averaged quantities, dashed yellow lines to disk-averaged quantities and the red dashed lines correspond to power-law profiles predicted in analytical works. Also, the jet sheath is close to isothermality.

We computed maps of κ [based on Eq. (18)] and w (based on Eq. (16), for both ϵ = 0 and ϵ = 0.015, these are shown in Fig. 6. The maps show that κ is low in the sheath, and w peaks in the sheath, therefore, most of the non-thermal emission will be produced in this region.

thumbnail Fig. 6.

κ and w parameter. Left panel: map of the κ parameter based on Eq. (18), where the jet spine is excluded (σ <  5). The map shows that the κ parameter is small in the jet sheath and large in the disk. Middle panel: map of the w parameter when ϵ = 0, showing large values in the sheath compared to the disk. Right panel: same as middle panel but now for ϵ = 0.015, the sheath has slightly larger values compared to the ϵ = 0 case.

3.2. Spectra and synchrotron images: dependency on electron distribution function

In this section, we discuss the spectral energy distributions (SEDs) of our thermal-jet and κ-jet models. The SEDs are calculated at an inclination of i = 160°, which ensures that the emitting region is in the south, as suggested by the EHT results (EHT Collaboration 2019a,e). Furthermore, the field-of-view of the camera is set to be 1000 rg in both the x and y-directions, while the resolution is set to be 2000 × 2000 pixels.

3.2.1. Fitting the SED

After averaging in time the SEDs from our models between t = 5000 rg/c and t = 104rg/c), these have been fit to non-simultaneous observations by Doeleman et al. (2012), Akiyama et al. (2015), Prieto et al. (2016), Walker et al. (2018), and Kim et al. (2018). The fit parameters are shown in Table 2, which highlights that the thermal-jet and κ-jet models differ in the accretion rate by a factor ≈2. The corresponding SEDs are shown in Fig. 7, which shows that κ-jet models recovers well the NIR flux. In particular, when comparing the ϵ = 0.0 and the ϵ = 0.015 models (the latter uses an injection radius of rinj = 10 rg and has a slightly lower accretion rate), it is possible to appreciate that the ϵ = 0.015 model has a larger and flatter radio spectrum at frequencies below ν = 228 GHz.

Table 2.

Model parameters.

thumbnail Fig. 7.

Spectral energy distributions. Shown are the thermal-jet (orange) and κ-jet with ϵ = 0 (orange) and ϵ = 0.015 (black) with their corresponding rms, overplotted with observational data points by Doeleman et al. (2012), Akiyama et al. (2015), Prieto et al. (2016), and Kim et al. (2018).

After 228 GHz both κ-jet models recover a power-law with an index of α ≈ −0.7, where α = −(p − 1)/2 for a power law distribution of non-thermal electrons Fν ∝ να. Furthermore, when compared to the thermal model, the flux in the κ-jet models is higher at lower frequencies (ν <  1011 Hz) and at the higher frequencies (5 × 1012 Hz <  ν).

When considering the various cut-off models, the cooling cut-off turned out to be unimportant, in agreement with the findings of Mościbrodzka et al. (2016), Broderick et al. (2015). When using the synchrotron burn-off, the correct cut-off is obtained if E/B ≈ 10−6, but no physical model is possible that recovers such a ratio. The only criterion that recovers the cut-off frequency is the Hillas criterion, which is obtained when the plasmoid size is set to L ≈ 105 − 107 cm, depending on the local magnetic field strength.

3.2.2. Synchrotron maps

The synthetic synchrotron maps are computed at three frequencies: 43, 86, and 228 GHz. The same inclination used for the SEDs is employed here and the images for the thermal case are shown in the top rows of Fig. 8, with the κ-jet models shown in the second and third rows. The maps shown are computed with a single GRMHD snapshot at t = 104rg/c. The forward jet at 43 GHz is aligned with the observed jet position angle at 43 GHz VLBI observations (Janssen et al. 2019), namely, 250°. The assumed mass and distance are MBH = 6.2 × 109M (Gebhardt et al. 2011) and d = 16.7 Mpc (Mei et al. 2007), which results in a field of view of: 0.744, 0.372 and 0.186 mas for the 43, 86, and 228 GHz maps, respectively.

thumbnail Fig. 8.

Synthetic synchrotron maps. Shown are all our models at three astronomically relevant frequencies. From left to right: 43, 86, and 228 GHz. Top row: synthetic images at a single snapshot of the thermal-jet at an inclination of i = 160°. Second row: same as top row but for the ϵ = 0.0 κ-jet. Bottom row: same as the first and second row but for the ϵ = 0.015 κ-jet.

The thermal-jet and ϵ = 0.0 κ-jet model show a similar source morphology at 43 GHz and 86 GHz, and ϵ = 0.015 κ-jet model is more extended in jet length. At 228 GHz both κ-jet models deviate from the thermal-jet model, the width of the ring around the shadow decreases when particle acceleration is present. In all 228 GHz images two rings are visible, the outer ring is the photon ring and marks the shadow of the black hole, the fainter smaller ring is emission originating from the jet facing the observer, see Appendix A for more details.

The logarithmic optical-depth maps at 228 GHz are shown in Fig. 9, where the size of the optically thick region (in blue) decreases when particle acceleration is present. This is in agreement with the less extended structure visible in the intensity-maps of Fig. 8. The reason behind this behavior is that lower mass-accretion rates decrease both the density and the magnetic field strength, hence decreasing the optical thickness of the jet base. As a result, for any given frequency, accelerated particles at lower mass-accretion rates contribute more than their thermal counterpart.

thumbnail Fig. 9.

Logarithmic optical-depth maps. Single snapshot at 228 GHz of the models at an inclination of i = 160°.

3.2.3. Origin of the jet emission

To obtain a quantitative understanding of how much flux originates either from the forward or the counter-jet, the jet facing towards or away from the observer, we computed synthetic images where the emission coefficient was set to zero either in the southern or northern hemisphere, while keeping the absorption coefficients in place. We computed the time-averaged ratios and the spread of the flux originating from the southern to flux from the northern hemisphere of both our models for all slices between 5000 and 104rg/c at 43, 86, and 228 GHz and have reported them in Table 3. When electron acceleration is present, the overall trend is that at 43 and 86 GHz, the ratio shifts to the counter-jet, while at 228 GHz no large shifts are seen. We therefore conclude that the counter jet at 43 and 86 GHz is more dominant in the κ-jet models compared to the thermal models. Overall, the forward jet is dominant at 43 GHz and 86 GHz, while at 230 GHz the counter jet dominates, this in agreement with earlier findings by Mościbrodzka et al. (2016). Appendix A provides a simple phenomenological model that is capable of reproducing this effect, where it is caused by a combination of gravitational lensing and the blocking of light by the black-hole’s event horizon.

Table 3.

Forward jet to counter jet ratio.

3.2.4. Core size and shift

We computed the source size of our models at 43, 86, and 228 GHz by using image moments (Hu 1962). The sources sizes are computed over a range of 5000–104rg/c and in Table 4 we report the time-averaged major and minor full-width half maxima (FWHM) and their corresponding spread.

Table 4.

Source size.

We computed the core shift with respect to the black-hole’s center at the following observational frequencies; 2.3, 5, 8.4, 15.4, 23.8, 43, 86, and 228 GHz. The core shift was calculated by computing the first-order image moments of time-averaged images and the comparison of the values obtained with the observational fit of Hada et al. (2011), meaning rRA(ν) = (1.4 ± 0.16)ν−0.94 ± 0.09, is shown in Fig. 10. The observed core shift is in agreement with the analytical predictions for which the core position should shift for a conical jet as a function of frequency as rcore  ∝  ν−1 (Blandford & Königl 1979; Falcke & Biermann 1995; Davelaar et al. 2018a), and in agreement with simulations of collimated jets (Porth et al. 2011). The κ-jet models show smaller core shifts with respect to the thermal-jet model, probably because the counter-jet is more dominant.

thumbnail Fig. 10.

Core shift. RA offset from the 43 GHz core as a function of frequency. Orange triangles correspond to a thermal-jet, black dots to a κ-jet, grey line represent the observational fit rRA(ν) = (1.4 ± 0.16)ν−0.94 ± 0.09 to the M 87 core by Hada et al. (2011).

3.2.5. Comparison with 43 GHz data

Finally, we compared our thermal-jet and κ-jet models with the 43 GHz VLBI observations, where M 87 was tracked for 8 hours with all VLBA stations1. The data was recorded with a bandwidth of 256 MHz, with the calibration and imaging of the data having been described by Janssen et al. (2019).

To compare with this observational data, we re-computed synthetic images with a large field of view of 3.7 mas and convolved them with a 0.3 × 0.1 mas2 beamsize by using the eht-imaging library (Chael et al. 2016, 2018a). The result of the comparison can be seen in Fig. 11 and highlights that the κ-jet models show more extended structure with respect to the thermal-jet model. At 43 GHz all models deviate from the VLBI observations at larger scales. Furthermore, in the observed image the flux levels upstream of the jet are higher and the jet-opening angle is wider.

thumbnail Fig. 11.

43 GHz radio map. Comparison of our thermal and κ-jet model with a VLBI map of M 87 at 43 GHz. First panel from the left: 43 GHz radio map of M 87 (Janssen et al. 2019). Second panel: synchrotron map of the thermal-jet model, convolved with a 2D Gaussian beam. Third panel: same as the second panel but now for a κ-jet model with ϵ = 0.0. Fourth panel: same as the third panel but now with ϵ = 0.015. The white ellipse indicates the beam used to convolve the images. All models produce a jet that is too narrow compared to the VLBI map. The extent of the jet increases when electron acceleration is present, and is maximum for ϵ = 0.015.

4. Discussion

4.1. CKS GRMHD simulations

Current models of the radio emission near the supermassive black hole in M 87 are based on GRMHD simulations using spherical polar coordinates. In this work, we used instead Cartesian coordinates, which do not require specialized treatment of the polar axis which represents a coordinate singularity (of the inverse metric) in spherical coordinates. The addition of AMR resulted in a highly resolved jet region, the jet diameter at z = 40 rg is resolved by 160 cells and at z = 1000 rg is resolved by 32 cells. The obtained jet resolution is well above the values of 20–26 cells per jet radius reported in convergence studies of jets by Anjiri et al. (2014) and Massaglia et al. (2016). We computed the mass-accretion rate and radial profiles of density, magnetic-field strength and temperature which are consistent with their spherical counterparts, the comparison of this can be found in Porth et al. (2019). The downside of the Cartesian grid is that in spherical grids it is possible to use a logarithmic grid in the radial direction, which results in higher resolutions close to the horizon. To ensure that the MRI is well resolved, we have computed the relativistic MRI quality factors (Noble et al. 2010) finding they are normally above ten in the bulk of the accretion disk, thus satisfying the requirements for a sufficiently resolved MRI found by Sano et al. (2004). Furthermore, in parallel works (Olivares et al. 2019; Porth et al. 2019) it was shown that the Cartesian CKS simulations show a behavior in the nonlinear phase that is similar to that of the spherical simulations. A short summary of the relevant quantifications reported in these works can be found in Appendix B.

4.2. The effect of electron acceleration on the SED

We computed spectra for our thermal-jet model and κ-jet models, wherein the latter case we parametrize the power-law index of the eDF based on sub-grid PIC simulations by Ball et al. (2018). The addition of accelerated electrons in the jet sheath leads to a fit to the observational data from radio up to the NIR. Our κ-jet model is an extension of the model presented by Mościbrodzka et al. (2016), which only studied the thermal-jet case. Our models use their best-fit inclination angle of i = 160°, such that the emitting region is in the south and the orientation of the asymmetry is in agreement with the image in EHT Collaboration (2019a). The radio SED shows a flat spectrum in both the thermal and κ-jet models. This is consistent with more recent work by Ryan et al. (2018) and Chael et al. (2019), who have evolved the thermal electron population as a separate fluid in the GRMHD simulation.

In contrast to previous works, our κ-jet models also recover the observed NIR flux by extending the optically thin emission with a power-law. The results are similar to the ones presented by Dexter et al. (2012), who also modeled accelerated electrons based on the amount of available magnetic energy. The κ-jet models yield a jet power of the order of 1043 ergs s−1, which is in agreement with observations of the jet core power by Reynolds et al. (1996), and is approximately two times lower than the thermal-jet models. This is probably due to the fact that in the κ-jet models there is a larger contribution of electrons in the tail of the distribution functions with respect to the thermal-models. Since these electrons emit at higher γ values, this results in a higher flux contribution per unit mass.

After defining the radiative efficiency as ϵrad = L/c2, we found that the thermal-jet has ϵrad = 0.003. This is to be contrasted with ϵrad = 0.013 and ϵrad = 0.020 for the κ-jet models with ϵ = 0.0 and ϵ = 0.015, respectively. An important note is that we do no include X-ray emission in this work. Although, the obtained values are well below the thin disk efficiency, thus justifying our assumption that the radiation can be decoupled from the evolution of the dynamics of the plasma.

4.3. The effect of electron acceleration on synchrotron maps

At 43 and 86 GHz, both κ-jet models show a more dominant counter-jet when compared to the thermal-jet model, hinting to a behavior that could be observable by future GMVA-ALMA observations. There is also a clear difference in the extent of the emission of the forward jet in the ϵ = 0.015 κ-jet model when compared to the ϵ = 0.0 κ-jet and to the thermal-jet model, with the emitting region being more compact in the κ-jet models at 228 GHz. The reason for this is that there is more energy available at higher γ in the eDF, which results in a higher flux contribution per unit mass. Indeed, to obtain a fit to the data, a lower mass-accretion rate is needed. Since our mass-accretion rate sets the scaling of the densities and magnetic fields, it also changes the optical thickness of the source. As a result, a more optically thin model will show a more compact emission region.

A comparison with the result from Mościbrodzka et al. (2016) shows that similar source morphologies at all frequencies for the thermal model. However, at 228 GHz our images show a more optically thick inner ring feature that partially blocks the view to the shadow. The reason for this is that our initial conditions differ from those of Mościbrodzka et al. (2016), as they used a disk with a pressure maximum at 24 rg, resulting in an outer radius of r = 240 rg, while we used a pressure maximum at 12 rg and outer radius of r = 40 rg. A larger disk is initially seeded with larger toroidal magnetic-field loops, and a larger loop increases the magnetic flux at the horizon at later times. These stronger magnetic fields will affect the overall source morphology, resulting in wider opening angles which lead to less obscuration of the shadow by the forward jet.

4.4. Core size, shift, and jet-opening angle

The obtained core sizes for our models are close to the observed values: θ43 GHz  =  0.13 ± 0.01, θ86 GHz  =  0.079 ± 0.021 (Hada et al. 2013), and θ228 GHz  =  0.040 ± 0.002 (Doeleman et al. 2012). If we compare these to values reported in Tables 4, we find that our models at 43 and 86 GHz are within the error margins of the observations. At 228 GHz, the ϵ = 0.0 κ-jet recovers the observational value. The thermal-jet model is slightly larger, this is probably caused by the larger emission region around the shadow. In the ϵ = 0.015 case, the deviation is caused by a more pronounced jet feature.

We obtain core-shift relations for both our models by calculating the core position that follows the trend found by Hada et al. (2011). They computed the core shift with respect to the 43 GHz core. Their obtained fit is then extrapolated to higher frequencies where they find an offset of 40 μas at 228 GHz. At frequencies below 10 GHz, deviations with the fit from Hada et al. (2011) are present. A possible explanation for this is the limited simulation domain of 1000 rg and the de-collimation of the jet after r ≈ 300 rg.

An important remark to make is that we have here considered a Standard And Normal Evolution (SANE) simulation. This results in a low magnetic flux at the event horizon when compared MAD simulations that result in the maximum amount of flux that can penetrate the event horizon (Narayan et al. 2003; Tchekhovskoy et al. 2011). If we compare our results with the MAD simulation from Chael et al. (2019), our jet-opening angle is smaller and our models are inconsistent with the observational constraints on the jet-opening angle at 43 GHz (55° Walker et al. 2018; Janssen et al. 2019); by contrast, Chael et al. (2019) showed that their thermal MAD simulations do match the observed opening angle.

4.5. Reconnection as the source of particle acceleration

The electrons’ energy-distribution function is one of the key open questions in modeling the appearance of jets launched by supermassive black holes. Simulations of these acceleration mechanisms rely on non-ideal effects, which are not captured in GRMHD-based simulations. Fully resistive treatments of the plasma using non-ideal GRMHD simulations (Palenzuela et al. 2009; Ohsuga et al. 2009; Dionysopoulou et al. 2013; Bucciantini & Del Zanna 2013; Del Zanna et al. 2016; Qian et al. 2017, 2018; Del Zanna & Bucciantini 2018; Mignone et al. 2019; Ripperda et al. 2019) or general-relativistic PIC simulations (Watson & Nishikawa 2010; Levinson & Cerutti 2018; Parfrey et al. 2019) are being developed and will help to provide detailed answers to these questions in the future. In principle, alternative acceleration mechanisms could be at work, such as shocks. In our model, the main region of emission is where the magnetization σ is around unity, where shocks are known to be less efficient (see e.g., Sironi et al. 2015).

4.6. The Event Horizon Telescope results

In EHT Collaboration (2019e), GRMHD models were used to interpret the first image of a black hole. In the post-processing of the GMRHD models, only a thermal distribution function were considered. In this work, we show the effect of electron acceleration by performing a comparison with a purely thermal model. The overall trend is that the emission region is optically thinner and smaller in size. Also the accretion rates and jet-power drop, which could have implication for some of the models reported in EHT Collaboration (2019e). In this work we performed a model-to-model comparison where we use best fit parameters (such as Rhigh = 100) from previous works by e.g. Mościbrodzka et al. (2016). From these works we know that when Rhigh is decreased the spectral slope of the radio spectrum will increase. If we vary our newly introduced parameters we see that decreasing the rinj parameter will result in overproducing the NIR, while increasing η would result in overproducing the radio in the jet, but it must be noted that in this work we did not perform a detailed parameter scan over the full parameter domain as is done in for example EHT Collaboration (2019e). This scan combined with a detailed comparison with respect to the EHT data is beyond the scope of this work, and will be addressed in future works. A first step to generate realistic synthetic data based on the models presented here will be discussed by Roelofs et al. (2019).

5. Conclusion

We have presented a κ-jet model for the accreting black hole in M 87 based on an AMR GRMHD simulation in Cartesian–Kerr–Schild coordinates, coupled to radiative-transfer calculations that include sub-grid models for electron acceleration based on reconnection in the magnetized jet. The use of a Cartesian grid with AMR resulted in a high-resolution jet simulation that we used to model the jet launching point in M 87. We have demonstrated that we can obtain a fit to the M 87 SED from radio up to NIR if there is an accelerated electron population present in the jet. The model does not evolve the electron distribution function in time and does not include cooling; both of these aspects will be considered in future works. The jet-opening angle at 43GHz is too narrow, Chael et al. (2019) showed that a MAD type accretion disk can recover this opening angle, and we plan to explore this setup in future works with the addition of particle acceleration. The model reproduces the broadband SED from radio up to NIR, observed source sizes, core shifts and recovers a jet power that is consistent with observations.


1

PI: R. Craig Walker, project code: BW0106.

Acknowledgments

The authors thank M. Moscibrodzka, C. Gammie, A. Philippov, Z. Younsi, and B. Ripperda for valuable discussions and feedback during the project. We thank the referee Luca Del Zanna for thoughtful review and constructive comments to improve our article. This work was funded by the ERC Synergy Grant “BlackHoleCam-Imaging the Event Horizon of Black Holes” (Grant 610058, Goddi et al. 2017). The GRMHD simulations were performed on the Dutch National Supercomputing cluster Cartesius and are funded by the NWO computing grant 16431. The VLBA data shown in Fig. 11 is from project code: BW0106 PI: R. Craig Walker. The authors used python (Oliphant 2007; Millman & Aivazis 2011), scipy (Jones et al. 2001), numpy (van der Walt et al. 2011), and matplotlib (Hunter 2007) for generating the figures shown in this manuscript and for performing parts of the analyses. This research has made use of NASA’s Astrophysics Data System.

References

  1. Akiyama, K., Lu, R.-S., Fish, V. L., et al. 2015, ApJ, 807, 150 [NASA ADS] [CrossRef] [Google Scholar]
  2. Anjiri, M., Mignone, A., Bodo, G., & Rossi, P. 2014, MNRAS, 442, 2228 [NASA ADS] [CrossRef] [Google Scholar]
  3. Baade, W., & Minkowski, R. 1954, ApJ, 119, 215 [NASA ADS] [CrossRef] [Google Scholar]
  4. Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
  5. Ball, D., Sironi, L., & Özel, F. 2018, ApJ, 862, 80 [NASA ADS] [CrossRef] [Google Scholar]
  6. Blandford, R. D., & Königl, A. 1979, ApJ, 232, 34 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bolton, J. G., Stanley, G. J., & Slee, O. B. 1949, Nature, 164, 101 [NASA ADS] [CrossRef] [Google Scholar]
  8. Broderick, A. E., & Loeb, A. 2009, ApJ, 697, 1164 [NASA ADS] [CrossRef] [Google Scholar]
  9. Broderick, A. E., Narayan, R., Kormendy, J., et al. 2015, ApJ, 805, 179 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bronzwaer, T., Davelaar, J., Younsi, Z., et al. 2018, A&A, 613, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Bucciantini, N., & Del Zanna, L. 2013, MNRAS, 428, 71 [NASA ADS] [CrossRef] [Google Scholar]
  12. Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016, ApJ, 829, 11 [NASA ADS] [CrossRef] [Google Scholar]
  13. Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018a, ApJ, 857, 23 [Google Scholar]
  14. Chael, A., Rowan, M., Narayan, R., Johnson, M., & Sironi, L. 2018b, MNRAS, 478, 5209 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chael, A., Narayan, R., & Johnson, M. D. 2019, MNRAS, 486, 2873 [NASA ADS] [CrossRef] [Google Scholar]
  16. Christie, I. M., Petropoulou, M., Sironi, L., & Giannios, D. 2019, MNRAS, 482, 65 [NASA ADS] [CrossRef] [Google Scholar]
  17. Curtis, H. D. 1918, Publ. Lick Observatory, 13, 9 [Google Scholar]
  18. Davelaar, J., Mościbrodzka, M., Bronzwaer, T., & Falcke, H. 2018a, A&A, 612, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Davelaar, J., Bronzwaer, T., Kok, D., et al. 2018b, Comput. Astrophys. Cosmol., 5, 1 [NASA ADS] [CrossRef] [Google Scholar]
  20. de Gasperin, F., Orrú, E., Murgia, M., et al. 2012, A&A, 547, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Del Zanna, L., & Bucciantini, N. 2018, MNRAS, 479, 657 [NASA ADS] [Google Scholar]
  22. Del Zanna, L., Papini, E., Landi, S., Bugli, M., & Bucciantini, N. 2016, MNRAS, 460, 3753 [NASA ADS] [CrossRef] [Google Scholar]
  23. Dexter, J., McKinney, J. C., & Agol, E. 2012, MNRAS, 421, 1517 [NASA ADS] [CrossRef] [Google Scholar]
  24. Dionysopoulou, K., Alic, D., Palenzuela, C., Rezzolla, L., & Giacomazzo, B. 2013, Phys. Rev. D, 88, 044020 [NASA ADS] [CrossRef] [Google Scholar]
  25. Doeleman, S. S., Fish, V. L., Schenck, D. E., et al. 2012, Science, 338, 355 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  26. EHT Collaboration, et al. 2019a, ApJ, 875, L1 (Paper I) [NASA ADS] [CrossRef] [Google Scholar]
  27. EHT Collaboration, et al. 2019b, ApJ, 875, L2 (Paper II) [NASA ADS] [CrossRef] [Google Scholar]
  28. EHT Collaboration, et al. 2019c, ApJ, 875, L3 (Paper III) [NASA ADS] [CrossRef] [Google Scholar]
  29. EHT Collaboration, et al. 2019d, ApJ, 875, L4 (Paper IV) [NASA ADS] [CrossRef] [Google Scholar]
  30. EHT Collaboration, et al. 2019e, ApJ, 875, L5 (Paper V) [NASA ADS] [CrossRef] [Google Scholar]
  31. EHT Collaboration, et al. 2019f, ApJ, 875, L6 (Paper VI) [Google Scholar]
  32. Falcke, H., & Biermann, P. L. 1995, A&A, 293, 665 [NASA ADS] [Google Scholar]
  33. Falcke, H., Melia, F., & Agol, E. 2000, ApJ, 528, L13 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  34. Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962 [NASA ADS] [CrossRef] [Google Scholar]
  35. Gebhardt, K., Adams, J., Richstone, D., et al. 2011, ApJ, 729, 119 [NASA ADS] [CrossRef] [Google Scholar]
  36. Goddi, C., Falcke, H., Kramer, M., et al. 2017, Int. J. Mod. Phys. D, 26, 1730001 [NASA ADS] [CrossRef] [Google Scholar]
  37. Guo, F., Li, H., Daughton, W., & Liu, Y.-H. 2014, Phys. Rev. Lett., 113, 155005 [NASA ADS] [CrossRef] [Google Scholar]
  38. Hada, K., Doi, A., Kino, M., et al. 2011, Nature, 477, 185 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  39. Hada, K., Kino, M., Doi, A., et al. 2013, ApJ, 775, 70 [NASA ADS] [CrossRef] [Google Scholar]
  40. Hada, K., Kino, M., Doi, A., et al. 2016, ApJ, 817, 131 [NASA ADS] [CrossRef] [Google Scholar]
  41. Hillas, A. M. 1984, ARA&A, 22, 425 [NASA ADS] [CrossRef] [Google Scholar]
  42. Hu, M. 1962, IRE Trans. Inf. Theory, 8, 179 [Google Scholar]
  43. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  44. Janssen, M., Goddi, C., van Bemmel, I. M., et al. 2019, A&A, 626, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python [Online] [Google Scholar]
  46. Junor, W., Biretta, J. A., & Livio, M. 1999, Nature, 401, 891 [NASA ADS] [CrossRef] [Google Scholar]
  47. Kerr, R. P. 1963, Phys. Rev. Lett., 11, 237 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  48. Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018, A&A, 616, A188 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Krichbaum, T. P., Graham, D. A., Bremer, M., et al. 2006, in J. Phys. Conf. Ser., eds. R. Schödel, G. C. Bower, M. P. Muno, S. Nayakshin, & T. Ott, 328 [NASA ADS] [CrossRef] [Google Scholar]
  50. Leung, P. K., Gammie, C. F., & Noble, S. C. 2011, ApJ, 737, 21 [NASA ADS] [CrossRef] [Google Scholar]
  51. Levinson, A., & Cerutti, B. 2018, A&A, 616, A184 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Löhner, R. 1987, Comput. Methods Appl. Mech. Eng., 61, 323 [Google Scholar]
  53. Ly, C., Walker, R. C., & Wrobel, J. M. 2004, AJ, 127, 119 [NASA ADS] [CrossRef] [Google Scholar]
  54. Markoff, S., Nowak, M., Young, A., et al. 2008, ApJ, 681, 905 [NASA ADS] [CrossRef] [Google Scholar]
  55. Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2016, A&A, 596, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Mei, S., Blakeslee, J. P., Côté, P., et al. 2007, ApJ, 655, 144 [NASA ADS] [CrossRef] [Google Scholar]
  57. Mignone, A., Mattia, G., Bodo, G., & Del Zanna, L. 2019, MNRAS, 486, 4252 [NASA ADS] [CrossRef] [Google Scholar]
  58. Millman, K. J., & Aivazis, M. 2011, Comput. Sci. Eng., 13, 9 [CrossRef] [Google Scholar]
  59. Mills, B. Y. 1952, Nature, 170, 1063 [NASA ADS] [CrossRef] [Google Scholar]
  60. Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, A&A, 586, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Mościbrodzka, M., Dexter, J., Davelaar, J., & Falcke, H. 2017, MNRAS, 468, 2214 [NASA ADS] [CrossRef] [Google Scholar]
  62. Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69 [NASA ADS] [Google Scholar]
  63. Noble, S. C., Leung, P. K., Gammie, C. F., & Book, L. G. 2007, Classical Quantum Gravity, 24, S259 [NASA ADS] [CrossRef] [Google Scholar]
  64. Noble, S. C., Krolik, J. H., & Hawley, J. F. 2010, ApJ, 711, 959 [NASA ADS] [CrossRef] [Google Scholar]
  65. Ohsuga, K., Mineshige, S., Mori, M., & Kato, Y. 2009, PASJ, 61, L7 [NASA ADS] [Google Scholar]
  66. Oliphant, T. E. 2007, Comput. Sci. Eng., 9, 10 [CrossRef] [PubMed] [Google Scholar]
  67. Olivares, H., Porth, O., Davelaar, J., et al. 2019, A&A, 629, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611 [NASA ADS] [CrossRef] [Google Scholar]
  69. Palenzuela, C., Lehner, L., Reula, O., & Rezzolla, L. 2009, MNRAS, 394, 1727 [NASA ADS] [CrossRef] [Google Scholar]
  70. Pandya, A., Zhang, Z., Chandra, M., & Gammie, C. F. 2016, ApJ, 822, 34 [NASA ADS] [CrossRef] [Google Scholar]
  71. Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Lett., 122, 035101 [NASA ADS] [CrossRef] [Google Scholar]
  72. Petropoulou, M., Giannios, D., & Sironi, L. 2016, MNRAS, 462, 3325 [NASA ADS] [CrossRef] [Google Scholar]
  73. Pierrard, V., & Lazar, M. 2010, Sol. Phys., 267, 153 [NASA ADS] [CrossRef] [Google Scholar]
  74. Porth, O., Fendt, C., Meliani, Z., & Vaidya, B. 2011, ApJ, 737, 42 [Google Scholar]
  75. Porth, O., Olivares, H., Mizuno, Y., et al. 2017, Comput. Astrophys. Cosmol., 4, 1 [Google Scholar]
  76. Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26 [NASA ADS] [CrossRef] [Google Scholar]
  77. Prieto, M. A., Fernández-Ontiveros, J. A., Markoff, S., Espada, D., & González-Martín, O. 2016, MNRAS, 457, 3801 [NASA ADS] [CrossRef] [Google Scholar]
  78. Qian, Q., Fendt, C., Noble, S., & Bugli, M. 2017, ApJ, 834, 29 [NASA ADS] [CrossRef] [Google Scholar]
  79. Qian, Q., Fendt, C., & Vourellis, C. 2018, ApJ, 859, 28 [NASA ADS] [CrossRef] [Google Scholar]
  80. Reynolds, C. S., Fabian, A. C., Celotti, A., & Rees, M. J. 1996, MNRAS, 283, 873 [NASA ADS] [Google Scholar]
  81. Rezzolla, L., & Zanotti, O. 2013, Relativistic Hydrodynamics (Oxford: Oxford University Press) [CrossRef] [Google Scholar]
  82. Ripperda, B., Porth, O., Sironi, L., & Keppens, R. 2019, MNRAS, 485, 299 [NASA ADS] [CrossRef] [Google Scholar]
  83. Roelofs, F., Janssen, M., Natarajan, I., et al. 2019, A&A, submitted [Google Scholar]
  84. Ryan, B. R., Ressler, S. M., Dolence, J. C., Gammie, C., & Quataert, E. 2018, ApJ, 864, 126 [NASA ADS] [CrossRef] [Google Scholar]
  85. Sano, T., Inutsuka, S.-I., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321 [NASA ADS] [CrossRef] [Google Scholar]
  86. Sironi, L., & Spitkovsky, A. 2014, ApJ, 783, L21 [NASA ADS] [CrossRef] [Google Scholar]
  87. Sironi, L., Petropoulou, M., & Giannios, D. 2015, MNRAS, 450, 183 [NASA ADS] [CrossRef] [Google Scholar]
  88. Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79 [NASA ADS] [CrossRef] [Google Scholar]
  89. Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2012, J. Phys.: Conf. Ser., 372, 012040 [CrossRef] [Google Scholar]
  90. Tsallis, C. 1988, J. Stat. Phys., 52, 479 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  91. Turland, B. D. 1975, MNRAS, 170, 281 [NASA ADS] [CrossRef] [Google Scholar]
  92. van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
  93. Walker, R. C., Ly, C., Junor, W., & Hardee, P. J. 2008, J. Phys. Conf. Ser., 131, 012053 [NASA ADS] [CrossRef] [Google Scholar]
  94. Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128 [NASA ADS] [CrossRef] [Google Scholar]
  95. Walsh, J. L., Barth, A. J., Ho, L. C., & Sarzi, M. 2013, ApJ, 770, 86 [NASA ADS] [CrossRef] [Google Scholar]
  96. Watson, M., & Nishikawa, K.-I. 2010, Comput. Phys. Commun., 181, 1750 [NASA ADS] [CrossRef] [Google Scholar]
  97. Werner, G. R., & Uzdensky, D. A. 2017, ApJ, 843, L27 [NASA ADS] [CrossRef] [Google Scholar]
  98. Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJ, 816, L8 [NASA ADS] [CrossRef] [Google Scholar]
  99. Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840 [NASA ADS] [CrossRef] [Google Scholar]
  100. Xiao, F. 2006, Plasma Phys. Controlled Fusion, 48, 203 [Google Scholar]

Appendix A: Phenomenological model explaining the dominance of the counter-jet at 228 GHz

In certain GRMHD-based models of M 87, when imaged at 228 GHz at low (∼20 deg) inclination angles, it is observed that most of the emission that reaches the observer originates from the “counter-jet”, the jet facing away from the observer. Here we describe a simple phenomenological model that is capable of reproducing this effect.

Figure A.1 shows a schematic overview of our model, which consists of two rings of luminous material. The model is symmetric with respect to the equatorial plane, and the black-hole’s rotation axis passes through the center of the rings. The rings are meant to be an approximation of the “jet base” which appears on both sides of the equatorial plane in many GRMHD simulations. We make the assumption that the luminous rings are perfectly optically thin (equivalently, we ignore absorption in this model).

thumbnail Fig. A.1.

Schematic diagram of the model. A diagram of our phenomenological model of the dominant counter-jet. The central circle represents the black-hole’s event horizon. Two parameters control the placement of the rings (the model is symmetric with respect to the equatorial plane). The ring’s thickness is set to 1 rg.

Two main parameters define the geometry of the rings: θoffset and rring determine their principal diameter and distance from the equatorial plane. A step-function-based emissivity profile is then used, which relates the emissivity at a location in space to the distance between that location and the nearest point on the ring; it is equal to unity if that distance is smaller than 0.25 rg, and zero otherwise. In other words, the cross-sectional thickness of the ring is 0.5 rg.

We also assign a velocity vector to the material in the ring; this is done using a simple Keplerian model for the orbital velocity of a particle with orbital radius rring. The effect of this velocity vector is to cause the characteristic relativistic-boosting effect seen in most of our simulations; the ring is slightly brighter on the approaching side. This effect is minor in the present case, due to our low inclination angle.

Figure A.2 shows a typical image of our model, with its key features annotated for clarity. Figure A.3 illustrates the effect of varying rring. Figure A.4 compares two images that show only the upper (lower) ring. The flux observed from the lower ring is 30% higher than that of the upper ring. The lensed image(s) of the lower ring always appear “outside” (but near) the black-hole’s photon ring, potentially causing the observer to overestimate the black-hole shadow size. Although the black-hole mass can be derived from the size of the black-hole shadow, such an estimate should be seen as an upper limit in the present context.

thumbnail Fig. A.2.

Example image. Shown is the case of a Schwarzschild black hole (a* = 0). The observer inclination i is 20 deg. θoffset = 1 rad, rring = 6 rg. Note the double image of the lower ring, which appears larger in size than the upper ring due to lensing. The doubled image of the lower ring appears close to the photon ring, but is slightly larger.

thumbnail Fig. A.3.

Varying rring. Illustration of the effect of changing rring, for a Schwarzschild black hole (a* = 0) imaged with i = 20 deg and θoffset = 1 rad. Note how the two images of the lower ring coincide for the case rring = 4 rg.

thumbnail Fig. A.4.

Upper and lower ring. Comparison of images which show only the upper (left) or lower (right) ring, omitting the other (as before, a* = 0, i = 20 deg, θoffset = 1 rad, rring = 6 rg). The integrated flux density received from the lower ring is about 30% higher than that from the upper ring, due to gravitational lensing.

In this model, gravitational lensing causes most of the radiation emitted in the “polar regions” of a black hole to be redirected toward the opposite side of the black hole with respect to its origin. This effect could help to explain why the counter-jet in optically thin simulations of M 87 dominates over the observer-facing jet.

Figure A.5 shows an illustration of the lensing effect that causes an observer to record two images of the lower ring, and only one of the upper ring. We can understand this feature by considering the rays shown in Fig. A.5 from left to right (i.e., decreasing the rays’ impact parameter): the deflection angle increases as the rays curve more and more. The first object with which the rays intersect is the lower ring – hence we see that as the widest object on the observer’s image. The next rays, curving even more, intersect the lower ring again, but now they are moving back toward the observer (having traveled around the black hole). Moving on to rays with still smaller impact parameters, the rays now come very close to the photon ring. At this point, the deflection angle begins to diverge, causing rays to orbit the black hole an arbitrary number of times. These rays image the entire sky infinitely many times, producing a multitude of images of the environment. All of these images, however, are very small (and thus they don’t contribute much flux), and they are confined to a thin ring, which is infinitesimally close to the photon ring.

thumbnail Fig. A.5.

Visualization of the “doubling effect”. The observer in this image is directly above the black hole (i.e., the inclination angle is zero). The black-hole’s event horizon is marked by the dark-grey circle, while the unstable-photon region is marked by the light-grey circle. Gravitational lensing enhances the overall size of the lower ring, although the divergence of the rays near the lower ring causes its images to have a reduced thickness. Two images of the lower ring appear; one due to rays that intersect the ring while moving away from the observer, the other due to rays that curve around the black hole and intersect the ring while moving toward the observer. This causes most of the flux that reaches the observer to originate in the lower ring, on the far side of the black hole. Adapted from an image by Alessandro Roussel.

The doubling effect is only visible at low inclination angles, when the system is viewed in a face-on manner; the effect vanishes entirely at inclination angles near 90° (symmetry demands that both rings then contribute equally to the integrated flux density of the image). Complications also arise when absorption is taken into account; an optically thick accretion disk may absorb much of the lensed radiation originating from the lower ring.

As a final comment, we note that the doubling of the lower ring due to gravitational lensing occurs everywhere along the lower ring. Therefore, a partial ring or even a very compact structure (e.g., a plasmoid or “Gaussian hot-spot”) will show the same behavior; most of the radiation will come from the hot-spot on the opposite side of the black hole, away from the observer.

Appendix B: Comparison with spherical grids

We benchmarked our CKS simulation to a spherical simulation in modified Kerr–Schild (MKS) coordinates evolved from the same initial condition. The MKS simulation used an effective resolution of 512 × 192 × 192 in log(r), θ and ϕ, where the radial domain extended up to rout = 3333 rg. We compared the accretion rate and magnetic flux through the horizon ΦBH from both simulations, which are shown in Fig. B.1. In the nonlinear phase both the accretion rate and magnetic flux are of the same order of magnitude. At the end of the simulation, the normalized magnetic flux is, in both cases, Φ BH / M ˙ 1.0 $ \Phi_{\mathrm{BH}}/ \sqrt{\dot{M}}\approx1.0 $. Therefore, the CKS model belongs to the category of Standard and Normal Evolution (SANE; Tchekhovskoy et al. 2012).

thumbnail Fig. B.1.

Comparing Cartesian with a spherical simulation. Top panel: accretion rate through the event horizon as a function of time for the CKS simulations used in this work and a reference simulation in spherical MKS coordinates. Bottom panel: normalized magnetic flux through the horizon for the same pair of simulations. Both quantities show a consistent behavior.

The MRI quality factor Q MRI ( θ ) $ Q^{(\theta)}_{\mathrm{MRI}} $ is shown in Fig. B.2. In the majority of the disk the values of the quality factor are sufficiently above six as recommended in Sano et al. (2004), although close to the horizon a drop is visible. However, this does not reflect into significant differences in the reported accretion rates and magnetic fluxes from Fig. B.1. Throughout the jet, the resolution of the Cartesian grid exceeds that of the spherical grid significantly.

thumbnail Fig. B.2.

MRI quality factor Q MRI ( θ ) $ Q^{(\theta)}_{\mathrm{MRI}} $. Shown is the xy plane (left) and on the xz plane (right). Q MRI ( θ ) $ Q^{(\theta)}_{\mathrm{MRI}} $ is sufficient in most of the disk, with a drop visible only very close to the event horizon.

All Tables

Table 1.

Maximum AMR refinement radii in rg for the different AMR levels.

Table 2.

Model parameters.

Table 3.

Forward jet to counter jet ratio.

Table 4.

Source size.

All Figures

thumbnail Fig. 1.

Example of the κ-distribution (orange) for a κ value of 3.5, and for comparison a Maxwell–Jüttner distribution (black), and a power-law distribution (yellow) with p = 2.5. The κ-distribution function is a combination of a powerlaw with a thermal core.

In the text
thumbnail Fig. 2.

Parameterization of the κ parameter. κ as function of β and σ as found by Ball et al. (2018). A high value of κ corresponds to steep particle spectra with power-law index p = κ − 1. We overplotted contours of constant κ in black.

In the text
thumbnail Fig. 3.

GRMHD simulation snapshot. Left panels: slice in the (x, z) and (x, y) planes of the density in code units. Middle panels: slice along the (x, z) and (x, y) planes of the magnetization parameter b2/ρ, over-plotted with the grid block sizes. Right panels: slice along the (x, z) and (x, y) planes of the dimensionless ion temperature. Shown with a black circle is the location of the event horizon.

In the text
thumbnail Fig. 4.

Mass accretion rate and jet and wind power. We computed global quantities of the Cartesian GRMHD simulation, accretion rate in code units as a function of time (left panel), and jet and wind power in code units as a function of time (right panel).

In the text
thumbnail Fig. 5.

Lineprofiles. To compare our simulation with analytical prediction for plasma quantities by Blandford & Königl (1979) and Falcke & Biermann (1995) in the jet we computed radial profiles. We show radial profiles of the dimensionless electron density (left), magnetic-field strength (middle) and electron temperature (right). Black lines correspond to jet averaged quantities, dashed yellow lines to disk-averaged quantities and the red dashed lines correspond to power-law profiles predicted in analytical works. Also, the jet sheath is close to isothermality.

In the text
thumbnail Fig. 6.

κ and w parameter. Left panel: map of the κ parameter based on Eq. (18), where the jet spine is excluded (σ <  5). The map shows that the κ parameter is small in the jet sheath and large in the disk. Middle panel: map of the w parameter when ϵ = 0, showing large values in the sheath compared to the disk. Right panel: same as middle panel but now for ϵ = 0.015, the sheath has slightly larger values compared to the ϵ = 0 case.

In the text
thumbnail Fig. 7.

Spectral energy distributions. Shown are the thermal-jet (orange) and κ-jet with ϵ = 0 (orange) and ϵ = 0.015 (black) with their corresponding rms, overplotted with observational data points by Doeleman et al. (2012), Akiyama et al. (2015), Prieto et al. (2016), and Kim et al. (2018).

In the text
thumbnail Fig. 8.

Synthetic synchrotron maps. Shown are all our models at three astronomically relevant frequencies. From left to right: 43, 86, and 228 GHz. Top row: synthetic images at a single snapshot of the thermal-jet at an inclination of i = 160°. Second row: same as top row but for the ϵ = 0.0 κ-jet. Bottom row: same as the first and second row but for the ϵ = 0.015 κ-jet.

In the text
thumbnail Fig. 9.

Logarithmic optical-depth maps. Single snapshot at 228 GHz of the models at an inclination of i = 160°.

In the text
thumbnail Fig. 10.

Core shift. RA offset from the 43 GHz core as a function of frequency. Orange triangles correspond to a thermal-jet, black dots to a κ-jet, grey line represent the observational fit rRA(ν) = (1.4 ± 0.16)ν−0.94 ± 0.09 to the M 87 core by Hada et al. (2011).

In the text
thumbnail Fig. 11.

43 GHz radio map. Comparison of our thermal and κ-jet model with a VLBI map of M 87 at 43 GHz. First panel from the left: 43 GHz radio map of M 87 (Janssen et al. 2019). Second panel: synchrotron map of the thermal-jet model, convolved with a 2D Gaussian beam. Third panel: same as the second panel but now for a κ-jet model with ϵ = 0.0. Fourth panel: same as the third panel but now with ϵ = 0.015. The white ellipse indicates the beam used to convolve the images. All models produce a jet that is too narrow compared to the VLBI map. The extent of the jet increases when electron acceleration is present, and is maximum for ϵ = 0.015.

In the text
thumbnail Fig. A.1.

Schematic diagram of the model. A diagram of our phenomenological model of the dominant counter-jet. The central circle represents the black-hole’s event horizon. Two parameters control the placement of the rings (the model is symmetric with respect to the equatorial plane). The ring’s thickness is set to 1 rg.

In the text
thumbnail Fig. A.2.

Example image. Shown is the case of a Schwarzschild black hole (a* = 0). The observer inclination i is 20 deg. θoffset = 1 rad, rring = 6 rg. Note the double image of the lower ring, which appears larger in size than the upper ring due to lensing. The doubled image of the lower ring appears close to the photon ring, but is slightly larger.

In the text
thumbnail Fig. A.3.

Varying rring. Illustration of the effect of changing rring, for a Schwarzschild black hole (a* = 0) imaged with i = 20 deg and θoffset = 1 rad. Note how the two images of the lower ring coincide for the case rring = 4 rg.

In the text
thumbnail Fig. A.4.

Upper and lower ring. Comparison of images which show only the upper (left) or lower (right) ring, omitting the other (as before, a* = 0, i = 20 deg, θoffset = 1 rad, rring = 6 rg). The integrated flux density received from the lower ring is about 30% higher than that from the upper ring, due to gravitational lensing.

In the text
thumbnail Fig. A.5.

Visualization of the “doubling effect”. The observer in this image is directly above the black hole (i.e., the inclination angle is zero). The black-hole’s event horizon is marked by the dark-grey circle, while the unstable-photon region is marked by the light-grey circle. Gravitational lensing enhances the overall size of the lower ring, although the divergence of the rays near the lower ring causes its images to have a reduced thickness. Two images of the lower ring appear; one due to rays that intersect the ring while moving away from the observer, the other due to rays that curve around the black hole and intersect the ring while moving toward the observer. This causes most of the flux that reaches the observer to originate in the lower ring, on the far side of the black hole. Adapted from an image by Alessandro Roussel.

In the text
thumbnail Fig. B.1.

Comparing Cartesian with a spherical simulation. Top panel: accretion rate through the event horizon as a function of time for the CKS simulations used in this work and a reference simulation in spherical MKS coordinates. Bottom panel: normalized magnetic flux through the horizon for the same pair of simulations. Both quantities show a consistent behavior.

In the text
thumbnail Fig. B.2.

MRI quality factor Q MRI ( θ ) $ Q^{(\theta)}_{\mathrm{MRI}} $. Shown is the xy plane (left) and on the xz plane (right). Q MRI ( θ ) $ Q^{(\theta)}_{\mathrm{MRI}} $ is sufficient in most of the disk, with a drop visible only very close to the event horizon.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.