Issue 
A&A
Volume 660, April 2022



Article Number  A107  
Number of page(s)  22  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202142295  
Published online  20 April 2022 
Impact of nonthermal particles on the spectral and structural properties of M87
^{1}
Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA
^{2}
Institut für Theoretische Physik, Goethe Universität, MaxvonLaueStr. 1, 60438 Frankfurt, Germany
email: cfromm@th.physik.unifrankfurt.de, osorio@itp.unifrankfurt.de
^{3}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
^{4}
TsungDao Lee Institute and School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, PR China
^{5}
Department of Physics, National and Kapodistrian University of Athens, Panepistimiopolis, 15783 Zografos, Greece
^{6}
Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey RH5 6NT, UK
^{7}
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
^{8}
Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010 6500 GL Nijmegen, The Netherlands
^{9}
Department of Astronomy and Columbia Astrophysics Laboratory, Columbia University, 550 W 120th St, New York, NY 10027, USA
^{10}
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
^{11}
Jodrell Bank Centre for Astrophysics, University of Manchester, Machester M13 9PL, UK
^{12}
Frankfurt Institute for Advanced Studies, RuthMoufangStrasse 1, 60438 Frankfurt, Germany
^{13}
School of Mathematics, Trinity College, Dublin 2, Ireland
Received:
24
September
2021
Accepted:
15
January
2022
Context. The recent 230 GHz observations of the Event Horizon Telescope are able to image the innermost structure of M 87 and show a ringlike structure that agrees with thermal synchrotron emission generated in a torus surrounding a supermassive black hole. However, at lower frequencies, M 87 is characterised by a largescale and edgebrightened jet with clear signatures of nonthermal emission. In order to bridge the gap between these scales and to provide a theoretical interpretation of these observations, we perform general relativistic magnetohydrodynamic simulations of accretion onto black holes and jet launching.
Aims. M 87 has been the target for multiple observations across the entire electromagnetic spectrum. Among these, very large baseline interferometry (VLBI) observations provide unique details of the collimation profile of the jet down to several gravitational radii. We aim to model the observed broadband spectrum of M 87 from the radio to the nearIR regime and at the same time, fit the jet structure as observed with global millimeterVLBI at 86 GHz.
Methods. We used general relativistic magnetohydrodynamics and simulated the accretion of the magnetised plasma onto Kerr black holes in 3D. The radiative signatures of these simulations were computed taking different electron distribution functions into account, and a detailed parameter survey was performed in order to match the observations.
Results. The results of our simulations show that magnetically arrested disks around fastspinning black holes (a_{⋆} ≥ 0.5) together with a mixture of thermal and nonthermal particle distributions are able to simultaneously model the broadband spectrum and the innermost jet structure of M 87.
Key words: black hole physics / magnetohydrodynamics (MHD) / accretion / accretion disks / radiative transfer / radiation mechanisms: nonthermal / globular clusters: individual: M87
© ESO 2022
1. Introduction
The recent observations of the elliptical galaxy M87 at 230 GHz performed by the Event Horizon Telescope (EHT) using the very long baseline interferometry (VLBI) technique are consistent with a supermassive black hole (SMBH) with mass ${M}_{\mathrm{BH}}=6.{5}_{0.7}^{+0.7}\times {10}^{9}\phantom{\rule{0.166667em}{0ex}}{M}_{\odot}$ at the centre of the galaxy (Event Horizon Telescope Collaboration 2019a). The compact radio source that was captured by the EHT revealed an asymmetric ring with a diameter d = 42 ± 3 μas, which is assumed to be located below the actual base of a relativistic jet.
The M 87 jet has been systematically monitored from radio to γrays (see e.g., Reid et al. 1982; Hada et al. 2013; Kim et al. 2018a; Snios et al. 2019; MAGIC Collaboration 2020). During the EHT 2017 observations, several international facilities in space and on the ground partnered for a quasisimultaneous multiwavelength campaign of the M 87 nuclei (Algaba & Anczarski 2021). These efforts can reveal the underlying jetlaunching and particle acceleration physics. Moreover, multiwavelength observations can shed light on theoretical expectations of a black hole jet and a magnetised disk wind (Blandford & Znajek 1977; Blandford & Payne 1982; Blandford et al. 2019; LyndenBell 2006).
More specifically, the jet collimation profile can be revealed from multifrequency (10–230 GHz) VLBI observations. These observations have shown that up to hundreds of Schwarzschild radii^{1}, the M87 jet has a parabolic shape (Asada & Nakamura 2012; Doeleman et al. 2012; Hada et al. 2013). The shape and width of a jet can be best described by its opening angle, and for the M87 case, the jet base has revealed a rather wide structure (Kim et al. 2018a). The EHT observations of M 87 show a ringlike emission structure in agreement with theoretical models of accretion onto an SMBH (Event Horizon Telescope Collaboration 2019b). The key physical ingredients are the spin of the SMBH, the flow properties, and the magnetic field, which is thought to be responsible for jetlaunching. A theoretical model needs to capture the compact ringlike image (on horizon scales) and fit the spectrum, but also account for the morphology and kinematics of the observed jet (Mertens et al. 2016; Walker et al. 2018). To solve this multiparametric problem, works have developed semianalytical techniques to fit the spectrum and model the basic jet characteristics (Di Matteo et al. 2003; Broderick & Loeb 2009). To draw a selfconsistent picture of accretion onto a black hole, general relativistic magnetohydrodynamic (GRMHD) simulations together with general relativistic radiative transfer (GRRT) calculations (see e.g., Dexter et al. 2012) are employed. The characteristics of the M87 radio core, a flat spectrum and an increasing size with wavelength, can be reproduced by a twotemperature accretion flow and a hot singletemperature jet (Mościbrodzka et al. 2016). The addition of nonthermal electrons is expected to better reproduce the nearIR (NIR) flux and the flat radio spectrum, but also acquire a more extended jet structure at 43 GHz and 86 GHz (Davelaar et al. 2019).
We study the impact of a nonthermal emission model on the structure and morphology of the M87 jet and its spectrum. We extended the study presented in CruzOsorio et al. (2022) and included additional models. Furthermore, we perform a detailed parameter survey to understand the impact of the various emission parameters on the resulting broadband spectrum and image structures. To this end, we perform longterm, highresolution 3D GRMHD simulations. The two important parameters of these simulations are the black hole spin, a_{⋆} = cJ/GM^{2} with angular momentum J, and the magnetic flux across the horizon, ϕ_{BH}. Depending on the magnetic flux, the models can be divided into standard and normal evolution (SANE) ($\mathrm{\Psi}={\varphi}_{\mathrm{BH}}/\sqrt{\dot{M}}<10$, where Ṁ is the mass accretion rate) and magnetically arrested disks (MAD) (Ψ≥15). Our definition of the MAD parameter Ψ differs by a factor $\sqrt{4\pi}$ from the one defined in Tchekhovskoy et al. (2011). We use SANE and MAD models with spins a_{⋆} = −0.9375, −0.50, 0, 0.50, and 0.9375. The output of these simulations is coupled with GRRT calculations, where we employ a mixture of thermal and nonthermal electrons using the kappa distribution (see Pierrard & Lazar 2010 for a review and Davelaar et al. 2019 for the application to M 87). We present a comprehensive analysis focusing on the impact of the different parameters of the emission model on the spectrum and image structure of M 87. More specifically, the free parameters in the emission model are the local ratio of the proton to electron temperature, the maximum magnetisation in the jet spine, and the percentage of the magnetic energy in the kappa model.
The paper is organised as follows: In Sect. 2 we discuss the details of the GRMHD simulations, and in Sect. 3 we describe the GRRT postprocess calculations. In Sect. 4 we analyse the impact of the free parameters of our model on the broadband spectrum and image structure of M 87, followed by a summary and discussion in Sect. 5. Lastly, we present our conclusions and an outlook in Sect. 6. Throughout this work, we use a black hole mass of 6.5 × 10^{9} M_{Å} for M 87 and a distance of 16.8 Mpc.
2. General Relativistic magnetohydrodynamic simulations
We used the latest version of the GRMHD code called black hole accretion code, BHAC (Porth et al. 2017; Olivares et al. 2019) to simulate accretion onto black holes and jet launching. BHAC solves the 3D GRMHD equations written in a coordinate basis (t, x^{i}) on a generic spacetime fourmetric (g_{μν}) with metric determinant g, while keeping the magnetic field divergence free. In this work, Greek indices run through [0, 1, 2, 3], while Roman indices cover [1, 2, 3]. For the purpose of this work, we chose Kerr black holes in spherical coordinates. The equations solved in BHAC in geometric units (GM = c = 1 and $1/\sqrt{4\pi}$ is absorbed in the magnetic field) are the conservation of mass, the local conservation of energymomentum, and covariant Maxwell equations,
$$\begin{array}{c}\hfill {\mathrm{\nabla}}_{\mu}(\rho {u}^{\mu})=0,\phantom{\rule{1em}{0ex}}{\mathrm{\nabla}}_{\mu}{T}^{\mu \nu}=0,\phantom{\rule{0.166667em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}{\mathrm{\nabla}}_{\mu}^{\ast}\phantom{\rule{0.166667em}{0ex}}{F}^{\mu \nu}=0,\end{array}$$(1)
where ρ is the restmass density, and u^{μ} is the fourvelocity. The stressenergy tensor for ideal MHD is defined as
$$\begin{array}{c}\hfill {T}^{\mu \nu}=(\rho {h}_{\mathrm{tot}}){u}^{\mu}{u}^{\nu}+(p+\frac{1}{2}{b}^{2}){g}^{\mu \nu}{b}^{\mu}{b}^{\nu},\end{array}$$(2)
where h_{tot} = h + b^{2}/ρ is the total specific enthalpy, and p is the fluid pressure. The magnetic field fourvector b^{μ} is given by
$$\begin{array}{c}\hfill {b}^{t}=\frac{\mathrm{\Gamma}\left({B}^{i}{u}_{i}\right)}{\alpha}\phantom{\rule{1em}{0ex}}{b}^{i}=\frac{{B}^{i}+\alpha {b}^{t}{u}^{i}}{\mathrm{\Gamma}},\end{array}$$(3)
where B^{i} is the magnetic field measured by an Eulerian observer, α is the lapse function, and Γ is the Lorentz factor. Using the definitions above, the diagnostics used in this paper, that is, plasma beta β = p/b^{2} and magnetisation σ = b^{2}/ρ, are computed using
$$\begin{array}{c}\hfill {b}^{2}={b}^{\mu}{b}_{\mu}=\frac{{B}^{2}+{\alpha}^{2}{({b}^{0})}^{2}}{{\mathrm{\Gamma}}^{2}}=\frac{{B}^{2}}{{\mathrm{\Gamma}}^{2}}+{({B}^{i}{v}_{i})}^{2},\end{array}$$(4)
where B^{2} = B^{i}B_{i}, where B^{i} is the magnetic field threevector in the Eulerian frame. For the closure of the system of conservation laws, we assumed an equation of state for an ideal gas, connecting the specific enthalpy, h, to the pressure, p, and density, ρ,
$$\begin{array}{c}\hfill h=1+\frac{\widehat{\gamma}p}{(\widehat{\gamma}1)\rho},\end{array}$$(5)
where $\widehat{\gamma}$ is the adiabatic index (see Antón et al. 2006; LoraClavijo et al. 2015; Porth et al. 2017; CruzOsorio et al. 2020 for more details of the evolution equations). We solved the GRMHD equations on a spherical polar grid (r, θ, ϕ) where the grid spacing is logarithmic in radial and linear in polar and azimuthal direction. The dimensions and resolution of the grid are presented in Table 1.
Dimension and resolution of the numerical grid used for the GRMHD simulations.
We used inflow boundary condition and outflow boundary conditions in a radial coordinate, solid reflective wall along the polar boundaries (Olivares et al. 2019), and in the azimuthal direction, we employed periodic boundary conditions to all physical quantities across the cells at ϕ = 0. To solve the GRMHD equations, we employed finitevolume and highresolution shockcapturing methods, the Laxvan Leer Friedrichs (LLF) flux formula in combination with piecewise parabolic method (PPM) reconstruction, and a twostep predictorcorrector scheme. At the same time, the magnetic fields were evolved with the staggered upwind constrained transport scheme (Olivares et al. 2019). We initialised our 3D GRMHD simulations with a magnetised torus in hydrodynamic equilibrium with the central rotating Kerr BH. The torus has a constant specific angular momentum and a weak poloidal magnetic field, to which a single loop is added at the top of the torus, defined by the vector potential
$$\begin{array}{c}\hfill {A}_{\varphi}\propto \mathrm{max}(q0.2,0),\end{array}$$(6)
where
$$\begin{array}{cc}& q=\frac{\rho}{{\rho}_{\mathrm{max}}}{\left(\frac{r}{{r}_{\mathrm{in}}}\right)}^{3}{sin}^{3}\theta exp\left(\frac{r}{400}\right)\phantom{\rule{1em}{0ex}}\mathrm{for}\phantom{\rule{0.166667em}{0ex}}\mathrm{MAD}\hfill \end{array}$$(7)
$$\begin{array}{cc}\hfill & q=\frac{\rho}{{\rho}_{\mathrm{max}}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\mathrm{for}\phantom{\rule{0.166667em}{0ex}}\mathrm{SANE}\hfill \end{array}$$(8)
(see e.g., Fishbone & Moncrief 1976; Font & Daigne 2002; Rezzolla & Zanotti 2013).
The initial conditions of the GRMHD simulations were the constant specific angular momentum, l, adiabatic index of the fluid $\widehat{\gamma}$, plasma beta, β = 2p/b^{2}, inner r_{in} radius and radius of the density maximum in the torus, r_{c}, as well as the dimensionless spin parameter a_{⋆}. In Table 2 we list the values we used, where model names MT.M.* refer to MAD states and MT.S.* to SANE states. To trigger the magnetorotational instability (MRI), we perturbed the equilibrium torus by a white noise to the fluid pressure so that its resulting value was p = p(1 + X_{p}) with a random number X_{p}< 0.04.
Physical parameters of the equilibriummagnetised torus for the MAD and SANE accretion models (for details, see text).
During the course of the GRMHD simulations, we monitored the mass accretion rate, Ṁ, and the accreted magnetic flux across the horizon, ϕ_{BH}. In Fig. 1 we present Ṁ, ϕ_{BH} and the MAD flux parameter $\mathrm{\Psi}={\varphi}_{\mathrm{BH}}/\sqrt{\dot{M}}$ for the MAD models (top three panels) and for the SANE models (bottom three panels). Since the inner edge of the MAD torus is located at a larger radii from the black hole than in the SANE configurations (see Table 2), the saturation times are different between the models (t_{sat} ∼ 10 000 M for MADs and t_{sat} ∼ 7000 M for SANEs, see Fig. 1). When the simulations obtain a quasistable mass accretion rate (t ∈ [8000 M, 10 000 M] for the SANE models and t ∈ [13 000 M, 15 000 M] for the MAD models), we computed the average mass accretion rate, ⟨Ṁ⟩, average magnetic flux across the horizon, ⟨ϕ_{BH}⟩, and the average MAD parameter ⟨Ψ⟩ (see Table 3). The average mass accretion rate approaches a constant value ⟨Ṁ_{MAD}⟩ ∼ 5 and ⟨Ṁ_{SANE}⟩ ∼ 0.3 independently of the BH spin. On the other hand, the average magnetic flux and average MAD parameter increase with spin and exhibit a maximum for a BH spin of a = +1/2 for the SANE and MAD models. All MAD models except for model MT.M.1 reach the MAD state, according to the threshold value of Ψ ∼ 10 (Tchekhovskoy et al. 2011). Models MT.S.* are well within the lowmagnetisation regime ⟨Ψ⟩≤1 and agree with results presented in Porth et al. (2019).
Fig. 1. Mass accretion rates Ṁ, magnetic flux ϕ, and MAD flux parameter Ψ in code units for black holes with different dimensionless spins. The first three panels correspond to the MAD state and the bottom three panels to the SANE simulations. The scales for the MAD and SANE models are different. 
Averaged accretion rates and their standard deviations for different spins a_{⋆} and accretion models.
In Fig. 2 we show the time^{2} and azimuthally averaged largescale morphology of the plasma magnetisation, σ = b^{2}/ρ, plasma beta, β, gas temperature, T_{code} = p/ρ, and Lorentz factor, Γ, for a BH with spin a = +15/16. The top row corresponds to a SANE model and the bottom row to a MAD model. The morphology of our simulations can be divided into two distinct regions: the jet and the disk. We used the Bernoulli parameter −hu_{t} < 1.02 to separate disk (bound plasma) from the jet region (unbound plasma). The jet region can be further divided into a jet spine and a jet sheath region. The former is characterised by a high magnetisation σ > σ_{cut}^{3} and low plasma beta (β ≪ 1). In Fig. 2 the solid black lines indicate log_{10}σ = −1.0, 0, and 0.5. Within the jet spine, the hot plasma is accelerated along the poles, reaching Lorentz factors Γ ∼ 10 in the MAD case and Γ ∼ 5 for the SANE models. We define as the jet sheath the region where σ < σ_{cut}, while the Bernoulli parameter Be = −hu_{t} > 1.02, indicating outflowing plasma (dashed white line in Fig. 2). In the jet sheath, we find typical Lorentz factors of Γ ∼ 8 for MAD models and Γ ∼ 5 for the SANE models. The regions outside of these boundaries are attributed to the disk and disk wind. In Fig. 2 the main differences between MAD and SANE models are clearly visible: (i) MAD models show larger jet opening angles than SANE models (compare the area contained within the dashed white lines and jet axis (x = 0) in Fig. 2) (ii) MAD models show higher magnetisation in the jet spine and sheath than SANE models (see the first and fifth panel in Fig. 2) (iii) MAD models produce hotter jets, that is, the temperature in the jet region is higher than in SANE models (see the third and seventh panel in Fig. 2) (iv) MAD models generate faster jets than their SANE counterparts (see the fourth and eight panel in Fig. 2)
Fig. 2. From left to right, the panels show the time and azimuthally averaged magnetisation, σ, plasma β, gas temperature (proton) in code units, T_{code}, and Lorentz factor Γ for a BH with a = +15/16. All parameters are plotted in logarithmic scale. The top row shows these quantities for the SANE accretion, and the bottom row shows them for the MAD case. The dashed white line represents the boundary between bounded (Be < 1) and unbounded (Be > 1.02) plasma defined through the Bernoulli parameter Be, and the black contour lines log_{10}σ = −1.0, 0, and 0.5 correspond to jet sheath and jet spine boundaries. 
3. General relativistic radiative transfer calculations
In order to compare our GRMHD models with singledish and VLBI observations of M 87, we need to compute their radiative signatures across the electromagnetic spectrum. The covariant radiative transfer equation can be written as
$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathcal{I}}{\mathrm{d}{\tau}_{\nu}}=\mathcal{I}+\frac{\eta}{\chi},\end{array}$$(9)
where ℐ is the Lorentzinvariant specific intensity and is connected to the specific intensity, I_{ν}, by
$$\begin{array}{c}\hfill \mathcal{I}=\frac{{I}_{\nu}}{{\nu}^{3}}.\end{array}$$(10)
The Lorentzinvariant emissivity, η, and absorptivity, χ, are related to the emission, j_{ν} and absorption coefficients, α_{ν} evaluated at frequency ν via
$$\begin{array}{c}\hfill \eta ={j}_{0,\nu}/{\nu}^{2}\phantom{\rule{2em}{0ex}}\chi ={\alpha}_{0,\nu}\nu ,\end{array}$$(11)
where the subscript 0 indicates quantities measured in the local rest frame of the plasma. Using the definition of the optical depth dτ_{ν} = ∫α_{ν}ds, the equation for the covariant radiative transport can be decoupled into two differential equations,
$$\begin{array}{c}\hfill \frac{\mathrm{d}{\tau}_{\nu}}{\mathrm{d}\lambda}={\gamma}^{1}{\alpha}_{0,\nu},\phantom{\rule{1em}{0ex}}\frac{\mathrm{d}\mathcal{I}}{\mathrm{d}\lambda}={\gamma}^{1}\left(\frac{{j}_{0,\nu}}{{\nu}^{3}}\right)exp({\tau}_{\nu}),\end{array}$$(12)
with affine parameter λ and energy shift between the oberver’s and the comoving frame γ^{−1} = ν_{0}/ν = −k_{α}u^{α}_{λ}/k_{β}u^{ β}_{∞}, where k_{α} is the wave vector of the photon. For more details about the radiative transfer scheme, see Younsi et al. (2012). The emission and absorption coefficients depend on the assumed emission process and electron distribution function (see e.g., Pandya et al. 2016).
The GRRT equations (Eq. (12)) were solved along the geodesic trajectories in the fastlight approximation, assuming that the dynamical time of the GRMHD simulations is greater than the lightcrossing time using the code called black hole observations in stationary spacetimes (BHOSS) of Younsi et al. (2020). Within BHOSS, the null geodesics (propagation of the electromagnetic radiation) are solved using a RungeKuttaFehlberg (RKF45) integrator with fourthorder adaptive step sizing and fifthorder error control. The radiative transfer equations were evolved using an Eulerian method where the step size is provided from the geodesic integrator. Throughout our GRRT calculations, we used a field of view of 10^{3} M (which corresponds to 4 mas for the black hole mass and distance of M 87) with a resolution of 800 × 800 pixels. The accuracy of the geodesic integration within the RKF45 method was set to Δϵ = 10^{−14}, and an optical depth cut of τ_{ν, cut} = 5 was used to speed up the GRRT calculations. The GRMHD simulations described in the previous section only provide the proton temperature in the plasma. Therefore, the properties of the radiating electrons, especially the electron temperature, T_{e}, need to be reconstructed from plasma quantities. Here, we followed the work of Mościbrodzka et al. (2016) and Event Horizon Telescope Collaboration (2019b) and computed the electron temperature via the socalled Rβ model, which assumes a plasmaβ depending temperature ratio between protons and electrons,
$$\begin{array}{c}\hfill {\mathrm{\Theta}}_{\mathrm{e}}=\frac{p{m}_{\mathrm{p}}/{m}_{\mathrm{e}}}{\rho {T}_{\mathrm{ratio}}},\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{1em}{0ex}}{T}_{\mathrm{ratio}}\equiv \frac{{T}_{\mathrm{p}}}{{T}_{\mathrm{e}}}=\frac{{R}_{\mathrm{low}}+{R}_{\mathrm{high}}{\beta}^{2}}{1+{\beta}^{2}},\end{array}$$(13)
where m_{p} and m_{e} the proton and electron masses, respectively. In addition to the plasmaβ, which is provided from our GRMHD simulations, two additional free parameters are involved in the calculation of the temperature ratio: R_{low} and R_{high}.
To better understand the impact of R_{low} and R_{high} on the electron temperature, we evaluated Eq. (13) for a wide range of plasmaβ (10^{−3} − 10^{3}, a typical range found in our GRMHD simulations) for several selected values for R_{low} and R_{high}. The result of this calculation can be found in Fig. 3. We plot the electrontoproton temperature ratio, that is, (T_{p}/T_{e})^{−1}. The solid lines correspond to R_{low} = 1.0 and the dotted lines to R_{low} = 0.1. The colour of the lines indicates different R_{high} values from black (R_{high} = 1) to green (R_{high} = 160). The vertical dashed line marks the transition between jet region (β < 1) and disk regions (β ≥ 1) ( see also the second and sixth panel in Fig. 2). The variation in the electrontoproton temperature ratio with plasmaβ resembles the shape of a sigmoidcurve, that is, it saturates for low plasmaβ to the value of 1/R_{low} and for high plasmaβ to the value given by 1/R_{high}. In other words, increasing the R_{high} value while keeping the R_{low} value fixed reduces the electron temperature in the disk regions, but does not alter the electron temperature in the jet.
Fig. 3. Electrontoproton temperature ratio T_{e}/T_{p} as a function of plasma β varying the R_{low} and R_{high} (definition given below in Eq. (13)). The vertical dashed black line marks the transition between the jet sheath region (β < 1) and the disk region (β > 1 ). The solid lines correspond to models with R_{low} = 1, and dotted lines show models with R_{low} = 0.1. The different colours indicate different R_{high} values. 
After the electron temperature, Θ_{e}, was computed, we needed to assign an emission process and the electron distribution function (eDF). These three parameters together with the plasma properties (e.g., the magnetic field) define the emission, j_{0, ν}, and absorption, α_{0, ν}, coefficients required for radiative transfer (see. Eq. (12)). As mentioned earlier, we assumed synchrotron emission, that is, the gyration of relativistic electrons around the magnetic field line, as the main radiative process responsible for the generation of the observed emission in M 87 from the radio to NIR range (see e.g., Yuan & Narayan 2014). Additional radiation processes such as Bremsstrahlung can provide additional seed photons that are important during inverse Compton scattering and thus might be responsible for the highenergy emission in M 87 (see e.g., Yarza et al. 2020). However, given their contributions to the frequency regime under investigation in this paper (10^{9} Hz ≤ 10^{16} Hz), these processes can be neglected.
The two eDFs we explore in this work are the MaxwellJüttner distribution function for thermal electrons and the kappa distribution function for nonthermal electrons. The application of the former is well justified by the modelling of the horizon scale emission of M 87 which is in good agreement with a thermal eDF (Event Horizon Telescope Collaboration 2019b). The use of the kappa eDF is motivated from the distribution of the spectral index, α, computed from VLBI observations of M 87 at 22 GHz, 43 GHz, and 86 GHz, which exhibit values between −1 ≤ α ≤ 0 within a distance of r < 4 mas from the core (see Fig. 15 in Algaba & Anczarski 2021). These values indicate a powerlaw electron distribution with 2 ≤ p ≤ 4. In addition, Davelaar et al. (2019) showed that the M 87 spectrum between (10^{9} Hz ≤ 10^{16} Hz), consisting of a flat radio spectrum and a steep NIR part, can be well reproduced by a kappa eDF. In the following paragraph, we provide details about the two eDFs and how we connect them to our GRMHD simulations (see Xiao 2006; Pandya et al. 2016; Davelaar et al. 2019 for further details). The MaxwellJüttner distribution function is given by
$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{e}}}{\mathrm{d}{\gamma}_{\mathrm{e}}}=\frac{{n}_{\mathrm{e}}}{4\pi {\mathrm{\Theta}}_{\mathrm{e}}}\frac{{\gamma}_{\mathrm{e}}\sqrt{{\gamma}_{\mathrm{e}}^{2}1}}{{K}_{2}(1/{\mathrm{\Theta}}_{\mathrm{e}})}exp(\frac{{\gamma}_{\mathrm{e}}}{{\mathrm{\Theta}}_{\mathrm{e}}}),\end{array}$$(14)
where n_{e} is the electron number density, γ_{e} is the electron Lorentz factor, and K_{2} is the Bessel functions of second kind. The second eDF is the kappa distribution, which smoothly connects a thermal distribution for small γ_{e} with a powerlaw tail for larger γ_{e}. The kappa distribution can be written as
$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{e}}}{\mathrm{d}{\gamma}_{\mathrm{e}}}=\frac{N}{4\pi}{\gamma}_{\mathrm{e}}\sqrt{{\gamma}_{\mathrm{e}}^{2}1}{(1+\frac{{\gamma}_{\mathrm{e}}1}{\kappa w})}^{(\kappa +1)},\end{array}$$(15)
where N is a normalisation factor (see Pandya et al. 2016, for details). The powerlaw exponent of a nonthermal particle distribution $\mathrm{d}{n}_{\mathrm{e}}/\mathrm{d}{\gamma}_{\mathrm{e}}\propto {\gamma}_{\mathrm{e}}^{s}$ is related to the κ value by s = κ − 1. The energy of the kappa distribution is set by its width, w. We followed Davelaar et al. (2019) and attributed in addition to the thermal energy also a fraction ε of the magnetic energy to the width of the kappa distribution,
$$\begin{array}{c}\hfill w:=\frac{\kappa 3}{\kappa}{\mathrm{\Theta}}_{\mathrm{e}}+\frac{\epsilon}{2}[1+tanh(r{r}_{\mathrm{inj}})]\phantom{\rule{0.166667em}{0ex}}\frac{\kappa 3}{6\kappa}\frac{{m}_{\mathrm{p}}}{{m}_{\mathrm{e}}}\sigma .\end{array}$$(16)
In the equation above, r_{inj} corresponds to the injection radius, that is, the distance from which we start injecting electrons with a magnetic energy contribution (see also Davelaar et al. 2019). Throughout this work, we set r_{inj} = 10 M, which is in agreement with the jet stagnation surface (u^{r}=0) typically located between 5 M and 10 M (Nakamura et al. 2018). We assumed magnetic reconnection as the main particle acceleration mechanism. Because this acceleration mechanism acts on scales that we cannot resolve in our global GRMHD simulations, we again followed Davelaar et al. (2019) and employed the particleincell simulationbased subgrid model from Ball et al. (2018). More precisely, we used the magnetisation and plasmaβ dependent powerlaw slope of the electron distribution and connected this to the κ value,
$$\begin{array}{c}\hfill \kappa :=2.8+0.7{\sigma}^{1/2}+3.7{\sigma}^{0.19}tanh(23.4{\sigma}^{0.26}\beta ).\end{array}$$(17)
Given the width, w, of the kappa distribution (see Eq. (16)) together with the numerical approximations for the emission, j_{ν}, and absorption coefficients, α_{ν} (see Pandya et al. 2016), the values for κ are restricted to the interval 3 < κ ≤ 8.
In Fig. 4 we compare the MaxwellJüttner and the kappa electron distribution functions. The solid black lines corresponds to a MaxwellJüttner eDF with an electron temperature of Θ_{e} = 10. As expected from Eq. (14), the MaxwellJüttner eDF decreases exponentially for large electron Lorentz factors, γ_{e}. The kappa eDF for κ = 3.5 and a width of w = 10 is plotted as the dashed red curve. The powerlaw tail of the kappa eDF is clearly visible for γ_{e} > 10^{3} and is well matched by a powerlaw with exponent s = 2.5 shown as the solid blue curve. When the width of the kappa eDF is fixed while the κ value is increased, the kappa eDF approximates a MaxwellJüttner eDF. This behaviour is clearly illustrated by the dashed orange curve, corresponding to a kappa eDF with w = 10 and κ = 10^{6}. Using Eq. (17) together with typical values for the magnetisation, σ, and plasmaβ in the jet sheath from our GRMHD simulations (see Fig. 2), we find 3 ≤ κ ≤ 7. The corresponding kappa eDFs using a fixed width of w = 10 are shown as the cyan region in Fig. 4.
Fig. 4. Comparison of eDFs: MaxwellJüttner electron distribution (black), powerlaw electron distribution (blue) with s = 2.5, κ electron distribution for κ = 3.5 and w = 10 (red), and for κ = 10^{6} and w = 10 (orange). For more details, see text. 
In order to interpret the results of the GRRT calculations, it is important to understand effects of the adjustable parameters of the temperature description, for instance R_{high}, in combination with the plasma properties (σ, plasmaβ) on the distribution of the defining parameters of the eDFs, that is, the electron temperature, Θ_{e}, the width, w, and the κ value. In Fig. 5 we show the distribution of the electron temperature, Θ_{e}, for models MT.S.5 (top) and MT.M.5 (bottom) for three different values of R_{high} while keeping R_{low} = 1 fixed. The distribution of the electron temperature can be understood in the following way: According to Eq. (13), the electron temperature depends on the code temperature T_{code} ∝ p/ρ and on the temperature ratio, T_{ratio} = 1 + R_{high}β^{2}/(1 + β^{2}) setting R_{low} = 1. The distribution of the code temperature and the plasmaβ are shown for models MT.S.5 and MT.M.5 in Fig. 2. From the distribution of the code temperature, it is clear that the highest electron temperatures will be located in the jet region and that the electron temperature in the jet will be higher for the MAD model than for the SANE model (compare the top and bottom panels in Fig. 5). The variation of the electron temperature with the R_{high} value can be explained by the distribution of the plasmaβ (see second and sixth panel in Fig. 2) and the behaviour of the electrontoproton temperature ratio plotted in Fig. 3. Independent of the accretion model, low plasmaβ values are found in the jet region and high values in the wind and disk regions. Thus, increasing the R_{high} value will decrease the electrontoproton temperature ratio in disk and wind regions while keeping the ratio in the jet region nearly unchanged (see the solid lines in Fig. 3). This behaviour can be best seen in the bottom row in Fig. 5 by the decrease in electron temperature in the disk region (100r_{g} ≤ x ≤ 300r_{g} and −100r_{g} ≤ z ≤ 100r_{g}).
Fig. 5. Logarithm of the electron temperature, Θ_{e}, for models MT.S.5 (top) and MT.M.5 (bottom) using different R_{high} values while keeping R_{low} = 1 fixed. 
In Fig. 6 we analyse the variation and distribution of the κ value and the width, w, similarly to the electron temperature. The κ value depends on the magnetisation, σ, and on the plasmaβ (see Eq. (17)) and their distributions for models MT.S.5 and MT.M.5 are shown in Fig. 2. As mentioned earlier, MAD models exhibit higher magnetisation in the jet region than their SANE counterparts. In combination with Eq. (17), this leads to lower κ values in the jet for the MAD models than in the SANE models. Lower κ values correspond to flatter slopes in the highenergy part of the kappa eDF and thus to more particles with large electron Lorentz factors, γ_{e} (see Fig. 4). The second important parameter for the kappa eDF is its width, w. The width is computed from the electron temperature, Θ_{e}, the κ value together with some fraction ε of the magnetic energy (see Eq. (16)). Since the fraction of the magnetic energy is typically < 1, the distribution of w mainly follows the distribution of the electron temperature (compare Figs. 5 and 6). It is important to mention again that we only applied the kappa eDF in regions where 3 < κ ≤ 8 due to the limitation of the numerical approximations for the emission and absorption coefficients. Outside of this range, we switched back to a MaxwellJüttner distribution and to the corresponding emission and absorption coefficients (Pandya et al. 2016). This implies that the kappa eDF is mainly used in the jet sheath region.
Fig. 6. Logarithm of the κ value (first and fifth panel) and the width of the kappa eDF, w, for models MT.S.5 (top) and MT.M.5 (bottom) using different R_{high} values while keeping R_{low} = 1, r_{inj} = 10 M, and ε = 0.5 fixed. 
4. Parameter space study
Given our large set of GRMHD models, that is, accretion type and spin (see Table 2), together with the different choices of the electron microphysics and hyperparameters, for example, σ_{cut}, during the GRRT, it is necessary to explore the parameter space in a robust and agnostic manner. However, given the computational cost, we restricted ourselves to a small set of parameters that allowed us to explore the hypersurfaces of the parameter space in adequate and scientifically sufficient way while keeping the computational cost within an affordable range. In addition to the different spins and accretion models (MAD and SANE), we list in Table 4 the free GRRT parameters of our study and their explored range.
GRRT parameters and their values explored during the parameter space study.
During the GRRT, we computed both the broadband spectrum from 10^{9} Hz to 10^{16} Hz and raytraced 86 GHz images. This allowed us to compare the spectral properties of our models with multifrequency observations of M 87 (see Table B.1) and at the same time, contrast the 86 GHz jet structure with the structure obtained via global millimetre VLBI (GMVA) observations of M 87 (Kim et al. 2018a). As mentioned earlier, the GRMHD simulations are scale free, and we applied a black hole mass of 6.5 × 10^{9} M_{⊙} and a distance of 16.8 Mpc (Event Horizon Telescope Collaboration 2019a) during the radiative transfer to adjust our simulations to M 87. In addition, we used a fixed inclination angle of 160° and iterated the mass accretion rate, ṁ, for each model until we obtained an average flux density of 1.0 Jy at 230 GHz (Doeleman et al. 2012; Akiyama et al. 2015) within a time window of 2000 M. This time window was taken between 13 000 M ≤ t ≤ 15 000 M for the MAD models and between 8 000 M ≤ t ≤ 10 000 M for the SANE models, ensuring that the GRMHD simulations are well within the quasistable state (see Fig. 1).
In the following, we investigate the influence of (i) the electron temperature via R_{low, high} (ii) the electron distribution function, that is, the MaxwellJüttner (thermal) and kappa distribution (nonthermal) (iii) the jet spine via the magnetisation cutoff σ_{cut} (iv) the fraction of the magnetic energy contributing to the total energy contained in the kappa distribution via ε (v) the black hole spin a_{⋆} on the spectral and the structural properties of our GRRT calculations. For the comparison between the theoretical and observed structure of M 87 at 86 GHz, we used the superresolved stacked 86 GHz GMVA image (presented in Kim et al. 2018a). In order to mimic the GMVA observations, we rotated the average GRRT image to a position angle^{4} of 288° (Walker et al. 2018) and convolved the image with a beam of 0.123 mas × 0.051 mas at a position angle of 0°. After the convolution, we rotated the image by −18°, that is, the jet axis coincided with the xaxis, and sliced the jet perpendicular to the xaxis. During the slicing, we limited the dynamical range to 10^{4} and only considered flux profiles with peak flux five times above the noise level for the jet diameter analysis. We defined as the jet diameter the location where the flux density profile is above 50% of peak flux. The obtained profiles for the jet diameter, D_{jet} and the apparent opening angle, ϕ_{app}, along the jet axis were compared to those from the GMVA observations.
4.1. Influence of the electron temperature
In our first scan of the parameter space, we analysed the influence of the electron temperature, Θ_{e}, on the spectrum and the 86 GHz image structure. The two free parameters defining the electron temperature, Θ_{e}, according to Eq. (13) are R_{high} and R_{low}. In the following section we analyse their influence on the images and spectra in detail.
4.1.1. Influence of the R_{high} parameter
We first focus on the influence of the R_{high} parameter using values of 10, 80, and 160 while keeping R_{low} = 1 fixed. The modelling of 230 GHz EHT observations shown at R_{high} = 1 for the MAD and SANE models disagrees with observational constraints (see Table 3 in Event Horizon Telescope Collaboration 2019b). This statement was obtained for thermal models. However, within our GRRT setup, we only applied the kappaeDF in the jet region. Thus, the disk region includes a thermal eDF as in the case of Event Horizon Telescope Collaboration (2019b). For this reason, we excluded R_{high} = 1 in our parameter sweep. Furthermore, we used σ_{cut} = 1 and three different eDFs, namely a thermal and a kappa eDF with ε = 0 and with ε = 1. Together with the three R_{high} values and the two different accretion models (MAD and SANE), this leads to 18 models explored here, and the results for a black hole with spin a_{⋆} = 0.9375 are presented in Fig. 7.
Fig. 7. Influence of the R_{high} parameter on the spectral and structural properties of our M 87 simulations for a black hole with spin a_{⋆} = 0.9375. The left column corresponds to a MAD model, and the right column shows a SANE model. Top panels: broadband spectrum averaged for 2000 M for R_{high} = 10 (dotted), R_{high} = 80 (dashed), and R_{high} = 160 (solid) while keeping R_{low} = 1 fixed. The different eDFs are marked by different colours: thermal (blue), kappa with ε = 0 (green), and kappa with ε = 1 (red). The middle image shows the raytraced 86 GHz images for three different eDFs (columns) and different R_{high} (rows). Bottom panels: jet diameter and opening angles profiles using the same colour and line styles as for the spectrum. 
The influence of the R_{high} parameter on the broadband spectrum can be best seen for the thermal models (MaxwellJüttner eDF) in the MAD case (top left panel, blue curves in Fig. 7). The highfrequency part of the spectrum (ν > 10^{12} Hz) steepens significantly with increasing R_{high}, while the lowfrequency part (ν ≤ 10^{12} Hz) slightly increases in flux without a change in the spectral slope (see inset in top left panel). Because the highenergy emission is mainly emitted from the disk region predominated by high plasmaβ values, the electron temperature can be written as Θ_{e} ∝ 1/R_{high}^{5} (see Eq. (13)). This means that increasing the R_{high} value decreases the electron temperature, Θ_{e}, in the disk (see Figs. 3 and 5). Using the relation between the electron temperature and the R_{high} parameter, the thermal emissivity in the highfrequency regime is given by
$$\begin{array}{c}\hfill {j}_{\nu ,\mathrm{th}}\propto exp({R}_{\mathrm{high}}^{2/3}{\nu}^{1/3}),\end{array}$$(18)
which finally explains the exponential decrease with frequency and the variation with R_{high} for the highfrequency part of the spectrum.
In a similar manner, we can explain the lowfrequency behaviour of the spectrum. Assuming that the lowfrequency emission is generated in the jet and disk wind regions, which are characterised by low plasmaβ values, the temperature ratio is not affected by the choice of R_{high} and thus the electron temperature is not altered (see Eq. (13) and Figs. 3 and 5).
Next, we compare the eDFs. The kappa eDF was only applied within the jet, and elsewhere, we used a thermal eDF. Thus, we need to consider the emissivities and absorptivities from both eDFs during the radiative transfer depending on the regions that are crossed by the ray. Again, we first focused on the highenergy part of the spectrum and used the approximations of Pandya et al. (2016) for the kappa emissivities. Thus we can write the kappa emissivity as
$$\begin{array}{c}\hfill {j}_{\nu ,\phantom{\rule{0.166667em}{0ex}}\kappa}\propto {\nu}^{(\kappa 2)/2}{(w\kappa )}^{\kappa 2}.\end{array}$$(19)
Inserting the definition of the width of the kappa eDF (Eq. (16)), the kappa emissivity can be approximated by^{6}
$$\begin{array}{c}\hfill {j}_{\nu ,\phantom{\rule{0.166667em}{0ex}}\kappa}\propto {\nu}^{(\kappa 2)/2}{[{\mathrm{\Theta}}_{\mathrm{e}}+\epsilon \sigma ]}^{\kappa 2}.\end{array}$$(20)
Together with the thermal emissivitiy, the total emissivity has the simplified form
$$\begin{array}{c}\hfill {j}_{\nu ,\mathrm{tot}}\propto exp({\mathrm{\Theta}}_{\mathrm{e}}^{2/3}{\nu}^{1/3})+{\nu}^{(\kappa 2)/2}{[{\mathrm{\Theta}}_{\mathrm{e}}+\epsilon \sigma ]}^{\kappa 2}.\end{array}$$(21)
If we consider that the highenergy thermal emission originates mainly from the disk (high plasmaβ region) and that the jet sheath^{7} has intermediate plasmaβ values, the equation above can be further simplified and the dependence on the R_{high} can be included,
$$\begin{array}{c}\hfill {j}_{\nu ,\mathrm{tot}}\propto exp({R}_{\mathrm{high}}^{2/3}{\nu}^{1/3})+{\nu}^{(\kappa 2)/2}{[1/{R}_{\mathrm{high}}+\epsilon \sigma ]}^{\kappa 2}.\end{array}$$(22)
With this simplification of the emissivity at hand, we can understand the spectral behaviour in the highfrequency regime for the thermal–nonthermal eDFs. First we concentrate on the ε = 0 models (green curves in the top panels of Fig. 7). Including the kappa eDF leads to a powerlaw tail in the highfrequency emission as compared to exponential decay in the thermal model, which is set by the second term in Eq. (22), that is, the contribution of the nonthermal electrons with large γ_{e} within the tail of the kappa eDF. For the ε = 0 models, the R_{high} dependence of the spectrum is still visible, that is, a steepening of the spectrum with increasing R_{high}. This effect is introduced by the R_{high}dependent terms in Eq. (22). For the ε = 1 models, there is no clear R_{high} dependence on the highfrequency spectrum (see the red curves in the top left panel of Fig. 7. This can be explained by the second term in Eq. (22), which includes the contribution from the magnetic field, that is, the nonthermal particles gain additional energy from the magnetic field. With ε = 1, the second term in Eq. (22) is the dominating term in the emissivity, and thus the dependence on R_{high} is no longer visible in the highfrequency emission.
The relations and interpretation of the spectral behaviour above also apply to the SANE models. However, the trends are less clear than in the MAD case. The reason for this is connected to the less magnetised structure of the SANE simulations (see Fig. 2). In addition, the required mass accretion rate during the flux normalisation process varies strongly for the SANE models across the different eDFs and R_{high} values. The turnover frequency, ν_{t,th}, for thermal emission can be approximated as ${\nu}_{t,\mathrm{th}}\propto {\mathrm{\Theta}}_{\mathrm{e}}^{2}B$ (see e.g., Zdziarski et al. 1998). Inserting the scaling relations from code units to cgs units together with the R_{high}dependent electron temperature leads to
$$\begin{array}{c}\hfill {\nu}_{t,\mathrm{th}}\propto \sqrt{{b}^{2}\dot{m}}/{R}_{\mathrm{high}}^{2}.\end{array}$$(23)
Here, $\sqrt{{b}^{2}}$ is the magnetic field strength. Thus, changes in the mass accretion rate and R_{high} value lead to shifts in the turnover position. The strong variations in the mass accretion rate are not found for the MAD models, which explains the nearly identical turnover positions in the MAD cases. The turnover frequency for the SANE models with R_{high} = 10 is located around 10^{11} Hz, compared to ∼10^{12} Hz for higher values of R_{high}. Taking the dependence of the turnover frequency on ṁ and R_{high} in account together with Eqs. (18) and (22), the spectral behaviour seen in the SANE models can be explained in a similar manner as for the MAD models (top right panel in Fig. 7).
So far, our discussion was focused on the spectrum, and in the following, we elaborate on the influence of the R_{high} parameter and the eDF on the 86 GHz image structure. In the middle panel of Fig. 7 we present the 86 GHz images from our models for M 87 for different R_{high} values (columns) and different eDFs (rows). The overall structure of the 86 GHz images can be described by a bright innermost region and a fainter jet. The innermost emission (r < 0.1 mas for the MAD models) is mainly generated in the disk region and is nearly independent of the choice of R_{high}. In case of the SANE models, the innermost emission for low R_{high} values is produced in the disk, and for high R_{high} values, the counterjet contributes the majority of the emission (see also Fig. 11 in Event Horizon Telescope Collaboration 2019b). On the other hand, the extent of the jet depends on the eDF and the R_{high} parameter. From the spectrum, we can see that all MAD models are optically thin at 86 GHz (see the inset in the left top panel), and we can use the approximation of the thermal emissivity (Eq. (18)) to qualitatively understand the obtained emission structure. The jet emission could be understood as a blend of radiation generated in the jet sheath and disk wind. These two regions are characterised by plasmaβ values ∼0.1 (see Fig. 2), and emission roughly follows Eq. (18). Thus, the shortening of the jet and the steepening of the broadband spectrum with increased R_{high} values can be explained.
Including nonthermal electrons via the kappa eDF in the jet sheath leads to wider and more extended jets (see the first row in the middle panel of Fig. 7). Similar to the change in the spectral behaviour between thermal and nonthermal eDF, we can explain this behaviour by Eq. (22). The addition of nonthermal electrons with large electron Lorentz factors, γ_{e}, increases the emission in the jet sheath at large distances from the black hole. For the model with ε = 0, increasing the R_{high} value leads to a shorter and narrower jet. This is expected from Eq. (22) and is in agreement with our explanation for the spectral properties for this model (green curves in the top left panel of Fig. 7). Adding the magnetic energy of plasma as additional energy source for the nonthermal particles further enhances the emission of the jet on large scales and increases the jet width. Analogously to the discussion of the spectral behaviour for the ε = 1 models, there is no significant dependence on the R_{high} parameter (see the red curves in the top left panel and the right column in the left middle panel in Fig. 7). The 86 GHz image structure of the SANE models (right middle panel in Fig. 7) follows the same discussion as above for the MAD models. The only exception is the 86 GHz thermal image for R_{high} = 10. In this case, it is optically thick at 86 GHz at the source as compared to thermal R_{high} > 10 (see the dotted blue curve in the inset in the top right panel of Fig. 7) and characterised by a lower total flux density than the other models.
In the bottom panels of Fig. 7 we compare the jet diameter and opening angle, ϕ_{app}, from our models with the results from the GMVA observation of M 87. The dotted faint red and black curves correspond to jet diameter fits to the GMVA observations (Kim et al. 2018a). Independent of the electron temperature and eDF, all MAD models exhibit opening angle profiles that agree with the observations, but the opening angle for the SANE models decreases faster than the observations with distance from the black hole. As mentioned above, the inclusion of nonthermal particles slightly increases the jet opening angle and significantly enhances the emission on large scales.
4.1.2. Influence of the R_{low} parameter
The second variable parameter of the electron temperature description (Eq. (13)) is the R_{low} value. In the previous analysis, we used R_{low} = 1. In this section we investigate the changes in the spectrum and the 86 GHz image structure if R_{low} = 0.1 is used and compare the results to the result from the R_{low} = 1 models. The resulting broadband spectrum and 86 GHz images are presented in Fig. 8. The comparison between the R_{low} = 1 (Fig. 7) and the R_{low} = 0.1 models show that the latter ones are brighter at high frequencies ν > 10^{14} Hz and show more extended jets for the thermal and ε = 0 models. In order to understand this behaviour, we need to analyse the dependence of the electron temperature on R_{low} (see Eq. (13)). Reduction in the R_{low} value increases the electron temperature in regions where plasmaβ < 1 while keeping the electron temperature unchanged in high plasmaβ regions (see Fig. 3). This implies that electrons in the jet sheath are hotter than those in R_{low} = 1 models. The heating of the jet sheath is also clearly visible in the 86 GHz images. The jets are more extended and wider than for the models with R_{low} = 1 (compare the bottom panels in Fig. 7 and Fig. 8). The steepening of the spectrum at high frequencies (ν > 10^{14} Hz) with increased R_{high} can be explained by the inverse proportionality of the electron temperature and the R_{high} values (T_{e} ∝ 1/R_{high}) in the disk region ( high plasmaβ region). However, the middle panel in Fig. 8 also shows that the jets are shortened with R_{high}. In case of R_{low} < 1, Eq. (13) can be used and determines for which combination of R_{high} and plasmaβ the electrons in the jet will be cooled, that is, T_{ratio} > 1,
$$\begin{array}{c}\hfill \beta >\sqrt{\frac{1{R}_{\mathrm{low}}}{{R}_{\mathrm{high}}1}}.\end{array}$$(24)
Fig. 8. Influence of the R_{low} parameter on the spectral and structural properties of our M 87 simulations for a black hole with spin a_{⋆} = 0.9375. The left column corresponds to a MAD model, and the right column shows a SANE model. Top panels: broadband spectrum averaged for 2000 M for R_{high} = 10 (dotted), R_{high} = 80 (dashed), and R_{high} = 160 (solid) while keeping R_{low} = 0.1 fixed. The different eDFs are marked by different colours: thermal (blue), kappa with ε = 0 (green), and kappa with ε = 1 (red). The middle image shows the raytraced 86 GHz images for three different eDFs (columns) and different R_{high} (rows). Bottom panels: jet diameter and opening angle profiles using the same colour and line styles as for the spectrum. 
Using our values of R_{high} = 10, 80, and 160, this leads to plasmaβ values of 0.3, 0.1, and 0.07. Because the typical plasmaβ in the jet sheath ∼0.1, we can finally understand the shortening of the jet in the 86 GHz images with increased R_{high}: In order to cool the electrons in the jet sheath if R_{low} = 0.1 and R_{high} = 10, typical plasmaβ values around 0.3 are required. Since the typical values are ∼0.1, the sheath electrons will be heated instead of being cooled, and the jet will be more extended. When we increase the R_{high} values, the required plasmaβ for cooling the sheath electrons drops below the typical plasmaβ found in this region, and the electrons will be cooled. As a result, the jet will be shortened with increased R_{high}. This effect also contributes to the steepening of the broadband spectrum. In the SANE case, only the thermal model for R_{high} = 10, which is optically thick at 86 GHz, does not follow this explanation (see the dotted blue curves in the inset of the top right panel in Fig. 8). Similar as in the previous analysis, the spectrum and image structure of the ε = 1 models are independent of the choice of the R_{high} value.
4.2. Influence of the jet spine
In our next parameter sweep, we explored the influence of the position of the jet spine on the broadband spectrum and the 86 GHz image structure. Because the jet spine is not well defined in the literature, the magnetisation of the plasma above a certain value, σ_{cut}, is commonly used to define the boundary between the jet sheath and jet spine. Here we used the parameters from the first analysis, R_{low} = 1.0, R_{high} = 160, and we applied a nonthermal eDF with ε = 1 in the jet sheath while varying the jet spine location via σ_{cut} = 0.1, 0.5, 1.0, 5.0, and 10.
Figure 9 shows the broadband spectrum (top), the 86 GHz jet structure (middle), and the jet width and opening angle profiles (bottom) for five different values of the σ_{cut} for a MAD model (left) and for a SANE model with black hole spin a_{⋆} = 0.9375. In the broadband spectrum the increase in the σ_{cut} value leads to higher flux density and flatter spectra in the highfrequency regime (ν > 10^{12} Hz). This behaviour could be explained in the following way: Increasing the σ_{cut} will include more magnetised plasma from the jet spine in the jet sheath (see the contour lines in the first and fifth panel Fig. 2). At the same time, these regions have lower plasmaβ, which according to Eq. (13) leads to higher electron temperatures in the jet sheath (see the second and sixth panel in Fig. 2). In addition to the increase in electron temperature in the jet sheath with increased σ_{cut}, the powerlaw slope of the nonthermal eDF will be flattened (see Eq. (17) and Fig. 6). These two effects combined, namely the increased electron temperature and the decrease in the κ value, will lead to a larger width of the kappaeDF (see Eq. (16)), that is, to more energetic nonthermal particles. Because changes in the σ_{cut} value will not alter the temperature of the disk region, the changes in the emission are entirely due to the jet sheath emission (second term in Eq. (21)). From this equation, the increase of the flux density and the flattening of the spectrum with increased σ_{cut} can be explained. For the SANE model, there is nearly no difference between the spectrum for σ_{cut} = 5 and σ_{cut} = 10, which is an indication that the largest magnetisation in this model is around σ ∼ 3.
Fig. 9. Influence of the σ_{cut} parameter on the spectral and structural properties of our M 87 simulations for a black hole with spin a_{⋆} = 0.9375. The left column corresponds to a MAD model, and the right column shows a SANE model. The top panels show the broadband spectrum averaged for 2000 M for six different values of σ_{cut} (colourcoded) while keeping R_{low} = 1, R_{high} = 160, and ε = 1.0 fixed. The middle image shows the raytraced 86 GHz images for different σ_{cut} values. The bottom panels show the jet diameter and the opening angles profiles using the same colourcoding as for the spectrum. 
With increasing σ_{cut} value, the 86 GHz jet becomes brighter at larger distances from the black hole. For the MAD models, the jet width increases (see the bottom left panels) and the emission signature of the jet spine becomes visible as a central ridge for σ_{cut} ≥ 3. In contrast to the MAD models, the SANE models show the largest jet width for σ_{cut} = 0.5, and higher values of σ_{cut} do not affect the jet width and opening angle (see the bottom right panels). As mentioned above, with increased σ_{cut} value, more magnetised and hotter plasma is included in the jet sheath, and at the same time, the powerlaw slope of the nonthermal eDF flattens. Thus, more radiating electrons with large γ_{e} are included in the jet sheath, and the jet becomes brighter on large distances. Similar as in the broadband spectrum, the independence of the jet width and opening angle for σ_{cut} ≥ 3 for the SANE models is a clear signal that we include the entire jet (jet spine and jet sheath) in the GRRT.
4.3. Influence of the magnetic energy
The last free parameter we considered in our radiative transfer calculation is the fraction of magnetic energy, ε, which is added to the nonthermal electrons (see Eq. (16)). So far, we considered two extreme values, namely ε = 0 and ε = 1.0, and we included here three additional values for ε = 0.25, 0.5, and 0.75 and for comparison reasons, a purely thermal model. During this sweep, we used R_{low} = 1.0, R_{high} = 160 and σ_{cut} = 1 in the parameter space, and the results are presented in Fig. 10. The broadband spectrum (top panels) shows that increasing the ε value increases the flux density at high frequencies (ν > 10^{12} Hz) and flattens the spectrum. By increasing the ε value, more energy from the magnetic field is included in the nonthermal particles (see Eq. (16)) and the nonthermal term in the combined emission coefficient becomes larger and dominates for high values of ε (see Eq. (21)). As a result, the flux density at high frequencies (ν > 10^{12} Hz) increases, and the exponential decay (due to the thermal disk emission) is transitioning into a power law. This behaviour is independent of the accretion model. Because of the different turnover positions in the SANE models, this trend is not as clear as in the MAD models. The shift in the turnover position is connected to the variation in mass accretion rate during the flux normalisation (see Eq. (23)), as mentioned before.
Fig. 10. Influence of the ε parameter (fraction of magnetic energy used for the acceleration of nonthermal particles) on the spectral and structural properties of our M 87 simulations for a black hole with spin a_{⋆} = 0.9375. The left column corresponds to a MAD model, and the right column shows a SANE model. Top panels: broadband spectrum averaged for 2000 M for thermal and five different values of ε (colourcoded) while keeping R_{low} = 1, R_{high} = 160, and σ_{cut} = 1.0 fixed. The middle image shows the raytraced 86 GHz images for different ε values. Bottom panels: jet diameter and opening angles profiles using the same colourcoding as for the spectrum. 
Increasing the fraction of magnetic energy in the kappa eDF leads to more extended jets (see the middle panel in Fig. 10). For ε ≥ 0.5, no further changes in the jet width and opening angle are obtained (see the bottom panels in Fig. 10). Increasing ε will increase the number of electrons with large γ_{e}, that is, more highenergy electrons in the tail of the kappa distribution (see Eq. (16)). Because we applied the kappa eDF in the jet sheath region, this region will become brighter. For ε ∼ 0.5, the entire jet sheath region includes electrons with large γ_{e}. The geometry of the jet sheath is defined by the spin of the black hole and the accretion model (see Fig. A.1). Thus, increasing the ε value above 0.5 has no effect on the jet structure and will only increase the emission (best seen at high frequencies in the broadband spectrum).
4.4. Influence of the black hole spin, a_{⋆}
The black hole spin parameter a_{⋆} is a fundamental parameter of our GRMHD simulations because it modifies the event horizon size, the gravitational potential, and characteristic radii such as the last marginally stable Keplerian orbit and innermost stable circular orbit (Font & Daigne 2002; Daigne & Font 2004; Rezzolla & Zanotti 2013). In Table 2 we summarise the black hole spin parameters and the corresponding initial data for the magnetized torus. Figure 1 shows the evolution of the accretion rates and the magnetic flux. The differences between MAD and SANE models are clearly visible, but differences between corotating and counterrotating black holes within the same accretion model are visible as well. This behaviour is best seen in the magnetic flux (third and sixth panel in Fig. 1) and in the average values presented in Table 3.
As mentioned in the last paragraph, the spin of the black hole has major impact on the geometry of the jet. Therefore, we analyse in this section the influence of the black hole spin on the broadband spectrum and the 86 GHz image structure. Figure 11 shows the average spectra and averaged 86 GHz images using R_{low} = 1, R_{high} = 160, and σ_{cut} = 3.0 together with a kappa electron distribution with ε = 0.5 for five different spin values a_{⋆} = −0.9375, −0.5, 0.0, 0.5, and 0.9375. The values for σ_{cut} and ε are motivated from the previous slices through the parameter space and lead to an improved fit to the broadband spectrum (see the top panels in Fig. 11). Interestingly, changes in the black hole spin for this set of bestfit parameter introduce only a minor scatter in the broadband spectrum, especially in the MAD case. The lowfrequency part of the spectrum (ν < 10^{11} Hz) is flatter for the MAD models than for their SANE counterparts (see also the total flux at 86 GHz given in the second panel of Fig. 11). This can be explained the lack of a highly magnetised jet in the SANE models (see Fig. A.1), which leads to higher κ values (see Eq. (17) and Fig. A.2). Therefore, less energetic particles are in the tail of the eDF, leading to lower emissivities (see Eq. (22)), resulting in shorter jets (see the second panel in Fig. 11).
Fig. 11. Influence of the black hole spin, a_{⋆}, on the spectral and structural properties of our M 87 simulations. The left column corresponds to a MAD model, and the right column shows a SANE model. Top panels: broadband spectrum averaged for 2000 M for thermal and five different spin values a_{star} = −0.9375, −0.5, 0, 0.5, and 0.9375 (colourcoded) while keeping R_{low} = 1, R_{high} = 160 σ_{cut} = 3.0, and ε = 0.5 fixed. The middle image shows the raytraced 86 GHz images for different spin values. Bottom panels: jet diameter and opening angles profiles using the same colourcoding as for the spectrum. 
The morphology of the 86 GHz images presented in the middle panel of Fig. 11 shows the already mentioned differences between SANE (top) and MAD models (bottom): MAD models produce wide and highly magnetised jets that are lit up by the inclusion of nonthermal particles via the kappaeDF. Independent of the accretion model, Schwarzschild black holes (a_{⋆} = 0) do not launch a wellestablished jet, and the bulk of the emission is produced in the disk with some minor contribution from the disk wind. In the MAD case, the differences between the jet width and the jet length for the different spins can be related to the variability of the accretion flow and the accreted magnetic flux. Following Tchekhovskoy et al. (2011) and Tchekhovskoy & McKinney (2012), the jet efficiency, η, for MAD models is proportional to ${\varphi}_{\mathrm{BH}}^{2}$. Together with the average values for ϕ_{BH} given in Table 3, the differences in the jet morphology between pro and retrograde models can be explained.
5. Summary and discussion
M87 is the perfect laboratory to challenge our current understanding of jet launching and plasma physics. Therefore, we investigated the models that are able to simultaneously reproduce the observed broadband spectrum and image structure of M 87. The workflow for this analysis can be divided into three main branches: (i) 3D GRMHD simulations, (ii) multifrequency GRRT calculations and (iii) comparison with observations. Because our models are based on initial assumptions and require a set of adjustable parameters, we performed a detailed parameter space survey to better understand their impact on the spectrum and image structure.
5.1. GRMHD simulations
Given the apriori unknown spin and accretion model for M 87, our 3D GRMHD simulations cover five different spin values, including pro and retrograde black holes and two accretion models. This set of numerical simulations allowed us to cover a fraction of the possible configuration of the diskjet system in M 87. The simulations were performed in modified KerrSchild coordinates together with three AMR levels that capture the accretion process close to the black hole with high resolution and ensures that we resolve the fastest growing MRI mode throughout the course of the simulation (see Porth et al. 2019). All of our GRMHD simulations were evolved until a quasi steadystate in the mass accretion is achieved. In the case of the MAD models, the all models, expect for the fastest spinning retrograde black hole, are above the MAD parameter of Ψ > 10. In agreement with the literature, our MAD models show wider and more highly magnetised jets than the SANE models (see e.g., Tchekhovskoy et al. 2011; Chael et al. 2019). These highresolution 3D GRMHD simulation provide accurate plasma dynamics for our radiative transfer calculations and our investigation of the physical conditions in M 87.
5.2. GRRT calculations
We performed a detailed parameter survey to understand the influence of the different parameters, including the eDF, on the broadband spectrum and the 86 GHz image structure. To ensure that the compact flux of our models at 230 GHz agreed with the EHT observations, we iterated the mass accretion rate ṁ to obtain an average flux of 1 Jy at 230 GHz within a time window of 2000 M. In addition to the black hole spin a_{⋆} and accretion model (MAD or SANE), there are five additional parameters, namely, R_{high}, R_{low}, the eDF (thermal or nonthermal), σ_{cut}, and ε. During the parameter survey, we allowed one parameter to vary while keeping the other parameters fixed. In Fig. 12 we present a qualitative summary of the variation in the spectrum and image structure for the different parameters under investigation. In each of the panels, we show the evolution of the spectrum (top) and the image structure for the indicated parameter, that is, outermost contour line (bottom). The black, red, and orange lines indicate a low, medium, and high value for the parameter under investigation and its impact on the spectrum and on the 86 GHz image structure. The separation between the contour lines indicates the sensitivity of the jet structure to changes in a certain parameter. For example, switching from a thermal to a nonthermal eDF significantly alters the jet structure (see panel c), while changes in the R_{high} value for a given eDF modify the jet structure in a less strong way.
Fig. 12. Schematic variation of the broadband spectrum and image structure as a function of various model parameters. 
The distribution of the plasmaβ in the GRMHD simulations together with the R − β description for the temperature ratio (Eq. (13)) leads during the GRRT to lower NIR fluxes for increasing R_{high} (see panel a and Fig 12), while the NIR fluxes increase if R_{low} is decreased. In addition, our parameter space survey shows independently of the spin and accretion model that the eDF has the largest impact on the broadband spectrum (best shown in the slope at high frequencies (ν > 10^{12} Hz) in panel c in Figs. 12 and 7), in agreement with the findings of Davelaar et al. (2019) for a fastspinning SANE model. The second important outcome of our parameter survey is the influence of the σ_{cut} on the spectrum. For both accretion models, excluding most of the magnetised jet region via σ_{cut} < 1 during the GRRT highly underproduces the NIR fluxes, while σ_{cut} > 3 overproduces the NIR flux and leads to too flat spectral slopes (see panel d in Figs. 12 and 9). A similar behaviour was found by Chael et al. (2019) for a MAD model with a_{⋆} = 0.9375 including electron heating. It is interesting to mention that a combination of emission parameters can be found that leads to very similar broadband spectra for different spins (see panel f in Figs. 12 and 11). Thus, we conclude that fitting the broadband spectrum alone will not provide constraints on the black hole spin.
These variations in the broadband spectrum are mainly independent of the black hole spin and accretion model. However, the structure of the jet strongly depends on these two parameters. As mentioned earlier, the black hole spin and the accretion model are initial conditions of the GRMHD simulations and determine the plasma dynamics, including the properties of the jets, for example, the width and its magnetisation. MAD models produce wider and more highly magnetised jets than their SANE counterparts. Thus, in the numerical simulations used in this work, the jet diameter and jet opening angle of MADs will always be larger than for SANE models, regardless of the adjustable parameters during the GRRT (see the bottom panels in Figs. 7–11). Similar, fastrotating black holes generate more powerful jets and more extended jets than slowly rotating black holes (see Fig. 11).
5.3. Comparison with observations
After understanding the impact of the initial conditions of the GRMHD and the free parameters during the GRRT on the broadband spectrum and the image structure, we can finally compare our models to the observations of M 87. The broadband spectrum for M 87 is remarkably well fitted for the MAD models using nonthermal particles via the kappa eDF with ε = 0.5, R_{low} = 1, R_{high} = 160, and σ_{cut} = 3 (see the top left panel in Fig. 11). However, as mentioned above, fitting the SED alone does not provide constraints on the black hole spin. Therefore, we extracted the jet collimation profiles from our 86 GHz images and compared them to those obtained from 86 GHz GMVA observations (Kim et al. 2018a). This analysis shows that the highspin prograde MAD models produce collimation profiles that agree with the currently highest resolution VLBI images^{8} of the jet in M 87 (see the bottom panels in Fig. 11). In Fig. 13 we present our bestfit images convolved with a 0.123 mas × 0.051 mas beam and two slices through the jet at r ∼ 0.6 mas and at r ∼ 0.8 mas, which can be compared to Fig. 4 in Kim et al. (2018a). Both models show a large jet opening angle and a significant signature of a counter jet. Furthermore, the jet slices show edgebrightening, together with a third flux density maximum enclosed by the outer ridges. This is best seen for the MAD a_{⋆} = 0.9375 model at r = 0.8 mas (dotted blue line in the right panel of Fig. 13).
Fig. 13. Best jet models (left panel for MAD a_{⋆} = 0.50 and middle panel for MAD a_{⋆} = 0.9375) convolved with a beam of 0.123 mas × 0.051 mas (plotted as the grey ellipse in the bottom left corner of the panels) in order to mimic the averaged and superresolved 86 GHz GMVA image of Kim et al. (2018a). The right panel presents slices through the jet a r ∼ 0.6 mas (red curves) and at r ∼ 0.8 mas. 
In addition to the convolved images presented in Fig. 13, we list in Table 5 the average magnetic field and bulk Lorentz factors for the bestfit models. Our favoured accretion model (MAD) together with the estimates for the magnetic field and mass accretion rate agree with the results obtained from the modelling of the horizon scale polarisation structure of M 87 (1−30 G and ṁ = 3−20 × 10^{−4} M_{Å} yr^{−1}) (Event Horizon Telescope Collaboration 2021). It is important to mention that these two sets of values are obtained from independent studies using different methods and underlying observations, confirming the physical conditions in the jetlaunching region of M 87 presented here. Furthermore, our values for the bulk Lorentz factors, Γ, roughly agree with the values derived from 43 GHz VLBA and 86 GHz GMVA observations (Mertens et al. 2016; Walker et al. 2018; Kim et al. 2018a). We note the difference in the nomenclature between our work and the observations: because we excluded the jet spine from our GRRT calculations, our definition of the jet sheath is similar to definition of the jet spine in the observations. With this in mind, the observations suggest Γ > 6 in the jet spine (our jet sheath) depending on the viewing angle and Γ ≤ 2 in the jet wind region (Mertens et al. 2016; Walker et al. 2018; Kim et al. 2018a), similar to the listed values in Table 5.
Averaged quantities for the bestfit MAD models.
5.4. Limitations of the models
We used ideal GRMHD, which does not include resistivity or heating or cooling of the plasma. Including resistivity could lead to the production of hotter plasmoids in the disk and jet sheath region, which could lead to more pronounced edgebrightening effects (Ripperda et al. 2020) and stronger disk winds (Vourellis et al. 2019). By using electron heating models such as turbulent and reconnection heating, the electron temperature can be directly obtained from the GRMHD simulations and the postprocessing step via Eq. (13) is not required (see e.g., Chael et al. 2019).
Despite these limitations, we do not expect large impacts on the results presented here: Nathanail et al. (2020), Dihingia et al. (2021), and Chashkina et al. (2021) showed that the formation of plasmoids using ideal GRMHD is similar to those from resistive GRMHD. In addition, the recent work by Mizuno et al. (2021) revealed that the R − β model for the electron temperature is well matched to heating models for R_{low} = 1 and R_{high} ≤ 160. The numerical setup applied in this work consists of a torus in hydrodynamical equilibrium seeded with a weak poloidal magnetic field, β ∼ 10^{2}, (see Table 2 for details). However, recent works of nonideal GRMHD including a mean field dynamo showed that poloidal magnetic fields can be selfsimilarly generated from an almost zero initial magnetic field (plasma beta 10^{6} − 10^{9}) (Tomei et al. 2020; Vourellis & Fendt 2021). Future studies should overcome the current limitations by including resistivity, electron heating, and dynamos in the disks, which will provide deeper insights into the jetlaunching process while reducing the number of free parameters in the modelling.
6. Conclusion and outlook
We presented longterm and highresolution 3D GRMHD simulations of standard and normal evolution disks and magnetically arrested disk accretion around a spinning black hole to investigate the launching of the M 87 jet. From the 3D GRMHD simulations, we computed the broadband spectrum (10^{9} ≤ ν ≤ 10^{16} Hz) and the jet collimation profiles. During the GRRT, we investigated the influence of the electron temperature via the R − β model, the electron distribution function (thermal and nonthermal particles), and the jet sheath on the radiative signatures. From our detailed parameter survey, we found that the best fit to the broadband spectrum (see Table B.1 for the used observational data) is obtained if nonthermal particles energetically supported by the magnetic field (ε = 0.5) are included in a magnetised jet sheath (σ_{cut} = 3). Interestingly, the broadband spectrum could equally well be fitted independent of the spin of the black hole and the accretion model. To break this degeneracy in the models, we included the jet collimation profile in our analysis and compared it to the observed profile extracted from 86 GHz GMVA observations (Kim et al. 2018a). The combined constraint of broadband spectrum and collimation profiles favoured the MAD accretion model together with fastspinning prograde black holes (a_{⋆} ≥ 0.5). We showed that our results agree with the values obtained from polarisation modelling and from estimates extracted from VLBI observations of M 87. In addition to the provided bestfit models, one of the main results of this work is that the combined modelling of broadband spectrum and additional constraints such as the image structure can help to further constrain possible models for M 87.
Future work of modelling the jetlaunching and propagation in M 87 should include more advanced, nonstandard GRMHD simulations including resistivity and electron heating models. In order to further investigate the physical conditions in M 87, additional observational constraints such as the polarisation structure and highenergy radiation (Xrays) should be computed during the radiative transfer. Finally, to reduce the influence of the hyperparameter included in the image reconstructions and to account for realistic observational conditions (e.g., sparse sampling of the u − v plane and antenna calibration uncertainties), a direct comparison between synthetic and measured visibilities should be performed. However, such a detailed direct fitting of the GRMHD simulations to the observations is computationally challenging and requires specialised numerical methods and codes that currently under development (see e.g., Fromm et al. 2019).
Commonly, a value of σ_{cut} = 1 (see e.g., Event Horizon Telescope Collaboration 2019b) is used.
Acknowledgments
We thank George Wong for careful reading and useful comments to the manuscript. C.M.F. is supported by the Black Hole Initiative at Harvard University, which is supported by a grant from the John Templeton Foundation. This research is supported by the DFG research grant “Jet physics on horizon scales and beyond” (Grant No. FR 4069/21), ERC synergy grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (Grant No. 610058) and ERC advanced grant JETSET: Launching, propagation and emission of relativistic jets from binary mergers and across mass scales (Grant No. 884631) A.N. was supported by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the 2nd Call for H.F.R.I. Research Projects to support PostDoctoral Researchers (Project Number: 00634). Z.Y. is supported by a UKRI Stephen Hawking Fellowship and acknowledges support from a Leverhulme Trust Early Career Fellowship. J.D. is supported by NASA grant NNX17AL82 and a Joint Columbia/Flatiron Postdoctoral Fellowship. Research at the Flatiron Institute is supported by the Simons Foundation. The simulations were performed on GOETHEHLR LOEWE at the CSCFrankfurt, Iboga at ITP Frankfurt and Pi2.0 at Shanghai Jiao Tong University. Software: BHAC (https://bhac.science/) (Porth et al. 2017), BHOSS (Younsi et al. 2020), ehtim (https://achael.github.io/ehtimaging/) (Chael et al. 2018)
References
 Akiyama, K., Lu, R.S., Fish, V. L., et al. 2015, ApJ, 807, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Algaba, J. C., Anczarski, J., et al. 2021, ApJ, 911, L11 [NASA ADS] [CrossRef] [Google Scholar]
 An, T., Sohn, B. W., & Imai, H. 2018, Nat. Astron., 2, 118 [Google Scholar]
 Antón, L., Zanotti, O., Miralles, J. A., et al. 2006, ApJ, 637, 296 [Google Scholar]
 Asada, K., & Nakamura, M. 2012, ApJ, 745, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Ball, D., Sironi, L., & Özel, F. 2018, ApJ, 862, 80 [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [CrossRef] [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A. E., & Loeb, A. 2009, ApJ, 697, 1164 [NASA ADS] [CrossRef] [Google Scholar]
 Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23 [Google Scholar]
 Chael, A., Narayan, R., & Johnson, M. D. 2019, MNRAS, 486, 2873 [NASA ADS] [CrossRef] [Google Scholar]
 Chashkina, A., Bromberg, O., & Levinson, A. 2021, MNRAS, 508, 1241 [NASA ADS] [CrossRef] [Google Scholar]
 CruzOsorio, A., GimenoSoler, S., & Font, J. A. 2020, MNRAS, 492, 5730 [NASA ADS] [CrossRef] [Google Scholar]
 CruzOsorio, A., Fromm, C. M., Mizuno, Y., et al. 2022, Nature, 6, 103 [NASA ADS] [Google Scholar]
 Daigne, F., & Font, J. A. 2004, MNRAS, 349, 841 [NASA ADS] [CrossRef] [Google Scholar]
 Davelaar, J., Olivares, H., Porth, O., et al. 2019, A&A, 632, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dexter, J., McKinney, J. C., & Agol, E. 2012, MNRAS, 421, 1517 [NASA ADS] [CrossRef] [Google Scholar]
 Dihingia, I. K., Vaidya, B., & Fendt, C. 2021, MNRAS, 505, 3596 [NASA ADS] [CrossRef] [Google Scholar]
 Di Matteo, T., Allen, S. W., Fabian, A. C., Wilson, A. S., & Young, A. J. 2003, ApJ, 582, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Doeleman, S. S., Fish, V. L., Schenck, D. E., et al. 2012, Science, 338, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019a, ApJ, 875, L1 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019b, ApJ, 875, L5 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2021, ApJ, 910, L13 [Google Scholar]
 Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962 [NASA ADS] [CrossRef] [Google Scholar]
 Font, J. A., & Daigne, F. 2002, ApJ, 581, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Fromm, C. M., Younsi, Z., Baczko, A., et al. 2019, A&A, 629, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hada, K., Kino, M., Doi, A., et al. 2013, ApJ, 775, 70 [Google Scholar]
 Hada, K., Park, J. H., Kino, M., et al. 2017, PASJ, 69, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018a, A&A, 616, A188 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kim, J.Y., Lee, S.S., Hodgson, J. A., et al. 2018b, A&A, 610, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12 [CrossRef] [Google Scholar]
 LoraClavijo, F. D., CruzOsorio, A., & Guzmán, F. S. 2015, ApJS, 218, 24 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D. 2006, MNRAS, 369, 1167 [NASA ADS] [CrossRef] [Google Scholar]
 MAGIC Collaboration (Acciari, V. A., et al.) 2020, MNRAS, 492, 5354 [NASA ADS] [CrossRef] [Google Scholar]
 Mertens, F., Lobanov, A. P., Walker, R. C., & Hardee, P. E. 2016, A&A, 595, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mizuno, Y., Fromm, C. M., Younsi, Z., et al. 2021, MNRAS, 506, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, A&A, 586, A38 [Google Scholar]
 Nakamura, M., Asada, K., Hada, K., et al. 2018, ApJ, 868, 146 [Google Scholar]
 Nathanail, A., Fromm, C. M., Porth, O., et al. 2020, MNRAS, 495, 1549 [Google Scholar]
 Olivares, H., Porth, O., Davelaar, J., et al. 2019, A&A, 629, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pandya, A., Zhang, Z., Chandra, M., & Gammie, C. F. 2016, ApJ, 822, 34 [Google Scholar]
 Perlman, E. S., Sparks, W. B., Radomski, J., et al. 2001, ApJ, 561, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Pierrard, V., & Lazar, M. 2010, Sol. Phys., 267, 153 [NASA ADS] [CrossRef] [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 [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]
 Reid, M. J., Schmitt, J. H. M. M., Owen, F. N., et al. 1982, ApJ, 263, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Rezzolla, L., & Zanotti, O. 2013, Relativistic Hydrodynamics (Oxford, UK: Oxford University Press) [CrossRef] [Google Scholar]
 Ripperda, B., Bacchini, F., & Philippov, A. A. 2020, ApJ, 900, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Snios, B., Nulsen, P. E. J., Kraft, R. P., et al. 2019, ApJ, 879, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Tchekhovskoy, A., & McKinney, J. C. 2012, MNRAS, 423, L55 [NASA ADS] [CrossRef] [Google Scholar]
 Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Tomei, N., Del Zanna, L., Bugli, M., & Bucciantini, N. 2020, MNRAS, 491, 2346 [Google Scholar]
 Vourellis, C., & Fendt, C. 2021, ApJ, 911, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Vourellis, C., Fendt, C., Qian, Q., & Noble, S. C. 2019, ApJ, 882, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128 [Google Scholar]
 Whysong, D., & Antonucci, R. 2004, ApJ, 602, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Xiao, F. 2006, Plasma Phys. Controlled Fusion, 48, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Yarza, R., Wong, G. N., Ryan, B. R., & Gammie, C. F. 2020, ApJ, 898, 50 [Google Scholar]
 Younsi, Z., Wu, K., & Fuerst, S. V. 2012, A&A, 545, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Younsi, Z., Porth, O., Mizuno, Y., Fromm, C. M., & Olivares, H. 2020, in Perseus in Sicily: From Black Hole to Cluster Outskirts, eds. K. Asada, E. de Gouveia Dal Pino, M. Giroletti, H. Nagai, & R. Nemmen, 342, 9 [Google Scholar]
 Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Zdziarski, A. A., Poutanen, J., Mikolajewska, J., et al. 1998, MNRAS, 301, 435 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Averaged jet morphologies
In Fig. A.1 we show the azimuthal and time averaged distribution of the logarithm of the magnetisation, σ, for five different spins (from left to right: 0.9375, 0.50, 0, 0.50, and 0.9375) for the SANE (top) and MAD (bottom) models. The dashed white line corresponds to −hu_{t} = 1.02 and separates the outflowing plasma (region enclosed by the jet axis and the hu_{t} contour line) and bound plasma, and the black lines show log_{10}(σ) = −1.0, 0.5, and 0 contours. As mentioned in Sect. 2, the jet regions in the MAD models exhibit higher magnetisation than the SANE models.
Fig. A.1. Azimuthal and timeaveraged distribution of the logarithm of the magnetisation, σ, for different spins (left to right) and for the SANE (top) and MAD (bottom) models. The dashed white line corresponds to the Bernoulli parameter, −hu_{t} = 1.02, and the black lines show log_{10}(σ) = −1.0, 0.5, and 0 contours. 
The averaged distribution of the κ parameter for different spins and two different accretion models are presented in Fig. A.2. Similar to Fig. A.1, the dashed white line corresponds to the Bernoulli parameter, −hu_{t} = 1.02, and the black lines show log_{10}(σ) = −1.0, 0.5, and 0 contours. Since the κ value depends on the magnetisation and the plasmaβ, the distribution of the κ values roughly follows the distribution of the magnetisation (see Eq. 17). We find low κ values, corresponding to flat particle distributions in the MAD models and steep particle distributions in the SANE case (see Eq. 15). Thus, the MAD models have more electrons with large Lorentz factors, γ_{e} in the jet region than their SANE counterparts. The emission and absorption coefficients derived by Pandya et al. (2016) for the kappa eDF are only valid for κ ≤ 8, and we switched to a thermal eDF in these regions during the GRRT.
Fig. A.2. Azimuthal and timeaveraged distribution of the κ value for different spins (left to right) and for the SANE (top) and MAD (bottom) models. The dashed white line corresponds to the Bernoulli parameter, −hu_{t} = 1.02, and the black lines show log_{10}(σ) = −1.0, 0.5, and 0 contours. 
Appendix B: Multifrequency observations of M 87
In Table B.1 we provide the observational data with the references we used for the broadband spectrum in Figs. 7 11.
Observed flux densities for M 87. Columns show the observing frequency in Hz, the total flux density, and the uncertainty, both in Jansky, details of the observations, and references. The first part of the table corresponds to M 87 during a quiescent phase, and the second part lists values during an active phase of M 87. The third group lists values taken during 2005 and the last one for 2017.
All Tables
Physical parameters of the equilibriummagnetised torus for the MAD and SANE accretion models (for details, see text).
Averaged accretion rates and their standard deviations for different spins a_{⋆} and accretion models.
Observed flux densities for M 87. Columns show the observing frequency in Hz, the total flux density, and the uncertainty, both in Jansky, details of the observations, and references. The first part of the table corresponds to M 87 during a quiescent phase, and the second part lists values during an active phase of M 87. The third group lists values taken during 2005 and the last one for 2017.
All Figures
Fig. 1. Mass accretion rates Ṁ, magnetic flux ϕ, and MAD flux parameter Ψ in code units for black holes with different dimensionless spins. The first three panels correspond to the MAD state and the bottom three panels to the SANE simulations. The scales for the MAD and SANE models are different. 

In the text 
Fig. 2. From left to right, the panels show the time and azimuthally averaged magnetisation, σ, plasma β, gas temperature (proton) in code units, T_{code}, and Lorentz factor Γ for a BH with a = +15/16. All parameters are plotted in logarithmic scale. The top row shows these quantities for the SANE accretion, and the bottom row shows them for the MAD case. The dashed white line represents the boundary between bounded (Be < 1) and unbounded (Be > 1.02) plasma defined through the Bernoulli parameter Be, and the black contour lines log_{10}σ = −1.0, 0, and 0.5 correspond to jet sheath and jet spine boundaries. 

In the text 
Fig. 3. Electrontoproton temperature ratio T_{e}/T_{p} as a function of plasma β varying the R_{low} and R_{high} (definition given below in Eq. (13)). The vertical dashed black line marks the transition between the jet sheath region (β < 1) and the disk region (β > 1 ). The solid lines correspond to models with R_{low} = 1, and dotted lines show models with R_{low} = 0.1. The different colours indicate different R_{high} values. 

In the text 
Fig. 4. Comparison of eDFs: MaxwellJüttner electron distribution (black), powerlaw electron distribution (blue) with s = 2.5, κ electron distribution for κ = 3.5 and w = 10 (red), and for κ = 10^{6} and w = 10 (orange). For more details, see text. 

In the text 
Fig. 5. Logarithm of the electron temperature, Θ_{e}, for models MT.S.5 (top) and MT.M.5 (bottom) using different R_{high} values while keeping R_{low} = 1 fixed. 

In the text 
Fig. 6. Logarithm of the κ value (first and fifth panel) and the width of the kappa eDF, w, for models MT.S.5 (top) and MT.M.5 (bottom) using different R_{high} values while keeping R_{low} = 1, r_{inj} = 10 M, and ε = 0.5 fixed. 

In the text 
Fig. 7. Influence of the R_{high} parameter on the spectral and structural properties of our M 87 simulations for a black hole with spin a_{⋆} = 0.9375. The left column corresponds to a MAD model, and the right column shows a SANE model. Top panels: broadband spectrum averaged for 2000 M for R_{high} = 10 (dotted), R_{high} = 80 (dashed), and R_{high} = 160 (solid) while keeping R_{low} = 1 fixed. The different eDFs are marked by different colours: thermal (blue), kappa with ε = 0 (green), and kappa with ε = 1 (red). The middle image shows the raytraced 86 GHz images for three different eDFs (columns) and different R_{high} (rows). Bottom panels: jet diameter and opening angles profiles using the same colour and line styles as for the spectrum. 

In the text 
Fig. 8. Influence of the R_{low} parameter on the spectral and structural properties of our M 87 simulations for a black hole with spin a_{⋆} = 0.9375. The left column corresponds to a MAD model, and the right column shows a SANE model. Top panels: broadband spectrum averaged for 2000 M for R_{high} = 10 (dotted), R_{high} = 80 (dashed), and R_{high} = 160 (solid) while keeping R_{low} = 0.1 fixed. The different eDFs are marked by different colours: thermal (blue), kappa with ε = 0 (green), and kappa with ε = 1 (red). The middle image shows the raytraced 86 GHz images for three different eDFs (columns) and different R_{high} (rows). Bottom panels: jet diameter and opening angle profiles using the same colour and line styles as for the spectrum. 

In the text 
Fig. 9. Influence of the σ_{cut} parameter on the spectral and structural properties of our M 87 simulations for a black hole with spin a_{⋆} = 0.9375. The left column corresponds to a MAD model, and the right column shows a SANE model. The top panels show the broadband spectrum averaged for 2000 M for six different values of σ_{cut} (colourcoded) while keeping R_{low} = 1, R_{high} = 160, and ε = 1.0 fixed. The middle image shows the raytraced 86 GHz images for different σ_{cut} values. The bottom panels show the jet diameter and the opening angles profiles using the same colourcoding as for the spectrum. 

In the text 
Fig. 10. Influence of the ε parameter (fraction of magnetic energy used for the acceleration of nonthermal particles) on the spectral and structural properties of our M 87 simulations for a black hole with spin a_{⋆} = 0.9375. The left column corresponds to a MAD model, and the right column shows a SANE model. Top panels: broadband spectrum averaged for 2000 M for thermal and five different values of ε (colourcoded) while keeping R_{low} = 1, R_{high} = 160, and σ_{cut} = 1.0 fixed. The middle image shows the raytraced 86 GHz images for different ε values. Bottom panels: jet diameter and opening angles profiles using the same colourcoding as for the spectrum. 

In the text 
Fig. 11. Influence of the black hole spin, a_{⋆}, on the spectral and structural properties of our M 87 simulations. The left column corresponds to a MAD model, and the right column shows a SANE model. Top panels: broadband spectrum averaged for 2000 M for thermal and five different spin values a_{star} = −0.9375, −0.5, 0, 0.5, and 0.9375 (colourcoded) while keeping R_{low} = 1, R_{high} = 160 σ_{cut} = 3.0, and ε = 0.5 fixed. The middle image shows the raytraced 86 GHz images for different spin values. Bottom panels: jet diameter and opening angles profiles using the same colourcoding as for the spectrum. 

In the text 
Fig. 12. Schematic variation of the broadband spectrum and image structure as a function of various model parameters. 

In the text 
Fig. 13. Best jet models (left panel for MAD a_{⋆} = 0.50 and middle panel for MAD a_{⋆} = 0.9375) convolved with a beam of 0.123 mas × 0.051 mas (plotted as the grey ellipse in the bottom left corner of the panels) in order to mimic the averaged and superresolved 86 GHz GMVA image of Kim et al. (2018a). The right panel presents slices through the jet a r ∼ 0.6 mas (red curves) and at r ∼ 0.8 mas. 

In the text 
Fig. A.1. Azimuthal and timeaveraged distribution of the logarithm of the magnetisation, σ, for different spins (left to right) and for the SANE (top) and MAD (bottom) models. The dashed white line corresponds to the Bernoulli parameter, −hu_{t} = 1.02, and the black lines show log_{10}(σ) = −1.0, 0.5, and 0 contours. 

In the text 
Fig. A.2. Azimuthal and timeaveraged distribution of the κ value for different spins (left to right) and for the SANE (top) and MAD (bottom) models. The dashed white line corresponds to the Bernoulli parameter, −hu_{t} = 1.02, and the black lines show log_{10}(σ) = −1.0, 0.5, and 0 contours. 

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.