© The Authors 2026

1. Introduction

Messier 87 (M87), a giant elliptical galaxy in the Virgo cluster, has provided an important laboratory in the field of extragalactic jet astrophysics for many decades. The jet, originating from a compact radio core in the galaxy center, where the active galactic nucleus (AGN) is located, has been studied extensively across the electromagnetic spectrum, ranging from the radio regime up to γ-rays (e.g. Schmitt & Reid 1985; Hada et al. 2013; Kim et al. 2018; Snios et al. 2019; Acciari et al. 2020; Lu et al. 2023). The Event Horizon Telescope Collaboration (EHTC; EHTC 2019c) resolved the radio core of M87* at 1.3 mm employing the very long baseline interferometry (VLBI) technique and interpreted the ring-like structure as the gravitationally lensed synchrotron emission, originating from the accretion disk of a supermassive black hole (SMBH). In more recent observations at 86 GHz, Lu et al. (2023) revealed that the jet and the radio core are indeed linked, which was interpreted as the jet base connecting to the accretion flow of relativistic plasma. Super-resolution images of mentioned observations with the Global Millimetre VLBI Array (GMVA) and Very Long Baseline Array (VLBA) at 43 GHz were obtained, employing the resolve algorithm (see Kim et al. 2025, 2024). The polarimetric EHT observations obtained in 2017 provide viable information on the dynamically important magnetic-field structure. From these observations, a linear polarization fraction of ∼15% (EHTC 2021a) and an upper bound on the circularly polarized emission fraction of ⟨|ν|⟩ < 3.7% (EHTC 2023) was found. These results constrain theoretical models such as general-relativistic magnetohydrodynamic (GRMHD) simulations of the plasma flow in the vicinity of the black hole (BH). These simulations can be compared to observations by means of general-relativistic radiative-transfer (GRRT) post-processing calculations on the GRMHD data. Yang et al. (2024) found the 86 GHz and 43 GHz jet morphology of similar GRMHD simulations to fit the observations very well. Cruz-Osorio et al. (2022) and Fromm et al. (2022) found high spin parameters of a ∼ 0.5 − 0.94 to be favorable to match the broad-band spectrum and jet morphology of M87* in wide parameter-space studies of GRMHD models. In this work, we focused on the comparison of five 3D GRMHD models and their impact on the electromagnetic footprint of the M87 jet. All share similar initial conditions, but differ with respect to their initial magnetic-field configuration. We used the Fishbone-Moncrief (FM) torus (Fishbone & Moncrief 1976) and initialized it with a weak poloidal magnetic field. In most studies, the poloidal extent of the initial magnetic field is set with the offset of the initial vector potential from a zero cutoff to 0.80. We compared this model with four simulations that fill the initial torus with varying spatial extents of a weak poloidal magnetic field. These models exhibit drastically altered energetic-flux, angular-momentum transport properties and jet-generation efficiencies. Chatterjee & Narayan (2022) studied the behavior and drivers of angular-momentum transport in magnetically arrested disk (MAD) simulations compared to standard- and normal-evolution (SANE) simulations in non-spinning BH systems to isolate the side effects of frame dragging. We complemented this study with the case of BHs with high spin and compared the impact of the magnetic filling of the torus. Recent studies by Jacquemin-Ide et al. (2025) complement explanations for the underlying mechanism of MAD states by studying the inward advective and outward diffusive cycles these simulations exhibit, and, based on that, they derived a first analytical expression for the recurrence timescale of MAD states. Further parametrical investigation of the MAD setup by Cho & Narayan (2025) found that the overall variability in the simulations is sensitive to the ratio of rotational and magnetic energy in the initial setup, which connects well to our study, where we tuned the initial magnetic-energy content.

This paper contains two major parts. In Sect. 2, we introduce the underlying GRMHD equations and simulation setups. We describe the tests we conducted on the resolution of the magneto-rotational instability (MRI). Thereafter, we explain our study of the impact of the initial magnetic energy in the torus on the jet energetics and on the efficiency of energy extraction from the BH and disk to the energy flux in the jet. We connect this analysis further with the investigation of the angular momentum transport of the five simulations. In connecting it to the morphological study of the radiative propertities of a subset of three simulations in Sect. 3, we considered the impact of magnetic-flux eruptions from the BH on the hydrodynamical jet structure of these simulations. In Sect. 3, we describe how we ray-traced the polarized synchrotron emission from a population of thermal and nonthermal electrons asserted by the GRMHD simulations. After the ray-tracing setup is introduced, we compare the jet morphology at 86 GHz and present spatially resolved, as well as time-dependent, linear polarization fractions. Finally, we compare the jet edge and width profiles of the three different models with observations from Lu et al. (2023) and discuss the impact of flux eruptions on the synchrotron emission from the jet and disk.

2. General-relativistic magnetohydrodynamics

We employed the 3D GRMHD code Kokkos-based High-Accuracy Relativistic Magnetohydrodynamics with AMR (KHARMA; Prather et al. 2021; Prather 2024), which simulates accretion onto and jet launching from BHs. KHARMA solves the ideal magneto-hydrodynamics (MHD) equations on curved space times with a four-metric g_μν and a metric-determining g in geometric units (G = c = 1, and a factor of $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{4\pi }$}\right.$ is absorbed in the magnetic field). Throughout this paper, we use code units for the time and length scales; i.e., the light-crossing time is redefined as t_g = 1 G_NM_BH/(c³)≡1 M, and the gravitational radius is r_g = 1 G_NM_BH/(c²)≡1 M, where G_N is Newton’s constant, c is the vacuum speed of light, M_BH is the BH mass, and M is the mass unit that scales the units in the simulations. Greek indices run through (0, 1, 2, 3) and Roman indices through (1, 2, 3), and we made use of the Einstein summation convention. We used a Kerr black-hole solution on a modified spherical coordinate system including funky-modified–Kerr–Schild coordinates (FMKS). These coordinates exhibit a decreasing grid cell size in the θ-direction toward the equator ($\theta =\raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$2$}\right.$) that becomes increasingly cylindrical along the jet-axis (θ = 0 and θ = π). The radial extent of the cells increases exponentially; therefore, high resolution is given at small radii, where finely resolving the accretion process is vital. (r, θ, ϕ) refer to the ordinary spherical coordinates. The simulation grids are equal in all simulations, in which we employed N_r = 256, N_ϕ = 128, and N_θ = 128 cells in the respective directions. The radial boundaries range from r = (1.1659 − 1000) M and otherwise use the complete map in ϕ = (0 − 2π) rad and θ = (0 − π) rad. The underlying equations solved by KHARMA are the covariant conservation of mass, local conservation of energy-momentum, and the covariant Maxwell equations (see Mizuno & Rezzolla 2024 for more details on the derivation):

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mu }(\rho u^{\mu })=0;{\mathrm{\nabla }}_{\mu }{T}^{\mu \nu }=0;{\mathrm{\nabla }}_{\mu }^{\ast }{F}^{\mu \nu }.\end{array}$$(1)

ρ denotes the rest-mass density and u^μ the four-velocity of the fluid. ^*F^μν is the dual Faraday tensor, and the stress-energy tensor is defined in the case of ideal MHD as follows:

$$\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)

with a total specific enthalpy of $h_{\mathrm{tot}}=h+\raisebox{1ex}{$b^2$}\!\left/ \!\raisebox{-1ex}{$\rho $}\right.$ and a fluid pressure of p. The magnetic-field four-vector b^μ in its temporal component and spatial components reads

$$\begin{array}{c}\hfill b^t=\frac{\mathrm{\Gamma }(B^i{u}^{\mu })}{\alpha };b^i=\frac{(B^i+\alpha b^t{u}^i)}{\mathrm{\Gamma }},\end{array}$$(3)

where Bⁱ are the spatial components of the magnetic-field three-vector measured by an Eulerian observer, α is the lapse function, and Γ is the bulk Lorentz factor. Note that with the definition of b^μ in Eq. (3), we made use of the notation b² = b^μb_μ in Eq. (2). We used the plasma beta, $\beta =\raisebox{1ex}{$2p$}\!\left/ \!\raisebox{-1ex}{$b^2$}\right.$, and magnetization, $\sigma =\raisebox{1ex}{$b^2$}\!\left/ \!\raisebox{-1ex}{$\rho $}\right.$, parameters for diagnostics with

$$\begin{array}{c}\hfill b^2=\frac{B^2+{\alpha }^2{(b^t)}^2}{{\mathrm{\Gamma }}^2}=\frac{B^2}{{\mathrm{\Gamma }}^2}+{(B^i{v}_i)}^2.\end{array}$$(4)

The system of equations is closed by an equation of state. We assumed an ideal gas that relates the 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 the adiabatic index is $\widehat{\gamma }$ (Antón et al. 2006).

