Two-temperature treatments in magnetically arrested disc GRMHD simulations more accurately predict light curves of Sagittarius A*

ABSTRACT

The Event Horizon Telescope Collaboration (EHTC) observed the Galactic centre source Sagittarius A|$^\ast$| (Sgr A|$^{\ast }$|⁠) and used emission models primarily based on single ion temperature (1T) general relativistic magnetohydrodynamic (GRMHD) simulations. This predicted emission is strongly dependent on a modelled prescription of the ion-to-electron temperature ratio. The most promising models are magnetically arrested disc (MAD) states. However, nearly all MAD models exhibit larger temporal fluctuations in radiative 230 GHz emission compared to observations. This limitation possibly stems from the fact that the actual temperature ratio depends on microphysical dissipation, radiative processes, and other effects not captured in ideal fluid simulations. Therefore, we investigate the effects of two-temperature (2T) thermodynamics in MAD GRMHD simulations of Sgr A|$^{\ast }$|⁠, where the temperatures of both species are evolved. We find that the 230 GHz synchrotron flux variability more closely matches historical observations when we include the 2T treatment compared to 1T simulations. For the low accretion rates of Sgr A|$^{\ast }$|⁠, a common assumption is to neglect radiative cooling. However, we find that the radiative cooling of electrons – via synchrotron, inverse Compton, and bremsstrahlung processes – reduces the electron temperature in the inner disc, where the EHT observes, by about 10 per cent, which, in turn, decreases both the (sub)millimetre synchrotron flux and its temporal fluctuations compared to uncooled simulations.

1 INTRODUCTION

Sgr A|$^{\ast }$| was discovered as a bright compact radio source in the centre of the Milky Way galaxy (Balick & Brown 1974; Ekers et al. 1975; Lo et al. 1975). Subsequent observations provided compelling evidence that this object is a supermassive black hole (SMBH) by analysing its proper motion and the dynamics of individual stars in orbit around it (Schödel et al. 2002; Ghez et al. 2003, 2008; Gillessen et al. 2009; GRAVITY Collaboration 2018, 2019; Do et al. 2019). More recently, the Event Horizon Telescope Collaboration (EHTC) presented direct evidence for the presence of this SMBH via imaging of near-event horizon regions (Event Horizon Telescope Collaboration et al. 2022a, 2024a). The accretion flow of Sgr A|$^{\ast }$| is collisionless based on the density inferred by e.g. Yuan, Quataert & Narayan 2003, Dexter et al. 2010, and Bower et al. 2019; Event Horizon Telescope Collaboration et al. (2022c, 2024b). In collisionless plasmas, the electron–ion collision time-scale is much longer than the accretion time-scale; electrons and ions are decoupled and not in thermal equilibrium, such that they may have two different temperatures (2T) (Shapiro, Lightman & Eardley 1976; Rees et al. 1982; Mahadevan & Quataert 1997; Quataert 1998). These accretion flows are typically modelled with ideal general relativistic magnetohydrodynamics (GRMHD), which does not capture its collisionless and 2T nature. The most common GRMHD models treat the fluid as composed of single temperature ions (1T), where the electron density and temperature are not considered in the evolution equations (e.g. Gammie, McKinney & Tóth 2003; Tchekhovskoy, Narayan & McKinney 2011; McKinney, Tchekhovskoy & Blandford 2012; Narayan et al. 2012). The (sub)millimetre emission is dominated by synchrotron radiation from relativistic electrons, making it crucial to accurately model electron thermodynamics. Ressler et al. (2015, 2017) introduced a method to evolve the GRMHD equations describing a gas consisting of ions and electrons that share the same dynamics but have independent thermodynamical evolution. In this so-called 2T treatment, there is only an additional electron entropy equation, while the particle number and energy-momentum equations continue to assume a single fluid (ions). Therefore, we investigate how incorporating this 2T treatment and the commonly ignored electron radiative cooling impacts the predicted (sub)millimetre variability from simulations of Sgr A|$^{\ast }$|⁠.

Decades of observations of Sgr A|$^{\ast }$| give very strong constraints on the (sub)millimetre variability, which is usually quantified with the modulation index |$M_{3} \equiv \sigma _{3}/\mu _{3}$|⁠, where the standard deviation |$\sigma _{3}$| and the mean |$\mu _{3}$| are measured over three hour time intervals of the light curve. Wielgus et al. (2022) used the Atacama Large Millimeter/submillimeter Array (ALMA) and the Submillimeter Array (SMA) as individual interferometers during the EHT observations on 2017 April 5–11. They studied the light curves of Sgr A|$^{\ast }$| at four frequencies bands between 213 and 229 GHz, with a minimum cadence of approximately 10 s. They reported that Sgr A|$^{\ast }$| exhibited a low flux density of |$2.4 \pm 0.2\, \mathrm{Jy}$| and overall low variability, with non overlapping values |$M_{3}=[0.024\!-\!0.051]$| across April 5–10. On April 11, the ALMA observations immediately followed an X-ray flare, with a corresponding enhanced variability |$M_{3}=[0.084\!-\!0.117]$|⁠. The modulation index is consistent with other observations in 2005–2019 at frequencies around 230 GHz documented in various published works with longer cadence and lower number of collected data points (Marrone et al. 2006, 2008; Yusef-Zadeh et al. 2009; Dexter et al. 2014; Bower et al. 2018; Fazio et al. 2018; Iwata et al. 2020; Murchikova & Witzel 2021; Witzel et al. 2021).

1T GRMHD simulations include spherical accretion models (e.g. Ressler et al. 2021; Lalakos et al. 2024; Galishnikova et al. 2025), wind-fed models (Ressler et al. 2020; Ressler, White & Quataert 2023), and torus-initialized models in the weakly magnetized (SANE) and magnetically arrested (MAD; Bisnovatyi-Kogan & Ruzmaikin 1974; Narayan, Igumenshchev & Abramowicz 2003) regimes. In the MAD regime, the accretion is choked by the strong horizon penetrating magnetic field (e.g. Tchekhovskoy et al. 2011; Porth et al. 2021). In these 1T simulations, the ion-to-electron temperature¹ ratio |$T_i/T_e$| is determined using the so-called |$R(\beta)$| prescriptions (e.g. Mościbrodzka, Falcke & Shiokawa 2016; Anantua, Ressler & Quataert 2020), which is the main uncertainty in EHT modelling. None of the EHTC models of Sgr A|$^{\ast }$| fully satisfy all the constraints drawn from multiwavelength observations at 86 GHz, 230 GHz, 2.2 |$\mathrm{\mu m}$|⁠, and in the X-ray (for a detailed explanation, see Event Horizon Telescope Collaboration et al. (2022b,c)). 230 GHz light curve variability at poses a particularly stringent challenge, as nearly all MAD models – and a significant portion of SANE models – exhibit greater variability (⁠|$M_3 \lesssim 0.5$|⁠) than seen in historical observations (Wielgus et al. 2022 and references therein). On the other hand, the more realistic stellar wind-fed accretion models better predict the submillimetre variability in Sgr A|$^{\ast }$| due to the comparatively lower levels of small-scale turbulence compared to SANE and MAD models (Murchikova, White & Ressler 2022). This variability problem could potentially be attributed to not modelling the evolution of |$T_e$| when using the |$R(\beta)$| prescriptions. In reality, |$T_e$| is fundamentally influenced by microphysical plasma and radiation interactions, and does not depend trivially on |$T_i$|⁠. A first-principles kinetic approach is required to model these collisionless effects (Parfrey, Philippov & Cerutti 2019; Crinquand et al. 2022; Galishnikova et al. 2023a). None the less, it is possible to effectively model the electron thermodynamics of thermal electrons with 2T treatments in extended GRMHD simulations (e.g. Ressler et al. 2015) and by including radiative effects (e.g. Sądowski et al. 2017; Chael et al. 2018).

