© The Authors 2026

1. Introduction

The nearby radio galaxy M87 is an ideal target for studying black hole accretion and jet formation (Yuan & Narayan 2014; Blandford et al. 2019). It is believed to have a low accretion rate consisting of a geometrically thick, optically thin accretion disk in a radiatively inefficient accretion flow state (RIAF; Narayan & Yi 1995a; Yuan & Narayan 2014). The emission at millimeter wavelengths of such low-luminosity accretion flows mainly comes from the synchrotron radiation of electrons (Ichimaru 1977; Narayan & Yi 1994; Yuan & Narayan 2014). The Event Horizon Telescope (EHT) collaboration, utilizing very-long-baseline interferometry (VLBI) techniques, successfully observed horizon-scale shadow images of M87^* at millimeter wavelengths (Event Horizon Telescope Collaboration 2019a, 2024).

General relativistic magnetohydrodynamic (GRMHD) simulation is a valuable tool for improving our understanding of accretion physics, and it can be used to simulate shadow images on horizon scales (e.g., Gammie et al. 2003; De Villiers & Hawley 2003; Noble et al. 2007; Mościbrodzka et al. 2009, 2014, 2016; Dexter et al. 2010; Shcherbakov et al. 2012; Davelaar et al. 2018, 2019; Mizuno et al. 2018; Dihingia & Fendt 2025). A ring-like structure is the most direct observational signature of a black hole shadow, as seen in synthetic images of the GRMHD simulations (Event Horizon Telescope Collaboration 2019b, 2022, 2025). Two distinct types of accretion flows have been identified. The first is the standard and normal evolution (SANE; e.g., Narayan et al. 2012; Sądowski et al. 2013) state, in which the magnetic field strength is weak, the angular momentum is transferred through magnetorotational instability (MRI), and the accreting gas does not become a magnetically arrested situation during the simulation. Therefore, the accretion process is smoother compared to the second type. The second type of accretion flow is the magnetically arrested disk (MAD; e.g., Narayan et al. 2003; Tchekhovskoy et al. 2011). In the MAD state, the magnetic flux accumulating near the horizon reaches the upper limits, and the accretion flow becomes magnetically arrested. However, this accumulated magnetic flux can be released substantially due to the magnetic flux eruptions, and the reduction of magnetic flux in the inner disk results in a temporary increase in the accretion rate until sufficient flux is advected again (e.g., McKinney et al. 2012; Sądowski et al. 2014). The intermediate state (INSANE) is also a possible third type of accretion flow. For instance, it is suggested to explain the transitions observed in X-ray binaries by Raha et al. (2026).

General relativistic radiative transfer (GRRT) calculations compute black hole shadow images utilizing GRMHD simulation data. The important parameters to model the electromagnetic emission are the electron distribution function (eDF) and the electron temperature, which are related synchrotron emission power. In their M87 GRRT calculation, the EHT collaboration assumed that the eDF is the Maxwell-Jüttner (thermal) distribution (Event Horizon Telescope Collaboration 2019b). However, the electron distribution can be affected by energy dissipation, particle acceleration, and thermalization (Yuan & Narayan 2014). Considering only thermal distribution ignores key processes such as magnetic reconnection and turbulent dissipation, which drive electrons toward a nonthermal power-law distribution (Ding et al. 2010; Hoshino 2013). In addition, some features of M87 caused by electron acceleration have been observed in the near-infrared and optical bands (Prieto et al. 2016). Recently, the application of the nonthermal κ distribution to MAD models has reproduced the wide opening angle jet morphology at 86 GHz and fit the broadband spectrum of M87 from the radio to the near-infrared bands (Cruz-Osorio et al. 2022). Compared to SANE models, Fromm et al. (2022) indicate that the combined constraint from broadband spectrum and jet collimation profile favors MAD model coupling with a κ distribution. Furthermore, Davelaar et al. (2023) show the propagation of waves along the shear layer of the jet wind using the κ distribution in the MAD regime, which provides a possible source to accelerate the electrons through turbulence or reconnection if the observation could confirm the imprints of such waves. In addition, Zhang et al. (2024) find that the jet emits more when a κ distribution is considered in the MAD model. They suggest that future multifrequency observations, which simultaneously resolve horizon-scale structure and the jet base, could be used to investigate the existence of nonthermal electrons. Recently, Tsunetoe et al. (2025) find that anisotropic nonthermal distribution functions can help produce sufficiently bright and limb-brightened jets. Therefore, nonthermal electron impacts on the near-horizon behavior of accretion flows around M87^* play a crucial role and should be examined closely.

Traditionally, the electron temperature is estimated using the gas temperature to calculate the emission using a single-fluid GRMHD simulation. One common parametric prescription is the so-called “R − β” model (Mościbrodzka et al. 2016), which estimates electron temperature using plasma-β and two parameters, R_l and R_h, respectively. Although this treatment is much more flexible in allowing a wide exploration of parameter space, it requires extensive parameter searching, and the physical processes related to the optimal parameters are difficult to interpret. In addition, the electron temperature is not only determined by ion temperature but also depends on the microscopic physical processes, such as heating, cooling, and the advection of the electrons (Event Horizon Telescope Collaboration 2021). Thus, the physically driven two-temperature model is needed to explore electron thermodynamics self-consistently. This approach further eliminates the dependence on hyperparameters R_l and R_h when calculating electron temperature. Several studies have shown that the self-consistent two-temperature model is well matched with the parameterized R − β model (Mizuno et al. 2021; Zhang et al. 2024; Mościbrodzka 2025). Apart from that, Cruz-Osorio et al. (2026) recently used the electron temperature obtained through the simulations of turbulent collisionless plasmas on a microscopic scale, providing a better depiction of the jet in morphology and width at 86 GHz than using the parameterized R − β model.

Ressler et al. (2015) introduced the two-temperature model to effectively evolve electron thermodynamics separately from ions by extending the equations in GRMHD simulations, while the energy-momentum and particle number conservation equations still assumed a single fluid. This approach has been applied to model M87^* (Ryan et al. 2018; Chael et al. 2019) and Sgr A^* (Ressler et al. 2017; Chael et al. 2018; Dexter et al. 2020; Yoon et al. 2020).

While our understanding of accretion flows has improved considerably in the last few decades, neglecting radiative cooling can lead to incomplete models or inconsistent predictions. For instance, radiative cooling may result in the electron temperature being lower than the ion temperature even in the high magnetization regions where the electron is heated efficiently (e.g., Mościbrodzka et al. 2011; Ryan et al. 2018; Chael et al. 2019). Moreover, radiative cooling affects the disk structure via the pressure balance inside the disk (e.g., Dibi et al. 2012; Yoon et al. 2020; Singh et al. 2025, 2026). Recently, Salas et al. (2025) found that the two-temperature model with radiative cooling better matches the historical observations in flux variability at 230 GHz for Sgr A^*. This is achieved by decreasing the total flux and its fluctuations. However, most of the aforementioned studies adopted thermal eDF, and radiative cooling impacts on black hole shadows and extended jets have not yet been fully investigated. This is the purpose of our study. In addition, the measured radiative efficiency of M87^* is relatively high for a hot accretion flow model (Event Horizon Telescope Collaboration 2021). These indicate that considering the two-temperature model with radiative cooling using nonthermal eDF may be crucial for accurately modeling M87^*’s black hole shadow and large-scale jet structure.