We employed inflow and outflow boundary conditions in a radial coordinate. Along the polar boundaries, we used a solid reflective wall (Olivares et al. 2019). In the azimuthal direction, we employed periodic boundary conditions to all physical quantities across the cells at ϕ = 0. The GRMHD equations were solved with finite-volume and high-resolution shock-capturing methods including the Lax-van Leer Friedrichs (LLF) flux formula in combination with the WENO5 reconstruction scheme and in synchronization with the IMEX scheme. The density floor was set to ρ_min = 10⁻⁶ and the internal energy floor to u_min, geo = 10⁻⁸; the ceiling of the magnetization was σ_max = 100, and the ceiling of the internal energy was ${\left(\raisebox{1ex}{$u$}\!\left/ \!\raisebox{-1ex}{$\rho $}\right.\right)}_{\max }=2$. The initial pressure maximum in the torus was set to ρ_max = 1. The 3D GRMHD simulations were initialized with a magnetized torus in hydrodynamic equilibrium with a rotating Kerr BH in the center of the simulation domain. Here, we used a setup that allows the disk to become partially magnetically arrested. In the MAD setup, we added a single poloidal loop in the vector potential on top of the FM torus (Fishbone & Moncrief 1976) in the following way:

$$\begin{array}{c}\hfill A_{\varphi }\propto max(q-(1-F),0),\end{array}$$(6)

with

$$\begin{array}{c}\hfill q=\frac{\rho }{{\rho }_{\max }}{\left(\frac{r}{r_{\mathrm{in}}}\right)}^3{sin}^3\theta exp\left(\frac{-r}{400}\right)·\end{array}$$(7)

We solely varied the filling factor, F, throughout the simulations from 0.6 − 0.99, keeping other initial conditions of the GRMHD simulations constant. Those are a constant adiabatic index, $\widehat{\gamma }=\frac{5}{3}$; inner radius, r_in = 20; and radius of the density maximum, r_c = 41, of the torus as well as the dimensionless spin parameter, a = 15/16, of the BH. The plasma beta was initialized such that its minimal value was above ${\beta }_{\min }=2\raisebox{1ex}{$p$}\!\left/ \!\raisebox{-1ex}{$b^2$}\right.=100$. We refer to the five models as K.60, K.80, K.90, K.95, and K.99; these correspond to the simulations in that the filling factor is set to F = 0.60, F = 0.80, F = 0.90, F = 0.95, and F = 0.99, respectively. In Appendix A, we present more details on the correspondance of the filling factor to the volumes initialized with the magnetic field relative to the torus volume.

In Eq. (6), the effect of F is to proportionally scale the offset from the zero cutoff of A_ϕ, which is visible in the bottom row of Fig. 1. With increasing F, the offset of A_ϕ increases, which in turn increases the poloidal extent or, to be precise, the total volume containing cells that have nonzero magnetic-field values. For further details, we refer the reader to Appendix A, where we also address that the filling factor does not scale the poloidal-field extent linearly to the extent of the torus; it is rather close to the latter, and for F ≳ 0.80 it is above it (see Table A.1).

Fig. 1. Invariant magnetic-field strength, b; the spatial magnetic-field components, Br and Bθ; and four-vector-potential component, Aϕ, over the polar angle at the initialization step (t = 0 M). Colors correspond to simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta), and K.99 (red). In the middle row, the solid lines show Br and the dashed lines Bθ. We chose polar slices at radii close to the inner edge of the torus at rin = 20 M and at the pressure maximum in the tori at r = 50 M. Note that in columns going from left to right, i.e., radially outwards, the simulations with higher filling factors show increased magnetic-field strengths and extended magnetic field, both radially and in the polar-angle, especially at the inner edge of the torus.

All simulations were calculated up to a simulation time of 30 000 M. We show the mass accretion rate, Ṁ; the magnetic flux through the event horizon, ϕ_BH; and the MAD parameter, $\mathrm{\Psi }=\raisebox{1ex}{${\varphi }_{\mathrm{BH}}$}\!\left/ \!\raisebox{-1ex}{$\sqrt{\dot{M}}$}\right.$, over the full temporal range of the simulations in Fig. 2. In Appendix B, Fig. B.1, we present azimuthal and time-averaged cross-sections of the meridional plane of ρ, σ, magnetic field magnitude $b=\sqrt{|b_{\mu }{b}^{\mu }|}$, plasma-β, and the ion temperature, Θ_i.

Fig. 2. Mass-accretion rates, Ṁ, and magnetic flux through the event horizon of the BH, ϕBH, and the MAD-parameter, $\mathrm{\Psi }={\varphi }_{\mathrm{BH}}/\sqrt{\dot{M}}$, for all simulations over time.

2.1. Resolution and MRI suppression tests

In this section, we describe our investigation of the MRI resolution, which is thought to be the main driver of angular momentum and energy transport of plasma accretion in SMBH systems. Due to the limitations of finite resolution, the spectrum of MRI modes is truncated in simulations; hence, it is able to produce various numerical artifacts. Truncation of the MRI modes affects the wave-number range of turbulence modes caused by the MRI and consequently alters the rate at which nonlinear mode–mode couplings transfer energy from large- to small-scale dynamics. Insufficient resolution of the fastest growing MRI mode brings the turbulent field amplification to a halt, and numerical dissipation takes over; hence, the field decays. We tested wether underlying MRI modes were well resolved by the physical criterion and that the number of grid cells that are able to capture the MRI wavelength λ_MRI was sufficiently large. The so-called quality or Q factors (Sano et al. 2004; Noble et al. 2010; Porth et al. 2019) were defined to reflect the resolution of the underlying MRI mode with the grid spacing in a given spatial direction, Δx⁽ⁱ⁾:

$$\begin{array}{c}\hfill Q^{(i)}\equiv \frac{{\lambda }_{\mathrm{MRI}}^{(i)}}{\mathrm{\Delta }{x}^{(i)}}·\end{array}$$(8)

The wavelength of the fastest growing MRI mode and the cell extent in a given direction can be defined as

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{MRI}}^{(i)}\equiv \frac{2\pi }{\sqrt{(\rho h+b^2)\mathrm{\Omega }}}{b}^{\mu }{e}_{\mu }^{(i)},\mathrm{\Delta }{x}^{(i)}\equiv {(\mathrm{\Delta }{x}^{\mu })}^{(i)}{e}_{\mu }^{(i)}.\end{array}$$(9)

Ω = u^ϕ/u^t is the azimuthal angular velocity of the fluid measured in the locally non-rotating reference frame (LNRF) in Kerr space time, and e_μ⁽ⁱ⁾ is the i-th one-form of the co-vielbein, determining the tetrad basis, where the LNRF is boosted into the comoving frame of the fluid (see Takahashi 2008 for the basis transformations). Δx⁽ⁱ⁾ is the grid spacing in a fluid frame tetrad basis oriented along the LNRF. (Δx^μ)⁽ⁱ⁾ is a vector containing, for each i, one of the grid spacings in the Kerr–Schild coordinates. For the θ coordinate and grid spacing, Δθ, in Kerr–Schild this would look as follows: (Δx^μ)^(θ) = (0, 0, Δθ, 0)^T. This is analogously calculated for the other coordinates (r, ϕ). We calculated the quality factors Q^θ and Q^ϕ, and we averaged them over a region that includes the wind and disk regions, which is defined by the bounds π/3 ≥ θ_disk ≥ 2π/3 and the radial upper bound r_disk ≤ 50 M, where we always took the full azimuthal map of 0 ≥ ϕ > 2π. The spatial average is calculated as the expectation value over the covariant volume as follows:

$$\begin{array}{c}\hfill {⟨Q^{(i)}⟩}_w=\frac{\int \int \int \sqrt{-g}\mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi Q^{(i)}(r,\theta ,\varphi )w(r,\theta ,\varphi )}{\int \int \int \sqrt{-g}\mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi w(r,\theta ,\varphi )},\end{array}$$(10)

where w(r, θ, ϕ) is a weighting function that acts as a masking function. In the case of the disk averages, this function is 1 inside the poloidal bounds specified above and 0 outside of them.

We show the temporal evolution of the disk’s averaged quality-factors and their median values in Appendix D, Fig. D.1. Sano et al. (2004) suggested Q factors of Q > 6 to correspond to a sufficient resolution of the MRI. We matched these requirements in Q^(θ) and Q^(ϕ) with consistent values in all simulations with median values in time of $Q_{\mathrm{med}}^{(\theta )}\gtrsim 35$, $Q_{\mathrm{med}}^{(\varphi )}\gtrsim 12$ and time averages of ⟨Q^(θ)⟩_t ≳ 60, ⟨Q^(ϕ)⟩_t ≳ 16, which are valid for all simulations over the time span [10 000, 30 000] M (exact values are shown in Fig. D.1, and meridional snapshots are given in Fig. D.2).

2.2. Electromagnetic and kinetic energy flux in the jet

For the time-series analysis presented here, we estimated the error of all time-averaged quantities by means of the correlation error, σ_corr (see Appendix E for the formula and more details). This takes into account that the number of data points that are independent of each other is reduced if the data are correlated.

The total energy flux is determined via the energy momentum tensor from Eq. (2) by contraction with the metric tensor and reformulating the enthalpy in terms of the internal energy density $u=\raisebox{1ex}{$p$}\!\left/ \!\raisebox{-1ex}{$(\widehat{\gamma }-1)$}\right.$ to be