The mechanisms of heating in collisionless plasma remain largely unconstrained. Typically used heating models include prescriptions for weakly collisional turbulent cascades (Howes 2010, H10; Kawazura, Barnes & Schekochihin 2019, K19) and magnetic reconnection mechanisms (Rowan, Sironi & Narayan 2017; 2019, R17, R19). These heating models have limited applicability, as reconnection and turbulence likely occur in diverse conditions across the disc and jet, making it improbable for a single, universally applied model to capture their complexity. Ryan, Dolence & Gammie (2015), Ryan et al. (2017), Sądowski et al. (2017), Chael et al. (2018), Chael, Narayan & Johnson (2019), and Liska et al. (2024) implemented the 2T treatment in a GR radiation MHD (GRRMHD) scheme with coupling between gas, radiation, magnetic fields, and gravity. This 2T treatment does not capture realistic dissipative heating (see Section 2.1). When assuming that ions and electrons are heated through magnetic reconnection (R17 model), it was demonstrated that electrons are always cooler than ions and that electron heating is more uniform (Chael et al. 2018). Additionally, the R17 model was favoured in a polarimetric parameter survey of Sgr A|$^{\ast }$|⁠, comparing H10, K19, and R17 (Dexter et al. 2020). Because of the 2T nature of the accretion flow, a gas mixture consisting of relativistic electrons with an adiabatic index |$\gamma _{e} = 4/3$| and non-relativistic ions with |$\gamma _{i} = 5/3$| was found for the mass accretion rate of Sgr A|$^{\ast }$| and R17 model (Chael et al. 2018; Liska et al. 2024). As a result, assuming a single fixed value for the gas adiabatic index could lead to inconsistencies between the total and individual temperatures and pressures. The adiabatic index characterizes the fluid response to compression, relating gas pressure |$p_g$|⁠, and density |$\rho$| via |$p_g \propto \rho ^\gamma$|⁠. Distinct adiabatic indices are important because electron thermodynamics govern observable emissions, while ion dynamics influence the bulk flow. Therefore, it is crucial to conduct a detailed comparison between the evolved |$T_e$| from 2T GRMHD simulations with variable adiabatic indices and the |$T_e$| calculated using the |$R(\beta)$| prescription.

The estimated Eddington ratio²|$\dot{M}_{\mathrm{edd}}\approx 0.05\, {\rm M}_{\odot }\,\mathrm{yr}^{-1}$| of Sgr A|$^{\ast }$| is |$f_{\mathrm{edd}}\equiv \dot{M}/\dot{M}_{\mathrm{edd}}\sim [1\!-\!100]\times 10^{-8}$| (e.g. Agol 2000; Bower et al. 2003; Marrone et al. 2007). A common assumption for Sgr A|$^{\ast }$| is that radiative cooling and transport do not have a significant impact on the accretion dynamics and (sub)millimetre emission. In the advection-dominated accretion flow (ADAF) model, the radiative efficiency |$\eta _{\mathrm{rad}}$| is less than 0.001 at all radii, meaning that less than 0.1 per cent of the available accretion energy is radiated (Narayan & Yi 1994; Narayan, Yi & Mahadevan 1995). Fragile & Meier (2009) and Dibi et al. (2012) were pioneers in incorporating electron radiative cooling processes (Esin et al. 1996), including bremsstrahlung, synchrotron, and inverse Compton, into simulations. When |$f_{\mathrm{edd}} \gt 10^{-7}$| for a 1T GRMHD SANE state, electron cooling lowers the scale-height of the accretion disc and the overall flux of the spectra from the submillimetre to the far-UV (Drappeau et al. 2013; Yoon et al. 2020). Similarly, variations in |$T_i/T_e$| were identified when |$f_{\mathrm{edd}} \gtrsim 10^{-7}$| in a 2T treatment (Dihingia et al. 2023). For a SANE state, radiation has a negligible effect on either the dynamics or the thermodynamics of the accreting gas when |$f_{\mathrm{edd}}\sim 2\times 10^{-8}$| (Sądowski et al. 2017). Radiative cooling was found to lower |$T_e$| in the inner regions of the accretion flow for |$f_{\mathrm{edd}} \gt 10^{-6}$| in SANE states (Ryan et al. 2015, 2017). On the other hand, a MAD state has significantly higher level of magnetic flux and therefore higher |$\eta _{\mathrm{rad}}$| due to more efficient synchrotron emission. When |$f_{\mathrm{edd}} \approx 10^{-7}$|⁠, |$\eta _{\mathrm{rad}}\approx 0.03$| for a MAD state and |$\eta _{\mathrm{rad}}\approx 0.002$| for a SANE state (Liska et al. 2024). Therefore, radiative cooling may be important for simulating Sgr A|$^{\ast }$| if the accretion flow is in a MAD state.

In this paper, we study the variability of synchrotron radiative emission at [43–1360] GHz predicted from 2T GRRMHD simulations in a MAD state. We include variable adiabatic index (Sądowski et al. 2017), and radiative cooling of electrons (synchrotron, inverse Compton, and bremsstrahlung; Esin, McClintock & Narayan 1997) appropriate for Sgr A|$^{\ast }$|⁠. In our implementation, the cooling source terms subtract energy and momentum (see Section 2.1; Fragile & Meier 2009; Dibi et al. 2012). Synchrotron emission dominates the radiative cooling rates and flux at [43–1360] GHz (Yoon et al. 2020). We also include Coulomb coupling that characterizes the transfer of energy between ions and electrons that occurs through collisions in a kinetic manner (see appendix  A; Stepney & Guilbert 1983). Additionally, we assume that electron/ion heating occurs via magnetic reconnection (R17 model), and we do not account for anisotropic thermal conduction along magnetic field lines, as considered in Ressler et al. (2015). Recent studies have incorporated similar or more advanced physics at lower grid resolutions (Chael et al. 2018; Dexter et al. 2020; Ressler et al. 2023). However, an analysis of variability that systematically adds one layer of physics at a time has been lacking. We quantify the variability of the total synchrotron flux density using the three-hour modulation index |$M_3$|⁠. Our results show that the inclusion of 2T treatment, variable adiabatic index, and cooling systematically shifts the |$M_3$| distribution towards lower variability values at 86, 228, 345, and 1360 GHz.

The structure of this paper is as follows: Section 2 provides the 2T GRMHD equations and numerical set-up for our simulations. The results of electron radiative cooling, heating, variable adiabatic index, and ray tracing are discussed in Section 3. The conclusions of the study are presented in Section 4.

2 SIMULATIONS

2.1 2T GRMHD equations

We use geometrized units with gravitational constant, black hole mass, and speed of light |$G=M=c=1$|⁠, and a factor of |$1/\sqrt{4\pi }$| is absorbed in the definition of the magnetic field. The gravitational radius is |$r_\mathrm{ g} \equiv GM/c^2$|⁠. Greek indices run through [0,1,2,3] and Roman indices through [1,2,3]. The metric determinant is g. A semicolon (;) is used to represent a covariant derivative. In the M1 closure scheme implemented in Sądowski et al. (2017), Chael et al. (2018), and Liska et al. (2024), radiation (approximated as a fluid) and gas are coupled through the exchange of four-momentum, a process known as radiative transport. We do not include the radiation energy–momentum equation, i.e. we neglect radiative transport. Instead, we solve the equations of GRMHD that comprise the particle number conservation equation:

the energy–momentum conservation equations:

and the Maxwell’s equations:

