Issue 
A&A
Volume 632, December 2019



Article Number  A2  
Number of page(s)  16  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201936150  
Published online  21 November 2019 
Modeling nonthermal emission from the jetlaunching region of M 87 with adaptive mesh refinement
^{1}
Department of Astrophysics/IMAPP, Radboud University, PO Box 9010 6500 GL Nijmegen, The Netherlands
email: j.davelaar@astro.ru.nl
^{2}
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
^{3}
Institut für Theoretische Physik, MaxvonLaueStraße 1, 60438 Frankfurt am Main, Germany
^{4}
Anton Pannekoek Instituut, Universiteit van Amsterdam, PO Box 94249 1090 GE Amsterdam, The Netherlands
^{5}
MaxPlanck Institute for Radio Astronomy, Auf dem Huegel 69, 53115 Bonn, Germany
Received:
21
June
2019
Accepted:
21
August
2019
Context. The galaxy M 87 harbors a kiloparsecscale relativistic jet, whose origin coincides with a compact source thought to be a supermassive black hole. Observational millimeter very long baseline interferometry campaigns are capable of resolving the jetlaunching region at the scale of the event horizon. In order to provide a context for interpreting these observations, realistic generalrelativistic magnetohydrodynamical (GRMHD) models of the accretion flow are constructed.
Aims. Electrons in the jet are responsible for the observed synchrotron radiation, which is emitted in frequencies ranging from radio to nearinfrared (NIR) and optical. The characteristics of the emitted radiation depend on the shape of the electrons’ energydistribution function (eDF). The dependency on the eDF is omitted in the modeling of the first Event Horizon Telescope results. In this work, we aim to model the M 87 spectralenergy distribution from radio up to optical frequencies using a thermalrelativistic Maxwell–Jüttner distribution, as well as a relativistic κdistribution function. The powerlaw index of the eDF is modeled based on subgrid, particleincell parametrizations for subrelativistic reconnection.
Methods. A GRMHD simulation in Cartesian–Kerr–Schild coordinates, using eight levels of adaptive mesh refinement (AMR), forms the basis of our model. To obtain spectra and images, the GRMHD data was postprocessed with the raytracing code RAPTOR, which is capable of ray tracing through GRMHD simulation data that is stored in multilevel AMR grids. The resulting spectra and images maps are compared with observations.
Results. We obtain radio spectra in both the thermaljet and κjet models consistent with radio observations. Additionally, the κjet models also recover the NIR and optical emission. The images show a more extended structure at 43 GHz and 86 GHz and more compact emission at 228 GHz. The models recover the observed source sizes and core shifts and obtain a jet power of ≈10^{43} ergs s^{−1}. In the κjet models, both the accretion rates and jet powers are approximately two times lower than the thermaljet model. The frequency cutoff observed at ν ≈ 10^{15} Hz is recovered when the accelerator size is 10^{6} − 10^{8} cm, this could potentially point to an upper limit for plasmoid sizes in the jet of M 87.
Key words: black hole physics / accretion / accretion disks / radiation mechanisms: nonthermal / acceleration of particles / radiative transfer
© 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 millimeterwavelengths, 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 M_{BH} = 6.2 × 10^{9} M_{⊙} (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 (mmVLBI) observations by the Event Horizon Telescope Collaboration, which has the aim to spatially resolve blackhole 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 blackholes event horizon, for a nonrotating black hole its size on the sky is given by $2\sqrt{27}GM/({c}^{2}D)$, with G the Gravitational constant, M the blackhole 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 mmVLBI 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, forcefree jet model. Their bestfit model is consistent with 43 GHz observations. The model parameters include a blackhole spin of a_{*} = Jc/GM^{2} = 0.998, a viewing angle of i = 25°, and a jet footpoint at r = 10 r_{g}, where the gravitational radius r_{g} is defined as r_{g} = GM/c^{2}. 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 blackhole mass of M = 3.4 × 10^{9} M_{⊙} (Walsh et al. 2013).
Generalrelativistic 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 highspin GRMHD simulation. Their models included a thermal electron population in the disk and a powerlaw based electron population in the jet. Their bestfit model, at an inclination of 25°, showed counterjet 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 massaccretion rate of Ṁ ≈ 10^{−3} M_{⊙} yr^{−1}, and a powerlaw index of the nonthermal electron distribution function of p = 3.25 − 3.5, where they used a constant electrontoproton temperature ratio of T_{p}/T_{e} = 10.
Mościbrodzka et al. (2016) used GRMHD simulations and a Monte Carlobased radiativetransfer code to model the full spectral energy distribution (SED) of an accreting supermassive black hole from radio to Xray, as well as raytraced 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 iontoelectron temperature as a function of plasma β = P/P_{mag}, where P is the gas pressure and P_{mag} 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 iontoelectron temperature ratio in the disk of T_{i}/T_{e} = 100. Smaller values of the iontoelectron temperature ratio resulted in an excess of Xray emission. The 230 GHz images showed counterjet 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 bestfit jetdominated 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 2Daxisymmetric twotemperature GRMHD simulations that include radiative cooling. The authors conclude that radiative cooling is important in the inner region (r < 10 r_{g}) of the accretion flow, and that the blackhole mass of M = 6.2 × 10^{9} M_{⊙} and spin a_{*} = 0.9375 simulation recovers the observed radio and Xray emission and image size at 230 GHz. The jetopening angle in their model at lower frequencies is too narrow compared to the millimeterobservations 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 twotemperature radiative GRMHD model of a magnetically arrested disc (MAD; Narayan et al. 2003; Tchekhovskoy et al. 2011). The model recovers observables such as jetopening 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 ringlike structure in the radio core of M 87 at 228 GHz. This ringlike 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 oftenmade assumption that the electrons in the full simulation domain are in a thermalrelativistic Maxwell–Jüttner distribution potentially breaks down in regions where nonideal effects are important.
These nonideal 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 powerlaw with powerlaw 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).
Fig. 1. Example of the κdistribution (orange) for a κ value of 3.5, and for comparison a Maxwell–Jüttner distribution (black), and a powerlaw 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 nonthermal prescription used in this model was an electron distribution with a constant powerlaw in the outflow of the simulation domain with a fixed powerlaw index.
To improve the model we here connect the electronacceleration parameters to information from local kinetic plasma simulations. Kinetic plasma, or particleincell (PIC) simulations, are capable of resolving the microphysics scales that GRMHD simulations cannot reach. Although local, these type of simulations can provide firstprinciple parametrizations of particleacceleration processes. For our model, we consider a parametrization of the powerlaw index for transrelativistic 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 magnetorotational instability (MRI) in the disk. We use the results of this simulation to generate SEDs, synthetic synchrotron maps (images), and opticaldepth maps of the jetlaunching region in M 87. We extend the generalrelativisticraytracing (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 radiativetransfer 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 electronphysics 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
$$\begin{array}{c}\hfill {\mathrm{\nabla}}_{\mu}(\rho {u}^{\mu})=0,\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(1\mathrm{a})\end{array}$$
$$\begin{array}{c}\hfill {\mathrm{\nabla}}_{\mu}{T}^{\mu \nu}=0,\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(1\mathrm{b})\end{array}$$
$$\begin{array}{c}\hfill {\mathrm{\nabla}}_{\mu}{\phantom{\rule{0.166667em}{0ex}}}^{\ast}\phantom{\rule{0.166667em}{0ex}}{F}^{\mu \nu}=0,\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(1\mathrm{c})\end{array}$$
where ∇_{μ} denotes the covariant derivative, ρ the restmass density, u^{μ} the fluid 4velocity, T^{μν} the energymomentum tensor of the combined perfect fluid and electromagnetic fields, and ^{*}F^{μν} the dual of the Faraday tensor (F^{αβ}).
The system is closed by the idealMHD 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(\rho ,P)=1+\frac{\widehat{\gamma}}{\widehat{\gamma}1}\frac{P}{\rho}$, where h and P are the specific enthalpy and gas pressure in the fluid frame, and the adiabatic index $\widehat{\gamma}=4/3$. The simulation is initialized with a Fishbone–Moncrief torus (Fishbone & Moncrief 1976) with its inner radius at 6 r_{g} and its pressure maximum at 12 r_{g}. 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 P_{max} and maximum magnetic pressure P_{mag, max} is P_{max}/P_{mag, 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 blackhole’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 r_{g}.
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)
$$\begin{array}{c}\hfill {g}_{\mu \nu}={\eta}_{\mu \nu}+f{l}_{\mu}{l}_{\nu},\end{array}$$(2)
where η_{μν} = ( − 1, 1, 1, 1) is the Minkowski metric, and
$$\begin{array}{c}\hfill f=\frac{2{r}^{3}}{{r}^{4}+{a}_{\ast}^{2}{z}^{2}},\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(3\mathrm{a})\end{array}$$
$$\begin{array}{c}\hfill {l}_{\nu}=(1,\frac{rx+{a}_{\ast}y}{{r}^{2}+{a}^{2}},\frac{ry{a}_{\ast}x}{{r}^{2}+{a}_{\ast}^{2}},\frac{z}{r}),\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(3\mathrm{b})\end{array}$$
where r is given by
$$\begin{array}{c}\hfill {r}^{2}=\frac{{R}^{2}{a}_{\ast}^{2}+\sqrt{{({R}^{2}{a}^{2})}^{2}+4{a}_{\ast}^{2}{z}^{2}}}{2},\end{array}$$(4)
and
$$\begin{array}{c}\hfill {R}^{2}={x}^{2}+{y}^{2}+{z}^{2}.\end{array}$$(5)
All units of length are scaled by the gravitational radius r_{g} which is given by r_{g} = GM/c^{2}. In the limit of R ≫ a_{*}, the radius r → R. The contravariant metric is given by
$$\begin{array}{c}\hfill {g}^{\mu \nu}={\eta}^{\mu \nu}f{l}^{\mu}{l}^{\nu},\end{array}$$(6)
where l^{ν} is given by
$$\begin{array}{c}\hfill {l}^{\nu}=(1,\frac{rx+ay}{{r}^{2}+{a}_{\ast}^{2}},\frac{ry{a}_{\ast}x}{{r}^{2}+{a}_{\ast}^{2}},\frac{z}{r})\xb7\end{array}$$(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 σ = b^{2}/ρ, where $b=\sqrt{{b}^{\mu}{b}_{\mu}}$ is the magneticfield 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 zdirections, respectively. The simulation domain is x ∈ ( − 500 r_{g}, 500 r_{g}), y ∈ ( − 500 r_{g}, 500 r_{g}) and z ∈ ( − 1000 r_{g}, 1000 r_{g}). We simulate up to t_{f} = 10^{4} r_{g}/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.
Maximum AMR refinement radii in r_{g} for the different AMR levels.
2.3. Ray tracing in AMR CKS grid
In order to perform generalrelativistic raytracing calculations in Cartesian coordinates within the blockbased AMR data structure of BHAC, it has been necessary to extend our generalrelativistic raytracing 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 wavevectors 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 GrammSchmidt 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 (x_{c}, y_{c}, z_{c}) in space which is computed based on the following parameters: (i) the radial distance between the camera and the black hole r_{c}; (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
$$\begin{array}{c}\hfill {t}_{0}^{\mu}=(1,0,0,0),\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(8\mathrm{a})\end{array}$$
$$\begin{array}{c}\hfill {t}_{1}^{\mu}=(0,sin({\theta}_{\mathrm{c}})cos({\varphi}_{\mathrm{c}}),sin({\theta}_{\mathrm{c}})sin({\varphi}_{\mathrm{c}}),cos({\theta}_{\mathrm{c}})),\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(8\mathrm{b})\end{array}$$
$$\begin{array}{c}\hfill {t}_{2}^{\mu}=(0,sin({\varphi}_{\mathrm{c}}),cos({\varphi}_{\mathrm{c}})),\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(8\mathrm{c})\end{array}$$
$$\begin{array}{c}\hfill {t}_{3}^{\mu}=(0,cos({\theta}_{\mathrm{c}})cos({\varphi}_{\mathrm{c}}),cos({\theta}_{\mathrm{c}})sin({\varphi}_{\mathrm{c}}),sin({\theta}_{\mathrm{c}})).\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(8\mathrm{d})\end{array}$$
The choice of trial vectors results in a righthanded basis where the observer is facing the black hole.
The integration of the geodesic equations is done by solving the secondorder differential equation
$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^{2}{x}^{\alpha}}{\mathrm{d}{\lambda}^{2}}={\mathrm{\Gamma}}_{\phantom{\rule{4pt}{0ex}}\mu \nu}^{\alpha}\frac{\mathrm{d}{x}^{\mu}}{\mathrm{d}\lambda}\frac{\mathrm{d}{x}^{\nu}}{\mathrm{d}\lambda}\xb7\end{array}$$(9)
where ${\mathrm{\Gamma}}_{\phantom{\rule{4pt}{0ex}}\mu \nu}^{\alpha}$ are the connection coefficients, x^{α} is the geodesic position, and λ is the affine parameter. We use a fourthorder Runge–Kutta algorithm, where the connection coefficients are evaluated using a finitedifference derivative of the metric.
The stepsizing for the geodesic integration in RAPTOR was adopted since it relies on spherical logarithmic coordinates. First we compute a required stepsize based on the geodesic wavevector
$$\begin{array}{c}\hfill \mathrm{d}{\lambda}_{x}=\mathrm{\Delta}\phantom{\rule{4pt}{0ex}}/(\left{k}^{x}\right+\delta ),\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(10\mathrm{a})\end{array}$$
$$\begin{array}{c}\hfill \mathrm{d}{\lambda}_{y}=\mathrm{\Delta}\phantom{\rule{4pt}{0ex}}/(\left{k}^{y}\right+\delta ),\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(10\mathrm{b})\end{array}$$
$$\begin{array}{c}\hfill \mathrm{d}{\lambda}_{z}=\mathrm{\Delta}\phantom{\rule{4pt}{0ex}}/(\left{k}^{z}\right+\delta ),\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(10\mathrm{c})\end{array}$$
$$\begin{array}{c}\hfill \mathrm{d}{\lambda}_{\mathrm{geod}}=\frac{R}{{\left\mathrm{d}{\lambda}_{x}\right}^{1}+{\left\mathrm{d}{\lambda}_{y}\right}^{1}+{\left\mathrm{d}{\lambda}_{z}\right}^{1}},\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(10\mathrm{d})\end{array}$$
where k^{x, y, z} are the wavevector components in the x, y, z directions, δ is a small real number to prevent divisions by zero, and Δ is a scale factor for the stepsize (typically Δ ≈ 0.01). Then we compute a required stepsize based on the AMR cell size dx
$$\begin{array}{c}\hfill {k}_{\mathrm{max}}=max({k}^{x},max({k}^{y},{k}^{z})),\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(11\mathrm{a})\end{array}$$
$$\begin{array}{c}\hfill \mathrm{d}{\lambda}_{\mathrm{grid}}=\frac{\mathrm{d}x}{n{k}_{\mathrm{max}}},\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(11\mathrm{b})\end{array}$$
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 stepsizes and use the smallest of the two to ensure convergence; dλ = min(dλ_{geod}, dλ_{grid}).
For the radiativetransfer part of the raytracing 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 blockbased data structure is parsed by the code. When we integrate the geodesics we use a nearestneighbor approach to interpolate the gridbased plasma variables to the geodesics.
2.4. Electron model and radiativetransfer model parameters
Since GRMHD simulations are scalefree, we have to rescale the plasma variables from code units to c.g.s. units. Units of length are scaled with ℒ = r_{g}, while units of time are scaled with 𝒯 = r_{g}/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 × 10^{9} M_{⊙} (Gebhardt et al. 2011), the mass used in this work is slightly smaller than the mass of M = (6.5 ± 0.7)×10^{9} M_{⊙} 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}, u_{0} = ρ_{0}c^{2}, and ${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 postprocessing. The plasma is assumed to be chargeneutral, so that n_{e} = n_{p} throughout the domain. For the electron temperature we employ the parametrization of Mościbrodzka et al. (2016):
$$\begin{array}{c}\hfill {T}_{\mathrm{ratio}}={T}_{\mathrm{p}}/{T}_{\mathrm{e}}={R}_{\mathrm{low}}\frac{1}{1+{\beta}^{2}}+{R}_{\mathrm{high}}\frac{{\beta}^{2}}{1+{\beta}^{2}},\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(12\mathrm{a})\end{array}$$
$$\begin{array}{c}\hfill {\mathrm{\Theta}}_{\mathrm{e}}=\frac{U(\widehat{\gamma}1){m}_{\mathrm{p}}/{m}_{\mathrm{e}}}{\rho {T}_{\mathrm{ratio}}},\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(12\mathrm{b})\end{array}$$
where m_{p} is the proton mass, m_{e} is the electron mass, U is the internal energy, which is related to h and ρ via $h(\rho ,U)=1+\widehat{\gamma}\frac{U}{\rho}$, therefore, P = U(γ − 1), where P is the gass pressure. The dimensionless electron temperature Θ_{e} can be rescaled to c.g.s units via T = Θ_{e}m_{e}c^{2}/k_{b}, where k_{b} is the Boltzmann constant. The parameters R_{low} and R_{high} are free parameters of the model; R_{low} sets the temperature ratio in the jet, where β ≪ 1, and R_{high} sets the temperature ratio in the disk where β ≫ 1 (Fig. 2).
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 powerlaw index p = κ − 1. We overplotted contours of constant κ in black. 
For the electrons’ energydistribution 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)
$$\begin{array}{c}\hfill \frac{d{n}_{\mathrm{e}}}{d\gamma}=N\gamma \sqrt{{\gamma}^{2}1}{(1+\frac{\gamma 1}{\kappa w})}^{(\kappa +1)},\end{array}$$(13)
where γ is the Lorentz factor of the electrons, κ is the parameter that sets the powerlaw 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 n_{e} 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
$$\begin{array}{c}\hfill {E}_{\kappa}=\frac{3\kappa}{\kappa 3}{n}_{\mathrm{e}}w.\end{array}$$(14)
We couple this energy to the energy present in a thermal distribution (E_{thermal} = 3n_{e}Θ_{e}) and add a source term based on the magnetic energy
$$\begin{array}{c}\hfill {E}_{\kappa}=\frac{3\kappa}{\kappa 3}{n}_{\mathrm{e}}w=3{n}_{\mathrm{e}}{\mathrm{\Theta}}_{\mathrm{e}}+\stackrel{\sim}{\u03f5}\frac{{B}^{2}}{8\pi},\end{array}$$(15)
here $\stackrel{\sim}{\u03f5}$ is used to join smoothly between between the κdistribution and the magnetic energy. After a bit of algebra, we can rewrite the width as
$$\begin{array}{c}\hfill w=\frac{\kappa 3}{\kappa}{\mathrm{\Theta}}_{\mathrm{e}}+\stackrel{\sim}{\u03f5}\frac{\kappa 3}{6\kappa}\frac{{m}_{\mathrm{p}}}{{m}_{\mathrm{e}}}\sigma .\end{array}$$(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 $\stackrel{\sim}{\u03f5}$ parameter is set by
$$\begin{array}{c}\hfill \stackrel{\sim}{\u03f5}=\u03f5\frac{1}{2}(1+tanh(r{r}_{\mathrm{inj}})),\end{array}$$(17)
where r_{inj} is the injection radius from which the magnetic energy contributes to the w parameter, and ϵ is the base value for radii larger than r_{inj}; hereafter, we will consider two cases: where ϵ is zero or nonzero. 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 powerlaw index of the electrons distribution functions (eDFs) is based on subgrid particleincell (PIC) simulations of transrelativistic reconnection by Ball et al. (2018), who simulated twodimensional reconnection layer (Harris sheath) for an electronion 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 powerlaw index p as a function of β and σ as
$$\begin{array}{c}\hfill p={A}_{p}+{B}_{p}tanh\left({C}_{p}\beta \right)\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(18\mathrm{a})\end{array}$$
$$\begin{array}{c}\hfill {A}_{p}=1.8+0.7/\sqrt{\sigma}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(18\mathrm{b})\end{array}$$
$$\begin{array}{c}\hfill {B}_{p}=3.7\phantom{\rule{0.166667em}{0ex}}{\sigma}^{0.19}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(18\mathrm{c})\end{array}$$
$$\begin{array}{c}\hfill {C}_{p}=23.4\phantom{\rule{0.166667em}{0ex}}{\sigma}^{0.26}.\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(18\mathrm{d})\end{array}$$
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 jetmodels.
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 cutoff
The SED of M 87 shows a clear cutoff in flux around ν = 10^{15} Hz (Prieto et al. 2016). We will consider three potential sources for this cutoff.
First, we assume that the cutoff is caused by synchrotron cooling in the jet, which becomes important when the synchrotroncooling time of the electron is comparable with the typical dynamical time. Under these conditions, the cooling (cutoff) frequency is given by
$$\begin{array}{c}\hfill {\nu}_{\mathrm{cool}}=\frac{18\pi}{{\sigma}_{T}^{2}}\frac{{m}_{e}{c}^{2}e}{{B}^{3}{z}_{\mathrm{jet}}},\end{array}$$(19)
where σ_{T} is the Thomson crosssection, and z_{jet} the position along the jet.
Second, we assume that the break takes place at the synchrotron burnoff limit, that is, at the maximum energy that a particle can gain while emitting synchrotron radiation. The maximum Lorentz factor in this case is
$$\begin{array}{c}\hfill {\gamma}_{\mathrm{max}}=\sqrt{\frac{3{m}_{\mathrm{e}}^{2}{c}^{4}E}{4\pi {e}^{3}{B}^{2}}},\end{array}$$(20)
where E is the electric field, and the cutoff frequency is then given by
$$\begin{array}{c}\hfill {\nu}_{\phantom{\rule{0.333333em}{0ex}}\text{cutoff}}=\frac{3}{2}{\gamma}_{\mathrm{max}}^{2}{\nu}_{c},\end{array}$$(21)
with ν_{c} = eB/(2πm_{e}c).
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
$$\begin{array}{c}\hfill {\gamma}_{\mathrm{max}}=\frac{\mathit{eBL}}{{m}_{\mathrm{e}}{c}^{2}},\end{array}$$(22)
which results in a cutoff frequency of ν ≈ 10^{15} Hz after using Eq. (22) in (21). In this way, we can also we can estimate the typical size L of the accelerationregion
$$\begin{array}{c}\hfill L=\sqrt{\frac{4\pi {\nu}_{\mathrm{cut}\mathrm{off}}{m}_{\mathrm{e}}^{3}{c}^{5}}{{e}^{3}{B}^{3}}}\approx 4.5\times {10}^{7}\phantom{\rule{0.166667em}{0ex}}\mathrm{cm}\sqrt{\frac{({\nu}_{\mathrm{cut}\mathrm{off}}/{10}^{15}\phantom{\rule{0.166667em}{0ex}}\mathrm{Hz})}{{(B/1\phantom{\rule{3.33333pt}{0ex}}\mathrm{G})}^{3}}}\xb7\end{array}$$(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 wellresolved relativistic jet up to the edge of the simulation domain at 1000 r_{g} in the zdirection. At z = 40 r_{g} the jet diameter is resolved with 160 cells, and with 32 cells at z = 1000 r_{g}. 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 r_{g} surface
$$\begin{array}{c}\hfill \dot{E}={\displaystyle {\int}_{0}^{2\pi}{\int}_{0}^{\pi}({T}_{t}^{r}\rho {u}^{r})\phantom{\rule{0.166667em}{0ex}}{\chi}_{(\xb7)}\phantom{\rule{0.166667em}{0ex}}\sqrt{g}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\varphi}\end{array}$$(24)
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 b^{2}/ρ, overplotted 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. 
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
$$\begin{array}{c}\hfill {\chi}_{\mathrm{jet}}=({b}^{2}/\rho >1\phantom{\rule{3.33333pt}{0ex}}\mathrm{or}\phantom{\rule{3.33333pt}{0ex}}\mu >2)\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(25\mathrm{a})\end{array}$$
$$\begin{array}{c}\hfill {\chi}_{\mathrm{wind}}=(\mathrm{not}\phantom{\rule{3.33333pt}{0ex}}{\chi}_{\mathrm{jet}}\phantom{\rule{3.33333pt}{0ex}}\mathrm{and}\phantom{\rule{3.33333pt}{0ex}}h{u}_{t}>1)\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(25\mathrm{b})\end{array}$$
$$\begin{array}{c}\hfill {\chi}_{\mathrm{disk}}=(\mathrm{not}\phantom{\rule{3.33333pt}{0ex}}{\chi}_{\mathrm{jet}}\phantom{\rule{3.33333pt}{0ex}}\mathrm{and}\phantom{\rule{3.33333pt}{0ex}}\mathrm{not}\phantom{\rule{3.33333pt}{0ex}}{\chi}_{\mathrm{wind}}),\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}(25\mathrm{c})\end{array}$$
and μ denotes the energy flux normalized to the restmass energy in the radial direction $\mu =({T}_{t}^{r}\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 b^{2}/ρ > 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 jetmodels, 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 T_{e} ∝ 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}, magneticfield strength B, and electron temperature T_{e}. This is done by performing the following integral
$$\begin{array}{c}\hfill q(r)=\frac{1}{\mathrm{\Delta}t}{\displaystyle \int \left(\frac{\int {\int}_{0}^{2\pi}q(t,r,\theta ,\varphi )\sqrt{\gamma (r,\theta )}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\varphi}{\int {\int}_{0}^{2\pi}\sqrt{g(r,\theta )}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\varphi}\right)\phantom{\rule{0.166667em}{0ex}}\mathrm{d}t.}\end{array}$$(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 r_{g}/c and t = 10^{4} r_{g}/c, with a total of hundred snapshots. The computed radial profiles are shown in Fig. 5 and are overplotted 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 r_{g}, followed by is an increase of temperature that correlates with the break in the profile of the electron density. The break is caused by decollimation 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 r_{g}.
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), magneticfield strength (middle) and electron temperature (right). Black lines correspond to jet averaged quantities, dashed yellow lines to diskaveraged quantities and the red dashed lines correspond to powerlaw 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 nonthermal emission will be produced in this region.
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 thermaljet 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 fieldofview of the camera is set to be 1000 r_{g} in both the x and ydirections, 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 r_{g}/c and t = 10^{4} r_{g}/c), these have been fit to nonsimultaneous 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 thermaljet 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 r_{inj} = 10 r_{g} 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.
Model parameters.
Fig. 7. Spectral energy distributions. Shown are the thermaljet (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 powerlaw with an index of α ≈ −0.7, where α = −(p − 1)/2 for a power law distribution of nonthermal electrons F_{ν} ∝ ν^{α}. Furthermore, when compared to the thermal model, the flux in the κjet models is higher at lower frequencies (ν < 10^{11} Hz) and at the higher frequencies (5 × 10^{12} Hz < ν).
When considering the various cutoff models, the cooling cutoff turned out to be unimportant, in agreement with the findings of Mościbrodzka et al. (2016), Broderick et al. (2015). When using the synchrotron burnoff, the correct cutoff is obtained if E/B ≈ 10^{−6}, but no physical model is possible that recovers such a ratio. The only criterion that recovers the cutoff frequency is the Hillas criterion, which is obtained when the plasmoid size is set to L ≈ 10^{5} − 10^{7} 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 = 10^{4} r_{g}/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 M_{BH} = 6.2 × 10^{9} M_{⊙} (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.
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 thermaljet 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 thermaljet 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 thermaljet 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 opticaldepth 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 intensitymaps of Fig. 8. The reason behind this behavior is that lower massaccretion 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 massaccretion rates contribute more than their thermal counterpart.
Fig. 9. Logarithmic opticaldepth 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 counterjet, 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 timeaveraged 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 10^{4} r_{g}/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 counterjet, 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 blackhole’s event horizon.
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–10^{4} r_{g}/c and in Table 4 we report the timeaveraged major and minor fullwidth half maxima (FWHM) and their corresponding spread.
Source size.
We computed the core shift with respect to the blackhole’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 firstorder image moments of timeaveraged images and the comparison of the values obtained with the observational fit of Hada et al. (2011), meaning r_{RA}(ν) = (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 r_{core} ∝ ν^{−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 thermaljet model, probably because the counterjet is more dominant.
Fig. 10. Core shift. RA offset from the 43 GHz core as a function of frequency. Orange triangles correspond to a thermaljet, black dots to a κjet, grey line represent the observational fit r_{RA}(ν) = (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 thermaljet and κjet models with the 43 GHz VLBI observations, where M 87 was tracked for 8 hours with all VLBA stations^{1}. 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 recomputed synthetic images with a large field of view of 3.7 mas and convolved them with a 0.3 × 0.1 mas^{2} beamsize by using the ehtimaging 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 thermaljet 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 jetopening angle is wider.
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 thermaljet 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 r_{g} is resolved by 160 cells and at z = 1000 r_{g} 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 massaccretion rate and radial profiles of density, magneticfield 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 thermaljet model and κjet models, wherein the latter case we parametrize the powerlaw index of the eDF based on subgrid 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 thermaljet case. Our models use their bestfit 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 powerlaw. 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 10^{43} 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 thermaljet 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 thermalmodels. 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/Ṁc^{2}, we found that the thermaljet 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 Xray 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 counterjet when compared to the thermaljet model, hinting to a behavior that could be observable by future GMVAALMA 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 thermaljet 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 massaccretion rate is needed. Since our massaccretion 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 r_{g}, resulting in an outer radius of r = 240 r_{g}, while we used a pressure maximum at 12 r_{g} and outer radius of r = 40 r_{g}. A larger disk is initially seeded with larger toroidal magneticfield 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 jetopening 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 thermaljet 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 coreshift 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 r_{g} and the decollimation of the jet after r ≈ 300 r_{g}.
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 jetopening angle is smaller and our models are inconsistent with the observational constraints on the jetopening 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’ energydistribution 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 nonideal effects, which are not captured in GRMHDbased simulations. Fully resistive treatments of the plasma using nonideal 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 generalrelativistic 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 postprocessing 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 jetpower drop, which could have implication for some of the models reported in EHT Collaboration (2019e). In this work we performed a modeltomodel comparison where we use best fit parameters (such as R_{high} = 100) from previous works by e.g. Mościbrodzka et al. (2016). From these works we know that when R_{high} is decreased the spectral slope of the radio spectrum will increase. If we vary our newly introduced parameters we see that decreasing the r_{inj} 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 radiativetransfer calculations that include subgrid models for electron acceleration based on reconnection in the magnetized jet. The use of a Cartesian grid with AMR resulted in a highresolution 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 jetopening 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.
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 “BlackHoleCamImaging 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
 Akiyama, K., Lu, R.S., Fish, V. L., et al. 2015, ApJ, 807, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Anjiri, M., Mignone, A., Bodo, G., & Rossi, P. 2014, MNRAS, 442, 2228 [NASA ADS] [CrossRef] [Google Scholar]
 Baade, W., & Minkowski, R. 1954, ApJ, 119, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Ball, D., Sironi, L., & Özel, F. 2018, ApJ, 862, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Königl, A. 1979, ApJ, 232, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Bolton, J. G., Stanley, G. J., & Slee, O. B. 1949, Nature, 164, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A. E., & Loeb, A. 2009, ApJ, 697, 1164 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A. E., Narayan, R., Kormendy, J., et al. 2015, ApJ, 805, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Bronzwaer, T., Davelaar, J., Younsi, Z., et al. 2018, A&A, 613, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bucciantini, N., & Del Zanna, L. 2013, MNRAS, 428, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016, ApJ, 829, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018a, ApJ, 857, 23 [Google Scholar]
 Chael, A., Rowan, M., Narayan, R., Johnson, M., & Sironi, L. 2018b, MNRAS, 478, 5209 [NASA ADS] [CrossRef] [Google Scholar]
 Chael, A., Narayan, R., & Johnson, M. D. 2019, MNRAS, 486, 2873 [NASA ADS] [CrossRef] [Google Scholar]
 Christie, I. M., Petropoulou, M., Sironi, L., & Giannios, D. 2019, MNRAS, 482, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Curtis, H. D. 1918, Publ. Lick Observatory, 13, 9 [Google Scholar]
 Davelaar, J., Mościbrodzka, M., Bronzwaer, T., & Falcke, H. 2018a, A&A, 612, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Davelaar, J., Bronzwaer, T., Kok, D., et al. 2018b, Comput. Astrophys. Cosmol., 5, 1 [NASA ADS] [CrossRef] [Google Scholar]
 de Gasperin, F., Orrú, E., Murgia, M., et al. 2012, A&A, 547, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Del Zanna, L., & Bucciantini, N. 2018, MNRAS, 479, 657 [NASA ADS] [Google Scholar]
 Del Zanna, L., Papini, E., Landi, S., Bugli, M., & Bucciantini, N. 2016, MNRAS, 460, 3753 [NASA ADS] [CrossRef] [Google Scholar]
 Dexter, J., McKinney, J. C., & Agol, E. 2012, MNRAS, 421, 1517 [NASA ADS] [CrossRef] [Google Scholar]
 Dionysopoulou, K., Alic, D., Palenzuela, C., Rezzolla, L., & Giacomazzo, B. 2013, Phys. Rev. D, 88, 044020 [NASA ADS] [CrossRef] [Google Scholar]
 Doeleman, S. S., Fish, V. L., Schenck, D. E., et al. 2012, Science, 338, 355 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 EHT Collaboration, et al. 2019a, ApJ, 875, L1 (Paper I) [NASA ADS] [CrossRef] [Google Scholar]
 EHT Collaboration, et al. 2019b, ApJ, 875, L2 (Paper II) [NASA ADS] [CrossRef] [Google Scholar]
 EHT Collaboration, et al. 2019c, ApJ, 875, L3 (Paper III) [NASA ADS] [CrossRef] [Google Scholar]
 EHT Collaboration, et al. 2019d, ApJ, 875, L4 (Paper IV) [NASA ADS] [CrossRef] [Google Scholar]
 EHT Collaboration, et al. 2019e, ApJ, 875, L5 (Paper V) [NASA ADS] [CrossRef] [Google Scholar]
 EHT Collaboration, et al. 2019f, ApJ, 875, L6 (Paper VI) [Google Scholar]
 Falcke, H., & Biermann, P. L. 1995, A&A, 293, 665 [NASA ADS] [Google Scholar]
 Falcke, H., Melia, F., & Agol, E. 2000, ApJ, 528, L13 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962 [NASA ADS] [CrossRef] [Google Scholar]
 Gebhardt, K., Adams, J., Richstone, D., et al. 2011, ApJ, 729, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Goddi, C., Falcke, H., Kramer, M., et al. 2017, Int. J. Mod. Phys. D, 26, 1730001 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, F., Li, H., Daughton, W., & Liu, Y.H. 2014, Phys. Rev. Lett., 113, 155005 [NASA ADS] [CrossRef] [Google Scholar]
 Hada, K., Doi, A., Kino, M., et al. 2011, Nature, 477, 185 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Hada, K., Kino, M., Doi, A., et al. 2013, ApJ, 775, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Hada, K., Kino, M., Doi, A., et al. 2016, ApJ, 817, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Hillas, A. M. 1984, ARA&A, 22, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, M. 1962, IRE Trans. Inf. Theory, 8, 179 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Janssen, M., Goddi, C., van Bemmel, I. M., et al. 2019, A&A, 626, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python [Online] [Google Scholar]
 Junor, W., Biretta, J. A., & Livio, M. 1999, Nature, 401, 891 [NASA ADS] [CrossRef] [Google Scholar]
 Kerr, R. P. 1963, Phys. Rev. Lett., 11, 237 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018, A&A, 616, A188 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 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]
 Leung, P. K., Gammie, C. F., & Noble, S. C. 2011, ApJ, 737, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Levinson, A., & Cerutti, B. 2018, A&A, 616, A184 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Löhner, R. 1987, Comput. Methods Appl. Mech. Eng., 61, 323 [Google Scholar]
 Ly, C., Walker, R. C., & Wrobel, J. M. 2004, AJ, 127, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Markoff, S., Nowak, M., Young, A., et al. 2008, ApJ, 681, 905 [NASA ADS] [CrossRef] [Google Scholar]
 Massaglia, S., Bodo, G., Rossi, P., Capetti, S., & Mignone, A. 2016, A&A, 596, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mei, S., Blakeslee, J. P., Côté, P., et al. 2007, ApJ, 655, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Mattia, G., Bodo, G., & Del Zanna, L. 2019, MNRAS, 486, 4252 [NASA ADS] [CrossRef] [Google Scholar]
 Millman, K. J., & Aivazis, M. 2011, Comput. Sci. Eng., 13, 9 [CrossRef] [Google Scholar]
 Mills, B. Y. 1952, Nature, 170, 1063 [NASA ADS] [CrossRef] [Google Scholar]
 Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, A&A, 586, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mościbrodzka, M., Dexter, J., Davelaar, J., & Falcke, H. 2017, MNRAS, 468, 2214 [NASA ADS] [CrossRef] [Google Scholar]
 Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69 [NASA ADS] [Google Scholar]
 Noble, S. C., Leung, P. K., Gammie, C. F., & Book, L. G. 2007, Classical Quantum Gravity, 24, S259 [NASA ADS] [CrossRef] [Google Scholar]
 Noble, S. C., Krolik, J. H., & Hawley, J. F. 2010, ApJ, 711, 959 [NASA ADS] [CrossRef] [Google Scholar]
 Ohsuga, K., Mineshige, S., Mori, M., & Kato, Y. 2009, PASJ, 61, L7 [NASA ADS] [Google Scholar]
 Oliphant, T. E. 2007, Comput. Sci. Eng., 9, 10 [CrossRef] [PubMed] [Google Scholar]
 Olivares, H., Porth, O., Davelaar, J., et al. 2019, A&A, 629, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owen, F. N., Eilek, J. A., & Kassim, N. E. 2000, ApJ, 543, 611 [NASA ADS] [CrossRef] [Google Scholar]
 Palenzuela, C., Lehner, L., Reula, O., & Rezzolla, L. 2009, MNRAS, 394, 1727 [NASA ADS] [CrossRef] [Google Scholar]
 Pandya, A., Zhang, Z., Chandra, M., & Gammie, C. F. 2016, ApJ, 822, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Lett., 122, 035101 [NASA ADS] [CrossRef] [Google Scholar]
 Petropoulou, M., Giannios, D., & Sironi, L. 2016, MNRAS, 462, 3325 [NASA ADS] [CrossRef] [Google Scholar]
 Pierrard, V., & Lazar, M. 2010, Sol. Phys., 267, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Porth, O., Fendt, C., Meliani, Z., & Vaidya, B. 2011, ApJ, 737, 42 [Google Scholar]
 Porth, O., Olivares, H., Mizuno, Y., et al. 2017, Comput. Astrophys. Cosmol., 4, 1 [Google Scholar]
 Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Prieto, M. A., FernándezOntiveros, J. A., Markoff, S., Espada, D., & GonzálezMartín, O. 2016, MNRAS, 457, 3801 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, Q., Fendt, C., Noble, S., & Bugli, M. 2017, ApJ, 834, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, Q., Fendt, C., & Vourellis, C. 2018, ApJ, 859, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Reynolds, C. S., Fabian, A. C., Celotti, A., & Rees, M. J. 1996, MNRAS, 283, 873 [NASA ADS] [Google Scholar]
 Rezzolla, L., & Zanotti, O. 2013, Relativistic Hydrodynamics (Oxford: Oxford University Press) [CrossRef] [Google Scholar]
 Ripperda, B., Porth, O., Sironi, L., & Keppens, R. 2019, MNRAS, 485, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Roelofs, F., Janssen, M., Natarajan, I., et al. 2019, A&A, submitted [Google Scholar]
 Ryan, B. R., Ressler, S. M., Dolence, J. C., Gammie, C., & Quataert, E. 2018, ApJ, 864, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Sano, T., Inutsuka, S.I., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Sironi, L., & Spitkovsky, A. 2014, ApJ, 783, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Sironi, L., Petropoulou, M., & Giannios, D. 2015, MNRAS, 450, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2012, J. Phys.: Conf. Ser., 372, 012040 [CrossRef] [Google Scholar]
 Tsallis, C. 1988, J. Stat. Phys., 52, 479 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Turland, B. D. 1975, MNRAS, 170, 281 [NASA ADS] [CrossRef] [Google Scholar]
 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
 Walker, R. C., Ly, C., Junor, W., & Hardee, P. J. 2008, J. Phys. Conf. Ser., 131, 012053 [NASA ADS] [CrossRef] [Google Scholar]
 Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Walsh, J. L., Barth, A. J., Ho, L. C., & Sarzi, M. 2013, ApJ, 770, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Watson, M., & Nishikawa, K.I. 2010, Comput. Phys. Commun., 181, 1750 [NASA ADS] [CrossRef] [Google Scholar]
 Werner, G. R., & Uzdensky, D. A. 2017, ApJ, 843, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Werner, G. R., Uzdensky, D. A., Cerutti, B., Nalewajko, K., & Begelman, M. C. 2016, ApJ, 816, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840 [NASA ADS] [CrossRef] [Google Scholar]
 Xiao, F. 2006, Plasma Phys. Controlled Fusion, 48, 203 [Google Scholar]
Appendix A: Phenomenological model explaining the dominance of the counterjet at 228 GHz
In certain GRMHDbased 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 “counterjet”, 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 blackhole’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).
Fig. A.1. Schematic diagram of the model. A diagram of our phenomenological model of the dominant counterjet. The central circle represents the blackhole’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 r_{g}. 
Two main parameters define the geometry of the rings: θ_{offset} and r_{ring} determine their principal diameter and distance from the equatorial plane. A stepfunctionbased 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 r_{g}, and zero otherwise. In other words, the crosssectional thickness of the ring is 0.5 r_{g}.
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 r_{ring}. The effect of this velocity vector is to cause the characteristic relativisticboosting 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 r_{ring}. 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 blackhole’s photon ring, potentially causing the observer to overestimate the blackhole shadow size. Although the blackhole mass can be derived from the size of the blackhole shadow, such an estimate should be seen as an upper limit in the present context.
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, r_{ring} = 6 r_{g}. 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. 
Fig. A.3. Varying r_{ring}. Illustration of the effect of changing r_{ring}, 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 r_{ring} = 4 r_{g}. 
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, r_{ring} = 6 r_{g}). 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 counterjet in optically thin simulations of M 87 dominates over the observerfacing 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.
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 blackhole’s event horizon is marked by the darkgrey circle, while the unstablephoton region is marked by the lightgrey 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 faceon 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 hotspot”) will show the same behavior; most of the radiation will come from the hotspot 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 r_{out} = 3333 r_{g}. 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, ${\mathrm{\Phi}}_{\mathrm{BH}}/\sqrt{\dot{M}}\approx 1.0$. Therefore, the CKS model belongs to the category of Standard and Normal Evolution (SANE; Tchekhovskoy et al. 2012).
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}_{\mathrm{MRI}}^{(\theta )}$ 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.
Fig. B.2. MRI quality factor ${Q}_{\mathrm{MRI}}^{(\theta )}$. Shown is the xy plane (left) and on the xz plane (right). ${Q}_{\mathrm{MRI}}^{(\theta )}$ is sufficient in most of the disk, with a drop visible only very close to the event horizon. 
All Tables
All Figures
Fig. 1. Example of the κdistribution (orange) for a κ value of 3.5, and for comparison a Maxwell–Jüttner distribution (black), and a powerlaw distribution (yellow) with p = 2.5. The κdistribution function is a combination of a powerlaw with a thermal core. 

In the text 
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 powerlaw index p = κ − 1. We overplotted contours of constant κ in black. 

In the text 
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 b^{2}/ρ, overplotted 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 
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 
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), magneticfield strength (middle) and electron temperature (right). Black lines correspond to jet averaged quantities, dashed yellow lines to diskaveraged quantities and the red dashed lines correspond to powerlaw profiles predicted in analytical works. Also, the jet sheath is close to isothermality. 

