© 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 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 quasi-simultaneous multi-wavelength campaign of the M 87 nuclei (Algaba & Anczarski 2021). These efforts can reveal the underlying jet-launching and particle acceleration physics. Moreover, multi-wavelength 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; Lynden-Bell 2006).

More specifically, the jet collimation profile can be revealed from multi-frequency (10–230 GHz) VLBI observations. These observations have shown that up to hundreds of Schwarzschild radii¹, 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 ring-like 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 jet-launching. A theoretical model needs to capture the compact ring-like 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 multi-parametric problem, works have developed semi-analytical techniques to fit the spectrum and model the basic jet characteristics (Di Matteo et al. 2003; Broderick & Loeb 2009). To draw a self-consistent 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 two-temperature accretion flow and a hot single-temperature jet (Mościbrodzka et al. 2016). The addition of non-thermal electrons is expected to better reproduce the near-IR (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 non-thermal emission model on the structure and morphology of the M87 jet and its spectrum. We extended the study presented in Cruz-Osorio 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 broad-band spectrum and image structures. To this end, we perform long-term, high-resolution 3D GRMHD simulations. The two important parameters of these simulations are the black hole spin, a_⋆ = cJ/GM² 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) (, where Ṁ is the mass accretion rate) and magnetically arrested disks (MAD) (Ψ≥15). Our definition of the MAD parameter Ψ differs by a factor 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 non-thermal 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 post-process calculations. In Sect. 4 we analyse the impact of the free parameters of our model on the broad-band 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⁹ 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ⁱ) on a generic space-time four-metric (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 is absorbed in the magnetic field) are the conservation of mass, the local conservation of energy-momentum, and covariant Maxwell equations,

(1)

where ρ is the rest-mass density, and u^μ is the four-velocity. The stress-energy tensor for ideal MHD is defined as

(2)

where h_tot = h + b²/ρ is the total specific enthalpy, and p is the fluid pressure. The magnetic field four-vector b^μ is given by

(3)

where Bⁱ 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² and magnetisation σ = b²/ρ, are computed using

(4)

where B² = BⁱB_i, where Bⁱ is the magnetic field three-vector 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, ρ,

(5)