In ideal MHD, the dual of the electromagnetic field tensor is |$F^{\ast \mu \nu }=b^{\mu }u^{\nu }-b^{\nu }u^{\mu }$| and the stress–energy tensor is

here n is the particle number density, |$\rho =n m_\mathrm{ p}$| is the rest mass density³, |$m_\mathrm{ p}$| is the proton mass, |$u_\mathrm{ g}$| is the gas energy density, |$p_\mathrm{ g}$| is the gas pressure, and |$b^2/2$| is the magnetic energy density. The 4-velocity and the magnetic field 4-vector are |$u^{\mu }$| and |$b^{\mu }$|⁠, respectively. |${\delta ^{\mu }}_{\nu }$| is the Kronecker delta. The gas pressure is proportional to the gas temperature (⁠|$p_\mathrm{ g}\propto T_\mathrm{ g}$|⁠). The magnetic pressure is |$p_\mathrm{ b}=b^2/2$|⁠, where b is the magnetic field strength in the frame co-moving with the fluid. |$q^-$| is the radiative cooling rate of electrons (see Section 3.1). We assume that proton cooling is insignificant, as proton synchrotron emission is negligible in the (sub)millimetre regime (e.g. Event Horizon Telescope Collaboration 2022c). We consider a purely hydrogen fluid, so the relative mass abundances of hydrogen and helium are |$X=1$| and |$Y=0$|⁠. Therefore, the number densities are |$n_{e}=n_{i}=n$|⁠.

Single fluid GRMHD assumes that both electrons and ions move at the same bulk velocity |$u_\mathrm{ e}^{\mu }=u_\mathrm{ i}^{\mu }=u^{\mu }$|⁠, so their momenta are considered separately. The first law of thermodynamics determines the evolution of the species entropies,

where |$s_{\mathrm{ e},\mathrm{ i}}$| and |$T_{\mathrm{ e},\mathrm{ i}}$| are the entropy per particle and temperature of electrons and ions, respectively. |$q^{\textrm {v}}=q^\mathrm{v}_e+q^\mathrm{v}_i$| is total dissipative heating rate (or ‘viscous’ heating rate), |$\delta _e=q^\mathrm{v}_e/q^{\textrm {v}}$| and |$\delta _i=1-\delta _e$| the electron- and ion-to-overall heating ratio, and |$q^C$| is the Coulomb coupling rate (energy exchange between electrons and ions).

The entropy per particle of an ideal gas with fixed adiabatic index |$\gamma$| is |$s \propto \ln (p/\rho ^{\gamma })$|⁠. We follow the approach in Sądowski et al. (2017), and use an approximate relativistic entropy formula,

which is analytically invertible for the dimensionless temperature |$\Theta _{\mathrm{ e},\mathrm{ i}}=k_\mathrm{ B}T_{\mathrm{ e},\mathrm{ i}}/m_{\mathrm{ e},\mathrm{ i}}c^2$|⁠, where |$k_\mathrm{ B}$| is the Boltzmann constant. The ideal gas equation of state provides the means to determine the effective temperature of the gas |$T_\mathrm{ g}$|⁠, which is a mix of electrons and ions,

The electron temperature |$T_\mathrm{ e}$| is obtained from the entropy density equation (5). Both species satisfy the equation of state,

The gas pressure/energy is the sum of the electron and ion pressures/energies, |$p_\mathrm{ g}=p_\mathrm{ e}+p_\mathrm{ i}$| and |$u_\mathrm{ g}=u_\mathrm{ e}+u_\mathrm{ i}$|⁠. Therefore, the temperatures must satisfy, |$T_\mathrm{ g}= \left(T_\mathrm{ e} +T_\mathrm{ i}\right)/2$|⁠, and the effective adiabatic index of the mixture is,

where |$\gamma _{\mathrm{ e},\mathrm{ i}}$| are the electron and ion adiabatic indices (Narayan, Kumar & Tchekhovskoy 2011; Sądowski et al. 2017),

We follow the method of Ressler et al. (2015) to numerically identify the total dissipative heating |$q^{\textrm {v}}$| by evolving the thermal entropies adiabatically over a time-step |$\Delta \tau$|⁠. By comparing the sum of the adiabatically evolved energy densities, |$u_{\mathrm{ i}, \textrm {ad}}$| and |$u_{\mathrm{ e}, \textrm {ad}}$|⁠, to the separately evolved total gas energy |$u_\mathrm{ g}$|⁠, we estimate the dissipative heating in the total fluid as |$q^{\textrm {v}}=(u_\mathrm{ g}-u_{\mathrm{ i}, \textrm {ad}}-u_{e, \textrm {ad}})/\Delta \tau$|⁠. The numerical implementation of mixing finite-sized fluid parcels into a single homogenized fluid within a cell artificially increases the entropy of the gas, which should not be treated as dissipation. As a result, this method loses accuracy in regions of the accretion flows with large discontinuities, such as those found in MAD states (for more details, see Sądowski et al. 2017; Chael 2025). This method does not capture real dissipation. In reality, the physical processes that produce dissipation occur at scales much smaller than the simulation grid. We do not explicitly model the heating mechanisms via effective resistivity or viscosity. Instead, we assume that heating at the grid scale occurs through sub-grid magnetic reconnection (Rowan et al. 2017).

2.1.1 Heating

Previous studies on grid-scale dissipation have considered mainly two heating prescriptions: turbulent heating and magnetic reconnection. The non-relativistic turbulent heating model H10 was originally developed for solar wind observations and is based on calculations of a turbulent cascade in a weakly collisional plasma with |$\sigma _i \ll 1$| (Howes 2010). The relativistic turbulent heating model K19 exhibits quantitatively similar behaviour to the H10 model (Kawazura et al. 2019). For Sgr A|$^{\ast }$|⁠, assuming H10 results in a negligible difference in electron temperatures in MAD models compared to K19 (Dexter et al. 2020).

We adopt the heating of electrons and ions resulting from magnetic reconnection with zero guide field (Rowan et al. 2017; Chael et al. 2018). This model is based on particle-in-cell simulations with physical mass ratio |$m_\mathrm{ i}/m_\mathrm{ e} = 1836$| and total magnetization for relativistic particles, that properly accounts for relativistic inertia |$\sigma _{\textrm {w}}=b^2/\textrm {w}\geqslant 0.03$|⁠. The electron-to-overall heating ratio |$\delta _\mathrm{ e}$| is

where the enthalpy density per unit volume is |$\textrm {w} = \rho _\mathrm{ i} c^2+\gamma _{\mathrm{ e}}u_\mathrm{ e}+\gamma _{\mathrm{ i}}u_\mathrm{ i}$|⁠. The ratio of ion thermal pressure to magnetic pressure is |$\beta _i=2k\rho _i T_i/\left(m_ib^2\right)$|⁠, where |$\beta _i\leqslant \beta _{i,\:\mathrm{max}}=1/\left(4\sigma _{\textrm {w}}\right)$|⁠. The heating prescription including the effect of guide fields by Rowan et al. (2019) is quantitatively similar to zero guide field case equation (12). Unlike H10, equation (12) varies less rapidly with |$\beta _i$|⁠, is always non-zero and never allocates more than half of the heat to electrons.

2.2 Numerical set-up