We organize this paper as follows: Sect. 1 provides a brief overview of the background. In Sect. 2, we describe the dissipation and cooling processes included in the GRMHD simulations as well as the nonthermal eDF used in the GRRT calculation. The results are presented and discussed in detail in Sect. 3. In Sect. 4, we discuss the implications of our findings, address the limitations of the study, and outline possible directions for future work.

Throughout this paper, we adopt units in which the speed of light is c = 1 and the gravitational constant is G = 1. We absorb a factor of $\sqrt{4\pi }$ into the definition of the magnetic field four-vector, b^μ.

2. Numerical setup

2.1. General relativistic magnetohydrodynamic simulations

We performed GRMHD simulations of magnetized accretion flows onto a rotating black hole, considering Coulomb interaction and radiative cooling in a two-temperature framework following Dihingia et al. (2023). Here, we considered radiative cooling as a source term without solving the radiation fields explicitly. We utilized the BHAC code (Porth et al. 2017; Olivares et al. 2019) for this study. The metric adopted in the simulation consists of spherical modified Kerr-Schild (MKS) coordinates. The torus is initialized by a Fishbone-Moncrief hydrodynamic equilibrium solution (Fishbone & Moncrief 1976) with r_in = 20 r_g and r_max = 40 r_g, where r_g = GM_BH/c² and M_BH is the black hole mass. An ideal-gas equation of state with a constant relativistic adiabatic index of γ = 4/3 was used (Rezzolla & Zanotti 2013). We note that a value close to 5/3 is also a considerable choice (Chael 2025; Gammie 2025), and it is also possible to update the adiabatic index self-consistently (Sądowski et al. 2017; Salas et al. 2025). We put a weak single magnetic field loop in this equilibrium torus defined by the vector potential with only one nonzero component, A_ϕ ∝ max(q − 0.2, 0), where

$$\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}$$(1)

Here, ρ is the fluid rest-mass density, and ρ_max is maximum density in the torus. This field configuration supplies enough magnetic flux onto the black hole to reach the MAD state (e.g., Narayan et al. 2003; Tchekhovskoy et al. 2011). To excite the MRI inside the torus, a 4 percent random perturbation was applied to the gas pressure within the torus.

The electron temperature time evolution in the two-temperature GRMHD simulations is based on solving the electron entropy equation (Ressler et al. 2015; Mizuno et al. 2021). Electron heating is provided by grid-scale dissipation models (e.g., Ressler et al. 2015). The physical processes include turbulent heating, magnetic reconnection, shock waves, and Ohmic dissipation. Our tests used two heating prescriptions: turbulence (Kawazura et al. 2019) and magnetic reconnection (Rowan et al. 2017). Apart from that, the energy transfer from protons to electrons through Coulomb interaction (Spitzer 1965; Colpi et al. 1984) was also considered for radiative cooling cases. Electron energy loss was considered through radiative cooling processes, namely bremsstrahlung, cyclo-synchrotron radiation of the thermal electrons (Esin et al. 1996), and multiple inverse Compton scattering of the cyclo-synchrotron photons by the thermal electrons (Narayan & Yi 1995b). Note that the radiative cooling processes due to nonthermal electrons in the plasma were ignored for simplicity throughout the simulations. The detailed initial setup about Coulomb interaction and radiative cooling is described in Dihingia et al. (2023), Dihingia & Fendt (2025).

The outer radial boundary is located at r = 2 500 r_g. The inner radial position of the simulation domain is well inside the black hole horizon. The simulation domain was discretized using an effective grid resolution of 384 × 192 × 192 with three layers of static mesh refinement. In particular, Mościbrodzka (2025) demonstrate that the results are independent of the grid resolution for the turbulent heating model. Here, we considered black hole spin, a = 0.9375.

First, we ran the GRMHD simulations without Coulomb interactions and radiative cooling until t = 10 000 t_g, where t_g = GM_BH/c³. The simulations mostly reached a quasi-stationary MHD state at this time. Then, we applied radiative cooling and Coulomb interactions with different mass accretion rates. Here, based on the accretion rate estimated from Event Horizon Telescope Collaboration (2021) and the total observed flux S₂₃₀ = 0.5 Jy (Event Horizon Telescope Collaboration 2019a, 2024), we chose three different mass accretion rates normalized to the Eddington rate, ṁ = Ṁ_BH/Ṁ_Edd = 1 × 10⁻⁵, 5 × 10⁻⁶, and 1 × 10⁻⁶ at the horizon for the turbulent heating models as well as ṁ = 5 × 10⁻⁶ for the reconnection heating model as a reference. We carried out our simulations up to t = 15 000 t_g.

2.2. General relativistic radiative transfer calculations

To calculate black hole shadow images from GRMHD simulations, we used the GRRT code BHOSS (Younsi et al. 2012, 2020, 2023), which solves covariant radiative transfer equations via the ray-tracing method. Here, we considered electron synchrotron radiation as a mechanism for calculating the shadow. Additionally, we employed a hybrid distribution composed of thermal and variable κ components for our GRRT calculations. The thermal distribution used in this study follows the Maxwell-Jüttner distribution given by Eq. (2). The κ distribution can simultaneously represent both thermal eDFs and extended power-law characteristics by adjusting the κ parameter, given in Eq. (3). The Maxwell-Jüttner distribution is expressed as follows:

$$\begin{array}{c}\hfill \frac{d{n}_{\mathrm{e}}}{d{\gamma }_{\mathrm{e}}}=\frac{n_{\mathrm{e}}}{{\mathrm{\Theta }}_{\mathrm{e}}}\frac{{\gamma }_{\mathrm{e}}\sqrt{{\gamma }_{\mathrm{e}}^2-1}}{K_2(1/{\mathrm{\Theta }}_{\mathrm{e}})}exp(-\frac{{\gamma }_{\mathrm{e}}}{{\mathrm{\Theta }}_{\mathrm{e}}}),\end{array}$$(2)

where n_e is the electron number density, γ_e is the electron Lorentz factor, K₂ is the Bessel function of the second kind, and Θ_e is the dimensionless electron temperature (e.g., Mizuno et al. 2021).

The relativistic nonthermal κ distribution (Xiao 2006) is expressed as follows:

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

where N is the normalization factor (Pandya et al. 2016; Davelaar et al. 2018) and κ is related to the slope of the power-law distribution, s = κ − 1. When γ_e is large, particles satisfy $d{n}_{\mathrm{e}}/d{\gamma }_{\mathrm{e}}\propto {\gamma }_{\mathrm{e}}^{-s}$, and the nonthermal κ distribution approximates the power-law distribution. The parameter w specifies the width of the κ distribution. Considering the contribution of both thermal and magnetic energies to heating and accelerating electrons (Davelaar et al. 2019; Cruz-Osorio et al. 2022; Fromm et al. 2022), the specific expression of w is

$$\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}$$(4)

where r_inj is the injection radius, m_e is the electron mass, m_p is the proton mass, σ = b²/ρ is the magnetization, b² is the square of the four-magnetic field, ρ is the fluid rest-mass density, and ε is a tunable parameter for the region with a radius larger than r_inj. The energy is dominated by thermal energy with a limit of σ ≪ 1, while the magnetic energy contributes to highly magnetized regions. We set ε = 0.5 to account for the contribution of magnetic energy to the GRRT images and the spatial distribution of w. The jet stagnation surface is a potential injection site and defines the injection radius. The stagnation surface is located at u^r = 0, where the potential injection radius is usually between 5 and 10 r_g (e.g., Nakamura et al. 2018). We thus assumed r_inj = 10 r_g in this study.