$$\begin{array}{cc}\hfill T_{\mu }^{\nu }& =(\rho +\widehat{\gamma }u+b^2)u^{\nu }{u}_{\mu }-b^{\nu }{b}_{\mu }+(p+\raisebox{1ex}{$b^2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.){\delta }_{\mu }^{\nu },\hfill \end{array}$$(11)

$$\begin{array}{cc}& =(\rho +p+u+b^2)u^{\nu }{u}_{\mu }-b^{\nu }{b}_{\mu }+(p+\raisebox{1ex}{$b^2$}\!\left/ \!\raisebox{-1ex}{$2$}\right.){\delta }_{\mu }^{\nu },\hfill \end{array}$$(12)

where δ_μ^ν is the Kronecker delta. The sub-tensor responsible for the EM contribution to the r − t component is

$$\begin{array}{c}\hfill T_{\mathrm{EM}t}^r=b^2{u}^r{u}_t-b^r{b}_t.\end{array}$$(13)

Analogously, the sub-tensor for the kinetic contributions is

$$\begin{array}{c}\hfill T_{\mathrm{Kin}t}^r=(\rho +p+u)u^r{u}_t.\end{array}$$(14)

The total, electromagnetic (EM), and kinetic energy fluxes were then determined by

$$\begin{array}{cc}\hfill \dot{E}& =-{\displaystyle \int \int }\sqrt{-g}\mathrm{d}\theta \mathrm{d}\varphi T_t^r,\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill {\dot{E}}_{\mathrm{EM}}& =-{\displaystyle \int \int \sqrt{-g}}\mathrm{d}\theta \mathrm{d}\varphi T_{\mathrm{EM}t}^r,\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill {\dot{E}}_{\mathrm{Kin}}& =-{\displaystyle \int \int \sqrt{-g}}\mathrm{d}\theta \mathrm{d}\varphi T_{\mathrm{Kin}t.}^r\hfill \end{array}$$(17)

Before executing the integral, we imposed the condition that the Bernoulli parameter had to be > 1.02, corresponding to the outflowing part of the plasma shown by the area inside the contours in the bottom row of Fig. 10. We further calculate the energy fluxes at varying radial coordinates and times.

In Fig. 3, we show that the time dependence of the total and EM energy fluxes calculated at a fixed radius of r = 10 M. K.99 on average exhibits a significantly higher energy flux of ∼50 [a.u.], exceeding the other simulations by ≳3σ_corr; this directly boosts the efficiency of the jet launching due to the comparable mass-accretion rates over all simulations. The energy flux at the jet base is highly dominated by the EM part of the energy flux (∼87%) for K.99 and causes the total energy flux excess. The dependence of the time-averaged total energy flux through the jet on the distance from the BH along the jet’s parallel axis can be taken from Fig. 4. The EM and kinetic contributions are in equipartition at different distances from the BH (vertical dashed lines) for each simulation. This is especially evident for simulation K.99, which is significantly more electromagnetically dominated both in extent along the jet and magnitude. We note that the energy fluxes increase with BH separation along the jet, since plasma from the wind region becomes gravitationally unbound; i.e., it reaches hu_t > 1.02.

Fig. 3. Temporal evolution of total energy flux (solid lines) and EM energy flux (dashed lines) measured at r = 10 M in the jet for simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta), and K.99 (red). The horizontal lines show the average of the total energy fluxes for each run over the shown time interval.

Fig. 4. Time-averaged total energy fluxes (dash-dotted lines), EM energy fluxes (solid lines), and kinetic energy fluxes (dashed lines) in the time interval (25 000 − 30 000) M for varying BH separations along the jet for simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta), and K.99 (red). The vertical dotted lines indicate the radii at which the kinetic and EM energy fluxes are in equipartition.

2.3. Outflow efficiencies

We define the efficiency of energy extraction from the BH spin and disk to outflowing energy at the poles as follows:

$$\begin{array}{c}\hfill \eta =\frac{{\dot{M}}_{5\mathrm{M}}-\dot{E}}{⟨{\dot{M}}_{5\mathrm{M}}⟩},\end{array}$$(18)

where we made use of Eq. (15) and calculated the mass-accretion rate at r = 5 M, Ṁ_5M and its averaged ⟨Ṁ_5M⟩ over the time interval t = [25 000, 30 000] M. This time interval is the most suitable one, as the accretion disk is well converged into a steady state. This can be seen by the mass-accretion rate decelerating its decline toward the end of the simulations (see Fig. 2). For comparison, we also show Ṁ_5M in the aforementioned interval in Fig. C.1. The nearly identical average mass-accretion rates also aid the comparison of the efficiencies (see Eq. (18)).

We explicitly did not use the mass-accretion rate at the event horizon, Ṁ (compare Fig. 2); we instead calculated the flux at 5 M since this excludes the injected matter to set the density floor at the inner simulation boundary. The radial rest-mass energy density flux is defined as follows:

$$\begin{array}{c}\hfill {\dot{M}}_{5\mathrm{M}}={\displaystyle \int {\int }_{r=5\mathrm{M}}}\sqrt{-g}\mathrm{d}\theta \mathrm{d}\varphi \rho u^r.\end{array}$$(19)

The resulting efficiencies of the energy-extraction process in all three models are shown in Fig. 5. On average, the K.99 simulation exhibits a vastly improved efficiency of ∼60% excess compared to the commonly used MAD setup, K.80. The averages and peak values of all simulations can be found in Table 1.

Fig. 5. Jet efficiencies for simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta), and K.99 (red), where the dashed horizontal lines show the time-averaged efficiencies and the correlation-corrected 2σcorr error band over the shown interval (listed in Table 1).

2.4. Angular-momentum transport

In this section, we investigate the angular-momentum transport by linking it to the jet efficiency and to the impact of the magnetic energy content in the torus on the dynamics from the initialization of the simulations. Following the approach and definitions of the angular-momentum fluxes of Chatterjee & Narayan (2022), we define the radial component of the time-averaged and shell-integrated angular-momentum flux vectors, ${\dot{J}}^i$ (i ∈ (r, θ)), accordingly. The total flux is defined as

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{tot}}^i={⟨T_{\varphi }^i⟩}_{t,\varphi },\end{array}$$(20)

and the advective flux reads

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{adv}}^i={\langle (\rho +u_{\mathrm{g}}+b^2/2)u^i\rangle }_{t,\varphi }{⟨u_{\varphi }⟩}_{t,\varphi }.\end{array}$$(21)

The stress-induced flux is given by

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{stress}}^i={\dot{J}}_{\mathrm{tot}}^i-{\dot{J}}_{\mathrm{adv}}^i,\end{array}$$(22)

which can be further subdivided into the flux induced by Maxwell-stresses,

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{stress},\mathrm{M}}^i={⟨b^2{u}^i{u}_{\varphi }/2-b^r{b}_{\varphi }⟩}_{t,\varphi },\end{array}$$(23)

and by Reynolds stresses, which are subdominant and therefore do not add anything to the discussions in this paper; hence, ${\dot{J}}_{\mathrm{stress}}^i\sim {\dot{J}}_{\mathrm{stress},\mathrm{M}}^i$.

We time-averaged the quantities again over the t = [25 000, 30 000] M interval in order to compare it with the jet energetics and efficiencies from Sects. 2.2 and 2.3.

In Fig. 6, we present the polar and radial dependence of time- and azimuthally averaged angular-momentum transport. In the left panel, we show the time- and azimuthally averaged total angular-momentum vector field for the first time for simulations in a MAD setup with a high prograde BH spin (cf. Chatterjee & Narayan 2022 Fig. 10, showing the case of a non-rotating BH in the MAD-setup). The meridional cross-section of the angular-momentum vector field exhibits quite distinct characteristics from the case of the MAD setup with a non-rotating BH, in that it exhibits a monopole field geometry in the jet and wind region up to a polar angle of ∼70° and transitioning into a quadrupole field geometry in the disk and torus. Through this magnetospheric configuration, the polar outflow from the jet and wind feeds back into the equatorial inflow in the disk, as opposed to the non-rotating MADs, where there is a transition from an equatorially parallel inflow inside the disk scale height (similar opening angle to the transition to the loops) to the parabolically outflowing angular-momentum flux in the jet. The major difference here is the dominance of the angular momentum transport driven by the rotating BH, where the dominantly radial outflow is redirected toward the equatorial plane into the disk, where the torus has to facilitate the radially outward angular-momentum transport in the wind region.

Fig. 6. Left panel: Time- and azimuthally averaged density distribution in the meridional plane and total angular-momentum flux streamlines (color-coded by a log-scale of the absolute value; increasing from black to white; simulation: K.80). Middle column: Polar slices of the radial angular-momentum flux at r = 2 M, indicated by the dashed gray lines in the left panel and right column. Shown are the total, advective and Maxwell-stress-induced angular-momentum fluxes from top to bottom. Right column: Shell-integrated radial component of the angular-momentum fluxes. Solid lines indicate positive values, i.e., outward fluxes and dashed lines are the absolute values of inward (negative) flux.