We perform the GRMHD simulations using the H-AMR code (Liska et al. 2022), which builds on the HARM2D code (Gammie et al. 2003; Noble et al. 2006). H-AMR employs GPU acceleration within a hybrid CUDA-OpenMP-MPI framework, integrating adaptive mesh refinement (AMR) and local adaptive time-stepping (LAT) to facilitate efficient and scalable simulations. We use spherical Kerr–Schild coordinates, where t, r, |$\theta$|⁠, and |$\phi$| are the temporal, radial, polar, and azimuthal coordinates, respectively. The number of cells in the radial, polar, and azimuthal directions are |$N_{\mathrm{ r}}\times N_{\theta} \times N_{\phi } = 520 \times 224 \times 224$|⁠. The radial domain is |$r=\left[1.2\!-\!2000 \right]r_\mathrm{ g}$|⁠. We use outflow boundary conditions in r, transmissive boundary conditions in |$\theta$|⁠, and periodic boundary conditions in |$\phi$|⁠, as described in Liska et al. (2022). In all runs, the disc is initialized using a torus in hydrostatic equilibrium (Fishbone & Moncrief 1976) around a Kerr black hole with dimensionless spin |$a=0.9375$|⁠. The initial inner edge of the torus is located at |$r=20 r_\mathrm{ g}$| and the pressure maximum at |$r=41 r_\mathrm{ g}$|⁠. The torus is threaded with a single poloidal magnetic field loop, defined by the |$\phi$|-component of the vector potential |$A_{\phi } \propto {\mathrm{max}}[\rho /\rho _{\mathrm{max}}(r/r_{\rm{in}})^3 \sin ^3\theta \exp (-r/400)-0.2 , 0]$|⁠, and normalized to obtain |$\beta =p_{\mathrm{ g}}/p_{\mathrm{ b}}=100$| at pressure maxima, specifically to reach a MAD state. The following floor and ceiling values are employed; the rest-mass density floor is |$\rho _{\mbox{fl}}=\mathrm{max}[b^2/25,10^{-7}r^{-2},10^{-20}]$|⁠, the gas energy density floor is |$u _{\mathrm{ g},\mbox{fl}}=\mathrm{max}[b^2/750,10^{-9}r^{-26/9},10^{-20}]$| and the magnetization ceiling is |$\sigma _{\mbox{max}}=25$| where |$\sigma =b^2/\rho$|⁠. The disc is initialized with the same energy density perturbation |$u_\mathrm{ g}(1 + 0.04(\textrm {rand([0,~1])}-0.5))$|⁠, where |$\textrm {rand([0,~1])}$| is a random deviate between 0 and 1.

For the 1T simulation, we assume a non-relativistic gas with a constant adiabatic index of 5/3. Table 1 summarizes the main simulation set-up parameters, and results of 228 GHz flux density and Eddington ratio. The 1T simulation is scale-free, whereas the 2TC simulations require a density scale |$\rho _\mathrm{ s}$| to compute the physical radiative cooling rates. To obtain the appropriate density scale, we performed an iterative process: running multiple 2TC simulations, performing ray-tracing, and verifying the time-averaged flux density at 228 GHz. The 2TC simulation (without extra labels) matches the 228 GHz flux density 2.7 Jy (Event Horizon Telescope Collaboration et al. 2022b). The 2TC-x.x|$\mathrm{Jy}$| simulations were initialized with higher density scales, resulting in higher accretion rates and increased 228 GHz synchrotron flux densities, |$F(\nu)\gt 8\mathrm{Jy}$|⁠. A higher plasma density increases opacity, thereby reducing synchrotron flux variability. Therefore, we compare the variability of synchrotron radiative emission only between simulations 1T, 2T, and 2TC with |$2.7\mathrm{Jy}$|⁠. 2TC follows the same initial set-up as 2T but with cooling enabled. The values reported in Table 1 and Section 3 for simulation 2T are obtained using the same density scale as in simulation 2TC. The primary simulations (1T, 2T, 2TC) are evolved for a time of approximately |$3\times 10^4r_g/c$|⁠, with a cadence of |$10r_\mathrm{ g}/c$|⁠, equivalent to 3.5 min for Sgr A|$^{\ast }$|⁠.

3 RESULTS

In order to determine quantitative differences between the simulations, we compare various fluid parameters. The averaged profile of a variable X is calculated by integrating over |$\theta$| and |$\phi$|⁠. The accretion flow is taken to be the region that satisfies magnetization |$\sigma =b^2/\rho \lt 1$|⁠. Additionally, we include a density weight to give more relevance to regions of the disc with higher density and therefore higher thermal synchrotron emissivity |$j_{\nu } = j_{\nu }(\rho ,B,T_\mathrm{ e})$|⁠,⁴

The fluxes are defined as follows. The mass accretion rate |$\dot{M}$| is given by:

the energy flux |$\dot{E}$| is given by:

the magnetic flux is defined as:

For |$r\lesssim 100r_\mathrm{ g}$|⁠, there is an inflow equilibrium in the disc so |$\dot{M}_{\textrm {disc}}$| is approximately constant. The disc properties still depend on the initial conditions of the Fishbone–Moncrief torus for radii larger than |$100 r_\mathrm{ g}$|⁠. To study MAD simulations, it is convenient to use the normalized magnetic flux, |$\varphi \equiv \Phi /\sqrt{ \left\langle \dot{M} \right\rangle }$|⁠, known as the ‘MAD parameter’ which, for spin |$a=0.9375$| and torus scale height |$H/R \approx 0.3$|⁠, has the critical value |$\varphi _{max}\approx 40-50$| (when using cgs units in the ratio; Tchekhovskoy, McKinney & Narayan 2012). For SANE models, |$\varphi _{max}\approx 7$| (Chatterjee & Narayan 2022). H and R are the full height and cylindrical radius of the disc, respectively. The scale height is defined geometrically as |$H/R\equiv \left\langle \left| \theta -\pi /2 \right| \right\rangle _{\rho [\theta ,\phi ]}^{\textrm {disc}}$|⁠. In general, subscripts indicate where in the radial domain the variables are analysed, e.g. |$X_{\mathrm{ BH}}$| at the event horizon and |$X_{10r_g}$| at |$10r_g$|⁠.

3.1 Radiative cooling

The total cooling rate for an optically thin gas is computed from the cooling function,

where |$q_{\rm br}^{-}$| and |$q_{\rm Sync}^{-}$| are the bremsstrahlung and synchrotron cooling rates, respectively. We find that Bremsstrahlung cooling is negligible, consistent with Yoon et al. (2020). |$\eta _{\rm br,C}$| and |$\eta _{\rm Sync,C}$| are the Compton enhancement factors for bremsstrahlung and synchrotron, respectively. These factors are the average energy gain of the photon in an assumption of single scattering (Esin et al. 1996). A detailed description of the cooling rates is available in Yoon et al. (2020).

Fig. 1 shows a comparison of |$T_\mathrm{ e}$| for simulations 2T and 2TC. Even for the very low accretion rate relevant for Sgr A|$^{\ast }$|⁠, |$f_{\mathrm{edd}}\approx 10^{-7}$|⁠, we find that radiative cooling still lowers |$T_\mathrm{ e}$| in the inner disc and the average synchrotron flux by 10 per cent (see Table 1). The importance of the radiative losses is tested by computing the radiative cooling time-scale |$\tau = u_e/q^-$| and compare it to the dynamical accretion time-scale. We calculate the synchrotron time-scale as

the Comptonization time-scale as

and the accretion time-scale as

MAD states are characterized by large fluctuations caused by episodic magnetic flux eruptions. Once magnetic flux becomes oversaturated, magnetic reconnection occurs near the black hole, ejecting low-density, highly magnetized flux tubes into the disc (e.g. Chatterjee et al. 2021; Porth et al. 2021; Ripperda et al. 2022). Fig. 2 presents time series between |$[23-28]\times 10^3 r_\mathrm{ g}/c$| of the normalized magnetic flux and the previously defined time-scales. Across all simulations, the average flux |$\left\langle \varphi \right\rangle \approx 35$| reveals distinct signatures of several magnetic flux eruption events. Our analysis indicates that the time-scales, ranked by relevance, are as follows: synchrotron |$\left\langle \tau _{\mbox{Sync}} \right\rangle / \left\langle \tau _{\mbox{accr}}\right\rangle \sim 3$|⁠, and Comptonization |$\left\langle \tau _{\mbox{IC}}\right\rangle / \left\langle \tau _{\mbox{accr}}\right\rangle \sim 500$|⁠. Therefore, we find that synchrotron emission dominates over inverse Compton, consistent with Yoon et al. (2020) and Liska et al. (2024).