where is the adiabatic index (see Antón et al. 2006; Lora-Clavijo et al. 2015; Porth et al. 2017; Cruz-Osorio 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.

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 finite-volume and high-resolution shock-capturing methods, the Lax-van Leer Friedrichs (LLF) flux formula in combination with piecewise parabolic method (PPM) reconstruction, and a two-step predictor-corrector 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

(6)

where

(7)

(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 , plasma beta, β = 2p/b², 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 magneto-rotational 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.

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 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 quasi-stable 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 low-magnetisation 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.

In Fig. 2 we show the time² and azimuthally averaged large-scale morphology of the plasma magnetisation, σ = b²/ρ, 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³ and low plasma beta (β ≪ 1). In Fig. 2 the solid black lines indicate log₁₀σ = −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 out-flowing 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, Tcode, 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 log10σ = −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 single-dish 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

(9)

where ℐ is the Lorentz-invariant specific intensity and is connected to the specific intensity, I_ν, by

(10)

The Lorentz-invariant emissivity, η, and absorptivity, χ, are related to the emission, j_ν and absorption coefficients, α_ν evaluated at frequency ν via

(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,

(12)

with affine parameter λ and energy shift between the oberver’s and the co-moving frame γ⁻¹ = ν₀/ν = −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 fast-light approximation, assuming that the dynamical time of the GRMHD simulations is greater than the light-crossing 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 Runge-Kutta-Fehlberg (RKF45) integrator with fourth-order adaptive step sizing and fifth-order 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³ 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⁻¹⁴, 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 so-called R-β model, which assumes a plasma-β depending temperature ratio between protons and electrons,

(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⁻³ − 10³, 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 electron-to-proton temperature ratio, that is, (T_p/T_e)⁻¹. 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 electron-to-proton temperature ratio with plasma-β resembles the shape of a sigmoid-curve, 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. Electron-to-proton temperature ratio Te/Tp as a function of plasma β varying the Rlow and Rhigh (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 Rlow = 1, and dotted lines show models with Rlow = 0.1. The different colours indicate different Rhigh 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 high-energy 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⁹ Hz ≤ 10¹⁶ Hz), these processes can be neglected.

The two eDFs we explore in this work are the Maxwell-Jüttner distribution function for thermal electrons and the kappa distribution function for non-thermal 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 power-law electron distribution with 2 ≤ p ≤ 4. In addition, Davelaar et al. (2019) showed that the M 87 spectrum between (10⁹ Hz ≤ 10¹⁶ 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 Maxwell-Jüttner distribution function is given by

(14)

where n_e is the electron number density, γ_e is the electron Lorentz factor, and K₂ 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 power-law tail for larger γ_e. The kappa distribution can be written as

(15)

where N is a normalisation factor (see Pandya et al. 2016, for details). The power-law exponent of a non-thermal particle distribution 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,

(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 particle-in-cell simulation-based subgrid model from Ball et al. (2018). More precisely, we used the magnetisation and plasma-β dependent power-law slope of the electron distribution and connected this to the κ value,

(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 Maxwell-Jüttner and the kappa electron distribution functions. The solid black lines corresponds to a Maxwell-Jüttner eDF with an electron temperature of Θ_e = 10. As expected from Eq. (14), the Maxwell-Jü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 power-law tail of the kappa eDF is clearly visible for γ_e >  10³ and is well matched by a power-law 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 Maxwell-Jüttner eDF. This behaviour is clearly illustrated by the dashed orange curve, corresponding to a kappa eDF with w = 10 and κ = 10⁶. 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: Maxwell-Jüttner electron distribution (black), power-law electron distribution (blue) with s = 2.5, κ electron distribution for κ = 3.5 and w = 10 (red), and for κ = 106 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β²/(1 + β²) 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 electron-to-proton 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 electron-to-proton 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 Rhigh values while keeping Rlow = 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 high-energy 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 Maxwell-Jü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 Rhigh values while keeping Rlow = 1, rinj = 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 micro-physics and hyper-parameters, 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 hyper-surfaces 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.

During the GRRT, we computed both the broad-band spectrum from 10⁹ Hz to 10¹⁶ Hz and ray-traced 86 GHz images. This allowed us to compare the spectral properties of our models with multi-frequency 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⁹ 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 quasi-stable 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 Maxwell-Jüttner (thermal) and kappa distribution (non-thermal) (iii) the jet spine via the magnetisation cut-off σ_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 super-resolved 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⁴ 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 x-axis, and sliced the jet perpendicular to the x-axis. During the slicing, we limited the dynamical range to 10⁴ 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 kappa-eDF 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 Rhigh 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: broad-band spectrum averaged for 2000 M for Rhigh = 10 (dotted), Rhigh = 80 (dashed), and Rhigh = 160 (solid) while keeping Rlow = 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 ray-traced 86 GHz images for three different eDFs (columns) and different Rhigh (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 broad-band spectrum can be best seen for the thermal models (Maxwell-Jüttner eDF) in the MAD case (top left panel, blue curves in Fig. 7). The high-frequency part of the spectrum (ν >  10¹² Hz) steepens significantly with increasing R_high, while the low-frequency part (ν ≤ 10¹² Hz) slightly increases in flux without a change in the spectral slope (see inset in top left panel). Because the high-energy emission is mainly emitted from the disk region predominated by high plasma-β values, the electron temperature can be written as Θ_e ∝ 1/R_high⁵ (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 high-frequency regime is given by

(18)

which finally explains the exponential decrease with frequency and the variation with R_high for the high-frequency part of the spectrum.

In a similar manner, we can explain the low-frequency behaviour of the spectrum. Assuming that the low-frequency 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 high-energy 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

(19)

Inserting the definition of the width of the kappa eDF (Eq. (16)), the kappa emissivity can be approximated by⁶

(20)

Together with the thermal emissivitiy, the total emissivity has the simplified form

(21)

If we consider that the high-energy thermal emission originates mainly from the disk (high plasma-β region) and that the jet sheath⁷ has intermediate plasma-β values, the equation above can be further simplified and the dependence on the R_high can be included,

(22)

With this simplification of the emissivity at hand, we can understand the spectral behaviour in the high-frequency regime for the thermal–non-thermal eDFs. First we concentrate on the ε = 0 models (green curves in the top panels of Fig. 7). Including the kappa eDF leads to a power-law tail in the high-frequency 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 non-thermal 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 high-frequency 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 non-thermal 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 high-frequency 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 (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

(23)

Here, 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¹¹ Hz, compared to ∼10¹² 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 counter-jet 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 broad-band spectrum with increased R_high values can be explained.

Including non-thermal 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 non-thermal eDF, we can explain this behaviour by Eq. (22). The addition of non-thermal 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 non-thermal 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 non-thermal 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 broad-band 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¹⁴ 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¹⁴ 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,

(24)

Fig. 8. Influence of the Rlow 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: broad-band spectrum averaged for 2000 M for Rhigh = 10 (dotted), Rhigh = 80 (dashed), and Rhigh = 160 (solid) while keeping Rlow = 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 ray-traced 86 GHz images for three different eDFs (columns) and different Rhigh (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 broad-band 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 broad-band 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 non-thermal 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 broad-band 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 broad-band spectrum the increase in the σ_cut value leads to higher flux density and flatter spectra in the high-frequency regime (ν >  10¹² 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 power-law slope of the non-thermal 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 kappa-eDF (see Eq. (16)), that is, to more energetic non-thermal 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 broad-band spectrum averaged for 2000 M for six different values of σcut (colour-coded) while keeping Rlow = 1, Rhigh = 160, and ε = 1.0 fixed. The middle image shows the ray-traced 86 GHz images for different σcut values. The bottom panels show the jet diameter and the opening angles profiles using the same colour-coding 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 power-law slope of the non-thermal 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 broad-band 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 non-thermal 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 broad-band spectrum (top panels) shows that increasing the ε value increases the flux density at high frequencies (ν >  10¹² Hz) and flattens the spectrum. By increasing the ε value, more energy from the magnetic field is included in the non-thermal particles (see Eq. (16)) and the non-thermal 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¹² 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 non-thermal 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: broad-band spectrum averaged for 2000 M for thermal and five different values of ε (colour-coded) while keeping Rlow = 1, Rhigh = 160, and σcut = 1.0 fixed. The middle image shows the ray-traced 86 GHz images for different ε values. Bottom panels: jet diameter and opening angles profiles using the same colour-coding 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 high-energy 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 broad-band 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 co-rotating and counter-rotating 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 broad-band 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 broad-band spectrum (see the top panels in Fig. 11). Interestingly, changes in the black hole spin for this set of best-fit parameter introduce only a minor scatter in the broad-band spectrum, especially in the MAD case. The low-frequency part of the spectrum (ν <  10¹¹ 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: broad-band spectrum averaged for 2000 M for thermal and five different spin values astar = −0.9375, −0.5, 0, 0.5, and 0.9375 (colour-coded) while keeping Rlow = 1, Rhigh = 160 σcut = 3.0, and ε = 0.5 fixed. The middle image shows the ray-traced 86 GHz images for different spin values. Bottom panels: jet diameter and opening angles profiles using the same colour-coding 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 non-thermal particles via the kappa-eDF. Independent of the accretion model, Schwarzschild black holes (a_⋆ = 0) do not launch a well-established 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 . 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 broad-band spectrum and image structure of M 87. The workflow for this analysis can be divided into three main branches: (i) 3D GRMHD simulations, (ii) multi-frequency 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 a-priori 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 disk-jet system in M 87. The simulations were performed in modified Kerr-Schild 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 steady-state 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 high-resolution 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 broad-band 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 non-thermal 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 broad-band 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 broad-band spectrum (best shown in the slope at high frequencies (ν >  10¹² Hz) in panel c in Figs. 12 and 7), in agreement with the findings of Davelaar et al. (2019) for a fast-spinning 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 broad-band spectra for different spins (see panel f in Figs. 12 and 11). Thus, we conclude that fitting the broad-band spectrum alone will not provide constraints on the black hole spin.

These variations in the broad-band 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, fast-rotating 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 broad-band spectrum and the image structure, we can finally compare our models to the observations of M 87. The broad-band spectrum for M 87 is remarkably well fitted for the MAD models using non-thermal 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 high-spin prograde MAD models produce collimation profiles that agree with the currently highest resolution VLBI images⁸ of the jet in M 87 (see the bottom panels in Fig. 11). In Fig. 13 we present our best-fit 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 edge-brightening, 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 super-resolved 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 best-fit 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⁻⁴ M_Å yr⁻¹) (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 jet-launching 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.

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 edge-brightening 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 post-processing 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², (see Table 2 for details). However, recent works of non-ideal GRMHD including a mean field dynamo showed that poloidal magnetic fields can be self-similarly generated from an almost zero initial magnetic field (plasma beta 10⁶ − 10⁹) (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 jet-launching process while reducing the number of free parameters in the modelling.

6. Conclusion and outlook

We presented long-term and high-resolution 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 broad-band spectrum (10⁹ ≤ ν ≤ 10¹⁶ 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 non-thermal particles), and the jet sheath on the radiative signatures. From our detailed parameter survey, we found that the best fit to the broad-band spectrum (see Table B.1 for the used observational data) is obtained if non-thermal particles energetically supported by the magnetic field (ε = 0.5) are included in a magnetised jet sheath (σ_cut = 3). Interestingly, the broad-band 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 broad-band spectrum and collimation profiles favoured the MAD accretion model together with fast-spinning 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 best-fit models, one of the main results of this work is that the combined modelling of broad-band spectrum and additional constraints such as the image structure can help to further constrain possible models for M 87.

Future work of modelling the jet-launching and propagation in M 87 should include more advanced, non-standard 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 high-energy radiation (X-rays) should be computed during the radiative transfer. Finally, to reduce the influence of the hyper-parameter 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).

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/2-1), 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 Post-Doctoral 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 GOETHE-HLR LOEWE at the CSC-Frankfurt, 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/eht-imaging/) (Chael et al. 2018)

References

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₁₀(σ) = −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 time-averaged 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, −hut = 1.02, and the black lines show log10(σ) = −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₁₀(σ) = −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 time-averaged 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, −hut = 1.02, and the black lines show log10(σ) = −1.0, 0.5, and 0 contours.

Appendix B: Multi-frequency observations of M 87

In Table B.1 we provide the observational data with the references we used for the broad-band spectrum in Figs. 7 -11.

All Tables

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, Tcode, 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 log10σ = −1.0,  0,  and 0.5 correspond to jet sheath and jet spine boundaries.

In the text

Fig. 3. Electron-to-proton temperature ratio Te/Tp as a function of plasma β varying the Rlow and Rhigh (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 Rlow = 1, and dotted lines show models with Rlow = 0.1. The different colours indicate different Rhigh values.

In the text

	Fig. 4. Comparison of eDFs: Maxwell-Jüttner electron distribution (black), power-law electron distribution (blue) with s = 2.5, κ electron distribution for κ = 3.5 and w = 10 (red), and for κ = 106 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 Rhigh values while keeping Rlow = 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 Rhigh values while keeping Rlow = 1, rinj = 10 M, and ε = 0.5 fixed.
In the text

Fig. 7. Influence of the Rhigh 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: broad-band spectrum averaged for 2000 M for Rhigh = 10 (dotted), Rhigh = 80 (dashed), and Rhigh = 160 (solid) while keeping Rlow = 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 ray-traced 86 GHz images for three different eDFs (columns) and different Rhigh (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 Rlow 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: broad-band spectrum averaged for 2000 M for Rhigh = 10 (dotted), Rhigh = 80 (dashed), and Rhigh = 160 (solid) while keeping Rlow = 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 ray-traced 86 GHz images for three different eDFs (columns) and different Rhigh (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 broad-band spectrum averaged for 2000 M for six different values of σcut (colour-coded) while keeping Rlow = 1, Rhigh = 160, and ε = 1.0 fixed. The middle image shows the ray-traced 86 GHz images for different σcut values. The bottom panels show the jet diameter and the opening angles profiles using the same colour-coding as for the spectrum.

In the text

Fig. 10. Influence of the ε parameter (fraction of magnetic energy used for the acceleration of non-thermal 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: broad-band spectrum averaged for 2000 M for thermal and five different values of ε (colour-coded) while keeping Rlow = 1, Rhigh = 160, and σcut = 1.0 fixed. The middle image shows the ray-traced 86 GHz images for different ε values. Bottom panels: jet diameter and opening angles profiles using the same colour-coding 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: broad-band spectrum averaged for 2000 M for thermal and five different spin values astar = −0.9375, −0.5, 0, 0.5, and 0.9375 (colour-coded) while keeping Rlow = 1, Rhigh = 160 σcut = 3.0, and ε = 0.5 fixed. The middle image shows the ray-traced 86 GHz images for different spin values. Bottom panels: jet diameter and opening angles profiles using the same colour-coding as for the spectrum.

In the text

	Fig. 12. Schematic variation of the broad-band 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 super-resolved 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 time-averaged 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, −hut = 1.02, and the black lines show log10(σ) = −1.0, 0.5, and 0 contours.
In the text

	Fig. A.2. Azimuthal and time-averaged 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, −hut = 1.02, and the black lines show log10(σ) = −1.0, 0.5, and 0 contours.
In the text