In the middle and right columns of Fig. 6, we display the impact of the magnetic filling of the MAD torus on the angular momentum transport. The middle column shows the polar dependence of the time- and azimuthally averaged radial component of the angular-momentum vectors at r = 2 M, where there is a maximum in the shell-integrated fluxes. The right panel shows the shell-integrated (see Sect. 2.2 for the definition of shell-integration) and time-averaged radial angular-momentum fluxes. The polar slices (indicated by dashed gray lines in the left panel and right column) of ${\dot{J}}_{\mathrm{tot}}^r$ and ${\dot{J}}_{\mathrm{stress},\mathrm{M}}^r$ show a clear model trend. Higher initial magnetic energy from the torus leads to an increase in outward angular-momentum flux, which is also sustained at larger radii of r ≲ 50 M (presented in Fig. F.1).

In the right panel at the top, one can see that the total shell-integrated fluxes are positive; i.e., radially outgoing flux in the jet and wind are dominant, whereas for the case of a non-rotating BH in Chatterjee & Narayan (2022), Fig. 14, the shell-integrated radial component of the total angular-momentum flux is constant and negative up to r ∼ 100 M. Secondly, overall the angular-momentum fluxes peak close to the event horizon of the BH, specifically inside the ergosphere, and they drop outside of this region.

We note that the model trend is not identical for the shell-integrated total and Maxwell-induced fluxes, since these are affected by the inward advective flux, which is more sensitive to turbulent realizations and does not appear to exhibit any model trend. The model trend is the clearest in the jet sheath as opposed to the equatorial plane in both total and Maxwell-induced angular-momentum fluxes. Looking at the middle panel, one can clearly see that the advective equatorial flux for K.80 is the highest; hence, it decreases the total shell-integrated flux the most for this model. However, we observe that simulation K.99 overall exhibits the highest angular-momentum flux both inward, driven by advective stresses, and especially outward, where the electromagnetic stresses are dominant. Furthermore, we observe that the model trend in the jet efficiency and the shell-integrated radial angular-momentum fluxes is consistent with one another, which directly links the efficiencies to the angular-momentum transport of these simulations.

2.5. Flux eruptions

We averaged the data over a time span in which a MAD state occurs in each simulation from a subset of three simulations: K.80, K.90, and K.99. In Fig. 7, we use a solid line to emphasize the mass-accretion rate and the magnetic flux over the time spans under investigation for each simulation. The time spans of interest in this study were chosen so that the magnetic fluxes through the event horizon have a similar profile, especially for the maximum and minimum amplitudes to be comparable and the total time span to be Δt = 600 M. These choices affected the comparability of the different models, especially when interpreting averaged plasma parameters and ray-traced images of azimuthally asymmetric structures, created by the cyclic MAD states. The notion of better comparability in the indicated time intervals in Fig. 7 is based on the fact that the simulations undergo very similar increases and decreases in accreted magnetic field, which, at the peak of ϕ_BH, causes the arrestation of the disk. This arises from an advectively driven phase of accumulation of magnetic field in the vicinity of the BH, which then overreaches the threshold magnetic pressure to be accreted, and outward diffusion of the accreted plasma takes over (Jacquemin-Ide et al. 2025). This is observed to occur in each model at very similar evolution states of the simulations (cf. Fig. 7). The build up of magnetic flux in the beginning of each interval is very similar, as every simulation accretes in a non-arrested way over a very similar time span, as is given for the duration of the arrested states of the disks in all three simulations. The MAD state (or magnetic-flux eruption) is indicated by the decrease in the mass-accretion rate (cf. Fig. 7) after both magnetic fluxes and mass-accretion rates reach their peak in every simulation. We note that these time intervals are the only ones within the whole simulations, where we see such strong similarities in the magnetic flux profiles, especially at those similar and importantly late evolution times (∼17 000 M) in the simulations. Since the electromagnetic contribution of the jet power depends quadratically on ϕ_BH (Tchekhovskoy et al. 2011), the EM energy flux contributions to the jet can be ruled out from being a result of strongly differing ϕ_BH over the three simulations.

Fig. 7. Interval of investigation in order to compare the jet morphologies of all simulations. The magnetic flux, ϕBH, describes a similar curve in height amplitude and width in time intervals for all simulations (indicated by solid lines). The dashed lines solely indicate time intervals that are not of interest for the discussions of the jet morphology.

We wanted to characterize the magnetic flux eruptions or MAD states, since their impacts on the morphology are essential for consideration in the subsequent analysis. We show, respectively, simultaneous snapshots of the meridional and equatorial planes in the GRMHD simulations, where the flux eruptions are at their respective peaks in Figs. 8 and 9. We note that the comparison of the snapshots lies outside the time intervals of comparison in Table 2; nevertheless, the MAD states already impact the simulations inside these intervals, and we show the peak MAD states to emphasize these effects. We want to demonstrate the maximal impact, that the flux eruptions can have on the plasma parameters and on the radiative signatures. The timescales from the onset of these magnetically arrested states to their peak flux eruptions and, finally, termination of the arrested disk do not differ significantly across the three models. This in turn aids the comparability of the time intervals (compare Table 2) chosen for time-dependent and time-averaged quantities. We want to emphasize that although the amplitude of magnetic flux that is built up and released is very similar across all simulations (cf. Fig. 7), the flux tubes (cf. Porth et al. 2019) do vary strongly in terms of their electron temperature, magnetization, and kinetic-to-magnetic pressure ratio (see Figs. 8 and 9). In the simulation K.80, these magnetic eruptions occur two-sided in this instance, whereas both K.90 and K.99 are completely azimuthally asymmetric (cf. Fig. 8), and the one-sided flux eruptions are much more stable in the latter simulations. K.99 shows a significantly extended flux tube that is hot enough to imprint clear signatures of thermal synchrotron radiation in the ray tracings (see Sect. 3 and Fig. G.1). The jet structure is also highly dependent on the MAD states as these flux tubes are able to contribute to the jet morphology significantly (see Fig. 9; K.99 is again rotated azimuthally by 180° to compare the MAD states), and, especially, they introduce high degrees of azimuthal asymmetries, which is important to consider when we investigate time-averaged quantities and morphologies in the following sections.

Fig. 8. 2D slices of the equatorial plane showing plasma β, σ, and Te of the peak MAD state of the simulations K.80, K.90, and K.99. Note the extended, arrow-like tail of the flux tube in K.99 that is twice the length compared to that of K.90.

Fig. 9. Meridional 2D slices of peak MAD state showing plasma β, σ, and Te of the simulations K.80, K.90, and K.99 (equal times to those in Fig. 8 for each simulation, respectively). We show the 180° azimuthally rotated cross-section of K.99, for comparison of the MAD states.

2.6. Comparison of the jet structure

We investigated the jet structure and its variability by means of the magnetization, σ, and the Bernoulli parameter, −hu_t. The contour at σ = 1 ≡ σ_cut defines our cutoff of the jet spine in the ray tracings shown in Sect. 3 and can be interpreted as our inner boundary of the jet sheath. The Bernoulli parameter at −hu_t > 1.02 separates the outflowing plasma in the jet from inflowing or gravitationally bound plasma in the disk and wind, so we defined the −hu_t = 1.02 contour as the outer boundary of the jet sheath. In Fig. 10, we show meridional slices (x − z-plane) in 45° increments in the azimuthal angle, ϕ, of the σ_cut = 1 contour (top row) and the Bernoulli parameter contour at −hu_t = 1.02 (bottom row), which are averaged in the x-direction over the respective time intervals listed in Table 2 for each simulation. Moreover, to emphasize the structural variability we show the standard deviations in the x-direction, indicated by the colored areas. Note that by x-directions we mean every axis that is perpendicular to the BH spin vector (in our simulations aligned with the z-axis). We explicitly did not average over phi, since the MAD states that are undergone in the chosen time intervals (see Fig. 7) introduce strong azimuthal asymmetries that would otherwise smear out.

Fig. 10. Upper row: σ = 1 contours in xz-cross-section corresponding to K.80 (green lines), K.90 (blue lines), and K.99 (red lines), and the 1σ standard deviation in the x-directions over the time intervals listed in Table 2. Each column presents poloidal cross-sections at azimuthal angles in 45° increments. Bottom row: Bernoulli parameter −hut = 1.02 contours, where we otherwise show the exact analogs for it as in the σ = 1 contours.

In Fig. 10, the runs K.80 and K.90 trace similar jet boundaries, whereas K.99 manifests a consistently wider jet for ϕ = {0° ,45° } and a narrower jet for ϕ = {90° ,135° }. This hints at a rather elliptical paraboloid jet geometry in the K.99 run, likely induced by the stronger MAD state that causes the azimuthal axis asymmetry. We note that we rotated K.99 azimuthally by 180°, such that the flux eruptions of the simulations match in azimuthal angles. Interestingly, at ϕ = 135° all runs consistently show the strongest axis asymmetries, which coincides roughly with the foot point of the flux eruption (compare to Fig. 8, where the x-axis is at ϕ = 0° /180°).

The 1σ band in Fig. 10 shows consistent structural properties as the Bernoulli contour. The sigma contour is overall significantly narrower downstream of the jet and less variable for K.99, especially at distances of ∼100 M − 200 M from the center, compared to K.80 and K.90.