In 2T GRRMHD simulations, Liska et al. (2024) demonstrated that magnetic flux significantly influences radiative efficiency, primarily attributed to more efficient synchrotron emission in MAD models at low accretion rates (see also Ryan et al. 2017; Dexter, Scepi & Begelman 2021). At |$f_{\mathrm{edd}} \approx 10^{-7}$|⁠, |$\eta _{\mathrm{rad}} \approx 0.002$| and |$T_\mathrm{ i}/T_\mathrm{ e} \approx 10$| at |$10r_\mathrm{ g}$| in the disc for SANE models, while for MAD models, |$\eta _{\mathrm{rad}} \approx 0.03$| and |$T_\mathrm{ i}/T_\mathrm{ e} \approx 3$| at |$10r_g$| in the disc (Liska et al. 2024). Their more self-consistent calculation of radiative efficiency incorporates the radiation stress-energy tensor in the M1 scheme. In our simulations, which include only radiative cooling, we estimate the radiative efficiency |$\eta _{\mathrm{rad}} \approx \eta _{\mathrm{NT}} \left(\tau _{\mbox{accr}}/ \tau _{\mbox{cool}}\right)_{10r_g} \approx 0.05\pm 0.1$|⁠, where |$\tau _{\mbox{cool}}=(1/\tau _{\mbox{Sync}}+1/\tau _{\mbox{IC}})^{-1}$|⁠. Alternatively, using the cooling rate |$\eta _{\mathrm{rad}}\approx \sum _{5r_\mathrm{ g}}^{200r_\mathrm{ g}} q^{-}/\left\langle \dot{M}_{5r_\mathrm{ g}}\right\rangle \approx 0.02\pm 0.01$| (see Fig. 3).

3.2 Heating and variable adiabatic index

Fig. 4 shows the correlation between the magnetization (⁠|$\sigma _{\textrm {w}}$|⁠), ion plasma beta (⁠|$\beta _i$|⁠), and the electron-to-overall heating ratio |$\delta _e$| determined by equation (12). In general, |$\delta _e$| determines the temperature and adiabatic index of both species and the gas mixture (see equations 5, 6, 10, and 11). By assuming heating through magnetic reconnection, we observe the increase of |$T_\mathrm{ e}$| more prominently in regions with high |$\sigma _{\textrm {w}}$| and low |$\beta _i$|⁠. In those regions, the adiabatic index of the electron/ion mixture, determined by |$T_\mathrm{ e}$| and |$T_\mathrm{ i}$|⁠, approaches 4/3, consistent with the plasma becoming relativistic as a result of heating. The 1T simulation employs a fixed value of 5/3. During an eruption in MAD states, the expulsion of magnetic flux goes through the magnetic reconnection of field lines in a current sheet in the equatorial plane. Highly magnetized plasma from the jet spine region (⁠|$\sigma =b^2/\rho \gt 1$|⁠) supplies matter to the current sheet, with |$T \propto \sigma _{\mbox{max}}$|⁠, and the reconnection exhaust deposits this hot plasma in the jet-disc interface (Ripperda et al. 2022).

In the first-principles kinetic simulations of (Rowan et al. 2017), magnetic energy is transformed into heat or kinetic energy when magnetic fields reconnect. This process accelerates and heats particles, with electron heating being particularly prominent in regions where magnetic pressure is high, like in the disc–jet boundary and flux tubes during episodes of magnetic flux eruption events. In GRRMHD SANE simulations of Sgr A|$^{\ast }$|⁠, Chael et al. (2018) demonstrated that electrons are preferentially heated in the polar outflows with the H10 prescription (see also e.g. Ressler et al. 2017; Sądowski et al. 2017), whereas electrons are heated by nearly the same fraction everywhere in the accretion flow with the R17 prescription. Additionally, the models assuming the R17 heating model show variability in the 230 GHz light curve more consistent with the level observed from Sgr A|$^{\ast }$| (Wielgus et al. 2022 and references therein).

Fig. 5 shows radial profiles of the scale-height, adiabatic indices, electron and ion temperatures, and the synchrotron, Comptonization, and accretion time-scales (see Section 3.1). Our findings reveal that the scale-height in 1T simulations is larger compared to the 2T and 2TC simulations, but within the one standard deviation range. The higher adiabatic index in the 1T simulation results in a higher |$T_\mathrm{ i}$| and a larger scale-height. Comparing the 2T and 2TC simulations, we find that radiative cooling reduces |$T_\mathrm{ e}$| by approximately 15 per cent within |$50r_\mathrm{ g}$|⁠, with the most significant drop occurring between |$[30\!-\!50]r_\mathrm{ g}$|⁠. This cooling effect also lowers the values of |$\gamma _\mathrm{ e}$| and |$\gamma _\mathrm{ g}$|⁠. We find that while the radial profile of the time-scale ratio |$\tau _{\mbox{Sync}}/\tau _{\mbox{accr}}$| remains relatively constant, |$\tau _{\mbox{IC}}/\tau _{\mbox{accr}}$| increases significantly with increasing radius.

We observe a temperature ratio of |$\left\langle T_\mathrm{ i}/T_\mathrm{ e} \right\rangle \approx 3$| and |$\left\langle T_e \right\rangle \sim 10^{11} \text{K}$|⁠, both time averaged between |$16\times 10^3 r_\mathrm{ g}/c$| and |$29\times 10^3 r_\mathrm{ g}/c$|⁠, and measured at |$10r_g$|⁠, consistent with those reported in 2T simulations (e.g. Liska et al. 2024), and in spectral fitting modelling of 1T simulations where the ion-to-electron temperature ratio is a model parameter (e.g. Mościbrodzka et al. 2009; Drappeau et al. 2013). Mościbrodzka et al. (2009) found that models with |$T_\mathrm{ i}/T_\mathrm{ e} = 1$| were inconsistent with the submillimeter spectral slope, while those with |$T_\mathrm{ i}/T_\mathrm{ e} = 3$| and 10 provided a better match to spectral and VLBI observations. Furthermore, with the inclusion of radiative cooling (Dibi et al. 2012), Drappeau et al. (2013) identified an optimal spectral fit for models with |$T_\mathrm{ i}/T_\mathrm{ e} = 3$|⁠.

3.3 Ray tracing

The ray tracing of the GRMHD simulations is performed by solving the general-relativistic radiative transfer (GRRT) equations using the bhoss code (Younsi, Wu & Fuerst 2012; Younsi et al. 2016).Radiative processes are calculated using synchrotron emission and absorption, which are sufficient for imaging and light-curve analysis across radio and submillimetre frequencies. In the submillimetre regime, synchrotron emission dominates the spectrum (Yoon et al. 2020; Event Horizon Telescope Collaboration et al. (2022c)). The 2.2 |$\mu \mathrm{m}$| emission is usually dominated by synchrotron, and the X-ray can be dominated by either Compton scattering or bremsstrahlung (Event Horizon Telescope Collaboration et al. 2022c). To match the average synchrotron flux density of Sgr A|$^{\ast }$| at 2.7 Jy at 228 GHz (Event Horizon Telescope Collaboration et al. 2022b), the density scale is adjusted in the GRRT post-processing of the 1T simulation. For the 2T models, the density scale is directly obtained from the GRMHD simulations (see Table 1).