The κ value is variable in different locations, and it is defined to be parametrically dependent on magnetization, σ, and plasma beta, β = p_g/p_m, where p_g is the fluid pressure and p_m = b²/2 is the magnetic pressure. For the PIC-CS model of Ball et al. (2018), the κ function can be expressed as

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

which was empirically obtained from particle-in-cell (PIC) simulations of magnetic reconnection in a Harris current sheet. For the PIC-TURB model of Meringolo et al. (2023), κ function is expressed as

$$\begin{array}{c}\hfill \kappa :=2.8+0.2{\sigma }^{-1/2}+1.6{\sigma }^{-6/10}tanh\left(2.25{\sigma }^{1/3}\beta \right),\end{array}$$(6)

which was obtained from PIC simulations of decaying plasma turbulence (see also Imbrogno et al. 2024, 2025 for results from PIC simulations of turbulent plasma). To more self-consistently account for the nonthermal electron distribution for the reconnection heating model, we used the PIC-CS model for κ and the PIC-TURB model for the turbulent heating case.

Following (Event Horizon Telescope Collaboration 2022), we assumed that the proportion of nonthermal electrons depends on σ and β. The emission coefficients c_ν (emissivity and absorptivity) combine the thermal (Leung et al. 2011) and κ coefficients (Pandya et al. 2016) through

$$\begin{array}{c}\hfill c_{\nu ,\mathrm{tot}}=(1-\eta )c_{\nu ,\mathrm{thermal}}+\eta c_{\nu ,\kappa },\end{array}$$(7)

where the nonthermal efficiency is expressed as

$$\begin{array}{c}\hfill \eta (ϵ,\beta ,\sigma )=ϵ[1-e^{-{\beta }^{-2}}][1-e^{-{(\sigma /{\sigma }_{min})}^2}].\end{array}$$(8)

Here, η → 0 on the disk, and η → ϵ in the jet. Because the emission at highly magnetized regions (σ > σ_cut = 1) were removed, the nonthermal electrons are mostly confined to the jet sheath. We set σ_min = 0.01 and ϵ = 0.5. In this study, we directly calculated the electron temperature from the two-temperature GRMHD simulations (Mizuno et al. 2021; Dihingia et al. 2023).

For the GRRT calculation, we modeled M87^* as a target source, with a mass of M_BH = 6.5 × 10⁹ M_⊙ at a distance of D = 16.8 Mpc (Event Horizon Telescope Collaboration 2019c). The field of view (FoV) was set to 760 μas in both directions, with a resolution of 1520 × 1520 pixels. The GRRT calculations were performed over the time range t ∈ [12 000 t_g, 15 000 t_g], with a 10 t_g cadence, at 230 GHz and an inclination angle of 163°. We varied the mass accretion rates (ṁ) in the radiative cooling GRMHD simulations as follows: 1 × 10⁻⁶ Ṁ_Edd, 5 × 10⁻⁶ Ṁ_Edd, and 1 × 10⁻⁵ Ṁ_Edd, where Ṁ_Edd is the Eddington accretion rate defined as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{Edd}}=\frac{L_{\mathrm{Edd}}}{\zeta c^2}=1.4\times {10}^{17}\frac{M_{\mathrm{BH}}}{M_{\odot }}{\mathrm{gs}}^{-1}.\end{array}$$(9)

Here, L_Edd = 4πGM_BHcm_p/σ_T is the Eddington luminosity, and σ_T is the Thomson cross section. By setting the efficiency to ζ = 1 and adopting the black hole mass for M87^* as M_BH = 6.5 × 10⁹ M_⊙, the Eddington accretion rate becomes Ṁ_Edd ≈ 9.1 × 10²⁶ gs⁻¹ ≈ 14 M_⊙ yr⁻¹.

To exclude regions with strong magnetization, a ceiling in magnetization was set at σ_cut = 1 for all models. To assess the impact of this choice on the shadow images, we also explore different values of σ_cut = 2, 5, 10, and 25 in Appendix A.

3. Results

3.1. Evolution of mass accretion rate and magnetic flux

To analyze the temporal evolution of the simulation models, Fig. 1 shows the normalized mass accretion rate, ṁ_BH ≡ Ṁ_BH/Ṁ_Edd, measured at the event horizon in Eddington units and the normalized magnetic flux, ${\varphi }_{\mathrm{BH}}\equiv {\mathrm{\Phi }}_{\mathrm{BH}}/\sqrt{{\dot{M}}_{\mathrm{BH}}}$, at the horizon (see e.g., Porth et al. 2019, for the definition of the magnetic flux). As shown in Fig. 1, the mass accretion rate profile and the normalized magnetic flux settle to steady states at t ≳ 10 000 t_g, with the small oscillations in time due to flux eruption events. The averaged values (t = 12 000 t_g − 15 000 t_g) for the non-cooling case (NC), cooling with time-averaged mass accretion rates ṁ = 1 × 10⁻⁶ (C_KA1e-6), 5 × 10⁻⁶ (C_KA5e-6), 5 × 10⁻⁶ (C_MR5e-6), and 1 × 10⁻⁵ (C_KA1e-5) are ϕ_BH∼ 12.5, 11.0, 10.4, 11.1, and 10.6, respectively. The decrease in ϕ_BH with radiative cooling is consistent with the result found in Singh et al. (2025). Our simulations indicate that the decrease in ϕ_BH results from the decrease in Φ_BH and the increase in $\sqrt{{\dot{M}}_{\mathrm{BH}}}$. In radiative cooling, the local magnetic field near the horizon becomes weaker due to the lower MRI growth rate. From t = 12 000 t_g to 15 000 t_g, the ratios of the standard deviation to the average value of the normalized magnetic flux – used to characterize variability – are 0.091, 0.102, 0.113, 0.110, and 0.132 for the NC, C_KA1e-6, C_KA5e-6, C_MR5e-6, and C_KA1e-5 simulations, respectively. The variability in the normalized magnetic flux becomes greater with radiative cooling. Our simulations considered the accretion flow at low accretion rates. Under these conditions, the radiative-cooling processes are not strong enough to collapse the geometrically thick torus into a thin-disc structure. Consequently, the temporal evolution and general trends of the mass accretion rate and magnetic flux rate are mostly similar.

Fig. 1. Top: Accretion rates measured at the event horizon. Bottom: Normalized magnetic flux at the horizon. The different colored curves correspond to the various electron heating prescriptions, radiative cooling, and time-averaged accretion rates. The black curve corresponds to the model without cooling (ṁ = 1 × 10−6). The blue, green, and red curves depict the turbulent heating model with cooling for ṁ = 1 × 10−6, 5 × 10−6, and 1 × 10−5, respectively. The magenta curve represents reconnection heating with cooling (ṁ = 5 × 10−6).

3.2. Density distribution

