Accretion Flow Morphology in Numerical Simulations of Black Holes from the ngEHT Model Library: The Impact of Radiation Physics

In the past few years, the Event Horizon Telescope (EHT) has provided the first-ever event horizon-scale images of the supermassive black holes (BHs) (M87*) and Sagittarius A$^*$ (Sgr A*). The next-generation EHT project is an extension of the EHT array that promises larger angular resolution and higher sensitivity to the dim, extended flux around the central ring-like structure, possibly connecting the accretion flow and the jet. The ngEHT Analysis Challenges aim to understand the science extractability from synthetic images and movies to inform the ngEHT array design and analysis algorithm development. In this work, we compare the accretion flow structure and dynamics in numerical fluid simulations that specifically target M87* and Sgr A*, and were used to construct the source models in the challenge set. We consider (1) a steady-state axisymmetric radiatively inefficient accretion flow model with a time-dependent shearing hotspot, (2) two time-dependent single fluid general relativistic magnetohydrodynamic (GRMHD) simulations from the H-AMR code, (3) a two-temperature GRMHD simulation from the BHAC code, and (4) a two-temperature radiative GRMHD simulation from the KORAL code. We find that the different models exhibit remarkably similar temporal and spatial properties, except for the electron temperature, since radiative losses substantially cool down electrons near the BH and the jet sheath, signaling the importance of radiative cooling even for slowly accreting BHs such as M87*. We restrict ourselves to standard torus accretion flows, and leave larger explorations of alternate accretion models to future work.