Based on the mean values of Sgr A|$^{\ast }$| reported by Do et al. (2019) and GRAVITY Collaboration (2019), we assume the mass |$M=4.14\times 10^6\,{\rm M}_{\odot }$| (same as in the 2T simulations) and distance |$D=8.127\,\mathrm{kpc}$|⁠, where |${\rm M}_{\odot }$| is the solar mass (Event Horizon Telescope Collaboration et al. 2022c). Since there is no definitive evidence of a jet in Sgr A|$^{\ast }$| determining the exact orientation of the source relative to our line of sight remains challenging. Nevertheless, the model comparison of the EHTC suggests that high inclination angles (⁠|$i \gt 50^{\circ }$|⁠) are unlikely (Event Horizon Telescope Collaboration et al. 2022c). We analyse the 1T simulation considering six inclination angles, |$i=[10, 30, 50, 130, 150, 170]^{\circ }$|⁠, and the 2T simulations with the preferred value |$i=30^{\circ }$| that passes most total intensity and polarimetric constraints (Event Horizon Telescope Collaboration et al. 2022c, 2024b). We find that |$M_3$| distributions do not significantly depend on the inclination angle.

In the 1T simulation, the strength of electron–proton coupling is unknown, yet |$T_\mathrm{ e}$| is crucial to calculate the radio synchrotron spectra assuming a thermal electron distribution function. Similarly as in Mościbrodzka et al. (2016); Mościbrodzka et al. (2017), we calculate |$T_\mathrm{ e}$| using the parametrization with respect to the plasma beta (⁠|$\beta \equiv p_\mathrm{ g}/p_\mathrm{ b}$|⁠),

This |$R(\beta)$| model was initially developed phenomenologically for SANE simulations and may not be suitable for MAD states, particularly in predicting submillimetre light curves. In regions of the accretion flow where |$\beta \gtrsim 1$| or |$\lesssim 10^{-2}$|⁠, the electron temperature is approximately independent of |$\beta$|⁠. Consequently, |$T_\mathrm{ e}$| fluctuates approximately as much as |$T_i$| (see Section 3.4).

For Sgr A|$^{\ast }$|⁠, electron distribution functions with approximately equal proton and electron temperatures (⁠|$R_{\mbox{high}} = 1$|⁠) are unlikely based on the 2017 data observed by the EHT, and the two most promising models have |$R_{\mbox{high}} = 160$| (Event Horizon Telescope Collaboration et al. 2022c, 2024b), so that the electrons are much colder in the disc and hotter in the jet. However, 2T GRMHD simulations with H10, K19, and R17 heating prescriptions have consistently shown that time averaged |$T_\mathrm{ i}/T_\mathrm{ e}\approx [3\!-\!10]$| in the accretion disc (e.g. Chael et al. 2018; Mizuno et al. 2021; Dihingia et al. 2023; Liska et al. 2024). The sole value of |$R_{\mbox{low}}=1$| explored by the EHTC could potentially skew the high values of |$R_{\mbox{high}}$|⁠. When |$f_{\mathrm{edd}} \approx 10^{-6}$| and |$a=0.94$|⁠, Mizuno et al. (2021) found that MAD models using the R17 heating prescription show only a slight reduction in 230 GHz light-curve variability compared to predictions from 1T simulations using the |$R(\beta)$| prescription. Moscibrodzka (2024) demonstrated that 2T MAD models with K19 heating exhibit less variability at 228 GHz compared to |$R(\beta)$| models, even when using a constant gas adiabatic index of |$\gamma = 13/9$|⁠. Additionally, |$M_3$| increases with black hole spin and slightly decreases when physics of non-thermal electrons are included. Moreover, their resolved images most closely resemble 1T models with |$T_\mathrm{ i}/T_\mathrm{ e} = 10$| in both linear and circular polarization.

We explore five values |$R_{\mbox{high}}=[1,10,20,40,160]$|⁠. We find that the evolved |$T_\mathrm{ e}$| surpasses |$T_\mathrm{ i}$| in the jet–disc interface (see Fig. 4), whereas the |$T_\mathrm{ e}$| calculated using the |$R(\beta)$| prescription invariably remains lower than |$T_\mathrm{ i}$|⁠. Fig. 6 shows a visual comparison of the ray-traced images at inclination |$i=30^{\circ }$| for simulations 2TC and 1T at different frequencies. Images from the 2T simulation are not included in Fig. 6, as they are nearly identical to those from the 2TC simulation. The density scales to match the time-averaged flux density of 2.7 Jy are |$\rho _{\mathrm{scale}}=2.6$| and |$15.4 \, (\mathrm{g\, cm^{-3}})$| for |$R_{\mathrm{high}}=1$| and 160, respectively. As |$R_{\mathrm{high}}$| increases, emission from the disc (⁠|$\beta \gt 1$|⁠) diminishes while emission from the rapidly moving outflows (⁠|$\beta \lt 1$|⁠) intensifies. This effect results in increased synchrotron flux variability, as shown in Table 2.

Table 2 also presents the jet power |$P_{\mathrm{jet}}=\dot{M}_{\mathrm{jet}}-\dot{E}_{\mathrm{jet}}$|⁠, averaged over the region between 50 and |$500 r_\mathrm{ g}$|⁠, where the jet is identified based on the Bernoulli parameter |$Be = -\bar{h} u_t\gt 1.02$| where |$\bar{h}=(\rho +u_\mathrm{ g}+p_\mathrm{ g})/\rho$| is the specific gas enthalpy and |$u_t$| is the time component of the inverse four-velocity (Davelaar et al. 2018). The density scale derived from ray tracing the 1T simulations is used to calculate the jet power in cgs units. For comparison, the jet power in the 2TC and 2T simulations is |$(1.4 \pm 0.4) \times 10^{38}\ \mathrm{erg\,s}^{-1}$| and |$(1.3 \pm 0.3) \times 10^{38}\, \mathrm{erg\,s}^{-1}$|⁠, respectively, which matches with the jet power in the 1T simulation for |$R_{\mbox{high}}$| between 1 and 10.

We obtain distributions of the modulation index |$M_3$| measured over three-hour intervals. Fig. 7 shows the 228 GHz light curves and |$M_3$| distributions⁵ for |$i=30^{\circ }$|⁠. Values of |$M_3$| for the 1T simulation with different inclination angles and |$R_{\mbox{high}}$| are presented in Table 2. A comparison of |$M_3$| at different frequencies for 1T, 2T, and 2TC simulations is presented in Table 3. We find that the 2T treatment, variable adiabatic index and radiative cooling shifts the |$M_3$| distribution towards lower variability values. However, historical observations of Sgr A|$^{\ast }$| (Wielgus et al. 2022 and references therein) still show lower |$M_3$| values than 2T simulations.

From the 2T and 2TC simulations, we extract two light curve predictions per simulation: one using the evolved |$T_\mathrm{ e}$| directly from the simulation, and another where the evolved |$T_\mathrm{ e}$| is ignored and instead calculated in post-processing with the |$R(\beta)$| prescription. We find that |$M_3$| is lower when using the evolved |$T_\mathrm{ e}$|⁠, suggesting that the |$T_\mathrm{ e}$| derived from the |$R(\beta)$| prescription exhibits greater fluctuations (see Fig. 7).

3.4 Fluctuations in electron and ion temperatures