Figure 2 shows the time- and azimuthally averaged density distribution over the interval t = 12 000 t_g to 15 000 t_g for the electron heating prescriptions and radiative cooling conditions, along with the corresponding time-averaged normalized accretion rates. Specifically, panels (a), (b), (c), and (e) compare the density distribution: (a) without radiative cooling, (b), (c), and (e) with turbulent heating and radiative cooling at normalized mass accretion rates of ṁ = 1 × 10⁻⁶, 5 × 10⁻⁶, and 1 × 10⁻⁵, respectively. Panels (c) and (d) compare the density profile with radiative cooling at ṁ = 5 × 10⁻⁶ under turbulent heating (c) and reconnection heating (d). The dashed white and solid black curves represent σ = 0.1 and 1, respectively. All cases generally show a relatively dense disk around the equatorial plane, a low-density off-equatorial region, and a sparse funnel region around the bipolar directions. The disk structure is related to the mass accretion rate when radiative cooling is considered. As shown in panels (a), (b), (c), and (e), a thinner disk forms with a higher accretion rate, due to the increase in radiative cooling efficiency (Singh et al. 2025; Dihingia et al. 2025). Compared with panels (c) and (d), under the same mass accretion rate ṁ = 5 × 10⁻⁶, there is no significant difference in the density profile caused by the electron heating prescriptions (i.e., turbulent heating and reconnection heating). The low-density region contributes to the disk winds, and the sparse funnel region contributes to the relativistic Poynting-dominated jet (Vourellis et al. 2019; Dihingia et al. 2021). Due to the higher density region around the equatorial plane, the efficiency of Coulomb interactions and electron radiative cooling increase (Dihingia et al. 2023). Hence, we expect the electron temperature near the inner disk to be affected (see Sect. 3.3 for more details).

Fig. 2. Logarithmic density distribution averaged in time and azimuth over the interval t = 12 000 tg to 15 000 tg. From left to right: Without cooling (a); turbulent heating with cooling at ṁ = 1 × 10−6 (b), ṁ = 5 × 10−6 (c), and ṁ = 1 × 10−5 (e); and reconnection heating with cooling at ṁ = 5 × 10−6 (d). The dashed white and solid black curves represent magnetization for σ = 0.1 and 1, respectively.

3.3. Temperature distribution

The electron temperature is an important quantity for modeling the electromagnetic radiation of an accreting black hole. In this section, we discuss the dependence of the dimensionless electron temperature on radiative cooling and mass accretion rates.

Figure 3 shows the time- and azimuthally averaged distribution of the logarithmic dimensionless electron temperature Θ_e = k_BT_e/m_ec² (e.g., Mizuno et al. 2021) without (panels (a) and (f)) and with (panels (b)–(e)) radiative cooling at different mass accretion rates. The differences are shown in a linear scale in panels (g)–(j) and are compared with the corresponding non-cooling case. The red shaded region indicates the increase in density compared to the non-cooling case, while the blue shaded region highlights the decrease in density compared to the non-cooling case. A smooth transition appears from red to blue via white. The solid black curves represent σ = 1, and the dashed sky-blue thin to thick curves represent Θ_e= 10, 32, and 100, respectively. Figure 4 shows the angular distribution of the time- and azimuthally averaged dimensionless electron temperatures at the given radii. The different colored curves correspond to the non-cooling and radiative cooling models at different accretion rates. The black curves represents the turbulent heating model without cooling; and the red, blue, and magenta curves represents cooling at ṁ = 1 × 10⁻⁶, 5 × 10⁻⁶, and 1 × 10⁻⁵, respectively. The dash-dotted lines in the same color represent σ = 1 for each case. The vertical dashed lines in black on both sides correspond to the disk region in the non-cooling case. The lines for other cases are omitted for simplicity, as they are located at a similar position. Here, disk angular thickness at a given radius r is (h/r)_r = [∬_θ,  φ(θ−π/2)²ρ dA_θφ /∬_θ,  φρ dA_θφ]^1/2, where $\mathrm{d}{A}_{\theta \phi }=\sqrt{-g}\mathrm{d}\theta \mathrm{d}\phi$ is an area element in θ − φ plane, and g is the metric determinant. The integrals are over all θ,  φ on a sphere of radius r (see e.g., Tchekhovskoy & McKinney 2012, for more details about the disk angular thickness). The dotted green lines mark the boundary where the images are decomposed (see Sect. 3.4 for more details).

Fig. 3. Panels (a)–(f): Logarithm of the dimensionless electron temperature averaged in time and azimuth over the interval t = 12 000 tg to 15 000 tg. Panels (g)–(j): Differences in the linear scale, calculated by subtracting the dimensionless electron temperature in the corresponding non-cooling case. The solid black curves represent σ = 1. The dashed sky-blue thin to thick curves represent Θe= 10, 32, and 100, respectively.

From Fig. 3, we see that the electron temperature at the inner disk falls sharply for all cooling cases around the equatorial plane. This dramatic drop can be confirmed quantitatively from Fig. 4 at least within 10 r_g. In addition, in the slightly farther outer disk region (∼20 r_g), the electron temperature also drops, with an increase in the mass accretion rate. Interestingly, the temperature profiles look similar in both the C_KA5e-6 and C_MR5e-6 cases, even though the underlying heating functions are different. This may be because they evolve from the same fluid, so there is no significant difference induced by the MRI, which was excited by random perturbations at the beginning of the simulation. In addition, temperature differences remain in the disk and jets (see Appendix C for more details). In Fig. 4, the disk region for the non-cooling case is outlined by dashed black lines. At a smaller radius, we see a smaller angular thickness. This is because, at a smaller radius, the disk near the equatorial plane is compressed vertically by magnetic fields at the funnel region. The temperature in the jet regions (outside the midplane) slightly decreases when radiative cooling is considered, compared to the non-cooling case.

Fig. 4. Angular distribution of time- and azimuthally averaged dimensionless electron temperatures at the given radii, plotted on a logarithmic scale. From top to bottom: Radial increase from 7 rg to 20 rg. The different colored curves correspond to the non-cooling and radiative cooling models at the normalized mass accretion rates. The black curve represents the turbulent heating model without cooling. The red, blue, and magenta curves correspond to the models with cooling at ṁ = 1 × 10−6, 5 × 10−6, and 1 × 10−5, respectively. The dash-dotted lines in the same color represent σ = 1 for each case. The vertical dashed lines in black on both sides correspond to the disk region of the non-cooling case (see Sect. 3.3 for more details). The dotted lines in green mark the boundary where the images are decomposed (see Sect. 3.4 for more details).

In summary, when we consider radiative cooling, the electron temperature in the inner disk around the equatorial plane significantly decreases. Furthermore, the temperature slightly decreases farther from the equatorial plane. This may profoundly impact the shadow images of the black holes, which we discuss in Sect. 3.4.

3.4. Image decomposition

Decomposed images estimate the fraction of the total emission contributed by different regions. This allows us to better understand the emission sources and their morphology in the image. By analyzing the decomposed images, we can gain insights into the underlying physical processes that are either responsible for or impact the observed image. Following previous studies (Event Horizon Telescope Collaboration 2019b; Zhang et al. 2024), we divided the entire region into three parts: the midplane, the nearside jet, and the farside jet. Specifically, the polar angles per region range as follows: the farside jet spans from 0° to 57.3° (1.0 rad); the midplane lies between 57.3° (1.0 rad) and 122.7° (2.14 rad); and the nearside jet extends from 122.7° (2.14 rad) to 180°. Emissivity in our GRRT calculations are set at zero in regions outside the specific polar angles mentioned above.