3. General-relativistic radiative transfer

In order to investigate the radiative signatures of the 3D GRMHD simulations, we calculated the polarized synchrotron emission by means of a general-relativistic radiative-transfer (GRRT) post-processing step. We therefore employed IPOLE (Mościbrodzka & Gammie 2017): a semi-analytic code for relativistic polarized radiative transfer.

We assumed a model for the electron distribution function (eDF), since our GRMHD simulations evolve the electrons neither simultaneously nor separately; rather, they are assumed to tightly follow the evolution of the simulated plasma. The model chosen for this study combines the thermal Maxwell-Jüttner distribution,

$$\begin{array}{c}\hfill \frac{1}{n_{\mathrm{e}}}\frac{\mathrm{d}{n}_{\mathrm{e}}}{\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{e}}}=\frac{{\mathrm{\Gamma }}_{\mathrm{e}}^2\sqrt{1-1/{\mathrm{\Gamma }}_{\mathrm{e}}^2}}{{\mathrm{\Theta }}_{\mathrm{e}}{K}_2(1/\mathrm{\Theta })},\end{array}$$(24)

and the nonthermal κ-eDF (see Pandya et al. 2016), which smoothly converges to the Maxwell–Jüttner distribution for low Γ_e:

$$\begin{array}{c}\hfill \frac{1}{n_{\mathrm{e}}}\frac{\mathrm{d}{n}_{\mathrm{e}}}{\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{e}}}=N{\mathrm{\Gamma }}_{\mathrm{e}}\sqrt{{\mathrm{\Gamma }}_{\mathrm{e}}^2-1}{(1+\frac{{\mathrm{\Gamma }}_{\mathrm{e}}-1}{\kappa w})}^{-(\kappa +1)},\end{array}$$(25)

with the normalization factor N; the dimensionless electron temperature Θ_e = k_bT_e/m_ec²; the modified Bessel function of the second kind, K₂; the power-law index κ; and the width parameter w.

The thermal eDF is the relativistic generalization of the thermal Maxwell distribution, and the κ-eDF was inspired by observations of the solar wind and by results of collisionless plasma simulations (e.g. Kunz et al. 2015) and converges to the thermal distribution for κ → ∞. Both distributions need an assignment of the electron temperature, T_e, which we define as

$$\begin{array}{c}\hfill T_{\mathrm{e}}=\frac{p{m}_{\mathrm{p}}}{R_{\raisebox{1ex}{$\mathrm{p}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{e}$}\right.}\rho m_{\mathrm{e}}}=\frac{T_{\mathrm{code}}}{R_{\raisebox{1ex}{$\mathrm{p}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{e}$}\right.}}\frac{m_{\mathrm{p}}}{m_{\mathrm{e}}},\end{array}$$(26)

where k_B, m_p, and m_e are, respectively, the Boltzmann constant and the proton and electron masses. ρ and p are the rest-mass density and pressure from the GRMHD simulation, respectively. The code temperature, $T_{\mathrm{code}}=\raisebox{1ex}{$p$}\!\left/ \!\raisebox{-1ex}{$\rho $}\right.$, is thus defined by plasma parameters from the GRMHD code as well, and it describes solely the ion temperature in code units. The ratio of proton temperature T_p and electron temperature T_e, R, is defined as

$$\begin{array}{c}\hfill R_{\raisebox{1ex}{$\mathrm{p}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{e}$}\right.}=\frac{T_{\mathrm{p}}}{T_{\mathrm{e}}}\equiv R_{\mathrm{high}}\frac{{\beta }^2}{{\beta }^2+1}+R_{\mathrm{low}}\frac{1}{{\beta }^2+1},\end{array}$$(27)

with the plasma parameter $\beta =2\frac{p}{b^2}$ provided by the simulation. The so called R-β model then has two parameters in total, R_high and R_low (see Mościbrodzka et al. 2016), that together with the code temperature control the electron temperature. Notice, that we have to assert an electron temperature model as the GRMHD simulation does not intrinsically provide an electron temperature.

The kappa distribution in Eq. (25) is defined by its width, w, which is now attributed a fraction ε of the magnetic energy besides the thermal energy (see Davelaar et al. 2019; Fromm et al. 2022):

$$\begin{array}{c}\hfill w:=\frac{\kappa -3}{\kappa }{\mathrm{\Theta }}_{\mathrm{e}}+\frac{\epsilon }{2}[1+tanh(r-r_{\mathrm{inj}})]\frac{\kappa -3}{6\kappa }\frac{m_{\mathrm{p}}}{m_{\mathrm{e}}}\sigma .\end{array}$$(28)

The injection radius, r_inj, is the distance from the BH, where electrons with magnetic energy contribution are injected and $\sigma =\raisebox{1ex}{$b^2$}\!\left/ \!\raisebox{-1ex}{$\rho $}\right.$.

The non-thermal contributions to the electron energies are thought to be originated in magnetic reconnection events, in which current sheets accelerate the charged particles. These processes cannot be resolved in global GRMHD simulations, therefore, the κ-parameter depends on the magnetization and plasma β when employing a PIC simulation based sub-grid model from Ball et al. (2018):

$$\begin{array}{c}\hfill \kappa :=2.8+0.7{\sigma }^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}+3.7{\sigma }^{-0.19}tanh\left(23.4{\sigma }^{0.26}\beta \right).\end{array}$$(29)

With that we again follow Fromm et al. (2022) and Davelaar et al. (2019). (In the former study, a more detailed parameter space study on the R-β model can be found).

3.1. M87 ray-tracing setup

For the 86 GHz Ray-Tracings, we fixed the parameters of all runs to the mass of the BH of M_BH = 6.2 × 10⁹ M_⊙ (cf. Gebhardt et al. 2011; EHTC 2019c,a,b), the distance to the BH to D = 16.9 Mpc (cf. Blakeslee et al. 2009; Cantiello et al. 2018; EHTC 2019a,b) and the viewing angle i = 163° (as spanned by the BH spin axis and line of sight) of the BH and jet in M87. The parameters for the electron temperature model are set to R_high = 160 and R_low = 1 and for the κ-eDF in Eq. (25), we choose a non-thermal emission fraction of ε = 0.5 and a injection radius of non-thermal particles at r_inj = 10 M. We note that the numerical approximations of the emission and absorption coefficients are only valid in the range 3.1 ≥ κ ≤ 8, where solely the κ-eDF is applied, and the Maxwell-Jüttner eDF otherwise.

This means that we mainly calculate non-thermal emission in the the jet-sheath and it is noted that the Flux-tubes can also be sufficiently magnetized, such that they emit non-thermally.

A polar cut of θ_cut = 4° (opening angle of the cut-out cone is α_cone = 8°) is applied for the ray-tracings as shown before in the GRMHD data (cf. Appendix B, Figs. B.1 and A.1). Additionally, we omit the emission from the jet-spine in all regions where σ > σ_cut = 1. This is due to the fact that the highly dominant magnetic energy in the jet-spine is prone to introduce small errors in the total energy evolution, drastically altering the internal energy, hence the plasma temperature.

The time-averaged flux of each ray-traced simulation is adjusted such that the flux at horizon scale at 230 GHz matches the observed 0.9 Jy by the EHTC (2019c) by means of optimizing the mass-accretion rates in the GRMHD simulations.

3.2. Comparing jet structure and polarization at 86 GHz

In Fig. 11 we present polarized, time-averaged ray-tracings of the three simulations up to the 0.5 mas-scale, averaged over the time-intervals as used before (cf. Table 2). These images are the result of a line-integral convolution (LIC) of the electric vector position angle (EVPA) with a field of random generated numbers. This generates stream lines of the EVPA vector field, encoding directional information of the field lines. The intensity I is then super-imposed onto the LIC of the EVPAs, where the opacity of those is a function of the linear polarization (LP) fraction. The image is then rescaled so that the maximum of the convolved image max[f(I,  m,  EVPA_{field − lines})] equals max(I). The LP fraction is given by:

$$\begin{array}{c}\hfill m=\frac{\sqrt{Q^2+U^2}}{I}\times 100[\%].\end{array}$$(30)

Fig. 11. Time-averaged 86 GHz ray tracings of the M87 jet. Each simulation is averaged over the time intervals corresponding to each simulation that are listed in Table 2. Stokes I is superimposed with a line integral convolution of the EVPA vector field by scaling the opacity as a function of the linear polarization fraction and rescaling the image by the maximum of the convolved image max[f(I,  m,  EVPAvector − field)] to max(I).

Qualitatively, one can tell by the dominance of the stream lines over the smooth Stokes-I in Fig. 11, that K.99 exhibits the highest linear polarization of all runs. Quantitatively, we show the LP fraction and I in Fig. 12 as a function of time, offset by a time T_i, such that the time-intervals in Table 2 have a common starting time for comparison purposes. The LP fraction peaks at ∼16% for K.99 compared to ∼11% in the normal accretion states for both K.80 and K.90. The LP-pattern itself of K.99 differs as well, in that the EVPA longitudinal component (to the jet position angle) inside the jet is clearly dominant, as can be seen by the straight lines over a wide area inside the jet. Comparing this observation to K.80 and K.90, those on the other hand, exhibit visible transverse components in most parts of the EVPAs. At the jet edge, this behaves analogously for K.99, where the EVPAs now sharply transition from following the jet longitudinally in the inside to a dominant EVPA-field perpendicular to the jet edge. This hints at a highly ordered magnetic field in the jet spine.