In the |$R(\beta)$| prescriptions used in 1T simulations, |$T_\mathrm{ e}$| fluctuates as much as |$T_i$| in regions with |$\beta \gtrsim 1$| or |$\lesssim 10^{-2}$| where |$T_\mathrm{ e} \approx T_\mathrm{ i}/R$| (Mościbrodzka et al. 2016). However, that is not the case in a 2T treatment with different adiabatic indices for electrons and ions. Thermodynamically, a 2T treatment with |$\gamma _e\ne \gamma _i$| is expected to have an impact in the fluctuation of |$T_\mathrm{ e}$| (Gammie 2025). From the first law of thermodynamics, the convective derivatives of |$T_\mathrm{ e}$| and |$T_\mathrm{ i}$| are (Gammie 2025)

where the heating/cooling rate per unit volume |$q=q^{\textrm {v}}+q^{C}+q^{-}$| includes ion–electron energy exchange |$q^{C}$|⁠, dissipation heating rate |$q^{\textrm {v}}$| and cooling rate |$q^{-}$|⁠. For each species |$q^{\textrm {v}}_\mathrm{ e}=\delta _\mathrm{ e}q^{\textrm {v}}$| and |$q^{\textrm {v}}_\mathrm{ i}=\delta _\mathrm{ i}q^{\textrm {v}}$|⁠, with |$\delta _\mathrm{ e}+\delta _\mathrm{ i}=1$| and the total dissipation rate |$q^{\textrm {v}}=(u_\mathrm{ e}+u_\mathrm{ i})/\tau _{\textrm {diss}}$|⁠. We consider only radiative cooling of electrons, i.e. |$q^{-}_\mathrm{ e}=u_\mathrm{ e}/\tau _{\textrm {cool}}$| and |$q^{-}_\mathrm{ i}=0$|⁠. |$\tau _{\textrm {comp}}$|⁠, |$\tau _{\textrm {diss}}$|⁠, |$\tau _{\textrm {cool}}$| are the time-scales of compression, dissipation, and cooling, respectively.

From our simulations |$\gamma _\mathrm{ e}\approx 4/3$|⁠, |$\gamma _\mathrm{ i}\approx 5/3$|⁠, |$T_\mathrm{ i}/T_\mathrm{ e}\approx 3$|⁠, |$\delta _e\approx [0.2\!-\!0.4]$|⁠, and |$\delta _\mathrm{ i}\approx [0.8\!-\!0.6]$|⁠. Analysing the adiabatic process case |$q=0$|⁠, from equations (22) and (23), |$T_\mathrm{ e}$| fluctuates approximately half as much as |$T_\mathrm{ i}$|⁠,

in the incompressible case |$\tau _{\textrm {comp}}\rightarrow \infty$| with only heating,

and the addition of cooling can reduce the fluctuations of |$T_e$| as compared to |$T_i$|⁠,

Therefore, the difference in adiabatic indices between relativistic electrons and non-relativistic ions causes an effective reduction in the fluctuations of the electron temperature. This theoretical framework is in agreement with our results of |$M_3$| distributions presented in Fig. 7 and Table 3, and it aligns with our interpretation that the evolved |$T_\mathrm{ e}$| in 2T simulations with variable adiabatic indices fluctuates less than the |$T_\mathrm{ e}$| derived from the |$R(\beta)$| prescription.

3.5 Comparison to semi-analytical models

In MAD states, we find substantial heating and synchrotron radiative cooling in areas with significant gas compression and high magnetization. Our results for radiative efficiency align with those of Liska et al. (2024), and challenge the paradigm in ADAF models, with low |$\eta _{\mathrm{rad}}$| across all radii. In our 2T and 2TC simulations, we obtain |$\left\langle T_\mathrm{ e} \right\rangle \approx 10^{11}\, \text{K}$|⁠, time averaged between |$16\times 10^3 r_\mathrm{ g}/c$| and |$29\times 10^3 r_\mathrm{ g}/c$|⁠, and measured at |$10r_g$|⁠, consistent with the values reported in spectral fitting modelling (e.g. Markoff et al. 2001; Yuan, Markoff & Falcke 2002; Dexter et al. 2010) and EHT polarization modelling (Dexter 2016; Event Horizon Telescope Collaboration et al. (2022c, 2024b)).

In the ADAF model by Narayan et al. (1995), the radiative efficiency |$\eta _{\mathrm{rad}}$| is less than 0.001 at all radii, meaning that less than 0.1 per cent of the available accretion energy is radiated, with the flow being advection-dominated. For |$r \lt 200r_\mathrm{ g}$|⁠, ion–electron coupling weakens, causing |$T_\mathrm{ e}$| to saturate, while |$T_\mathrm{ i}$| continues to follow the virial temperature profile. They found that |$T_e$| rises up to |$\sim 10^{10}\,\text{K}$| in accreting black hole flows. For Sgr A|$^{\ast }$|⁠, Quataert & Narayan (1999) achieved a good spectral fit using radiatively inefficient accretion flow (RIAF) models with |$\delta _\mathrm{ e} = 0.55$| and required |$T_\mathrm{ e} \sim 10^{11}\, \text{K}$| close to the BH, which is larger than in ADAF models. In semi-analytical jet models of Sgr A|$^{\ast }$|⁠, Markoff et al. (2001) and Yuan et al. (2002) obtained a good spectral fit with electron temperatures reaching |$\sim 10^{11}\,\text{K}$|⁠.

Recently, the EHTC reported a large resolved polarization fraction of 24–28  per cent, with a peak around 40 per cent (EHTC Sgr A|$^{\ast }$|2017 VII). Synchrotron emission is intrinsically highly polarized, but as light traverses a magnetized medium, it undergoes Faraday depolarization (e.g. Ricarte et al. 2020). Faraday depth, which quantifies the rotation of polarized light, is inversely proportional to |$T_\mathrm{ e}^2$| (Jones & Hardee 1979). Consequently, higher |$T_e$| leads to reduced depolarization. In one-zone models, the total flux and optical depth constraints for Sgr A|$^{\ast }$| indicate small Faraday depths, which account for the high polarization fraction. Based on estimates of |$n_\mathrm{ e} \approx 10^6 \, \mathrm{cm}^{-3}$| and |$B \approx 29 \, \mathrm{G}$|⁠, it is expected that |$T_\mathrm{ e} \approx 10^{11} \, \mathrm{K}$| at |$\beta = 1$| (Dexter 2016; Event Horizon Telescope Collaboration et al. (2022c, 2024b)), consistent with our results. Locally, |$\beta _i$| can drop below unity, particularly at the disc–jet interface and within flux tubes. However, the time-averaged |$\left\langle \beta _i \right\rangle \approx 1$| in the disc at |$\lesssim 20r_\mathrm{ g}$|⁠, as shown in Fig. 4.

4 CONCLUSIONS

The EHTC has provided significant insights into the Galactic centre source Sgr A|$^{\ast }$|⁠. Models based on 1T GRMHD simulations have been able to explain aspects of observations across radio to X-ray wavelengths in quiescent (Event Horizon Telescope Collaboration et al. 2022c, 2024b), and flaring states (e.g. Chatterjee et al. 2021; Scepi, Dexter & Begelman 2022). However, none of the EHTC ray-traced models fully satisfy all the constraints drawn from multiwavelength observations at 86 GHz, 230 GHz, 2.2 |$\mathrm{\mu m}$|⁠, and in the X-ray (Event Horizon Telescope Collaboration et al. 2022b, 2024a). The 230 GHz variability constraint is particularly stringent, as nearly all EHTC models in a MAD state exhibit greater variability (Event Horizon Telescope Collaboration et al. 2022c) than historical observations (Wielgus et al. 2022 and references therein). This limitation is possibly attributed to the modelled prescription of |$T_\mathrm{ i}/T_\mathrm{ e}$|⁠, which is based on the local plasma magnetization (so called |$R(\beta)$| prescriptions, e.g. Mościbrodzka et al. 2016; Anantua et al. 2020). In reality, |$T_\mathrm{ e}$| is fundamentally influenced by microphysical plasma and radiation interactions, and does not depend trivially on |$T_i$|⁠. A first-principles kinetic approach is required to model these collisionless effects (Parfrey et al. 2019; Crinquand et al. 2022; Galishnikova et al. 2023a).