Figure 5 shows the time-averaged decomposed images without radiative cooling, along with their fractional emission contribution to the total image over the interval t ∈ [12 000 t_g,  15 000 t_g], at 230 GHz. The panels represent different heating prescriptions and various mass accretion rates with the eDF modeled as a hybrid of thermal and variable κ components. The percentage marked on the bottom right of each image represents the fraction of the total emission contributed by this region relative to the whole image. Morphologically, we can confirm that the extended structure is seen in cases with different heating prescriptions and mass accretion rates from the nearside jet. In the turbulent heating scenario, the emission fraction from the nearside jet increases from 4.7% to 10.3% as the mass accretion rates increase. Quantitatively, the increase in mass accretion rate enhances the optical thickness in the midplane (see Sect. 3.5 for more details). Consequently, the contribution from optically thin, extended (nearside) jet emission increases with the accretion rate in the non-cooling case. Comparing the non-cooling cases for turbulent and reconnection heating at the same accretion rate Ṁ_BH/Ṁ_Edd = 5 × 10⁻⁶, the decomposed images are similar. This is because the temperatures are comparable in the absence of cooling (see Appendix C for more details).

Fig. 5. Time-averaged GRRT decomposed images from MAD simulations from t = 12 000 tg to 15 000 tg, assuming a black hole spin of a = 0.9375, observed at 230 GHz with an inclination angle of 163°. From top to bottom: Accretion rates of ṀBH/ṀEdd = 1 × 10−6, 5 × 10−6, 5 × 10−6, and 1 × 10−5, respectively. The electron heating prescriptions are used for turbulent heating, turbulent heating, reconnection heating, and turbulent heating, respectively. From left to right: Emissions from every region, including the full, midplane, nearside jet, and farside jet regions. The eDF is modeled as a hybrid of thermal and variable κ components, with ε = 0.5 in κ width, w. Cooling is not included.

Similar to Fig. 5, Fig. 6 shows the decomposed images at different mass accretion rates, but including radiative cooling. Taking the non-cooling images as the reference, the effects of radiative cooling on the images can be studied. First and foremost, compared to the non-cooling case, including radiative cooling leads to a sharp decrease in the total flux for every region at the same mass accretion rate due to the reduced electron temperature. Furthermore, the images with radiative cooling clearly exhibit more extended jet structures, and the nearside jet emission contribution increases. Thus, radiative cooling leads to a dim disk, more extended and brighter jets, and reduced total flux. These are consistent with radiative cooling effects on electron temperature, as illustrated in Sect. 3.3: radiative cooling results in a cooler disk as well as a slight decrease in electron temperature in the jet sheath. Notably, under turbulent heating, as the mass accretion rates increase, the fraction of emission contributed by the nearside jet first increases from 11.8% to 16.3% and then decreases to 13.6%. This is because the viewing angle for the nearside jet is smaller for the scenario at a mass accretion rate of Ṁ_BH/Ṁ_Edd = 1 × 10⁻⁵ when we consider σ_cut = 1 (e.g., see the dash-dotted magenta lines at 15 r_g and at 20 r_g in Fig. 4). They are located further to the left than in the other cases, indicating that larger regions were excluded. Meanwhile, the viewing angle of the farside jet is similar across cases with different mass accretion rates. Hence, the farside jet contributes relatively more emission as disk becomes cooler with increasing mass accretion rate. Furthermore, as shown in Fig. A.1, when we consider σ_cut = 2 for the case with a mass accretion rate of Ṁ_BH/Ṁ_Edd = 1 × 10⁻⁵ (for which the viewing angle of the nearside jet is similar to that in the other cases with σ_cut = 1), the fraction of emission from the nearside jet increases to the level seen in the case with Ṁ_BH/Ṁ_Edd = 5 × 10⁻⁶. Nonetheless, the difference is very small. This may be because the Coulomb interactions dominate over radiative cooling in the inner dense disk (∼7 r_g) for the case with Ṁ_BH/Ṁ_Edd = 1 × 10⁻⁵. It leads to higher electron temperatures in the inner dense disk (see Fig. 4 for more details) and more emission from the midplane. Therefore, the increase in the fraction of emission from the nearside jet slows down. When comparing the radiative cooling case with turbulent and reconnection heating at the same mass accretion rate, Ṁ_BH/Ṁ_Edd = 5 × 10⁻⁶, we find that case C_KA5e-6 exhibits more emissions from the midplane. By contrast, case C_MR5e-6 exhibits more emissions from the nearside and farside jets. These differences using various electron heating prescriptions are significantly not shown in the non-cooling cases (see Appendix C for further details).

Fig. 6. As in Fig. 5 but with radiative cooling included.

3.5. Spectral energy distribution

Even though none of the models is precisely tuned to fit the M87^* data, we show the SEDs of different regions at various accretion rates without and with cooling in Fig. 7 to investigate how radiative cooling affects flux variation with frequencies at different mass accretion rates. The solid curves represent the average flux from t = 12 000 t_g to t = 15 000 t_g, the points correspond to the peaks, and the shaded regions denote the standard deviation resulting from the time variation relative to the average values. The black, red, blue, and green curves represent emission from the whole region, the midplane, the nearside jet, and the farside jet, respectively. To compare the effects of radiative cooling, non-cooling cases with various mass accretion rates are included in the left panel of Fig. 7.

Fig. 7. SED curves of different regions under turbulent heating at various accretion rates, without cooling (left) and with radiative cooling (right). From top to bottom: Time-averaged mass accretion rates of ṀBH/ṀEdd = 1 × 10−6, 5 × 10−6, and 1 × 10−5. All curves adopt the hybrid thermal and variable κ eDF. The solid curves represent average values. The points correspond to the peaks. The shaded regions denote the standard deviation relative to the average values. The vertical dash-dotted lines correspond to 86 GHz, 230 GHz, 340 GHz, and 136 THz.

First, it is useful to understand SED behavior as a whole. As shown in Fig. 7, SED curve initially increases with increasing frequency, followed by a decrease after reaching a turnover frequency. The turnover frequencies for the non-cooling cases with Ṁ_BH/Ṁ_Edd = 1 × 10⁻⁶, 5 × 10⁻⁶, and 1 × 10⁻⁵ are 78 GHz, 230 GHz, and 397 GHz, respectively. For the radiative cooling cases at the same accretion rates, the turnover frequencies are 43 GHz, 78 GHz, and 108 GHz, respectively. The peak frequency is related to the electron temperature, the magnetic field strength, and the optical depth (e.g., Zdziarski et al. 1998). Notably, radiative cooling reduces the temperature. The optical depth is related to the mass accretion rate and electron distribution. As a result, at the same accretion rate, the peak (turnover frequency) in the cooling case shifts to a lower frequency compared to the non-cooling case. Regardless of whether the case is non-cooling or cooling, the SED with a higher accretion rate reaches its peak at a higher frequency. Thus, the peak position is sensitive to accretion rates and radiative cooling.

Second, based on the position of the 230 GHz flux relative to the peak (turnover frequency) in the SED, synchrotron radiation is more self-absorbed at higher accretion rates in the non-cooling cases. This means that the emission at 230 GHz becomes optically thick at higher accretion rates. In the radiative cooling cases, however, the emission at 230 GHz remains in the optically thin regime, and the turnover frequencies maintain relatively similar positions. For the non-cooling cases, this self-absorption mainly occurs in the midplane, while the jet emission is exponentially cut off at 230 GHz. This explains why more emission originates from the jets as the accretion rate increases, as shown in Sect. 3.4.

Third, increasing the accretion rates leads to higher fluxes at higher frequencies. For a given accretion rate, radiative cooling reduces the flux at high frequencies compared to the non-cooling case. At low frequencies and at the same accretion rate, the non-cooling and radiative cooling cases behave similarly.