Fig. 12. Top row: 86 GHz light curve for simulations K.80 (green), K.90 (blue), and K.99 (red) in the time intervals shown in Table 2, but shifted by an offset for better comparability. Bottom row: Linear polarization fraction with the same color-code as before.

Notably, the jet structure of K.99 is very distinct from the usual paraboloidal shape that is associated with MAD-simulations. Beyond the jet base, about ∼30 μas from the core, the jet transitions into a cone-like shape, which will be discussed in more detail in the next section.

3.3. Jet edge and width profiles

In Fig. 13, we trace the profile of the jet analogously to the structural investigation of the jet sheath in Section 2.6. We extract the time-averaged contours at ⟨I⟩ = 5 × 10⁻⁶ Jy (solid lines) and illustrate with the corresponding colored areas the variability perpendicular to the jet position angle (PA), i.e. the 1σ-bands of the variation in jet-PA-perpendicular direction of the I = 5 × 10⁻⁶ Jy contour over time. Therefore, the images were de-rotated to a jet PA of −180° (South-North direction) such that the jet can be subdivided into jet PA perpendicular slices for each pixel along the jet direction. The closest pixel to the contour is determined at each time step. Eventually, the mean and standard deviation at each slice over time can be calculated. The chosen contour of I = 5 × 10⁻⁶ Jy is a sensible tracer of something that could be considered the edge of the jet, which is in itself loosely defined, since it traces the apparent jet edge visible from Fig. 11 at a dynamic range of about four orders of magnitude. In observations of the jet of M87* often the observed edge-brightening is the best indicator of the jet edge (often taken to be the HWHM towards the outside of the jet of the bright Gaussian edge features). In these simulations on most scales the jet perpendicular flux-profile does not exhibit bi- or multimodal distributions, especially not time-averaged, such that bright spots are not as pronounced on the edge of the jet.

Fig. 13. ⟨I⟩ = 5 × 10−6 [Jy] contours of time-averaged 86 GHz ray tracings (not the line integral convolved images shown in Fig. 11) for simulations K.80 (green), K.90 (blue), and K.99 (red). The colored areas are 1σ bands of the standard deviations in the toroidal direction.

The mean contours of K.80 and K.90 are very similar, whereas K.99 stands out by its partially cone-like shape as pointed out in the last section. The actual difference in the asymmetries of the models becomes very apparent. In a basis defined by coordinates (x, z) as in Fig. 13, we find that the asymmetry is strongly shifted to negative x-values in K.99. The variability of the jet structure provides confidence that the jet structures shown in Fig. 11 are not distorted because of strong time-dependent morphological variations. The variation of these contours in x-direction slowly increase in z-direction in the cases of K.80 and K.90, but increase rather abruptly and asymmetrically for K.99 on the positive x half-space at ∼250 M and on the negative one at ∼200 M. This can be explained by a helical arc propagating through the jet sheath, seen in Fig. 11 downstream the jet in K.99. Nevertheless, the jet shape of K.99 is completely altered in comparison to K.80 and K.90.

From these jet profiles of the 86 GHz ray-tracing we can now extract the jet widths as a function W(z) of the jet longitudinal axis z (as defined before) and compare these to the observations of Lu et al. (2023) shown in Fig. 14. We note that the data-points of the jet-width measurements do not show the error-bars in Lu et al. (2023), hence we cannot elaborate on the statistical certainty of this comparison. From Fig. 14, it is apparent that the theoretical predictions especially of K.80 and K.90 do match the observations well at the scales below 100 μas, contrarily to the claims of Lu et al. (2023) that MAD simulations may have difficulty explaining the observed jet width at these scales. At the scales shown above 100 μas the simulations K.80 and K.90 converge to the extrapolated W(z)∝z^0.58 power law fit from large scale observations (see Asada & Nakamura 2012; Hada et al. 2013, 2016). K.99 follows the power-law further up-stream and exhibits a sudden increase of the jet width at ∼40 μas, whilst even closer to the BH it becomes almost constant, but even wider than the other models. This sudden increase in the width of the jet can be explained by the counter-jet, which exhibits a wider opening angle and thus shifts the jet edge of the counter-jet to positive z-values up to ∼100 M (cf. Fig. 11, right image). For the common choice of the R-β electron temperature- and κ non-thermal electron prescriptions chosen here K.99 exhibits a steeper slope than this observation suggests, rather following the large-scale power-law. The new observational data points at the scales closer to the jet launching region suggest that K.90 satisfies these observations very well, as the data points are very well within the variability of K.90 (cf. green band in Fig. 14). We also see that this holds mostly for the common model K.80, as well.

Fig. 14. Jet width of M87 at 86 GHz along the longitudinal jet axis for simulations K.80 (green), K.90 (blue), and K.99 (red) with corresponding 1σ bands. The black dots are data points of the jet width of M87 from the observation of Lu et al. (2023). The dotted line is a power-law fit from observations of M87. The solid horizontal line is the estimated 86 GHz ring size, and the dashed horizontal line is the 230 GHz ring size with respective error band (gray areas). The vertical lines are the distances from the BH, where the intrinsic half-opening angles equals the jet viewing angle of θv = 17° (color-code is consistent).

3.4. Ray-tracing the flux eruptions

Since we mostly study the time-averaged quantities of the simulations presented in this work to extract general physical differences from the three models, we have to also shed light on the influence of the episodic magnetic flux eruptions on the radiative signatures (cf. in Sect. 2.5 we discussed the consequences of the flux eruptions on the GRMHD quantities). In Fig. G.1 we show an example of a ray-tracing of such a flux eruption in the case of model K.99, where on the left we only ray-traced the disk emission and on the right the full emission that includes emission from the jet, as well. By comparison of both images, the influence of the flux eruption on the jet morphology becomes clear. The jet base strongly deviates from a paraboloid as the flux eruption envelopes the jet in a azimuthally asymmetric manner (cf. the feature in the upper half of the right image in Fig. G.1). Since these flux eruptions are intermittent phenomena, they do not influence the radiative imprints on the jet morphology steadily and azimuthally symmetric, and are thus important to consider for the study of the long time-scale structural studies of MAD simulations. In order to be able to average out the strong azimuthal asymmetries introduced by MAD states, they have to be averaged over multiple ones. Since they occur quasi-periodically with a frequency of the order of ∼1000 M, this would make studies on time-scales of at least 10 000 M necessary. Hence, we limited our discussion to one MAD event per simulation.

4. Conclusions

We find that the initial magnetic energy content of the torus systematically alters the jet energetics and morphology of the developed simulations. It mostly affects the EM contribution in the jet, and from that it gains ∼60% efficiency in the conversion from rest-mass energy to the energy of the outflow by means of increasing the magnetic energy of the torus at the initialization. This in turn can make nonthermal particle acceleration from vacuum gaps and subsequent pair cascades more likely, due to the increased EM field in the jet spine (see Levinson & Rieger 2011; Wendel et al. 2021 for details on the vacuum-gap model and applications to low-luminosity AGNs such as M87 and Mrk501). We also find, that an increased magnetic filling of the FM torus elevates the outward angular-momentum transport, primarily in the jet, which is driven by Maxwell stresses. Comparing the angular-momentum transport of BHs with very high spin presented here to those of Chatterjee & Narayan (2022), in which the case of a non-rotating BH in a standard MAD-setup (F = 0.8) was investigated, the major differences are the following: (i) the total radial angular-momentum flux is dominated by radially outward-directed flux in the jet and wind, mostly originating from the jet base; and (ii) morphologically, the angular-momentum vector fields are distinct, in that the high-spin BH creates field lines that point radially outwards in a polar angle range from ∼(0 − 70)° or ∼(110 − 180)°. Toward the equatorial plane, the field lines transition from a monopole geometry to a quadrupole geometry in the disk and torus. The dominant, Maxwell-stress-induced outward polar flux is redirected to an equatorially inward, advectively driven flux. In contrast, the case of a non-rotating BH exhibits dominant, advectively driven inward angular-momentum flux in a larger portion of the disk and torus, whereas the outward-directed Maxwell-stress-induced flux is subdominant in the jet, such that the net angular momentum is transferred to the central BH. The outward angular momentum transport also covers a much smaller polar angle compared to simulations with high spin, contributing far less to the strength of outward winds. Our simulations of MADs with high-spin BHs demonstrate the effect of extracting the angular momentum of the BH to the angular momentum of the outflowing plasma in the jet and wind, as originally proposed by Blandford & Znajek (1977). We further find that there is a maximum of all net angular momentum fluxes inside the ergopsheric region of the BH, as opposed to the net-constant fluxes for MAD simulations with no spin.