Our investigation into 2T thermodynamics within MAD GRMHD simulations addresses some of these limitations by evolving both electron and ion temperatures (Ressler et al. 2015; Sądowski et al. 2017). We do not resolve the actual heating mechanisms; instead, we assume that heating at the grid scale occurs through sub-grid magnetic reconnection (Rowan et al. 2017). We do not model non-thermal electron distributions which are likely non-negligible (Moscibrodzka 2024). We perform an analysis of variability that systematically adds one layer of physics at a time (see Table 3). By incorporating the 2T treatment, variable adiabatic index, and radiative cooling of electrons, we achieve a closer match to historical 228 GHz variability compared to 1T simulations. We find an effective reduction of nearly 50 per cent in the values of the three-hour modulation index (⁠|$M_3$|⁠) distribution. Additionally, we find that |$M_3$| increases with frequency and does not significantly depend on the angle between the observer’s line of sight and the angular momentum vector of the accretion disc.

In a 2T GRMHD simulation, we extract two light curve predictions: one using the evolved |$T_\mathrm{ e}$| directly from the simulation, and another where the evolved |$T_\mathrm{ e}$| is ignored and instead calculated in post-processing with the |$R(\beta)$| prescription. We find that |$M_3$| is lower when using the evolved |$T_\mathrm{ e}$|⁠, suggesting that the |$T_\mathrm{ e}$| derived from the |$R(\beta)$| prescription exhibits greater fluctuations (see Fig. 7). This result is consistent with theoretical expectations for a 2T treatment, where the difference in adiabatic indices between relativistic electrons and non-relativistic ions causes an effective reduction in the fluctuations of the electron temperature (see Section 3.4).

Even for the very low accretion rates relevant for Sgr A|$^{\ast }$| (⁠|$f_{\mathrm{edd}} \approx 10^{-7}$|⁠), we find that radiative cooling of electrons – via synchrotron, inverse Compton, and bremsstrahlung processes – still affects the accretion flow, lowering |$T_\mathrm{ e}$| in the inner |$50r_\mathrm{ g}$| accretion disc (see Fig. 5), reducing average 228 GHz synchrotron flux, and shifting the |$M_3$| distribution towards lower variability values by roughly 10 per cent. We find that synchrotron emission dominates over inverse Compton, while bremsstrahlung is negligible, consistent with Yoon et al. (2020). Our estimated radiative efficiency, |$\eta _{\mathrm{rad}} \approx [0.02\!-\!0.05]$|⁠, is broadly consistent with the value |$\eta _{\mathrm{rad}} \approx 0.03$| obtained from GRRMHD MAD simulations that include radiative transport (Liska et al. 2024). These results challenge the paradigm in ADAF models, with |$\eta _{\mathrm{rad}}\lt 0.001$| across all radii. MAD states have higher radiative efficiency, while in SANE states, cooling is expected to be less significant especially at low accretion rates. Therefore, electron radiative cooling is non negligible if the accretion flow of Sgr A|$^{\ast }$| has a dynamically strong magnetic field, typical of a MAD state, as favoured by the current EHT observations.

Despite these improvements, further progress is required, as our 2T simulations still show more variability than historical observations of Sgr A|$^{\ast }$|⁠. For example, the more realistic stellar wind-fed accretion models better predict the submillimetre variability due to the comparatively lower levels of small-scale turbulence compared to SANE and MAD models (Murchikova et al. 2022). Additionally, our simulations do not resolve plasmoid-mediated magnetic reconnection, that requires higher resolution (Salas et al. 2024) or an effective resistivity or viscosity. Magnetic reconnection can cause the dissipation of magnetic energy into heat (Ripperda et al. 2019, 2022), thereby potentially influencing the radiative efficiency. Furthermore, ideal magnetohydrodynamics can never capture kinetic effects, like pressure anisotropy that can influence synchrotron emission and absorption (Galishnikova et al. 2023a; Galishnikova, Philippov & Quataert 2023b). An inherent uncertainty in our approach stems from not resolving the heating mechanisms and relying on a single heating prescription. In future work, we will explore heating prescriptions based on first-principles kinetic simulations in particular regions of the GRMHD domain, approximating a number of important subgrid effects of collisionless physics not captured by GRMHD simulations.

ACKNOWLEDGEMENTS

We thank Charles Gammie and Andrew Chael for stimulating discussions. We thank the referee for their valuable suggestions, which have improved the quality of this paper. LS thanks Martin Heemskerk for his valuable lessons on job arrays. LS and SM were supported by a Dutch Research Council (NWO) VICI award, grant no. 639.043.513 and by a European Research Council (ERC) Synergy Grant ‘BlackHolistic’ grant no. 101071643. In addition, LS was supported by a Colfuturo Scholarship, in partnership with the Colombian Ministry of Science. ML was supported by the John Harvard, ITC, and NASA Hubble Fellowship Program fellowships. KC was supported in part by grants from the Gordon and Betty Moore Foundation and the John Templeton Foundation to the Black Hole Initiative at Harvard University, and by NSF award OISE-1743747. GM was supported by a Canadian Institute of Theoretical Astrophysics (CITA) postdoctoral fellowship and by a Netherlands Research School for Astronomy (NOVA), Virtual Institute of Accretion (VIA) postdoctoral fellowship. GM acknowledges support from the Simons Collaboration on Extreme Electrodynamics of Compact Sources (SCEECS). OP acknowledges funding from VIA within NOVA. BR is supported by the Natural Sciences & Engineering Research Council of Canada (NSERC) and by a grant from the Simons Foundation (MP-SCMPS-00001470). Research at the Flatiron Institute is supported by the Simons Foundation. This research was enabled by using resources from Calcul Quebec (http://www.calculquebec.ca) and Compute Canada (http://www.computecanada.ca). This work used the Dutch national e-infrastructure with the support of the SURF Cooperative using grant no. EINF-3036, EINF-5383, and EINF-9222, which is (partly) financed by the Dutch Research Council (NWO), for post-processing of simulation data.

DATA AVAILABILITY

The simulation post-processed data used to plot the images in this work are available in Zenodo at http://doi.org/10.5281/zenodo.14793884.

Footnotes

REFERENCES

APPENDIX A: COULOMB COUPLING

The transfer of energy between ions and electrons is characterized by the Coulomb coupling rate in the comoving frame, which in units of |$\mathrm{erg\, cm^{-3}\, s^{-1}}$|is equivalent to,

with |$\theta _\mathrm{ m}=\left(1/\theta _\mathrm{ e}+1/\theta _\mathrm{ i} \right)^{-1}$| and |$m_\mathrm{ i}=m_\mathrm{ p}\left(X+4Y \right)$|⁠. |$K_i$| is the modified Bessel function of the ith order, |$\ln \Lambda \approx 20$| is the Coulomb logarithm, and |$\bar{n}=\left(X+Y \right)\rho /m_\mathrm{ p}$| is the number density (Stepney & Guilbert 1983; Sądowski et al. 2017). Ryan et al. (2017) suggested that Coulomb collisions become as important as viscous heating at |$f_{\mathrm{edd}} \approx 10^{-4}$|⁠. Dexter et al. (2021) and Liska et al. (2024) found Coulomb collisions are not important until |$f_{\mathrm{edd}} \approx 10^{-3}$|⁠. Conversely, semi-analytical models typically assume that Coulomb collisions dominate for |$f_{\mathrm{edd}} \gtrsim 10^{-2}$| (Esin et al. 1997).