3.6. Time variability

Figure 8 shows the 230 GHz light curves for the models with turbulent heating at various mass accretion rates both without and with radiative cooling. Generally, the average flux increases with the accretion rate. For a given accretion rate, the flux in the non-cooling case is higher than in the radiative cooling case. Additionally, the variability in both the non-cooling and radiative cooling cases decreases as the accretion rate increases.

Fig. 8. Light curves of the 230 GHz flux at a 163° inclination angle and spin a = 0.9375. The curves are plotted for accretion rates of ṀBH/ṀEdd = 1 × 10−6 (top), 5 × 10−6 (middle), and 1 × 10−5 (bottom). The solid curves represent the non-cooling cases, and the dashed curves represent the radiative cooling cases. All curves assume turbulent heating and adopt the eDF modeled as a hybrid of thermal and variable κ components.

To quantitatively compare the models at different accretion rates, without and with cooling, Fig. 9 shows the time-averaged flux at 230 GHz and its standard deviation. We fixed the accretion rates rather than the flux, so the average value is different from case to case. Overall, the average flux increases with the accretion rate for both non-cooling and cooling cases. At a given accretion rate, the non-cooling case exhibits higher flux than the radiative cooling case.

Fig. 9. Total flux variation in turbulent heating in the non-cooling (dots) and radiative cooling (squares) cases at 230 GHz at a 163° inclination angle and spin a = 0.9375. The different colors correspond to the different accretion rates: ṀBH/ṀEdd = 1 × 10−6 (black), 5 × 10−6 (red), and 1 × 10−5 (blue). The labels on the x-axis denote the emissions from each region: the full region is depicted first (Total), followed by the midplane (Mid), the nearside jet (Near), and the farside jet (Far). The colored dots and squares represent the time-averaged values, with the error bars indicating the standard deviation about the mean.

To understand the time variability for each case, we list the ratio of standard deviation relative to the average value in Table 1. For the non-cooling cases with emission from the entire region (total), the results show that the modulation index (the standard deviation divided by the average value) decreases from 0.21 to 0.15 as the mass accretion rate increases from Ṁ_BH/Ṁ_Edd = 1 × 10⁻⁶ to 1 × 10⁻⁵. This trend is the same for the midplane, nearside, and farside regions. It indicates that for non-cooling cases, the time variability of all regions decreases with increasing accretion rate. For radiative cooling cases, we see a similar dependence of time variability on the accretion rate for the total region and the midplane. However, no clear dependence of time variability on the accretion rate is evident for the nearside and farside jet regions in the cooling cases.

To determine the origin of the time variability in the total region, we also list the ratio of the standard deviation relative to the average value of the total case in Table 1. Our results show that, the time variability of the total region is primarily driven by the midplane. Furthermore, in the cooling cases, the jet regions also partially contribute to the variability. The proportion of emission from the jet regions is higher in the cooling than in the non-cooling cases (see Figs. 5 and 6 for these proportions).

4. Summary and discussion

In our previous study (Zhang et al. 2024), we ignored energy exchange impacts from Coulomb coupling and radiative cooling on black hole shadow and extended jet images. To address this issue, we adopted the electron heating prescriptions of Mizuno et al. (2021) and added both Coulomb interactions and radiative cooling, following Dihingia et al. (2023), in our GRMHD simulations. In our GRRT calculations, we adopted a hybrid thermal and variable nonthermal κ eDF, (e.g., Event Horizon Telescope Collaboration 2022; Cruz-Osorio et al. 2022; Davelaar et al. 2023; Zhang et al. 2024), with ε = 0.5. We also considered the non-cooling cases as a reference for exploring radiative cooling impacts on density distribution, electron temperature distribution, shadow images, jet morphology, SEDs, and total flux variation at 230 GHz at an inclination angle of i = 163°. Below we list our conclusions point-by-point.

1. We find that Coulomb coupling and radiative cooling significantly modify the electron temperature distribution, leading to a cooler accretion disk, a slightly cooler jet sheath, and changes in black hole shadow morphology.

2. Our results show that radiative cooling produces a dimmer ring, more extended and brighter jets, and a reduced total flux emission.

3. Our results show that radiative cooling mediates heating effects. At an accretion rate of Ṁ_BH/Ṁ_Edd = 5 × 10⁻⁶, the inclusion of radiative cooling results in a dimmer disk but brighter jets under reconnection heating compared to the turbulent heating case. Notably, the total flux remains higher under reconnection heating than under turbulent heating in both the non-cooling and radiative cooling scenarios.

4. We find that cooling reduces synchrotron self-absorption, shifting the turnover frequency to lower values at the same accretion rate.

5. Our results show that time variability primarily originates from the midplane in both the cooling and non-cooling cases, and higher accretion rates reduce time variability.

6. We find that although our simulations are scaled to M87^*, our results are generally applicable to other low-luminosity active galactic nuclei (LLAGNs) in the MAD state, where radiative cooling universally leads to a cooler disk, brighter jets, reduced total flux, and smoother light curves, independent of black hole mass and distance.

We note that, in this study, we only considered the σ_cut = 1, zeroing out 230 GHz emission from the regions with a higher magnetization. Several studies have examined the impact of different σ_cut assumptions. For example, the results are sensitive to the choice of σ_cut at frequencies ν > 230 GHz (Chael et al. 2019). Including contributions from regions with stronger magnetization enhances the jet (Chael et al. 2019; Zhang et al. 2024). Additionally, a lower σ_cut value requires a higher accretion rate to satisfy the observed flux density, resulting in more efficient radiative cooling and a greater Faraday rotation depth (Chael 2025). Hence, with radiative cooling, different treatments of σ_cut yield varying jet structures and polarization fractions, as these depend on the accretion rate (follow Appendix A for more details).

In addition, we used radiative cooling and heating to avoid overestimating electron temperature. We find that radiative cooling reduces the flux at high frequencies compared to the non-cooling case. However, future GRMHD simulations should incorporate cooling functions more self-consistently to account for different nonthermal distributions. This includes adding cooling terms for nonthermal electrons due to bremsstrahlung, synchrotron radiation, and inverse Compton scattering. These nonthermal electrons would radiate more efficiently than thermal electrons and might possibly cool the plasma further. Moreover, it would be interesting to find the proper mass accretion rates (likely between 1 × 10⁻⁶ Ṁ_Edd and 5 × 10⁻⁶ Ṁ_Edd) for different models targeting the same flux, and to assess which model best fits the observation. For instance, Cruz-Osorio et al. (2022) successfully reproduce the broadband spectrum of M87 from radio to near-infrared bands without considering radiative cooling. To reproduce the broadband spectrum while accounting for radiative cooling, it is necessary to include more nonthermal electrons. This ensures sufficient emission in simulations in the near-infrared and optical bands. This raises important questions about which electron heating prescriptions can accelerate the electrons and where this electron energization occurs.

Lastly, Zhang et al. (2024) show that nonthermal electrons result in more extended and brighter jets. Furthermore, as shown in this paper, radiative cooling still produces a faint large-scale jet emission on horizon scales, which appears brighter in GRRT images at 230 GHz. The total intensity of this emission varies with different treatments of σ_cut. The presence of such horizon-scale jet emission, as well as the explicit consideration of σ_cut, should be rigorously evaluated against observational data. Although currently below the dynamic range of EHT observations, both the horizon-scale shadow and the extended jet images may be resolved simultaneously at 230 GHz by the next-generation arrays, such as the ngEHT (Johnson et al. 2023; Ricarte et al. 2023; Ayzenberg et al. 2025). These arrays would have sufficient dynamic range to detect the features reported in this study.