It appears that the disk’s vertical flux tubes created by the magnetically arrested disk affect the jet structure more significantly for high filling factors. From our sampled MAD states, we find a tendency for more spatially extended flux tubes that affect the jet morphology by introducing significant azimuthal asymmetries in the case of K.99; these propagate along the jet at least up to 200 M. The jet’s sheath structure is more stable on these scales; i.e., the jet structure exhibits lower variability in the jet-perpendicular direction. However, since we restricted the jet-morphology analysis to a time interval, in that a single flux eruption occurs in each simulation, further detailed, statistical analyses of the impact of the filling factor on the quasi-periodic flux eruptions could provide valuable insights into the observable footprints that they have.

The radiative transport calculations also show that the jet’s sheath structure is more stable and collimated up to ∼200 μas, supported by the higher EM energy flux contribution in K.99. We find significantly elevated LP fractions of up to ∼15%, in agreement with observations of the EHTC (2021b). The jet width relation over the jet axial distance shows that GRRT calculations of GRMHD simulations of magnetically arrested accretion flow are able to match the observed jet width profile of Lu et al. (2023) very well for the reference models. For K.99, we find deviations from the jet width, as it follows the observed large-scale power-law relation W ∝ z^0.58 to shorter distances, whereas the observations follow a sub-power-law relation. The deviation of K.99 from the observations, however, is a result of the different techniques used for the determination of the jet width. In the observations of Lu et al. (2023), the HWHM of the outermost components perpendicular to the jet longitudinal axis determines the jet width, whereas we did not recreate the edge brightened jet morphology and thus used brightness contours to determine the jet width. This made the detailed jet–width relation difficult to compare to these observations; nevertheless, we show that all our MAD simulations provide a wide jet base substantially wider than the pure Blandford–Znajek (BZ) jet (see Blandford & Znajek 1977) for a high-spin BH (a → 1). It is evident from our analysis of the jet energy fluxes that the almost evacuated, highly magnetized jet spine (σ ≳ 1) is dominated by EM energy flux, and in the enveloping jet sheath the wind contributes to the outflowing kinetic energy flux that increases with distance. This jet-sheath region (σ ≲ 1 and the gravitationally unbound part of the plasma) is the origin of the jet emission, and the σ = 1 contour is thus not a good tracer for the outermost jet edge that could be observed. We remind the reader that σ = σ_cut = 1 is the value for the emission cutoff for the jet spine, where we do not expect significant synchrotron emission in the radio regime because of the extremely low electron densities, as well as efficient synchrotron cooling, due to the presence of high magnetic-field strengths in σ ≫ 1 regimes. We stress that a simple extrapolation and comparison of plasma parameters such as the magnetization from simulations to trace the edge or width of the jet observed at millimeter wavelengths (see Lu et al. 2023, Fig. 3; gray band shows σ contours from force-free electrodynamic simulations of a BZ jet for different BH spins) is insufficient to approximate the complex emission processes, which are yet to be fully understood. From the GRRT of our GRMHD simulations for the jet of M87*, it is evident that the apparent width of the jet base is not anchored at the ergosphere of the BH, as opposed to force-free electrodynamic simulations.

Overall, we find the jet of the K.99-model to exhibit more flux downstream of the jet compared to the reference models. This aids in solving the issue that GRMHD simulations produce jets that decline too much in flux downstream compared to observations. Hence, this model is well suited for modeling jets on larger scales, beyond the jet launching region. This is very helpful, since most SMBH systems in the Universe are of smaller angular size than M87*. With the enhanced jet power in this model, the jet remains brighter and more collimated on scales even beyond 1000 M.

The magnetic filling factor regulates the spatial extent of the initial vector potential and thus that of the magnetic field, as well as the total magnetic energy inside the initial torus configuration of our simulations, which are commonly used to model magnetically arrested accretion and also constitute the main portion of the simulation library used for M87 by the EHTC (2019b). In general, we conclude from this study that the evolution of these MAD simulations is very sensitive to minor changes in the initial magnetic field in the FM torus, which was set by the filling factor. In order to keep track of the vast parameter space of MAD simulations and consequences of the initial conditions, the specification of the detailed initial vector potential of simulations is therefore paramount.

Acknowledgments

We thank Charles F. Gammie for the helpful discussion regarding this work. This research is supported by the DFG research grant “Jet physics on horizon scales and beyond” (Grant No. 443220636) within the DFG research unit “Relativistic Jets in Active Galaxies” (FOR 5195). The numerical simulations and calculations have been performed on MISTRAL at the Chair of Astronomy at the JMU Wuerzburg. YM is supported by the National Key R&D Program of China (grant no. 2023YFE0101200), the National Natural Science Foundation of China (grant no. 12273022, 12511540053), and the Shanghai municipality orientation program of basic research for international scientists (grant no. 22JC1410600).

References

Appendix A: Initial magnetic-field configuration and filling-scaling of the tori

Here, we calculate the initial covariant volume fractions of the volume containing the cells in which the initial invariant magnetic-field strength b is defined and the initial volume enclosing the full mass density distribution of the torus. The initial magnetic field is normalized such that the minimum of plasma β in the torus is β_min = 100 in all simulations. We calculate the total and shell-integrated covariant volumes of the torus and enclosed magnetic-field as follows:

$$\begin{array}{cc}& V_{\mathrm{t}}={\displaystyle {\int }_{\{rho>8\times {10}^{-7}\}}}\sqrt{-g}drd\theta d\varphi ,\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& V_{\mathrm{b}}={\displaystyle {\int }_{\{b>1\times {10}^{-7}\}}}\sqrt{-g}drd\theta d\varphi .\hfill \end{array}$$(A.2)

We chose the masking conditions for the enclosed volume of the density distribution such that ρ is above the density distribution in the atmosphere (the space excluding the torus and BH), that should have its minimum at ρ_min = 10⁻⁶ (as stated in Sec. 2), however, increases up to ρ ∼ 8 × 10⁻⁷ in the vicinity of the BH. For b, the condition is not as strict, since its numeric value is only defined, where the four-vector potential A_ϕ contributes to the curl to generate the magnetic field. We choose a low threshold of b > 1 × 10⁻⁷ and set the non-defined values that are mostly outside the torus to 0 in order to capture the full volume of the initial magnetic field. The actual filling factor is then given by:

$$\begin{array}{c}\hfill f_{\mathrm{a}}=\frac{V_{\mathrm{b}}}{V_{\mathrm{t}}}·\end{array}$$(A.3)

We note that, f_a is a measure for the volumetric magnetic-field extent, or the magnetic filling of the initial magnetospheric configuration, that we calculate for each filling factor F and is not equivalent to the latter values as can be seen in Tab. A.1.

The total volume containing magnetic fields is generally quite close to the constant torus volumes throughout the simulations, much closer than the filling factors F suggest. However F, simply scales the offset of the initial vector potential A_ϕ (cf. Eq. (6)), essentially increasing its offset from the zero-cutoff maximum function in Eq. (6) as can also be seen in Fig. 1. Still, increasing F indeed leads to an increased volume that contains magnetic fields. Hence, we use the filling factor synonymously with the scaling of the magnetic-field extent and the magnetic energy in the FM-torus.

Fig. A.1. Initial disk-vertical magnetic flux permeating the disk. The colors correspond to simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta) and K.99 (red).

Appendix B: Average MHD parameters

Fig. B.1. Time and azimuthally averaged plasma (MHD) parameters of model K.99. The time is averaged over the interval indicated in Fig. 7 and listed in Tab. 2. We apply a polar cut of θcut = 4° in the post-processed ray tracings for all models, which is indicated by the black regions at the poles.

Appendix C: Mass-accretion rate

Fig. C.1. Mass-accretion rate Ṁ5M of simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta) and K.99 (red). The horizontal lines are the averages and the shaded area span the 2σcorr confidence interval of the correlated time series. The overlapping error-bands show that the mass-accretion rate is sufficiently similar in between all simulations, thus, indicating that these are in a steady-state.

Appendix D: MRI quality factors

Fig. D.1. Disk averaged MRI quality factors in θ and ϕ of each simulation, row-wise. Dashed lines are the median values taken over the time span shown in the graphs. Besides the median values we also give corresponding time averaged values. Due to the skewed nature of the temporal distriubtion of the disk averaged Q factors, the median is more appropriate for interpretation. It is noted that for Q(θ) ∼ (3 − 6)% of points lie outside the graph, i.e. are above 100 in the three simulations. For Q(ϕ) ∼ 1% are concerned.

Fig. D.2. Meridional cross-section of MRI Quality factors in θ and ϕ directions. Top row: Meridional cross-section. Bottom row: Volume averaged Q factors over the full azimuthal angle. The black lines indicate the regions where the disk is cut, poloidally, and subsequently averaged over the full azimuthal angle to obtain the disk averaged values given in the bottom left corners.

Appendix E: Error estimation of time series

Here we elaborate on the error estimation of the correlated time series, inspired by Lalakos et al. (2025). The reduction of the number of independent samples can be quantified by the integrated autocorrelation time τ_int of a time series f(t), which is defined as:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{auto}}=(1+2{\displaystyle \sum _{k=1}^N}\rho (k\mathrm{\Delta }t))\mathrm{\Delta }t,\end{array}$$(E.1)