In the text 
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 
Fig. 7. Spectral energy distributions. Shown are the thermaljet (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 
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 thermaljet 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 
Fig. 9. Logarithmic opticaldepth maps. Single snapshot at 228 GHz of the models at an inclination of i = 160°. 

In the text 
Fig. 10. Core shift. RA offset from the 43 GHz core as a function of frequency. Orange triangles correspond to a thermaljet, black dots to a κjet, grey line represent the observational fit r_{RA}(ν) = (1.4 ± 0.16)ν^{−0.94 ± 0.09} to the M 87 core by Hada et al. (2011). 

In the text 
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 thermaljet 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 
Fig. A.1. Schematic diagram of the model. A diagram of our phenomenological model of the dominant counterjet. The central circle represents the blackhole’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 r_{g}. 

In the text 
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, r_{ring} = 6 r_{g}. 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 
Fig. A.3. Varying r_{ring}. Illustration of the effect of changing r_{ring}, 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 r_{ring} = 4 r_{g}. 

In the text 
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, r_{ring} = 6 r_{g}). 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 
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 blackhole’s event horizon is marked by the darkgrey circle, while the unstablephoton region is marked by the lightgrey 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 
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 
Fig. B.2. MRI quality factor ${Q}_{\mathrm{MRI}}^{(\theta )}$. Shown is the xy plane (left) and on the xz plane (right). ${Q}_{\mathrm{MRI}}^{(\theta )}$ 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.