Acknowledgments

This research is supported by the National Key R&D Program of China (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). MZ is supported by the Doctoral Student Program of the Young S & T Talents Cultivation Project, CAST, and by the T.D. Lee scholarship. IKD acknowledges the TDLI postdoctoral fellowship for financial support. ACO acknowledges to DGAPA-UNAM (grant IN110522), the Ciencia Básica y de Frontera 2023–2024 program of SECIHTI México (projects CBF2023-2024-1102 and 257435), and the European Horizon Europe Staff Exchange (SE) programme HORIZON-MSCA2021-SE-01 under Grant No. NewFunFiCO-101086251. The simulations were performed on the Astro cluster at Tsung-Dao Lee Institute, Pi 2.0, and the Siyuan-1 cluster in the Center for High Performance Computing at Shanghai Jiao Tong University. This work has made use of NASA’s Astrophysics Data System (ADS).

References

Appendix A: Exclusion of magnetized region

Due to the density, pressure, and internal energy in simulations may reach the floor value in highly magnetized regions (where σ ≫ 1), we chose a conservative cut-off value σ_cut = 1 in the previous results in this study to avoid the potential numerical error in highly magnetized regions.

Here, we investigate the effects of various σ thresholds on flux and jet structure for a MAD state considering radiative cooling, under a mass accretion rate Ṁ_BH/Ṁ_Edd = 1 × 10⁻⁵. Figure A.1 shows the time-averaged (t = 12 000 t_g  −  15 000 t_g) decomposed GRRT images with various σ_cut at 230 GHz and an inclination angle of 163°, with hybrid thermal and variable κ eDF. Combing with decomposed GRRT images with σ_cut = 1 shown in the last row in Fig. 6 and the images in Fig. A.1, the results show that with σ_cut increasing to 25, the ratio of extended jet emission to total emission increases from 16.4 % to 25.3 %, compared to the case of σ_cut = 1. Meanwhile, the total flux increases from 12.5 Jy to 14.7 Jy. It is easy to confirm that the increase in flux mainly results from the brighter nearside jets, and the jet structure remains similar. Interestingly, the proportion of farside jet flux remains similar under the increase of σ_cut. That means the flux of the farside jet, though increasing, is moving away from us and is therefore beyond our reach. The similar diminishing rate of increase on flux of farside jet with the increase of σ_cut can also be found in Fig. 10 in Zhang et al. (2024). Therefore, consider radiative cooling, under a mass accretion rate Ṁ_BH/Ṁ_Edd = 1 × 10⁻⁵, there is no significant dependence of jet structure on σ_cut values, but the flux of nearside jet depends on σ_cut values. However, if one wants to find a proper mass accretion rate satisfying the total flux obtained from the observation, the calculated mass accretion rate and, thereafter, the radiative efficiency will be different for various choices of σ_cut value, which results in different theoretical predictions on jet emission, Faraday depth, and the polarization fraction.

Fig. A.1. Same as Fig. 5 but using σcut = 2, 5, 10, and 25, respectively, with radiative cooling, in turbulent heating, and the mass accretion rate ṀBH/ṀEdd = 1 × 10−5.

Appendix B: Difference maps of density

To better visualize the differences in density in different heating prescriptions, the radiative cooling, and the mass accretion rates, we plotted the pixel-by-pixel difference maps as shown in Fig. B.1. We plotted the time- (t = 12 000 t_g  −  15 000 t_g) and azimuthally averaged distribution of logarithmic density without radiative cooling in panel (a), and the differences of density in linear scale in panels (b) − (e) under various mass accretion rate and electron heating prescriptions with radiative cooling, compared with panel (a). The region shaded in red indicates the density exceeds that of the non-cooling case, while the region shaded in blue highlights the density is lower than that of the non-cooling case. There is a smooth transition from red to blue through white. The solid black curves represent σ = 1. As shown in panels (b), (c), and (e) in Fig. B.1, the density on the equatorial plane generally decreases, and a thinner disk is formed with a higher accretion rate, which is consistent with the results in Sec. 3.2. Compared with panels (c) and (d), the density is higher in the former case.

Fig. B.1. Same as Fig. 2 for panel (a), but panels (b) – (e) highlight the differences in linear scale by subtracting the density depicted in panel (a). The solid black curves represent σ = 1.

Appendix C: Angular distribution of electron temperature and spectrum in reconnection heating

To elucidate why there are more emissions from midplane in simulation C_KA5e-6 than C_MR5e-6 and more emissions from jets in latter case, and why the decomposed images become similar in both turbulent and reconnection cases without cooling as discussed in Sec. 3.4, we present the angular distribution of the time- and azimuthally averaged dimensionless electron temperature for reconnection heating at specified radii, and the spectral energy distribution curves of different regions in reconnection heating.

The temperature distribution curves under turbulent heating are also included in Fig. C.1 as a reference to compare the impacts of different electron heating prescriptions and radiative cooling. The curves in various colors correspond to different cases: the turbulent heating model without cooling (black), with cooling of Ṁ_BH/Ṁ_Edd = 5 × 10⁻⁶ (blue), and the reconnetion heating model without cooling (red), with cooling of Ṁ_BH/Ṁ_Edd = 5 × 10⁻⁶ (magenta).

Fig. C.1. Same as Fig. 4 but the curves in various colors correspond to different cases: the turbulent heating model without cooling (black), with cooling of ṀBH/ṀEdd = 5 × 10−6 (blue), and the reconnetion heating model without cooling (red), with cooling of ṀBH/ṀEdd = 5 × 10−6 (magenta).

As shown in Fig. C.1, the temperature without cooling is similar in both the turbulent and reconnection heating. This is why under the same accretion rate Ṁ_BH/Ṁ_Edd = 5 × 10⁻⁶ without cooling, the total flux, emission contribution from each divided region, and the jet structure are similar in both cases.

To understand why there are more emissions from midplane for simulation C_KA5e-6 than case C_MR5e-6 and more emissions from jets in the latter, we need to consider the differences in density profile, temperature distribution, and the emissions from nonthermal electrons. Even though Fig. C.1 indicates that the disk temperature in the C_MR5e-6 case is higher than the C_KA5e-6 case, the density is higher in the latter case (see panels (d) and (c) in Fig. B.1), which results in there being more emissions from the midplane in the case C_KA5e-6. Besides, the spectral energy distribution curves under turbulent heating are also included in Fig. C.2 as a reference to understand the impacts of nonthermal electrons in different heating prescriptions, at 230 GHz and on a broader range of frequencies. Compared with simulation C_MR5e-6, Fig. C.2 shows that the spectral curves of jets are flatter in case C_KA5e-6 and the flux at higher frequency is also higher. Furthermore, the emissions from jets at 230 GHz are lower in case C_KA5e-6. One reason is that the temperature of jets in case C_MR5e-6 is higher than that in simulation C_KA5e-6, as shown in Fig. C.1. The other reason is that the PIC-TURB model of Meringolo et al. (2023) shows smaller values of κ than the PIC-CS model of Ball et al. (2018) in the jet region (Cruz-Osorio et al. 2026), which indicates a broader distribution of nonthermal electrons. As a result, as shown in Fig. 6 and Fig. C.2, there are more emissions from jets at higher frequencies rather than 230 GHz in case C_KA5e-6, and more emissions are coming from jets in simulation C_MR5e-6 at 230 GHz. Furthermore, we plotted the time- (t = 12 000 t_g − 15 000 t_g) and azimuthally averaged distribution of κ in model PIC-TUBR for simulations NC and C_KA5e-6 in Fig. C.3. As shown in Fig. C.3, the κ values are smaller in the jet region in the cooling case C_KA5e-6 than those in the non-cooling case NC. Hence, the nonthermal impacts in the PIC-TURB model on the jets at 230 GHz and higher frequencies are more significant if the radiative cooling is considered (i.e., there are more emissions from jets in case C_MR5e-6 than C_KA5e-6 at 230 GHz as shown in Fig. 6, due to the nonthermal emission contributing more at higher frequencies in the PIC-TURB model with radiative cooling. And as shown in Fig. 5, the jet emission is similar in both the PIC-TURB model and PIC-CS model without cooling).