where ρ(t) is the normalized autocorrelation function of the unity normalized time series f(t) at lag t, N is the number of samples and Δt is the sampling interval (5M in our case) for normalization. The normalized autocorrelation function is defined as:

$$\begin{array}{c}\hfill \rho (k\mathrm{\Delta }t)=\frac{1}{{\sigma }^2}{\displaystyle \sum _{i=1}^{N-k\mathrm{\Delta }t}}(f(i\mathrm{\Delta }t)-\mu )(f(i\mathrm{\Delta }t+k\mathrm{\Delta }t)-\mu ).\end{array}$$(E.2)

The autocorrelation function is computed up to a lag of N/2 to avoid noise-dominated estimates at large lags and normalized to lie between values of −1 and 1 to enable comparability between the simulations.

The error of the mean μ of a correlated time series of duration T is then given by:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{corr}}=\sigma \sqrt{\frac{2{\tau }_{\mathrm{auto}}}{N}},\end{array}$$(E.3)

where σ is the standard deviation of the time series. We typically found values of τ_auto ≲ 10M corresponding to an error ≲5 higher than the standard error of the mean, when assuming the data are entirely uncorrelated. Therefore, we avoid underestimating the uncertainty of the mean.

Appendix F: Angular momentum fluxes at varying radii

Fig. F.1. Polar slices of the radial angular momentum flux at varying radii increasing to the right. Solid lines are the radially outward flux and dashed lines are inward flux. Note the varying scales on the y-axis.

Appendix G: Morphology of a flux eruption snapshot

Fig. G.1. 86 GHz ray-tracings of the flux eruption of the K.99 run (same procedure as in Fig. 11 for the line integral convolution). Left image: only ray-traced thermal disk emission. Right image: ray-traced non-thermal jet and thermal disk emission

All Tables

All Figures

Fig. 1. Invariant magnetic-field strength, b; the spatial magnetic-field components, Br and Bθ; and four-vector-potential component, Aϕ, over the polar angle at the initialization step (t = 0 M). Colors correspond to simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta), and K.99 (red). In the middle row, the solid lines show Br and the dashed lines Bθ. We chose polar slices at radii close to the inner edge of the torus at rin = 20 M and at the pressure maximum in the tori at r = 50 M. Note that in columns going from left to right, i.e., radially outwards, the simulations with higher filling factors show increased magnetic-field strengths and extended magnetic field, both radially and in the polar-angle, especially at the inner edge of the torus.

In the text

	Fig. 2. Mass-accretion rates, Ṁ, and magnetic flux through the event horizon of the BH, ϕBH, and the MAD-parameter, $\mathrm{\Psi }={\varphi }_{\mathrm{BH}}/\sqrt{\dot{M}}$, for all simulations over time.
In the text

	Fig. 3. Temporal evolution of total energy flux (solid lines) and EM energy flux (dashed lines) measured at r = 10 M in the jet for simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta), and K.99 (red). The horizontal lines show the average of the total energy fluxes for each run over the shown time interval.
In the text

Fig. 4. Time-averaged total energy fluxes (dash-dotted lines), EM energy fluxes (solid lines), and kinetic energy fluxes (dashed lines) in the time interval (25 000 − 30 000) M for varying BH separations along the jet for simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta), and K.99 (red). The vertical dotted lines indicate the radii at which the kinetic and EM energy fluxes are in equipartition.

In the text

	Fig. 5. Jet efficiencies for simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta), and K.99 (red), where the dashed horizontal lines show the time-averaged efficiencies and the correlation-corrected 2σcorr error band over the shown interval (listed in Table 1).
In the text

Fig. 6. Left panel: Time- and azimuthally averaged density distribution in the meridional plane and total angular-momentum flux streamlines (color-coded by a log-scale of the absolute value; increasing from black to white; simulation: K.80). Middle column: Polar slices of the radial angular-momentum flux at r = 2 M, indicated by the dashed gray lines in the left panel and right column. Shown are the total, advective and Maxwell-stress-induced angular-momentum fluxes from top to bottom. Right column: Shell-integrated radial component of the angular-momentum fluxes. Solid lines indicate positive values, i.e., outward fluxes and dashed lines are the absolute values of inward (negative) flux.

In the text

	Fig. 7. Interval of investigation in order to compare the jet morphologies of all simulations. The magnetic flux, ϕBH, describes a similar curve in height amplitude and width in time intervals for all simulations (indicated by solid lines). The dashed lines solely indicate time intervals that are not of interest for the discussions of the jet morphology.
In the text

	Fig. 8. 2D slices of the equatorial plane showing plasma β, σ, and Te of the peak MAD state of the simulations K.80, K.90, and K.99. Note the extended, arrow-like tail of the flux tube in K.99 that is twice the length compared to that of K.90.
In the text

	Fig. 9. Meridional 2D slices of peak MAD state showing plasma β, σ, and Te of the simulations K.80, K.90, and K.99 (equal times to those in Fig. 8 for each simulation, respectively). We show the 180° azimuthally rotated cross-section of K.99, for comparison of the MAD states.
In the text

Fig. 10. Upper row: σ = 1 contours in xz-cross-section corresponding to K.80 (green lines), K.90 (blue lines), and K.99 (red lines), and the 1σ standard deviation in the x-directions over the time intervals listed in Table 2. Each column presents poloidal cross-sections at azimuthal angles in 45° increments. Bottom row: Bernoulli parameter −hut = 1.02 contours, where we otherwise show the exact analogs for it as in the σ = 1 contours.

In the text

Fig. 11. Time-averaged 86 GHz ray tracings of the M87 jet. Each simulation is averaged over the time intervals corresponding to each simulation that are listed in Table 2. Stokes I is superimposed with a line integral convolution of the EVPA vector field by scaling the opacity as a function of the linear polarization fraction and rescaling the image by the maximum of the convolved image max[f(I,  m,  EVPAvector − field)] to max(I).

In the text

	Fig. 12. Top row: 86 GHz light curve for simulations K.80 (green), K.90 (blue), and K.99 (red) in the time intervals shown in Table 2, but shifted by an offset for better comparability. Bottom row: Linear polarization fraction with the same color-code as before.
In the text

	Fig. 13. ⟨I⟩ = 5 × 10−6 [Jy] contours of time-averaged 86 GHz ray tracings (not the line integral convolved images shown in Fig. 11) for simulations K.80 (green), K.90 (blue), and K.99 (red). The colored areas are 1σ bands of the standard deviations in the toroidal direction.
In the text

Fig. 14. Jet width of M87 at 86 GHz along the longitudinal jet axis for simulations K.80 (green), K.90 (blue), and K.99 (red) with corresponding 1σ bands. The black dots are data points of the jet width of M87 from the observation of Lu et al. (2023). The dotted line is a power-law fit from observations of M87. The solid horizontal line is the estimated 86 GHz ring size, and the dashed horizontal line is the 230 GHz ring size with respective error band (gray areas). The vertical lines are the distances from the BH, where the intrinsic half-opening angles equals the jet viewing angle of θv = 17° (color-code is consistent).

In the text

	Fig. A.1. Initial disk-vertical magnetic flux permeating the disk. The colors correspond to simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta) and K.99 (red).
In the text

	Fig. B.1. Time and azimuthally averaged plasma (MHD) parameters of model K.99. The time is averaged over the interval indicated in Fig. 7 and listed in Tab. 2. We apply a polar cut of θcut = 4° in the post-processed ray tracings for all models, which is indicated by the black regions at the poles.
In the text

Fig. C.1. Mass-accretion rate Ṁ5M of simulations K.60 (orange), K.80 (green), K.90 (blue), K.95 (magenta) and K.99 (red). The horizontal lines are the averages and the shaded area span the 2σcorr confidence interval of the correlated time series. The overlapping error-bands show that the mass-accretion rate is sufficiently similar in between all simulations, thus, indicating that these are in a steady-state.

In the text

Fig. D.1. Disk averaged MRI quality factors in θ and ϕ of each simulation, row-wise. Dashed lines are the median values taken over the time span shown in the graphs. Besides the median values we also give corresponding time averaged values. Due to the skewed nature of the temporal distriubtion of the disk averaged Q factors, the median is more appropriate for interpretation. It is noted that for Q(θ) ∼ (3 − 6)% of points lie outside the graph, i.e. are above 100 in the three simulations. For Q(ϕ) ∼ 1% are concerned.

In the text

	Fig. D.2. Meridional cross-section of MRI Quality factors in θ and ϕ directions. Top row: Meridional cross-section. Bottom row: Volume averaged Q factors over the full azimuthal angle. The black lines indicate the regions where the disk is cut, poloidally, and subsequently averaged over the full azimuthal angle to obtain the disk averaged values given in the bottom left corners.
In the text

	Fig. F.1. Polar slices of the radial angular momentum flux at varying radii increasing to the right. Solid lines are the radially outward flux and dashed lines are inward flux. Note the varying scales on the y-axis.
In the text

	Fig. G.1. 86 GHz ray-tracings of the flux eruption of the K.99 run (same procedure as in Fig. 11 for the line integral convolution). Left image: only ray-traced thermal disk emission. Right image: ray-traced non-thermal jet and thermal disk emission
In the text