Introduction
With the advent of the Event Horizon Telescope ([EHT; [1,2]), imaging the near-horizon structure of supermassive black holes (BHs) is now a reality. The primary targets of the EHT and the future next-generation EHT (or ngEHT 1 ) are M87 * (the supermassive BH in the elliptical galaxy M87; [1]) and Sgr A * in the Galactic Center [2], which are two of the most well-studied low-luminosity active galactic nuclei. Extracting information about the event horizon-scale accretion flows in these two sources using the EHT's enormous resolving power is an active area of research. With the ngEHT, we will achieve unprecedented levels of angular resolution and sensitivity to low-flux regions, with the dynamic range in flux expected to increase to ∼1000 compared to the EHT's current dynamic range of ∼10 (e.g., [3]). This would enable us to investigate the BH shadow shape with higher precision as well as provide a crucial connection between the accretion flow and the jet launching region. The expected advances in sensitivity require deeper investigations of feature extraction from simulated synthetic reconstructions of BH systems. Hence, we designed the ngEHT analysis challenges 2 [4] to test our ability to capture the complex dynamics of gas and magnetic fields around M87 * and Sgr A * using the ngEHT reference array (e.g., [5]) with various analysis methods.
Black hole accretion and jet physics have been intensively studied over the past few decades (e.g., [6][7][8][9][10][11][12][13][14][15]). In the context of M87 * and Sgr A * , we expect the accretion flow to be highly sub-Eddington, radiatively inefficient, and geometrically thick, popularly known as radiatively inefficient accretion flows (RIAFs). This accretion flow solution has been used to successfully model the multiwavelength spectrum of Sgr A * (e.g., [16]). On the other hand, semi-analytical models of jets are preferred to explain the spectrum of M87 * (e.g., [17]). Thus, these two sources already provide a means to probe two different components of BH accretion, namely, the inner accretion flow structure and turbulence in Sgr A * and the prominent jet feature in M87 * . The first three EHT numerical simulation papers [18-20] already provide us with important clues about the horizon-scale conditions of these BH systems based on numerical simulations: (1) these BHs probably have non-zero spin, (2) the accretion disk is expected to have colder electrons than the jet sheath, and (3) the observations favor the presence of dynamically important magnetic fields close to the BH. All of these results point us toward the magnetically arrested disk (MAD; [10,21]) state, an accretion mode where the BH magnetosphere becomes over-saturated with magnetic flux and exhibits quasi-periodic explosions of vertical magnetic field bundles. MAD flows around spinning BHs also have powerful relativistic jets, where the jet power can exceed the input accretion power [13], which is a definite signature of BH spin energy extraction via the Blandford and Znajek [22] process.
Building on the semi-analytical RIAF models, the time-dependent general relativistic magneto-hydrodynamic (GRMHD) simulations have become important tools for deciphering BH accretion physics in a variety of astrophysical systems (e.g., [11,13,[23][24][25][26][27]). Indeed, the EHT regularly uses the libraries of GRMHD simulations to model the observed horizon-scale BH images of M87 * and Sgr A * as well as larger-scale jet images (such as for Centaurus A; [28]) in order to constrain the time-variable plasma properties. In designing the ngEHT reference array, it is therefore crucial to use GRMHD simulations for understanding the attainability of specific science goals, such as resolving the photon ring and the disk-jet connection region as well as tracing out time-variable features via the ngEHT analysis challenges.
In this work, we discuss the numerical fluid simulations that were used as source models for the ngEHT analysis challenges. In particular, our objective is to compare between the models that incorporate increasingly complicated levels of accretion and electron physics, focusing on M87 * and Sgr A * . Our model set consists of a time-dependent shearing hotspot stitched to a steady-state RIAF solution, two standard GRMHD simulations of MAD accretion flows, a GRMHD MAD simulation with electron heating via incorporating two-temperature physics, and a fully radiative, two-temperature, GRMHD MAD simulation. We describe the equations and setup of our numerical models in Section 2, show our comparison results in Section 3, and, finally, conclude in Section 4.

Numerical Simulations
In this section, we provide a brief description of the semi-analytical stationary RIAF and shearing hotspot model as well as the (two-temperature/radiative) GRMHD simulations used for the ngEHT analysis challenges.

RIAF + Hotspot Solutions
RIAF models attempt to describe the time and azimuthally averaged appearance of accretion flows. This is performed using a set of building blocks. The specific RIAF models used in the challenges are based on Yuan et al. [16], Broderick and Loeb [29], Broderick and McKinney [30] approach. We decompose the accretion flow into a set of phenomenological models that describe the electron density, temperature, magnetic field, and velocity profile. We have a cylindrical coordinate system x µ ≡ (R cyl , ϕ, z). The electron density profile is defined in terms of the cylindrical radius R cyl = r| sin(θ)| and vertical displacement z = r cos(θ), and is provided by where X denotes the population of electrons. The disk height h is set to unity for this work. For the challenge dataset, we included both thermal synchrotron emitting (X ≡ th), and non-thermal synchrotron (X ≡ nth) emitting electrons. The thermal electrons have n e,th,0 = 1.3 × 10 8 and p th = −1.1, while the non-thermal electrons are provided by n e,nth,0 = 1.3 × 10 5 , p nth = −2.02. These numbers are from Tiede et al. [31] and are set to match the best-fit parameters for the Sgr A * spectrum from Broderick et al. [32] and Broderick et al. [33].
The temperature profile of the thermal electrons is also provided by a radial power law with a Gaussian envelope describing the height: where, for the challenge data, we set T e,0 = 6.3 × 10 10 K. Following Pu et al. [34], we define the gas pressure by assuming that the ultra-relavistic protons are in roughly virial equilibrium with the gravitational force providing: where m p , G, and M BH are the proton mass, gravitational constant, and the black hole mass. For the local magnetic field, we then assume a constant β = p gas /p mag = 10 plasma, which, combined with Equation (3), provide us the magnetic field strength. The orientation of the magnetic field is then provided by a purely toroidal configuration, relative to the plasma observer. Finally, we extract the emission and absorption coefficients from the synchrotron self-absorption model in Broderick and Blandford [35] for both the thermal and non-thermal synchrotron electrons. For the non-thermal electrons, we use a power law prescription with radiation coefficients from Jones and O'Dell [36], and a photon spectral power law index of 1.25 (see [31] for more details).
We follow the velocity field prescription from Pu et al. [34] to describe the accretion flow dynamics. Using the notation from Tiede et al. [31], the velocity field u µ = u t v µ is provided by where v µ K denotes the Keplerian velocity field, v µ f f is the free-fall velocity field, and α vel = 0.01 and κ vel = 1.0, i.e., the hotspot rotates with Keplerian angular velocity and is weakly free-falling. For α vel = 0, κ vel = 1, we have a Keplerian orbit, whereas, for α vel = 1, κ vel = 0, we have free-fall. The remaining component u t is provided by the normalization condition u µ u µ = −1. The radial component outside the inner stable circular orbit (ISCO) is v r K = 0 as the disk is in steady-state. However, inside the ISCO, we use plunging geodesics that are specified by matching the angular momentum and energy at the ISCO.
The hotspot evolution follows the model from Tiede et al. [31] (also see [29]), where we assume that the hotspot travels passively along the accretion flow velocity field Equation (4). This implies that the equation of motion is provided by the conservation of particle number Equation (6). For the emission, we assume a non-thermal synchrotron hotspot, with an initial Gaussian density profile n e (x µ ) = n 0 e −∆r µ r µ +(∆r where we have set n 0 = 6 × 10 6 , hotspot size R s = 0.5 r g = 0.5GM BH /c 2 , and ∆r µ is the displacement from the hotspot center.

GRMHD Simulations
Over the previous two decades, multiple GRMHD codes have been developed and utilized to model black hole accretion and jet launching physics over long dynamical timescales. The wide usage of GRMHD simulations is particularly encouraging since this allows for verification of code-specific numerical choices that users usually have to make even while solving the same base set of GRMHD equations. Indeed, recently, there was a community-wide effort to benchmark these codes against each other for a standard problem: evolving a weakly magnetized torus of gas around a spinning black hole [26]. It was found that these codes largely provide similar solutions, though some disk quantities remain unconverged with increasing grid resolutions, suggesting more investigation is required. For this work, we employ three different GRMHD codes to probe black hole accretion, increasing the complexity of the equations solved at each step: (1) single fluid GRMHD simulations from the H-AMR code [37], (2) a two-temperature single fluid GRMHD simulation from the BHAC code [38], and (3) a two-temperature radiative GRMHD simulation from the KORAL code [39].
First, we describe the set of GRMHD Equations (e.g., from [23,26]). We have the conservation of particle number and energy-momentum: Here, ρ is the rest-mass gas density and can also be written in the form of ρ = mn, where m is the mean rest-mass per particle and n is the particle number density. We also have the four-velocity u µ , stress-energy tensor T µ ν , metric determinant g ≡ det(g µν ), and the metric connection Γ β να . Note that the index t refers to the temporal component of the vector or tensor and i denotes the spatial indices.
The stress-energy tensor T µ ν is provided as where U gas and p gas are the gas internal energy and pressure, related by the ideal gas equation: p gas = (Γ gas − 1)U gas assuming a gas adiabatic index Γ gas . We also have the magnetic pressure p mag = b 2 /2 and the magnetic field 4-vector b µ , which can be defined in terms of the magnetic field 3-vector B i : Here, we included a factor of √ 4π into the definition of B i . We evolve the magnetic field B i using the spatial components of the induction equation, while the temporal component provides the no monopoles constraint, These equations are numerically integrated in the conservative form [23] to obtain the physically relevant quantities ρ, U gas , u µ , and B i . We refer the reader to the corresponding code papers for more information on the numerical techniques used to evolve these equations over space and time.
In this work, we use two GRMHD simulations performed with the H-AMR code, one targeting M87 * and the other Sgr A * . These simulations employ logarithmic Kerr-Schild coordinates and the grid resolutions are N r × N θ × N ϕ = 580 × 288 × 512 for the M87 * simulation and 348 × 192 × 192 for the Sgr A * simulation. All the simulations in this work adopt the geometrical unit convention, G = c = 1 and using M BH = 1, normalizing the length scale to the gravitational radius r g = GM BH /c 2 . The M87 * GRMHD simulation evolves a MAD flow around a black hole with spin a = 0.9375. The Sgr A * model also simulates a MAD flow but around a black hole with spin a = 1/2. H-AMR uses outflowing radial boundary conditions (BCs), transmissive polar BCs, and periodic azimuthal BCs (for more details, see [40]).
Since GRMHD simulations are scale-free, we determine the gas density code-to-CGS units conversion factor (hereafter, "density scaling") by raytracing the simulation at 230 GHz for a target source and flux. We use the general relativistic ray-tracing codes BHOSS [41] and IPOLE [42] to compute images at 230 GHz and set the compact flux to be approximately 0.5 Jy for M87 * [1] and 2.4 Jy for Sgr A * [2]. We use the black hole masses and distances of M BH = 6.2 × 10 9 M and D BH = 16.9 Mpc for M87 * (e.g., [18,43] and references therein) and M BH = 4.14 × 10 6 M and D BH = 8.127 kpc for Sgr A * [20, 44,45].
GRMHD simulations evolve a single temperature fluid. At the low accretion rates seen in M87 * and Sgr A * , Coulomb coupling between ions and electrons is inefficient and, therefore, the two particle species are not in thermal equilibrium. Since ions are much heavier than electrons, the ions dominate the single fluid thermodynamics evolved in GRMHD simulations. Hence, to calculate the radiative output from GRMHD simulations, we calculate the electron temperature T e using sub-grid models such as the R − β prescription [46] based on local gas plasma-β (≡ p gas /p mag ): For the ngEHT analysis challenges, we select R high values of 160 and 40 and source inclinations of 163 • and 50 • for the M87 * and Sgr A * simulations, respectively. We assume a thermal relativistic Maxwell-Jüttner distribution for describing the electron energy distribution in the Sgr A * model, and a hybrid thermal+non-thermal κ−distribution (e.g., [47,48]) for the M87 * model. The model images are shown in Roelofs et al. [4].

Two-Temperature Physics
Two-temperature GRMHD (or 2t-GRMHD) simulations (e.g., [49,50]) evolve the ion and electron entropy equation separately and, hence, provide the ion and electron temperature in a self-consistent manner. The main advantage of this method is that we remove the electron temperature as a free parameter when constructing images. However, we do have to make a choice about the sub-grid prescription that determines the fraction of local dissipative energy that heats the electrons. There are two heating mechanisms that are thought to be applicable to global accretion flows: turbulent heating [51,52], and magnetic reconnection [53,54]. For the ngEHT analysis challenges, we focus only on one simulation with reconnection heating, from Mizuno et al. [55], as heating via magnetic reconnection in equatorial current sheets formed in magnetically arrested flows, captured by the Rowan et al. [54] heating model, is arguably more important than heating via small-scale turbulent eddies that are more prevalent in weakly magnetized disks.
We assume that the number densities and velocities of ions and electrons are equal, i.e., n i = n e = n and u µ i = u µ e = u µ , thus maintaining charge neutrality. The electron entropy equation is provided as where the electron entropy is s e = p e /ρ Γ e , with p e and Γ e as the electron pressure and adiabatic index. The total heating rate Q is calculated by comparing the total internal energy of the gas and the internal energy obtained from the electron entropy conservation Equation (see [49] for more details). The fraction of dissipative heating that contributes to electron heating is provided by f e . For this particular simulation, f e is designed to capture electron/ion heating via magnetic reconnection from Rowan et al. [54]: where β max = 1/4σ h , defined using the hot gas magnetization σ h = b 2 /ρE th and the specific gas enthalpy E th = 1 + Γ gas p g /[ρ(Γ gas − 1)]. The 2t-GRMHD simulation from Mizuno et al. [55] assumes modified Kerr-Schild coordinates and a black hole spin of 0.9375. The grid resolution is 384 × 192 × 192. The accretion mode is of a magnetically arrested flow and the simulation is raytraced (using BHOSS) once the near-horizon flow has reached steady state. The target source is M87 * , assuming a black hole mass of M BH = 6.5 × 10 9 M and distance of 16.9 Mpc [1]. The accretion rate is normalized such that the 230 GHz compact flux density is 0.8 Jy. We assume a thermal electron distribution everywhere except in the jet sheath where we adopt a κ−distribution. More details about the image are provided in Roelofs et al. [4].

Radiative GRMHD
Two temperature GRMHD simulations do not include radiative cooling and, hence, are thought to be appropriate for low luminosity supermassive black holes such as M87 * and Sgr A * . To verify this assumption, we consider a two-temperature radiative GRMHD (2t-GRRMHD hereafter) simulations from Chael et al. [25]. This simulation accounts for selfconsistent radiation physics, incorporating both particle heating via magnetic reconnection (as in Section 2.3) and radiative cooling via bremsstrahlung, synchrotron, Compton, and Coulomb losses. The simulation is run using the 2t-GRRMHD code KORAL [39,56,57], which evolves a two-temperature magnetized fluid and treats the radiation field as a second fluid [58]. The conservation equations solved in 2t-GRRMHD are different from that of GRMHD: where R µ ν is the frequency-integrated radiation field, defined as, Here, the radiation field is described by its rest-frame energy density E and four-velocity u µ R following the M1 closure scheme. The ion and electron entropy equations are where q v is the dissipative heating rate and q C is the Coulomb coupling rate that captures the exchange of energy between ions and electrons. The heating fraction f e is from Rowan et al. [54] (see Chael et al. [57] for more details), same as the 2t-GRMHD simulation. Finally, G is the radiative cooling rate [39]. For further details about the equations, see Sądowski et al. [39], Chael et al. [57], Sądowski et al. [58]. The simulation assumes a black hole spin of a = 0.9375 and mass M BH = 6.2 × 10 9 M , targeting M87 * . The gas density is scaled to physical CGS units such that the compact emission at 230 GHz is roughly 0.98 Jy [59,60]. The simulation uses modified Kerr-Schild coordinates with a grid resolution of N r × N θ × N ϕ = 288 × 224 × 128. See Chael et al. [25] for more details about the simulation, while not utilized for the ngEHT analysis challenges, we included this simulation in this work since this model captures the coupling between gas and radiation, necessary for black holes accreting close the Eddington limit. Further, this model has been used in previous ngEHT reference array papers [3,5].

Results
We perform a series of comparisons focused on the time-evolution of horizon-scale quantities and radial dependence of disk and jet properties. The diagnostics are chosen such that any trends we find can inform EHT/ngEHT science applications, such as horizon-scale morphology and variability of the accretion flow. Further, the quantities are similar to those reported in the GRMHD code comparison project [26] and, thus, can be directly compared. There is a total of five models: three (2t-radiative) GRMHD simulations targeting M87 *, and one RIAF solution and one GRMHD simulation for Sgr A *. We further note that all three numerical simulations of M87 * have the same BH spin, favoring direct comparisons of the horizon-scale gas properties.

Temporal Behavior of Horizon Fluxes
We calculate the mass, magnetic, angular momentum and energy fluxes in the radial direction as follows: Ang.Mom. : where all quantities are calculated at the event horizon radius r hor = r g (1 + √ 1 − a 2 ). We note that there could be a substantial contribution of density floors when calculating the mass accretion rate for MAD systems. However, this radius was chosen for simplicity when comparing to previous simulations in the literature. Figure 1 shows the mass accretion rateṀ in units of solar masses per year (M /yr), the dimensionless magnetic flux φ = Φ/ Ṁ r 2 g c, the outflow power P out =Ṁc 2 −Ė, and the specific angular momentum fluxJ/Ṁ for simulations targeting M87 * and Sgr A *. The RIAF solution being a steady-state solution is excluded from this section (though the hotspot evolves with time). The quantities from the 2t-GRRMHD simulation are only shown for (11 − 16) × 10 3 r g /c, i.e., the time period over which the simulation was raytraced in Chael et al. [25]. Remarkably, despite the difference in electron the physics complexity, the simulations behave very similarly. The factor of 2 difference inṀ between the M87 * non-radiative simulations and the 2t-GRRMHD simulation can be explained by the lower electron temperatures in the near-horizon accretion flow due to radiative cooling (see Section 3.2) as well as the higher 230 GHz flux normalization used for the radiative model.
The accretion rate in all simulations show large variation with quasi-periodic sharp drops. These drops inṀ occur due to the emergence of magnetic flux eruptions, a characteristic feature of the magnetically arrested disks [61][62][63][64]. These eruptions also lower the value of φ since magnetic flux bundles escape from the vicinity of the BH, carrying away the magnetic flux accumulated in the BH magnetosphere. We see that φ often crosses the magnetic flux saturation value of 50 [13], overwhelming the BH magnetosphere with strong magnetic fields that eventually reconnect and trigger flux eruptions (see [63] for the detailed mechanism). Figure 2 shows a series of equatorial snapshots from the M87 * 2t-GRMHD simulation where we follow a particular magnetic flux eruption by the rotating gas that is hot (T e 10 12 K), has low density (ρ < −2 in code units), and primarily consists of vertical field lines. As these field line bundles move out and interact with the disk, they (1) hinder accretion, lowerinġ M, (2) remove magnetic flux from near the BH, lowering the jet power, and (3) push gas outwards, reducing the inward angular momentum flux. Curiously, we see larger drops in the specific angular momentum flux for the Sgr A * GRMHD model. This is possibly due to the smaller BH spin (a = 0.5 as opposed to 0.9375 for the M87 * models) as the weakly powered jet does not carry away angular momentum as efficiently as the higher BH spin models and flux eruptions play a bigger role in regulating disk angular momentum transport. Additionally, the reconnection events that trigger these eruptions accelerate electrons to higher energies, and are, thus, crucial for understanding flare activity in BH sources.
Values are calculated at the event horizon. Table 1. GRMHD simulations considered in this work. Simulation grid resolution, simulation time period over which raytracing was performed, modulation index (MI) of the mass accretion rateṀ, dimensionless magnetic flux φ, outflow efficiency P out /Ṁc 2 , and the specific angular momentum fluxJ/Ṁ for each GRMHD model. The MI is calculated over the final 5000 r g /c in runtime and at the event horizon (see Figure 1). To quantify the time-variability of the horizon-fluxes, we calculate the modulation index MI, which is defined as the ratio of the standard deviation and the mean of the quantity over time [20]. We show the MI for the different fluxes in Table 1. The MI(Ṁ) is usually a good proxy for the variability of the sub-millimeter (sub-mm) emission in these slowly accreting optically thin black hole sources (e.g., [65]). The MI(Ṁ) values we see from the simulations are ∼0.23-0.29 and are larger than expected from Sgr A * 230 GHz lightcurves (where MI ∼ 0.1; [66]). This suggests that careful analysis of the electron distribution function is needed to understand if we are substantially over-predicting the 230 GHz lightcurve variability. Further, in general, weakly magnetized accretion flows exhibit lower MI(Ṁ) values due to the absence of flux eruptions, which suggests that further study of the accretion mode in Sgr A * is also necessary. It is encouraging to note that our MI values forṀ and φ are consistent with the MI values from longer time-evolved GRMHD simulations of a = 0.9 BHs in Narayan et al. [27], indicating that our simulations are sufficiently converged with respect to horizon-scale quantities. Magnetic flux eruptions, a characteristic feature of magnetically arrested disks, eject bundles of vertical magnetic fields filled with relativistically hot, low-density plasma. These flux bundles spiral around the black hole and may be responsible for high-energy flares [63]. The time unit M is equivalent to r g /c.

Disk-Averaged Quantities
Here, we calculate the disk-averaged properties of each model, namely gas density ρ, thermal pressure p gas , magnetic pressure p mag , radial velocity |v r |, azimuthal velocity |v ϕ |, angular momentum u ϕ , disk scale height h/r, ion temperature T i , and the electron temperature T e . We define disk-averaging of a quantity q as where q ∈ {ρ, p gas , p mag , |v r |, |v ϕ |, u ϕ , h/r, T i [Kelvin], T e [Kelvin]}. Further definitions follow: where Γ ad and u are the adiabatic index and the internal energy of the gas.  The density profiles in the inner few tens of r g converge roughly to a n e ∝ r −1 profile, matching the RIAF density profile as well as the longer time-evolved MAD simulations [64]. The M87 * 2t-GRRMHD density is larger by a factor of ≈ 2 from the GRMHD/2t-GRMHD, as is expected from the difference in the mass accretion rate (Figure 1). The 2t-GRRMHD simulation exhibits a slightly more magnetized inflow within the inner 2 r g , but, overall, the GRMHD simulations have a similar plasma-β ≡ p gas /p mag disk profile. The stronger magnetic field seen in the 2t-GRRMHD model could explain the higher values of the horizon magnetic flux seen in Figure 1. The RIAF model assumes a constant disk plasmaβ = 10, (see Section 2.1), which is substantially higher when compared to the MAD GRMHD models. This value of plasma-β is chosen in order to match the observed 230 GHz flux density of Sgr A * . As we see from the disk scale height in Figure 4, the RIAF model has a much thicker disk than the GRMHD models, and, therefore, produces a lot more sub-mm emission even with a low electron temperature and weak magnetic field strength. Next, we see that the disk-averaged electron temperature T e in the 2t-GRRMHD M87 * model is more than an order of magnitude lower than the other GRMHD models within the inner 10 r g , but actually matches the Sgr A * RIAF T e profile and has a shallower profile T e ∝ r −1 instead of r −3/2 . It is also interesting to note that the disk ion temperatures T i are very similar in all the GRMHD simulations shown here. Therefore, despite the same reconnection-driven heating mechanism captured in both the 2t-GRMHD and the 2t-GRRMHD models, the radiative cooling of hot electrons plays a crucial role in determining the eventual T e . Due to the low T e , the required accretion rate normalization is higher in the 2t-GRRMHD, as we noted in the previous subsection.
In Figure 4, we show the average disk scale height h/r, the radial and angular velocities (v r and Ω and the specific angular momentum u ϕ ). The MAD simulations all show very similar disk properties. The h/r ≈ 0.1-0.3 with a sharp increase within 3 r g where the inflow becomes vertically supported by strong poloidal magnetic fields. The radial velocity has a profile of r −1/2 , similar to the scaling relation found in ADAF solutions assuming a constant viscosity parameter α [67,68]. The α parameter profile depends on how magnetized the accretion flow is, with α ∝ r −1 for the weakly magnetized flows and close to constant for the MAD-like flows (e.g., [64,69]). We also see highly sub-Keplerian angular velocity profiles in the GRMHD models, typical for magnetically supported disks. For the RIAF model, the RIAF disk is not infalling and has a constant Keplerian angular velocity. Instead, the hotspot, added to the RIAF solution, undergoes shearing and disappears into the BH with a radial velocity similar to the values found in the GRMHD MAD disks. This occurs because the hotspot is designed to travel along plunging geodesics (see Section 2.1), similar to rapid gas infall close to the BH in the GRMHD models. The angular momentum in the GRMHD models looks sub-Keplerian as expected for MADs. Figure 5. We show the jet radius R jet and the jet Lorentz factor γ from the M87 * GRMHD and 2t-GRRMHD models, and the Sgr A * GRMHD model. The gray circles indicates the deprojected jet radius of the M87 jet assuming a BH mass of 6.2 × 10 9 M and a source inclination of 14 • [70]. The data points are a compilation of various papers [59,60,71,[76][77][78].

Jet Properties
Here, we calculate the radial profiles of the jet half width R jet and Lorentz factor γ: where µ = −T r t /(ρu r ) is the specific radial energy density and α = 1/ −g tt . We calculate the average jet radius R jet by integrating over the total surface area covered by the two jets over a shell at a particular radius, and assume that the total surface area is approximately 2πR 2 jet . We define the jet boundary as µ > 2, i.e., the region over which the jet still remains highly energized. This definition of the jet boundary is quite similar to the condition σ = 1, which is widely used in the literature (e.g., [27]). Since µ = γ(σ + E th ), our condition µ > 2 also incorporates regions where the jet might not be magnetically dominated but is relativistically hot or fast. Since we restrict our jet profiles to within r 10 3 r g , the jet radius is primarily determined by the jet magnetization. 3 Figure 5 shows the jet radius R jet and Lorentz factor γ as a function of the radial distance from the BH for the M87 * GRMHD and 2t-GRRMHD models as well as the Sgr A * GRMHD model. The M87 * jet radius from our models matches the observed jet width from M87 (gray circles) quite well, with the radial profile roughly proportional to r 0.625 , which is the fitted powerlaw for the M87 jet [70], though the index value has also been reported to be slightly smaller in some works (0.57; [71,72]). The power law index of 0.625 is larger than that found using the σ = 1 condition from Narayan et al. [27], where the authors found a power law index of 0.428 for their MAD spin a = 0.9 GRMHD model. It is possible that we find larger jet radii as we incorporate a part of the hot jet sheath region within our definition of R jet (as suggested by Figure 7 in [73]). For the Sgr A * model, we also find a similar R jet profile. While there are no detections of an extended jet in Sgr A * (e.g., [74]), semi-analytical and GRMHD models of Sgr A * largely favor a jet component from a spinning BH (e.g., [20,75]).
We also show the Lorentz factor γ in Figure 5. Mostly, the jets accelerate to γ ≈ 3-4 by 10 3 r g in all of our GRMHD models. It is difficult to compare our γ profiles with values inferred from observations of the M87 jet (e.g., [79]) since our γ values are biased toward the jet spine while the observations generally capture the velocities of localized features in the sub-relativistic jet sheath/disk wind, especially at small distances from the BH. Indeed, both simulations and observations show that the jet Lorentz factor varies greatly as a function of jet radius (e.g., see [73]). We speculate that a better approach might be to calculate emissivityweighted Lorentz factors in order to compare to the measured γ from M87. Since our focus is on the comparison between GRMHD simulations, we leave direct comparisons to observational data to future work.

Axisymmetrized Profiles
In the previous sections, we found that the largest differences between the GRMHD models occur in electron temperature distribution. Figure 6 shows the time and azimuthally averaged 2D vertical plots of gas density n e and electron temperature T e . We show the normalized n e so as to capture the relative change in the disk/wind density distribution, which would provide us information about the disk structure. The large difference in disk scale height is immediately apparent between the RIAF and the MAD GRMHD models (also see Figure 4). The presence of a prominent wide jet component in MADs squeezes the inner disk and pushes against the accretion flow, a feature which is not captured in the constant h/r RIAF model. However, the RIAF model roughly reproduces the density profile of the disk midplane region, suggesting that the RIAF model could represent non/weakly jetted, quasi-spherical accretion flows quite well. For sources such as M87 * , where we see a prominent jet component, the density gradient in the vertical direction is expected to be steeper as strong magnetic stresses power winds carry gas away from the disk (e.g., [64]).
Overall, the disk/wind density distribution among the GRMHD models look similar with small differences in the lateral extension of the wind region and the steepness of the vertical gradient in density. For example, if we compare the 2t-GRRMHD model with the other two simulations, the density in the wind region is larger in the radiative model. The reason for the shallow vertical density profile in the 2t-GRRMHD model is unclear since the weakly magnetized thick disk simulations tell us that radiative cooling would lead to the loss of gas pressure in the disk and would result in the disk collapsing to a relatively dense structure in the midplane (e.g., [80,81]). However, in the presence of strong poloidal magnetic fields, i.e., in the MAD state, the plasma-β decreases to β ≈ 0.2-1 in the disk midplane (see Figure 3, third row, left panel), and can reach even lower values in the upper layers of the accretion flow. The high magnetic pressure could help support the disk against collapse while sufficiently strong magnetic stresses could power disk winds. Such behavior is also seen in recent GRMHD simulations of near-Eddington, geometrically thin, strongly magnetized disks, where the inner disk (or corona) has a larger h/r than the outer disk due to magnetic pressure support [82][83][84]. To verify how radiative cooling affects the inner disk/wind structure in highly sub-Eddington accretion flows such as M87 * and Sgr A * , we require longer 2t-GRRMHD simulations such that the disk is in an inflow-outflow equilibrium out to a radius of at least 50 r g . Figure 6. We show t-and ϕ-averaged data: electron number density n e (top row) and temperature T e (bottom row). We also denote the jet boundary with σ = 1 (black lines). The time-averaging is performed over the 5000 r g /c for each model. RIAF plots are for Sgr A * while the rest are for M87. The Sgr A * GRMHD model produces similar plots of n e and T e as the M87 * model, and, hence, we do not show it here.
The 2D temperature plot of the RIAF model also looks vastly different in the inner disk (r 20 r g ) when compared to the GRMHD and 2t-GRMHD simulations, but is similar to the temperature distribution in the 2t-GRRMHD disk midplane (also seen in the T e plot of Figure 3). The RIAF model does not capture gas heating in the jet sheath region (the region just outside of the jet boundary indicated by the σ = 1 dashed line) and, therefore, T e drops as we move away from the midplane toward the poles. In the GRMHD models, the jet sheath is as hot, if not hotter, than the inner accretion flow as temperatures reach T e > 10 11 K. For the GRMHD simulation, the electron temperature is provided as a fraction of the fluid temperature, where the fraction depends on how magnetized the gas is in the region, as per the R − β prescription from Equation (14). For the M87 * model, we chose a R high value of 160 to have a jet-dominated sub-mm image. This choice of R high suppresses the electron temperature in the disk, focusing higher temperatures in the jet sheath. Comparing the GRMHD model with the 2t-GRMHD model, the jet sheath region exhibits very similar T e values, but the disk midplane is hotter by a factor of a few in the 2t-GRMHD model. We note that this difference in T e in the midplane is more noticeable in the 2D plot rather than in the disk-averaged T e profile shown in Figure 3 as the upper layers of the disk become substantially hotter in the GRMHD model.
For the radiative 2t-GRRMHD model, the inner regions of the disk are cooler as electrons heated by magnetic reconnection quickly cool via synchrotron and Compton losses. From Figure 3, the drop in T e for the 2t-GRRMHD model is shown to be as large as an order of magnitude when compared to the (2t-) GRMHD models. Another interesting feature is that the hot region (T e > 10 11 K) in the jet sheath is much narrower in the 2t-GRRMHD model, which could have a significant bearing on the ray-traced image, possibly producing a thinner jet sheath. Finally, the difference in T e in the jet body between the GRMHD models is due to the different density/internal energy floor setups used by the corresponding codes. Since the gas in the jet sheath and the jet body undergo mixing due to boundary instabilities (e.g., [73,85]), it is possible that the choice of floors could affect the overall electron temperature in the jet sheath. Such a study is outside the scope of our paper and is left to future work.

Orbiting Hotspot in a RIAF Model
High-energy flares are commonly observed in AGNs, with GeV and MeV flares seen in M87 * (e.g., [86,87]) and quasi-daily nIR and X-ray flares in Sgr A * (e.g., [44,[88][89][90][91][92][93][94]). A number of attempts have been performed to explain the origin of flaring, such as magnetic reconnection in turbulent gas flows in the disk and the jet [65,73,[95][96][97] and magnetic flux eruptions [61,63,98,99]. For Sgr A * , the semi-analytical models found that we require highenergy electrons, assumed to be accelerated via an ad hoc process such as a large magnetic reconnection event or shocks, to describe the large flaring events [100][101][102]. The near-infrared observations from the GRAVITY telescope provided further evidence for orbiting hotspot features in the accretion flow [103] that may be linked to acceleration events. It has also been recently shown that orbiting hotspots can be used to model double-peaked X-ray flares [94,104] and prominent Stokes Q-U loops in the sub-mm emission of Sgr A * [105]. These results provide us with considerable motivation to test the capability of the ngEHT to detect hotspot formation in accretion flows around black holes.
Instead of isolating a particular magnetic flux eruption event in our simulations, we added a shearing hotspot to the RIAF solution as detailed in Section 2.1. Figure 7 shows the temporal evolution of the azimuthally averaged electron number density of the hotspot. We begin with a Gaussian distribution of gas that undergoes shearing as the gas falls in closer to the BH. The overall density normalization is much lower than in the RIAF disk since the optically thin hotspot gas produces a large enough non-thermal synchrotron emissivity. The hotspot is evolved over 800 r g /c, but the gas distribution comes to a near-steady-state profile within the first 200 r g /c, which is roughly one hour for Sgr A * . The shearing of the hotspot gas has a significant impact on the evolution of the 230 GHz image [4,31]. From Figure 4 (right column), we see that the radial velocity matches the disk-averaged gas velocity from the GRMHD model, showing nearly free-fall speeds, while the azimuthal velocity becomes highly sub-Keplerian. The velocity profiles show that our hotspot model should be able to reproduce the expected hotspot motion from the GRMHD models, and is ideal for investigating multiwavelength flare lightcurves. Emami et al. [106], a companion paper, goes into further details about how current dynamical reconstruction techniques can be used to trace out the motion and morphology of the shearing hotspot in the context of ngEHT observations. These hotspot models and reconstruction methods would be integral in deciphering the more complex gas dynamics of magnetic flux eruption events in MADs, which have been shown to produce a significant variation in the image structure of M87 * at 230 GHz (e.g., [107]).

Conclusions
In this work, for the first time, we compare a series of numerical solutions with increasing complexity that were specifically constructed to understand the accretion flow around the supermassive black holes M87 * and Sgr A * . We include a time-independent radiatively inefficient accretion flow (RIAF) model as well as fully 3D GRMHD simulations of accreting black holes, incorporating the effects of electron heating and cooling losses via two-temperature and radiation physics. In addition, each of our simulations are run with different GRMHD codes, which is similar to the approach of another community-wide code comparison effort [26]. We found that the simulations exhibit remarkably similar properties given that the simulations incorporate varying levels of complexity in electron physics. The notable exception is the electron temperature, where radiative cooling decreases the temperature by a factor of 10 within the inner 10 gravitational radii, the region that produces the bulk of the 230 GHz emission in M87 * , one of the two primary targets of the EHT and the ngEHT (the other being Sgr A * ). The main goal of this work is to understand the variation in the underlying accretion flow and jet properties in our models since synthetic ray-traced images constructed from these models are used as "truth" images of M87 * and Sgr A * for the ngEHT Analysis Challenges [4]. The ngEHT Analysis Challenges are an effort to determine how much information about the accretion flow and jet dynamics we can glean from the proposed ngEHT reference array, and what modifications in the image reconstruction tools are necessarily required to decode future ngEHT observational data.
Our paper deals with numerical models designed to investigate hotspot evolution, turbulent inspiralling gas flows, and extended powerful jets, targeting M87 * and Sgr A * . We restricted our model set to the community-standard setup: a rotating, geometrically thick, optically thin torus of magnetized gas around a spinning black hole, which is the fiducial model choice of the EHT [18][19][20]. This model choice leaves out the exploration of multiple new setups of black hole accretion, such as quasi-spherical wind-fed inflows (e.g., [108,109]), strongly wind-fed accretion (e.g., [110,111]), geometrically thin accretion disks (e.g., [83,84,112]), puffy radiation-dominated super-Eddington disks (e.g., [113,114]), and misaligned accretion disks (e.g., [24,40,115,116]). Apart from varying the accretion mode, the high resolution of the images from EHT and ngEHT could potentially help distinguish between different space-time metrics [117]. To date, only a limited number of non-Kerr GRMHD simulations have only been performed (e.g., [118][119][120]). The future of numerical studies is bright, given their rising popularity in the astrophysics community and the increase in computational resources. The breadth of the current investigations in accretion physics would result in a plethora of variable structures that should be thoroughly studied keeping the observational capabilities of the ngEHT in mind.  3 Specific enthalpy includes the rest-mass energy contribution in our definition from Section 2.3.