Fig. C.2. Same as Fig. 7 for simulation C_KA5e-6 but the case C_MR5e-6 is also added.

Fig. C.3. Time- and azimuthally averaged distributions of κ for the PIC-TURB obtained from decaying plasma turbulence [Eq. (6)] in cases NC and C_KA5e-6. The solid black and dashed white curves represent σ = 1 and the Bernoulli parameter, −hut = 1.02, respectively.

All Tables

All Figures

Fig. 1. Top: Accretion rates measured at the event horizon. Bottom: Normalized magnetic flux at the horizon. The different colored curves correspond to the various electron heating prescriptions, radiative cooling, and time-averaged accretion rates. The black curve corresponds to the model without cooling (ṁ = 1 × 10−6). The blue, green, and red curves depict the turbulent heating model with cooling for ṁ = 1 × 10−6, 5 × 10−6, and 1 × 10−5, respectively. The magenta curve represents reconnection heating with cooling (ṁ = 5 × 10−6).

In the text

Fig. 2. Logarithmic density distribution averaged in time and azimuth over the interval t = 12 000 tg to 15 000 tg. From left to right: Without cooling (a); turbulent heating with cooling at ṁ = 1 × 10−6 (b), ṁ = 5 × 10−6 (c), and ṁ = 1 × 10−5 (e); and reconnection heating with cooling at ṁ = 5 × 10−6 (d). The dashed white and solid black curves represent magnetization for σ = 0.1 and 1, respectively.

In the text

Fig. 3. Panels (a)–(f): Logarithm of the dimensionless electron temperature averaged in time and azimuth over the interval t = 12 000 tg to 15 000 tg. Panels (g)–(j): Differences in the linear scale, calculated by subtracting the dimensionless electron temperature in the corresponding non-cooling case. The solid black curves represent σ = 1. The dashed sky-blue thin to thick curves represent Θe= 10, 32, and 100, respectively.

In the text

Fig. 4. Angular distribution of time- and azimuthally averaged dimensionless electron temperatures at the given radii, plotted on a logarithmic scale. From top to bottom: Radial increase from 7 rg to 20 rg. The different colored curves correspond to the non-cooling and radiative cooling models at the normalized mass accretion rates. The black curve represents the turbulent heating model without cooling. The red, blue, and magenta curves correspond to the models with cooling at ṁ = 1 × 10−6, 5 × 10−6, and 1 × 10−5, respectively. The dash-dotted lines in the same color represent σ = 1 for each case. The vertical dashed lines in black on both sides correspond to the disk region of the non-cooling case (see Sect. 3.3 for more details). The dotted lines in green mark the boundary where the images are decomposed (see Sect. 3.4 for more details).

In the text

Fig. 5. Time-averaged GRRT decomposed images from MAD simulations from t = 12 000 tg to 15 000 tg, assuming a black hole spin of a = 0.9375, observed at 230 GHz with an inclination angle of 163°. From top to bottom: Accretion rates of ṀBH/ṀEdd = 1 × 10−6, 5 × 10−6, 5 × 10−6, and 1 × 10−5, respectively. The electron heating prescriptions are used for turbulent heating, turbulent heating, reconnection heating, and turbulent heating, respectively. From left to right: Emissions from every region, including the full, midplane, nearside jet, and farside jet regions. The eDF is modeled as a hybrid of thermal and variable κ components, with ε = 0.5 in κ width, w. Cooling is not included.

In the text

	Fig. 6. As in Fig. 5 but with radiative cooling included.
In the text

Fig. 7. SED curves of different regions under turbulent heating at various accretion rates, without cooling (left) and with radiative cooling (right). From top to bottom: Time-averaged mass accretion rates of ṀBH/ṀEdd = 1 × 10−6, 5 × 10−6, and 1 × 10−5. All curves adopt the hybrid thermal and variable κ eDF. The solid curves represent average values. The points correspond to the peaks. The shaded regions denote the standard deviation relative to the average values. The vertical dash-dotted lines correspond to 86 GHz, 230 GHz, 340 GHz, and 136 THz.

In the text

Fig. 8. Light curves of the 230 GHz flux at a 163° inclination angle and spin a = 0.9375. The curves are plotted for accretion rates of ṀBH/ṀEdd = 1 × 10−6 (top), 5 × 10−6 (middle), and 1 × 10−5 (bottom). The solid curves represent the non-cooling cases, and the dashed curves represent the radiative cooling cases. All curves assume turbulent heating and adopt the eDF modeled as a hybrid of thermal and variable κ components.

In the text

Fig. 9. Total flux variation in turbulent heating in the non-cooling (dots) and radiative cooling (squares) cases at 230 GHz at a 163° inclination angle and spin a = 0.9375. The different colors correspond to the different accretion rates: ṀBH/ṀEdd = 1 × 10−6 (black), 5 × 10−6 (red), and 1 × 10−5 (blue). The labels on the x-axis denote the emissions from each region: the full region is depicted first (Total), followed by the midplane (Mid), the nearside jet (Near), and the farside jet (Far). The colored dots and squares represent the time-averaged values, with the error bars indicating the standard deviation about the mean.

In the text

	Fig. A.1. Same as Fig. 5 but using σcut = 2, 5, 10, and 25, respectively, with radiative cooling, in turbulent heating, and the mass accretion rate ṀBH/ṀEdd = 1 × 10−5.
In the text

	Fig. B.1. Same as Fig. 2 for panel (a), but panels (b) – (e) highlight the differences in linear scale by subtracting the density depicted in panel (a). The solid black curves represent σ = 1.
In the text

	Fig. C.1. Same as Fig. 4 but the curves in various colors correspond to different cases: the turbulent heating model without cooling (black), with cooling of ṀBH/ṀEdd = 5 × 10−6 (blue), and the reconnetion heating model without cooling (red), with cooling of ṀBH/ṀEdd = 5 × 10−6 (magenta).
In the text

	Fig. C.2. Same as Fig. 7 for simulation C_KA5e-6 but the case C_MR5e-6 is also added.
In the text

	Fig. C.3. Time- and azimuthally averaged distributions of κ for the PIC-TURB obtained from decaying plasma turbulence [Eq. (6)] in cases NC and C_KA5e-6. The solid black and dashed white curves represent σ = 1 and the Bernoulli parameter, −hut = 1.02, respectively.
In the text