A Study of the Accretion State of Magnetically Arrested Disks across Black Hole Spins for Radiatively Inefficient Accretion Flows

The study of magnetically arrested disks (MAD) has attracted strong interest in recent years because these disk configurations were found to generate strong jets, as observed in many accreting systems. Here, we present the results of 14 general relativistic magnetohydrodynamic simulations of advection-dominated accretion flow in the MAD state across black hole (BH) spins, carried out with cuHARM. Our main findings are as follows. (i) The jets transport a significant amount of angular momentum to infinity in the form of Maxwell stresses. For positive, high spin, the rate of angular momentum transport is about five times higher than for negative spin. This contribution is nearly absent for a nonrotating BH. (ii) The mass accretion rate and the MAD parameter, both calculated at the horizon, are not correlated. However, their time derivatives are anticorrelated for every spin. (iii) For zero spin, the contribution of the toroidal component of the magnetic field to the magnetic pressure is negligible, while for a fast-spinning BH, it is on the same order as the contribution of the radial magnetic component. For high positive spin, the toroidal component even dominates. (iv) For negative spins, the jets are narrower than their positive-spin counterparts, while their fluctuations are stronger. The weak jet from the nonrotating BH is the widest with the weakest fluctuations. Our results highlight the complex nonlinear connection between the black hole spin and the resulting disk and jet properties in the MAD regime.


Introduction
Accretion disks are ubiquitous in many astronomical objects, such as active galactic nuclei (AGNs) and X-ray binaries.The structure of the accretion disk mainly depends on the accretion rate.At a high accretion rate, close to the Eddington limit, the disks are typically geometrically thin and optically thick, and the models of Novikov & Thorne (1973) and Shakura & Sunyaev (1973) are thought to accurately describe their physics.When the accretion rate is much lower than the Eddington accretion rate, the cooling time becomes longer than the accretion time, leading to a radiatively inefficient accretion flow (RIAF), and the disk becomes geometrically thick and optically thin.In this regime, there are several theoretical disk models, such as the advection-dominated accretion flow (ADAF; Narayan & Yi 1994, 1995;Abramowicz et al. 1995;Yuan & Narayan 2014).The estimated low luminosity of Sagittarius A å as well as the black hole (BH) at the center of the M87 galaxy compared to the Eddington luminosity suggests that these BHs accrete in the form of an ADAFs (Yuan et al. 2002).
It is widely believed that the structure of an accretion flow consists of a turbulent accretion disk, a bipolar jet, and a magnetized wind (see, e.g., De Villiers et al. 2003;McKinney & Gammie 2004).The details of this structure strongly depend on the configuration and the strength of the magnetic field inside and outside the disk.The magnetic fields in the disk, either advected from large distances or created in situ by the dynamo effect, are amplified by the magnetorotational instability (MRI; Balbus & Hawley 1991, 1998).They ultimately drive angular momentum transport, regulate accretion, and produce a bipolar, strongly magnetized jet.There are two distinct modes of accretion depending on the magnetic fields surrounding the BH, which ultimately lead to two different disk configurations.In the standard and normal evolution (hereinafter, SANE; Narayan et al. 2012;Sądowski et al. 2013), the magnetic field pressure is not strong and the accretion process is smooth.The accretion disk, although turbulent, extends nearly evenly up to the horizon.In this accretion mode, angular momentum is transported mostly radially inside the disk by MRI (Chatterjee & Narayan 2022).
The second type of disk is termed magnetically arrested disk (MAD; Bisnovatyi-Kogan & Ruzmaikin 1974; Narayan et al. 2003;Igumenshchev 2008).In this model, the magnetic flux accumulates near the horizon until it saturates.In fact, the accumulated magnetic field becomes so strong close to the BH that it can change the dynamics of the infalling matter, thereby regulating the accretion.It was found in 2D simulations that the accretion can be nearly fully stopped by the magnetic pressure and then resumes, following the reconnection of the magnetic field lines at the equator (see, e.g., Chashkina et al. 2021).This picture, however, only partially holds in 3D simulations: accretion proceeds continuously via the development of nonaxisymmetric instabilities (Spruit et al. 1995;Begelman et al. 2022), with the infalling gas being shaped into filaments by the strong magnetic field (see, e.g., Figure 1 of Xie &Zdziarski 2019, andWong et al. 2021).
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
The MAD state has attracted increasingly more attention in the past few years, following the observation of the closest region to the supermassive BH in M87 and Sagittarius A å by the Event Horizon Telescope (EHT) collaboration (Event Horizon Telescope Collaboration et al. 2019Collaboration et al. , 2022a)).By comparing the images taken in the radio band to postprocessed GRMHD simulations, it was determined that the accretion should operate in the MAD state for these two BHs (Event Horizon Telescope Collaboration et al. 2021Collaboration et al. , 2022b)).Yuan et al. (2022) independently arrived at the same conclusion for M87 by studying rotation measurements.Moreover, several studies, including Dexter et al. (2020b), Porth et al. (2021), Ripperda et al. (2022), and Scepi et al. (2022), proposed a model in which the flares observed in Sagittarius A * by the GRAVITY experiment (GRAVITY Collaboration et al. 2018, 2020) have their origin in the magnetic flux eruptions, which are characteristic of the MAD state (Igumenshchev 2008).It is clear that the MAD state is ubiquitous at low accretion rates, and a better understanding of its properties will shed light on key observations of BHs and their physics.
An important characteristic of the MAD state is the saturated magnetic flux at the horizon.Since numerical studies of disk evolution depend on assumed initial conditions, in order to numerically study these systems, one has to use an appropriate initial magnetic field configuration.Tchekhovskoy et al. (2011) found such a configuration that results in the transport of a large amount of magnetic flux to the horizon.Using the normalized BH spin parameter a = 0.99, they found that the disk would be in the MAD state when the MAD parameter reaches Φ B ∼ 50.Here and below, the MAD parameter is defined as the ratio of the magnetic flux to the square root of the mass accretion rate at the horizon.Later, Tchekhovskoy & McKinney (2012) performed two simulations with a = 0.9 and a = −0.9 and found that the MAD parameter of the retrograde disk is about 30, which is smaller than that of prograde disk.Narayan et al. (2022) studied the dependence on the spin of the MAD parameter, confirmed the results from Tchekhovskoy & McKinney (2012), and also found that the maximum value of the MAD parameter (for a given initial magnetic configuration) is, in fact, reached for a BH spin a ; 0.5.
The MAD parameter Φ B is an important quantity in quantifying both the disk structure and the emerging jet.It was found that disks in the MAD state around rotating BHs launch strong and powerful jets via the Blandford & Znajek (1977) mechanism (see, e.g., Tchekhovskoy et al. 2011).The emerging power of the magnetized jet is proportional to the square of the MAD parameter multiplied by the mass accretion rate (Tchekhovskoy et al. 2011).It is therefore important to constrain the accretion parameters and mechanisms that determine and regulate the value and the duty cycle of the MAD parameter.
The strong magnetic field during the MAD state pushes the gas out and stops, or at least regulates, the dynamics of the infalling of matter.Therefore, one may naively expect that an increase in the magnetic flux should result in a decrease in the mass accretion rate, namely, that they are anticorrelated.An anticorrelation is also expected if accretion proceeds by interchange instability, as matter replaces a highly magnetized region closer to the BH, resulting in a drop in the MAD parameter (Porth et al. 2021).However, Porth et al. (2021) did not find a correlation or an anticorrelation between the mass accretion rate  M and Φ B .Here, we expand their results for all BH spins: as we show below, none of our simulations show a correlation or an anticorrelation between  M and Φ B .On the other hand, we do find an anticorrelation between their time derivatives.
A complete understanding of the interplay between accretion and saturated magnetic field close to the horizon and how it impacts the structure and power of the disk, jet, and wind, as well as their evolution is still lacking.Narayan et al. (2022) found that prograde disks have wider jets and that the shape itself depends on the spin of the BH.The structure of the jet and the disk mainly depends on the magnetic field pressure and the gas pressure inside the disks and jets.Begelman et al. (2022), using the simulation in the MAD regime around a rotating BH with a = 0.9375 from Dexter et al. (2020bDexter et al. ( , 2020a)), proposed that the disk properties in this regime are due to a dynamically important toroidal field in the close vicinity of the BH.However, Chatterjee & Narayan (2022), studying accretion in the case of a nonrotating BH, showed that the radial magnetic field b r is actually stronger than the toroidal magnetic field at small radii.This apparent inconsistency needs to be resolved, even though its origin, namely the BH spin, is trivially identified.Here, we find that close to the BH horizon, the toroidal component of the magnetic field, b f , is increasing with the absolute value of the BH spin, |a|, being dynamically unimportant for a = 0, but dominant for |a| → 1.
It is generally thought that the strong magnetic field close to the BH suppresses the development of MRI because the disk height is smaller than the wavelength of the most unstable vertical mode (McKinney et al. 2012;Marshall et al. 2018;White et al. 2019).This in turn affects the rate of angular momentum transport inside the disk, the details of which are still poorly understood in the MAD state.Recently, Begelman et al. (2022) argued that MRI is actually not suppressed closest to the BH by considering nonaxisymmetric modes (Das et al. 2018).On the other hand, Chatterjee & Narayan (2022) demonstrated that a for nonrotating BH, angular momentum is transported predominantly by magnetic flux eruptions, which are characteristics of a disk in the MAD state.It remains to be understood whether the contribution of this process still dominates the transport of angular momentum over MRI in the case of a rotating BH.In addition, angular momentum from the disk is transported by the emerging winds above and below it.
An additional source of angular momentum from an accretion disk system is the jet, in the form of Maxwell stresses.This contribution to the angular momentum mainly originates from the BH and not from the disk, and as such acts to spin down the BH (Chatterjee & Narayan 2022).As we show here, this transport can significantly contribute to the amount of angular momentum deposited in the external medium.Since the power of the jet depends on the BH spin, this also constrains the cosmic period over which the system is active.The net rate of angular momentum strongly depends on the spin and is highest for prograde disks.
The shape of the jet itself also depends on the BH spin.It is determined by the balance between the internal and external stresses.As explained above, the main source of stresses are the magnetic field components, which in turn depend on the BH spin.Furthermore, the BH spin not only affects the shape of the time-averaged jet, but also fluctuations around it.As we show here, retrograde disks produce the narrowest jets with the widest fluctuations, while nonspinning BHs produce the widest jets with the narrowest fluctuations.
Given the complexity of accreting systems, many previous works used general relativistic (eventually, radiative and two temperatures) magnetohydrodynamic (GRMHD) simulations to investigate the properties and evolutionary process of MAD disks (Narayan et al. 2012;Sądowski et al. 2013;Porth et al. 2017;White et al. 2020;Porth et al. 2021;Begelman et al. 2022;Chatterjee & Narayan 2022;Narayan et al. 2022).In the past two decades, with the rapid increase in compute capability, these simulations have become increasingly popular and practical (Gammie et al. 2003;Anninos et al. 2005;Noble et al. 2006;Stone et al. 2008;Porth et al. 2017;Tchekhovskoy 2019;Liska et al. 2022;Bégué et al. 2023).In addition, general-purpose graphics processing units (GPU) started to be used to accelerate fluid simulations in recent years, as they are particularly well suited to be ran on GPUs.As a result, several GRMHD codes can now use GPU accelerators (see, e.g., Chandra et al. 2017;Liska et al. 2020;Bégué et al. 2023;Shankar et al. 2023).grim uses the library ArrowFire to achieve GPU compatibility (Chandra et al. 2017).Liska et al. (2020) developed H-AMR with openMP, MPI, and CUDA.Building on HARMPI, our group developed a new GPUaccelerated GRMHD code, cuHARM, which uses openMP and CUDA (Bégué et al. 2023).This code is thoroughly optimized to maximally harness the power available in NVidia GPUs, with more than 50% computation efficiency on NVidia A100 cards.For the results presented here, the simulations are made on a single multi-GPU workstation.
In this paper, using our GRMHD code cuHARM, we study the role played by the magnetic field in the MAD state for different BH spins and its effect on the structure of the disk and the jet.To do this, we present several simulations with different initial magnetic field strengths and BH spins.This paper is organized as follows.In Section 2 we present the setup of our simulations and introduce the numerical diagnostics used in our analysis.We discuss the dynamics of the accretion disk system in our simulations in Section 3. In particular, after specifying the inflow and outflow equilibrium radius and the time evolution of  M and f B in Sections 3 and 3.2, respectively, we study (i) the absence of a correlation between the mass accretion rate and the MAD parameter, and we introduce the anticorrelation between their time derivative in Section 3.3; (ii) the shape of the jet as a function of spin in Section 3.7, finding that retrograde disks are narrower than their corresponding prograde disks; and (iii) the component-wise contributions of the magnetic pressure to underline the differences between a spinning and a nonspinning BH in Section 3.8, where the toroidal component is found to be subdominant for a = 0, but similar to the radial component for large |a|.In Section 4 we discuss the transport of angular momentum for our simulations with spin a = − 0.94, a = 0, and a = 0.94, underlying the differences between each BH spins.The summary and conclusions of the paper are given in Section 5.

Simulations
We perform several simulations with cuHARM (Bégué et al. 2023), which uses the finite-volume method to numerically solve the conservative GRMHD equations (for reviews, see, e.g., Martí & Müller 2003;Font 2008;Rezzolla & Zanotti 2013).The code is written in CUDA-C and openMP, and all calculations of cuHARM are accelerated by GPU (only the data transfer and exports are powered by CPU).To perform the simulations, the results of which are presented in this article, we use an Nvidia DGX-V100 server with 8 Nvidia V100 GPUs.

Initial Setup
In this paper, we study the accretion flows in the MAD state around a spinning BH.Our simulations begin with the stationary axisymetric torus described by Fishbone & Moncrief (1976).We set the gas adiabatic index to Γ = 14/9 and consider an initially large disk with r in = 20r g and = r r 41 g max , where r in is the inner boundary of the disk, r max is the radius at which the pressure reaches its maximum, and r g is the gravitational radius.The matter and internal energy densities are normalized such that for the initial disk, the maximum matter density ρ in the entire disk is normalized to r = 1 max .The internal energy density is scaled accordingly.Since the initial torus is in equilibrium, it does not evolve spontaneously.We therefore add small random perturbations (set to 4%) to the internal energy density u as the seed of instabilities, which will promote accretion.
This initial torus is in full hydrodynamic equilibrium and therefore does not contain any magnetic field.We introduce a purely poloidal subdominant magnetic field defined by the vector potential A, such that A r = A θ = 0 and which has been previously employed in, e.g., Wong et al. (2021) and Narayan et al. (2022).Here, r, θ, and f are the horizon-penetrating spherical Kerr-Schild coordinates.The corresponding magnetic field is initially a single loop confined to the disk.The magnitude of the magnetic field is further normalized by the parameter , where p gas,max is the maximum gas pressure, = p b 2 b,max 2 is the maximum of the magnetic field pressure, and b = b μ b μ is the norm of the four-vector magnetic field (see Section 2.3 below).This expression of the magnetic field is designed to ensure that enough magnetic flux can be transported to the BH throughout the course of the simulation and saturates its magnetosphere (see further discussion in Tchekhovskoy et al. 2011).
We conduct a series of simulations with different BH spins, a ä {−0.985, −0.94, −0.85, −0.5, 0, 0.5, 0.85, 0.94, and 0.985} and an initial magnetization β 0 = 100.In the case of a retrograde disk with a = −0.94,we also varied the initial magnetic field strengths, with β 0 ä {100, 200, 400, and 800}.We evolve most of the simulations until t = 2 × 10 4 t g , where t g = r g /c, except for aM94b800, which is evolved to t = 2.5 × 10 4 t g due to the weak initial magnetic field and the longer time required to reach the MAD state for this setup.Additionally, simulation aM94b100h is evolved until t = 5 × 10 4 t g in order to study the long-time behavior of our accretion disk system.We use the spin and the initial β 0 to name the simulations: "a" stands for spin, "M" indicates a negative value, and "b" represents the initial β 0 , such that, for instance, aM94b100 stands for a simulation with the negative ("M") spin a = − 0.94 and an initial β 0 = 100.A summary of all the simulations used in this work is given in Table 1.

Numerical Aspects
Since we are studying accretion around rotating BH, we use the Kerr metric for our simulations.The Kerr-Schild (KS) coordinate system (t, r, θ, f) is used as the physical coordinates, while to both enhance the robustness of the calculation and focus the computation in the region of interest, namely close to the BH and at the equator, the modified Kerr-Schild (MKS; see, e.g., McKinney & Gammie 2004) coordinates (t, q 1 , q 2 , f) are used in the numerical calculation.The relation between these coordinates, as implemented in cuHARM, can be found in Section 4.1 of Bégué et al. (2023).
We use the inflow and outflow boundary conditions in the radial direction at small and large radii, respectively.In the θ direction, we use the reflective boundary condition, and the periodic boundary condition is used in the f direction.To address the potential numerical errors in an empty or strongly magnetized region, we adopt the same flooring model as in Bégué et al. (2023), which is used in many other papers, e.g., Porth et al. (2019).The density ρ and the internal energy u are limited using Matter and energy are added when needed to preserve the conditions b 2 /ρ < 50 and b 2 /u < 2.5 × 10 3 .The reference resolution of most simulations is 192,96,96), except for aM94b100h and a0b100h, which have a slightly higher resolution (N r × N θ × N f ) = (256, 128, 128).Here, the "h" appended to the name stands for "high resolution."The resolution of the simulations presented in this paper is somewhat lower than that used in some of the simulations presented in recent works.For example, Narayan et al. (2022) performed fairly similar simulations with a resolution of 288 × 192 × 144.White et al. (2020) examined the impact of different resolutions, and they argued that the accretion rate and the general disk structure agree across the simulations with different resolution.We use our higher-resolution simulations, aM94b100h and a0b100h, to check the solidity of our results to a change in the resolution.We did not find any significant difference between the low-and high-resolution simulations.

Diagnostics
Following Komissarov (1999), let b μ ≡ (åF) μν u ν represent the four-vector magnetic field and u μ be the four-velocity, which is orthogonal to b μ .In the ideal MHD limit, the dual to the Faraday tensor is given by In this limit of a perfectly magnetized fluid, the stress energy tensor T μ ν is given by Here, h = ρ + u + p g is the enthalpy, p g is the gas pressure, and b 2 = b μ b μ and g μ ν is the metric tensor with the determinant noted -g .Using cuHARM, we solve the general relativistic MHD equations, namely Note.The three last columns give the time-averaged value of the MAD parameter Φ B , of the accretion rate  M , and of the jet efficiency η for 10 4 t g < t < 2 × 10 4 t g (1.5 × 10 4 t g < t < 2.5 × 10 4 t g for aM94b800, and 10 4 t g < t < 5 × 10 4 t g for aM94b100h).The mass accretion rate, the MAD parameter, and the jet efficiency are also displayed in the top row of Figure 4 as a function of BH spin.
( ) which are the equations of mass conservation, the equation of energy and momentum conservation, and the homogeneous Maxwell equations, respectively.
The MAD state mainly depends on the magnetic flux through the horizon.Therefore, we define the following radial diagnostics, which can finally be evaluated at the horizon: 1.The mass accretion rate: The magnetic flux crossing the horizon (through one hemisphere): . 10 3. The energy flux through the horizon toward the BH: The angular momentum flux in the radial direction: The jet efficiency at the horizon: Note that our definitions of the angular momentum flux  J and energy flux  E are opposite to those employed by Narayan et al. (2022), but they agree with the definitions of, e.g., Porth et al. (2019).We are also interested in the structure of the disk and the jet.Therefore, we define the additional following diagnostics: 1.The disk height, denoted by (h/r): . 16 0 2 0 2 3.The disk average of a quantity q: , 0 2 0 0 2 0 4. The disk average of a quantity q but within a narrow θ range (used below in calculating the pressure): We will also present time-averaged quantities and time-and faveraged maps, which are computed for 10 4 t g < t < 2 × 10 4 t g , unless specified otherwise, and the full 2π range for f.

Inflow Equilibrium Radius and Radial Limit of Our Analysis
Most simulations are evolved until t = 2 × 10 4 t g , which is sufficient for the disk to be in the MAD state.We assume this state to be established when the MAD parameter reaches its average value at late times.We find that in all cases, the averaged MAD parameter is comparable to or greater than 15, the limiting value proposed by Tchekhovskoy et al. (2011) to define the MAD state.Only aM94b800 barely reaches the MAD state at t = 2 × 10 4 M because of the initial small magnetic field normalization.This simulation is therefore evolved further until t = 2.5 × 10 4 M, at which time the MAD state for this initial setup is well established.
We first study the radial profile of the mass accretion rate  M and of the angular momentum flux  J , as given by Equations ( 9) and (12), respectively.The results are displayed in the left and middle panel of Figure 1 for simulations aM94b100 and a94b100, serving as examples.The results are similar for all other simulations.For aM94b100 and a94b100, the radial profiles of  M and  J are averaged over three different time intervals, namely from 5 × 10 3 to 10 4 t g , from 10 4 to 1.5 × 10 4 t g , and from 1.5 × 10 4 to 2 × 10 4 t g .In the last time interval, the mass accretion rates of these two runs are independent of the radius r for r < 30r g .The angular momentum fluxes also exhibit a similar pattern, remaining constant at r < 30r g .This means that the inner region with r < 30r g is in inflow equilibrium state.
There are two key differences between the prograde and retrograde cases.First, the sign of the angular momentum flux is different.The prograde BH has a positive angular momentum flux, namely, the BH is losing angular momentum, as was previously reported in Narayan et al. (2022).Conversely, the angular momentum flux of the retrogade BH is negative, meaning that it accumulates positive angular momentum and spins down.The second difference is the existence of a radius at which the angular momentum flux changes sign for a retrograde BH.This indicates that there is a net flux of angular momentum away from the BH at large distances.This radius is also observed in the simulation presented in Narayan et al. (2022; see their Figure 3 with the sharp drop of angular momentum flux at r ∼ 10 2 r g ).We discuss in more detail the differences for the angular momentum transport between prograde and retrograde disks in Section 4.
The duration of most of our simulations is shorter than that of recent long-time simulations.For example, Narayan et al. (2022) evolved the accretion system until t = 10 5 t g , such that they obtained an inflow equilibrium radius at ∼100r g , larger than ours by a factor 3-4.To verify the reliability of our results over a longer time span, we extend the simulation aM94b100h to t = 5 × 10 4 t g .In the right panel of Figure 1, we present the radial profile of the mass accretion rate  M and the angular momentum flux  J for this extended period.The time-averaged quantities for 2 × 10 4 t g < t < 3 × 10 4 t g , 3 × 10 4 t g < t < 4 × 10 4 t g , and 4 × 10 4 t g < t < 5 × 10 4 t g are displayed.These profiles show the same behavior as aM94b100, but the inflow equilibrium radius now extends to r > 50r g .Therefore, the shorter simulation time does not affect the establishment of the equilibrium state, but only limits the inflow equilibrium to a smaller radius.In the following, we focus on the inner region characterized by r 30 g and demonstrate, where relevant, that we obtain results that are similar to those presented in Narayan et al. (2022).Considering that we are using (i) a different numerical scheme, (ii) a different numerical resolution, and (iii) analyzing the results earlier in the simulation, this shows the solidity of these results.

Mass Accretion Rate and Magnetic Flux at the Horizon
The time evolution of the mass accretion rate  M at the horizon, derived from Equation (9), is shown in Figure 2 for all simulations with different spins.Approximately the same behavior is observed in all cases.After an initial quiet period, the mass accretion rate steadily increases to reach its late-time average around t ∼ 5 × 10 3 t g .The left column of Figure 2 shows the time evolution of  M for retrograde disks, while the corresponding prograde disks are displayed in the right column.The y-axis scale is identical to facilitate comparison.There is a clear dependence of the mass accretion rate on the BH spin.First, as displayed by the dashed line, representing the time average of the mass accretion rate for 10 4 t g < t < 2 × 10 4 t g , prograde disks have a higher mass accretion rate than retrograde ones.The time-averaged accretion rate from t = 10 4 t g to t = 2 × 10 4 t g is listed in Table 1 alongside the 1σ temporal variation.Second, the variation in the mass accretion rate is more pronounced in prograde than in retrograde disks.
The time evolution of the MAD parameter can be described as follows.As the accretion proceeds, the magnetic flux at the horizon accumulates until it saturates.At this point, the MAD parameter, as well as the magnetic flux threading the horizon and the mass accretion rate, are regulated by eruptions of magnetic flux.These eruptions expel the magnetic field to far distances from the BH.In the MAD state, the pressure of the saturated magnetic field is balanced with the gas pressure (Bisnovatyi-Kogan & Ruzmaikin 1974; Narayan et al. 2003).
Although the flux eruptions cause the fluctuation in both the magnetic flux (Begelman et al. 2022) and the mass accretion rate, the MAD parameters generally remains stable around their time-averaged value.
To explore the differences in the MAD state for different spins, we present the time evolution of the MAD factors in Figure 3, where the left column shows negative-spin simulations and the right column shows positive-spin simulations, similar to Figure 2. The horizontal dashed lines represent the timeaveraged MAD parameters from t = 10 4 t g to t = 2 × 10 4 t g , at which interval the MAD state is well established.The timeaveraged MAD parameters during this period are also listed in Table 1.We find that the MAD parameters of prograde disks are higher than those of the corresponding retrograde disks.Additionally, similarly to the mass accretion rate, we find that the temporal fluctuations of the MAD parameters have a larger amplitude for a prograde disk.This agrees with the findings of Porth et al. (2021), who also found that their corotating simulation has weaker flux expulsions than the counter-rotating case.We also observe many more small flux eruptions for the corotating case, but they are accompanied by several strong eruptions, such as the eruption at ∼t = 9 × 10 3 for a94b100, as seen in Figure 3.
Both the mass accretion rate and the MAD parameter are strongly dependent on the BH spin, a.This relation is displayed in the top left panel of Figure 4. We find that the mass accretion rate is highest for low, positive spin, and it drops for higher values of the spin for both prograde and retrograde disks, with retrograde disks showing lower values than prograde disks.For the a = 0 spin, we show two results, one with the standard resolution, and the other with a higher resolution, which we believe is more accurate (see the discussion below).We also show in Figure 4 the 1σ temporal fluctuation of the mass accretion rate as error bars.Clearly, prograde disks have stronger fluctuations (represented by the larger error bars in Figure 4) than retrograde disks.
We present the relationship between the MAD parameters during the MAD state and the spin of BH in the top middle panel of Figure 4.The MAD parameter increases as the spin increases from a = −0.985 to a = 0.5, and then decreases as the spin increases from a = 0.5 to a = 0.985.Narayan et al. (2022) reported a similar trend and used a third-order polynomial to fit the relationship between a and the MAD Figure 1.The radial profile of the mass accretion rate  M (solid line) and of the angular momentum flux  J (dashed line) for aM94b100 (left panel), a94b100 (middle panel), and for the long-term evolution aM94b100h (right panel) at different time periods (see legend).As the time increases, the steady region extends to a larger radius.This suggests that aM94b100 and a94b100 have established inflow-outflow equilibrium at r < 30r g for 1.5 × 10 4 t g < t < 2 × 10 4 t g , while it is about 70r g at 4 × 10 4 < t g < 5 × 10 4 for aM94b100h.For the simulations with the negative spin in the left and right panels, the angular momentum fluxes change sign at r = 30-50r g and r = 50-100r g , respectively.
parameter Φ B .We include their fit result (displayed by the dashed red line) in the top middle panel of Figure 4 to compare with our results.Their definition of the MAD parameters differs slightly from ours, so we renormalized their fit formula by a factor of p 4 2 to account for this discrepancy.Our results agree well with theirs, which enhances the credibility of the results, considering that we use a different code, different resolutions, and a shorter simulation duration.
Our simulations have a somewhat lower resolution (by a factor of ≈4) than those of Chatterjee & Narayan (2022) and Narayan et al. (2022).To address this limitation, we conducted two simulations, aM94b100h and a0b100h, which have the same initial conditions as aM94b100 and a0b100, respectively, but differ by having a higher resolution.We show the mass accretion rates and the MAD parameters of these two simulations in Figures 2 and 3, alongside their counterpart with the lower resolution.In the case of aM94b100 and aM94b100h, there is no significant difference in terms of mass accretion rates and MAD parameters.However, for a0b100 and a0b100h, the mass accretion rates are different by a factor of M .The left column pertains to retrograde disks, while the right columns corresponds to prograde disks.For convenience, we set the y-axes to identical scales to facilitate direct comparison.The horizontal lines are the time-averaged mass accretion rate  M , where the averaged is taken from 10 4 t g to 2 × 10 4 t g .The last row shows a discrepancy for a = 0 between our low-(a0b100) and high-resolution (a0b100h) simulations.It shows, in this case, that the resolution is too low to resolve the mass accretion rate accurately.However, we found that the MAD parameters are nearly identical, indicating that the properties of the MAD state should be similar.We provide more details on this issue at the end of Section 3.2.nearly three.According to White et al. (2019), this indicates that the resolution is too low for our simulation with spin a = 0 (but not for our simulations with rotating BHs, as is suggested by aM94b100h).They demonstrate that resolving the mass accretion rate demands a better resolution than to resolve the MAD parameter, as is also shown here.Indeed, the saturation values of the MAD parameters are nearly identical.This indicates that the properties of the MAD state should be similar.
Most of our simulations end at t = 2 × 10 4 t g .As discussed in Section 3.1, we extend aM94b100h to t = 5 × 10 4 t g .The time evolution of the mass accretion rate  M and of the MAD parameter Φ B is shown in Figure 5.The MAD parameter remains relatively stable, oscillating around its average value after the MAD state is established.However, the mass accretion rate  M gradually decreases until t = 5 × 10 4 t g , which is caused by the spread of the disk as the simulation advances.

Characteristic Dynamical Timescale
In the MAD state, magnetic flux eruptions play a major role in the dynamics of the accretion: the strong magnetic pressure pushes the gas out, which leads to the fluctuations in the accretion rate and of the magnetic flux.Moreover, these magnetic flux eruptions have been proposed as the origin of the Sgr A * flares that were observed in infrared by the GRAVITY Figure 3.The time evolution of MAD parameters Φ B .The left column represents retrograde disks, while the right column is for prograde disks.We use the horizontal lines to denote the time-averaged Φ B from 10 4 t g to 2 × 10 4 t g , which is referred to as the typical MAD parameter during the MAD state.
collaboration (GRAVITY Collaboration et al. 2018, 2020).In particular, Dexter et al. (2020b) numerically estimated the recurrence time of the flare to be between 10 3 and 10 4 r g /c, corresponding to 5-50 hr for Sgr A * .The maximum of the intensity is correlated with the sharp drop of the MAD parameter at the onset of each eruption (Dexter et al. 2020b), lasting around 100 r g /c (Ripperda et al. 2022).The dynamic of the flux tube created in each eruptions was studied in detail by Porth et al. (2021).They found that the motion of the lowdensity/high-magnetization region is strongly radial at first because of the magnetic tension, and then tends to circularize when the field is nearly vertical.The typical lifetime of these magnetic flux tubes was found numerically to be around two orbits, depending on their size and magnetic energy.
We use the discrete Fourier transform to study the duty cycle of the MAD parameters.The analysis is performed considering data from t = 10 4 t g to t = 2 × 10 4 t g , for which, as a preparation step, the time series are detrended and their average removed, so that the fundamental Fourier coefficient is null.We show the power spectrum of the MAD parameter Φ B for our simulations in Figure 6.We focus on the power spectrum from f ; 4 × 10 −3 to f ; 2 × 10 −4 , which corresponds to the period from T = 250t g to T = 5000t g .All our simulations show the same trend: after a few dominant peaks, the power spectrum decays proportionally to f −1 , characteristic of pink noise.This is shown as dashed red lines in Figure 6.Janiuk & James (2022) performed a similar study, not of the MAD parameter, but rather of the mass accretion rate.Furthermore, their disks were smaller than ours, with r in = 6 or 12 and = r 12 max or 25 r g .They found a power spectrum of the mass accretion rate that was well approximated with a power-law of index ∼1.5.They further reported a dependence between the index and the BH spin.Our analysis, on the other hand, does not show such a dependence of the spin on the index for the MAD parameter.
For some simulations, such as aM94b100 and a5b100, the power spectrum shows clear peaks in the frequency window preceding the onset of the noise.The corresponding periods are t = 1, 111t g for aM94b100 and t = 1, 666t g for a5b100.For the other simulations, no dominant timescale can be identified unambiguously.Previous identifications of a cyclic behavior of M at the horizon for different spins.The blue points are the results of low-resolution simulations, while the red points are the high-resolution simulations.Top middle panel: normalized MAD parameters for different spins.Plotted on top (dashed green line) are the fit results of Narayan et al. (2022).Their formula is divided by p 4 2 to be consistent with our definition of the MAD parameter (see Section 3.2 for more details).Top right panel: jet efficiency as a function of spin.The green points are the predicted efficiencies, and the blue (red) points are the efficiencies derived from our simulation (red shows the higher resolution).The error bars represent the 1 sigma temporal variation of each quantity.η > 1 represents an efficiency higher than 100%.Bottom left: dependence of the slope on the spin a.The slope k decreases as the spin a increases until a = 0, at which point the slope k begins to increase as a continues to increases.However, a = 0.85 slightly deviates from this trend.Bottom right panel: time-averaged disk height h/r given by Equation (15) as a function of radius.
the mass accretion rate and magnetic flux were published by, e.g., Chashkina et al. (2021).By analyzing the data of their 2D GRMHD simulations, they found a period around t ; 500r g /c (see also, e.g., Igumenshchev 2008;Dihingia et al. 2023).However, as noted by Chashkina et al. (2021), this pattern is an artifact due to the nature of 2D simulations, preventing nonasymmetric instabilities and phenomena, which are inherently responsible for the continuous accretion in 3D, from developing.

Relation between the MAD Parameter and the Mass Accretion Rate
The strong magnetic field during the MAD state pushes out the gas and stops, or at least regulates the dynamics of the infalling of matter.Therefore, we naively expect that an increase in the magnetic flux should result in a decrease in the mass accretion rate, in other words, that they are anticorrelated.An anticorrelation is also expected if accretion proceeds by interchange instability, as matter replaces a highly magnetized region closer to the BH, resulting in a drop in the MAD parameter (Porth et al. 2021).However, Porth et al. (2021) did not find a correlation or an anticorrelation between  M and Φ B .We expand their results for all BH spins: none of our simulations shows a correlation or an anticorrelation between  M and Φ B .We further studied the relation between the time derivatives of the mass accretion rate  M and of the MAD parameter Φ B .We first remove their average and detrend the time series, then we use a Fourier transformation to filter out the high-(and weak) frequency components.We checked that this procedure does not change the results presented below.Instead, it allows for a better visualization of the result.The processed time derivatives are displayed in the left column of Figure 7 for aM94b100 and a94b100.The blue line is the normalized  dM dt, and the dashed red line represents the normalized dΦ B /dt.Here, by "normalize" we mean that the two time derivatives are stretched so that the maximum of their absolute value is 1.We show the resulting time series for aM94b100 and a94b100 as examples, but all other simulations exhibit a similar behavior.To the (positive) heights of  dM dt correspond the (negative) lows of dΦ B /dt and vice versa.
The time derivatives of  M and Φ B are clearly anticorrelated.In the right column of Figure 7, we show this anticorrelation.This means that when the mass accretion rate increases, the flux threading the horizon decreases.This can be understood as follows.Focusing on the region close to the equator, the horizon surface is separated into two regions: the accretion funnel from which the matter falls inward, and the highly magnetized, low-density regions of the magnetic flux eruption.As the MAD parameter increases, the magnetic pressure outside of the accretion funnel increases and it is compressed, resulting in a lower accretion rate.Similarly, as the magnetic flux eruption develops, the magnetic flux at the horizon drops, thereby reducing the magnetic pressure and allowing for a larger accretion funnel and a higher accretion rate.The latter is also enhanced by the interaction between the low-density magnetic flux tube and the inner region of the turbulent disk at r ∼ 10-15 r g .Indeed, the magnetic flux tube velocity is sub-Keplerian, reducing the velocity of the matter in the disk, which then falls radially onto the BH.
We further perform a linear fit to this anticorrelation and show the dependence of the slope on the spin a in the bottom left panel of Figure 4. We find that the nonspinning BH has the steepest slope k, which means that smaller variations in mass accretion rates correspond to a larger variation in the MAD parameter, compared to other BH spins.A smaller slope k is in general a characteristic of a high absolute value of the spin a, except for a = 0.5 or for a = 0.85, which would deserve further investigations.Finally, we point out that the filtering performed in preparation of the data changes the value of k by a factor of up to about 2, but the trend remains the same when this rescaling is accounted for.We further note that there seems to be a strong dependence of k on the resolution, as is shown in the bottom left panel of Figure 4 by the two red dots, which are the slopes obtained from simulations aM94b100 and a0b100.For aM94b100 and aM94b100h, the reason for this dependence may stem from the fact that while the resolution does not affect the amount of material and the magnetic field that is accreted, as shown by Figures 2 and 3, it might influence the relationship between the eruption and accretion.This in turn affects the spectrum of the eruptions.However, since the mass accretion rate is different for aM0b100 and aM0b100h, this analysis is inconclusive and would deserve further investigations.

Jet Efficiency
The jet efficiency, given by Equation ( 14), depends on the magnetic flux through the horizon normalized by the accretion rate  M. In agreement with many previous publications (see, e.g., Tchekhovskoy et al. 2011  efficiency of a positive-spin BH with a > 0.5 is higher than 100% on average, meaning that the BH is actually losing energy: the jet is powered by the Blandford and Znajek mechanism (Blandford & Znajek 1977).
We use the fitting formula proposed by Tchekhovskoy et al. (2010) and later used by Narayan et al. (2022) to interpret the efficiency, namely Here, k p = 0.08 is a numerical constant, and Ω H = ac/r H is the angular velocity at the horizon.This value of κ is chosen such that p accounts for the difference in the definition of Φ B used here and in Tchekhovskoy et al. (2010), and the numerical factor 0.08 is chosen to match the normalization of our numerical data.The predicted efficiency, i.e., using the measured value of the MAD parameter Φ B together with Equation (19), and the simulated efficiency, namely directly measured from our simulations as defined by Equation ( 14), are shown in the right panel of Figure 4 as a function spin.It is clear that they are consistent with each other.Narayan et al. (2022) reported the same behavior, but used a different value of κ = 0.05 that differs from our normalization by a factor smaller than 2. This discrepancy could be due to the time at which the measurements are performed.Indeed, when averaging the results of simulation aM94b100h at late times, 10 4 t g < t < 5 × 10 4 t g , we find a lower numerical factor of ∼0.06.

Effects of the Initial Magnetic Field Strength
Since the magnetic field plays a critical role in the accretion process and its regulation, we wish to assess the solidity of our result with respect to the initial magnetic field strength.To this Figure 6.Power spectrum of the MAD parameters in the time interval 10 4 t g < t < 2 × 10 4 t g (1.5 × 10 4 t g to 2.5 × 10 4 t g for aM94b800).In all simulations, the power is concentrated in the low-frequency region, which corresponds to periods about 500-2000t g .In the high-frequency region, the power spectrum decreases following ∼f −1 .The dashed red line shows this evolution.end, we also performed simulations with different values of β 0 for a = −0.94.In Table 1 we list three additional simulations with an initial magnetic field strength β 0 ä 100, 200, 400, and 800.The simulation aM94b100 has the strongest initial magnetic field strength, and aM94b800 has the weakest initial magnetic field strength.
In Figure 8 we present the mass accretion rates and MAD parameters for the additional simulations with a different initial Figure 7. Left: time evolution of the time derivatives of the mass accretion rate and of the MAD parameter.We show the results for aM94b100 and a94b100, but other simulations also exhibit similar behaviors.These two quantities,  dM dt and dΦ B /dt have opposite signs throughout the simulations.There is a clear anticorrelation between these two derivatives, shown in the right panels.
Figure 8. Mass accretion rates and MAD parameters of aM94b100, aM94b200, aM94b400, and aM94b800, which are simulations with the same initial conditions, but different initial magnetic field strengths.In the left panel, we show the mass accretion rates.The stronger the initial magnetization (the smaller β 0 ), the faster the increase in mass accretion rate.After t = 1.5 × 10 4 t g , the mass accretion rates agree across all simulations.In the right panel, the MAD parameters of these simulations are shown to be consistent with each other after t = 1.5 × 10 4 t g .However, the time required to reach the MAD state is clearly different.magnetic field β 0 .The left panel shows that at the beginning of the accretion process, the mass accretion rate  M increases with the initial magnetic field strength, which can be attributed to the fact that the transfer of angular momentum is initially driven by MRI, the development of which depends on the strength of the magnetic field.Therefore, the larger the initial magnetic field, the higher the mass accretion rate.However, at late time t > 15, 000t g , we see that the mass accretion rates achieve similar values across all simulations.
The right panel of Figure 8 shows that the MAD parameters also achieve similar values across all simulations during the MAD state, which is reached here after t > 15, 000t g specifically, for aM94b800.This suggests that the MAD state is independent of the initial magnetic field strength.However, the time required to reach the MAD state varies.To investigate the dependence of the required time on the initial β 0 , we use the following criteria to define t MAD , which is the time required to reach the MAD state: 1.The mean MAD parameter averaged from t = t MAD to t = t MAD + 1000 t g should be larger than the 1σ lower limit of the final MAD parameter.2. The MAD parameter in the MAD state should be highly variable.Therefore, we require the derivative of the MAD parameter to change sign several times between t = t MAD to t = t MAD + 1000 t g .
According to these two criteria, the time it takes to reach the MAD state t MAD for aM94b100, aM94b100h, aM94b200, aM94b400, and aM94b800 is 3650 t g , 3630 t g , 5220 t g , 8665 t g , and 15505 t g , respectively.In the left panel of Figure 9, we display the time required to reach the MAD state as a function of the initial magnetic field strength β 0 .We find that t MAD increases linearly with β 0 , namely, it increases as the initial magnetic field strength weakens.The best-fit result is a linear function, t MAD = 17β 0 + 1.9 × 10 3 .This result suggests that the accumulation of the magnetic flux linearly depends on the initial magnetic field strength.We also study the dependence of t MAD on the spin for the other simulations with the same initial β 0 .The right panel of Figure 9 displays the dependence of t MAD on the spin a.As the spin increases from a = −1 to a = 0, t MAD increases.It then decreases as a increases from 0 to 1.This dependence is fit with a third-order polynomial, which results in t MAD = 1.9 × 10 2 a 3 − 1.2 × 10 3 a 2 − 4.6 × 10 2 a + 4.4 × 10 3 .

Evolved Disk and Jet Structures
The strong magnetic field at the center of the accretion system shapes the disk and launches the bipolar jet at low radii.Therefore, the disk structure is an important characteristic of the MAD state.In the bottom right panel of Figure 4, we show the time-averaged disk height h/r, which is averaged from t = 10 4 t g to t = 2 × 10 4 t g .We find that retrograde disks are thicker than prograde disks.In the inner region (5r g < r < 15r g ), the disk height increases as the absolute value of the spin increases, which is consistent with the results of Narayan et al. (2022).However, in contrast to the findings in this aforementioned paper, we find that the nonspinning BH has the thinnest disks of all, while they argued that the disk of the nonspinning BH is thicker than that of prograde disks.A possible explanation of this discrepancy is the fact that our analysis is carried out at earlier times.
We show the time and azimuthal averages of the density ρ of the plasma parameter β and of the magnetization σ for the simulations aM94b100, a94b100, and a0b100 in Figure 10.These values are obtained using Equation ( 16), and they are then averaged from t = 10 4 t g to t = 2 × 10 4 t g .The difference between these three spins is clear.In the density ρ plot, the nonspinning BH has the thinnest disk, while the prograde disk is thickest, which is consistent with the bottom right panel of Figure 4.The β plots show similar behaviors as the density.
Figure 9. Left Panel: relation between the coasting time t MAD to reach the MAD state and the initial magnetic field.We show the results for aM94b100, aM94b200, aM94b100h, aM94b400, and aM94b800 in blue circles.The coasting time depends linearly on the initial magnetic field.We use a linear function to fit it, and the bestfit result is shown by the dashed red line.Right Panel: t MAD as a function of the value of the spin a.We use a third-order polynomial function to fit this relation, and the best-fit result is shown by the dashed red line.
Figure 10.Time and azimuthally averaged density ρ, ratio of gas to magnetic pressure β, and magnetization σ for aM94b100 (top), a0b100 (middle), and a94b100 (bottom).The black lines in the middle and right panels represent the β = 1 and σ = 1 conditions, respectively.The white wedges at the poles correspond to the first cell close to the pole.
There are clear vacuum regions inside the jet, especially for a0b100, which are attributed to the numerical flooring rather than to physical effects.
We show the polar angle of the jet boundary, which we associate with σ = 1, as a function of radius for aM94b100, a94b100, and a0b100 in the left panel of Figure 11.The vertical lines correspond to 1σ variations, and the points represent the median values of the polar angle θ for a specific radius r.These values are obtained for 10 4 t g < t < 2 × 10 4 t g .As shown in Figure 3, all three runs remain in the MAD state during this time period, and they maintain inflow-outflow equilibrium.It is seen that the prograde disk has a wider jet than the retrograde jet at r < 20r g , which is consistent with the findings of Narayan et al. (2022).It is also clearly seen that the nonspinning BH has the widest jet.This result is consistent with Figure 10.
In order to quantify the degree of variation compared to their mean, we use the standard error s q q , where σ θ is the 1σ variation of θ, and q is the median.The temporal and radial (between 2r g and 30r g ) averages of s q q for aM94b100, a94b100, and a0b100 are 0.11, 0.08, and 0.07, respectively.We further show the kernel density estimation (KDE) of θ σ=1 at r = 10r g in the right panel of Figure 11.The boundary of the retrograde disk has the largest variations, as expected from the strong shear stresses at the boundary between the jet and the disk, which have opposite toroidal velocity.Our results in this sense are consistent with those of Wong et al. (2021), who argued both analytically and numerically that retrograde disks have stronger shear across the jet-disk boundary than prograde disks.

Pressure
In Section 3.7 we demonstrated that prograde disks are thinner and have wider jets.We attribute this to the pressure distribution inside the disk and jets.In Figure 12 we show the different components of the magnetic pressure and of the gas pressure for aM94b100, a0b100, and a94b100.These components are calculated using Equation (18) and are further time averaged.This implies that they are biased toward highdensity regions.We find that the overall total pressure p + b 2 /2 is lowest for nonrotating BHs and highest for a prograde disk.At small radii, the magnetic pressure b 2 , represented by the blue line, dominates the gas pressure.As the radius increases, the gas pressure gradually becomes dominant.The transition radii for = a 0.94, 0, and 0.94 are r eq = 2.71r g , 3.13r g , and 2.35r g , respectively.
Begelman et al. (2022) analyzed a prograde disk, which is somewhat smaller than the disks we analyze here.They reported that the gas pressure evolves ∝r −2 , while the magnetic pressure was found to evolve faster closest to the BH.Here, we find that the gas-pressure evolution is somewhat steeper than r −2 .This could be due to the different disk structure or the different resolution we are using or to the shorter time interval over which the averaged is performed.We note that a similar radial evolution of the gas pressure, namely p g ∝ r −2 was reported by Tchekhovskoy et al. (2011) andMcKinney et al. (2012), who found that p g ∝ r −1.9 for their "thinner" disk models.
We further show in Figure 12 the radial dependence of all individual magnetic pressure components.The total magnetic pressure has a similar radial evolution in our work and in the work of Begelman et al. (2022), namely, its radial evolution is steeper than r −2 close to the BH.As is seen from Figure 12, the polar component q b 2 never dominates for any of the spins and is much lower than the radial and toroidal components.Second, at small radii inside the ISCO, the main difference between rotating and nonrotating BHs is the contribution of the toroidal field f b 2 : it is negligible compared to the radial component for a nonrotating BH, while both components are on the same order for rapidly rotating BHs.The toroidal field even dominates for the prograde disk very close to the BH.This explains the difference between the conclusions of Begelman et al. (2022), who studied disks around rapidly rotating BHs, and Chatterjee & Narayan (2022), who studied nonspinning BH systems, which is the importance of the toroidal component in a MAD Figure 11.Left Panel: 1σ region of the jet boundary for a = −0.94,0, and 0.94.We use the magnetization parameter σ = 1 condition as the jet-disk boundary and derive the 1σ variations and median of the polar angle θ for σ = 1 at a given radius r.The vertical lines are 1σ variations, and the dots are the median values.Right Panel: KDE of the polar angle θ such that σ = 1 at radius r = 10r g .The nonspinning BH has the widest jet, corresponding to the largest θ, while the prograde disk with a = 0.94 has a wider jet than the retrograde disk with a = −0.94.In addition, the nonspinning BH has the narrowest peak, which suggests that the jet boundary of a nonspinning BH is the most stable.accretion process.A similar radial evolution of the magnetic field components was reported by Tchekhovskoy et al. (2011) and McKinney et al. (2012), who found that b r ∝ r −1.5 and b f ∝ r −1 for their "thinner" disk models.Clearly, the radial dependence of b r does not seem to depend on the BH spin, while the toroidal component does.
We show the fand time-averaged maps of the pressure components of aM94b100, a0b100, and a94b100 in Figure 13.All panels use the same color scaling, allowing an easy comparison between them.At small radii close to the BH, the radial component b r 2 is large for all three simulations and decreases quickly as the radius increases.For a = ± 0.94, high values of b r 2 are also found at large radii r ∼ 10r g in the jet region.The polar component, q b 2 , is negligible in all three cases.Note that the blank regions in the color maps of q b 2 are caused by numerical truncating errors: this component is too small relative to the other two components.The distributions of the toroidal component f b 2 also shows a great variance.The nonspinning BH has the weakest f b 2 component, while it is stronger for the spining BHs.We note the drop in |b f | 2 at the equator: this is because b f changes sign at the equator due to the presence of a current sheet that separates the two hemispheres.This drop is more pronounced and visible for the spinning BHs.
The gas pressure p g of the nonspinning BH contains clear weak channels about  q p 3 and  q p 2 3 .These channels are also visible for the spinning BHs, but are less pronounced.As we previously discussed in Section 3.7, these channels are caused by the numerical floors applied in the regions closest to the pole, which are required to maintain numerical stability in our simulations.

Angular Momentum Flux
It has long been argued that magnetic fields drive the transfer of angular momentum during accretion via the development of the MRI (Balbus & Hawley 1991, 1998), but also by helping to launch a disk wind and producing a strong jet, both components capable of significantly transporting angular momentum.In this section, we investigate the dependence of angular momentum flux on the spin of the BH.We use the definitions in Chatterjee & Narayan (2022) to define the total angular momentum flux  J total , the advected angular momentum flux  J adv , and the angular momentum flux due to stresses  J stress , where 〈X〉 j represent the average of X with respect to variable j.Note that for this section, there is no weighting by the density ρ.We further decompose the stressed induced angular momentum flux into its Maxwell  J stress,M and Reynolds  J stress,R components, Following these definitions, the time-and f-averaged angular momentum flux in the radial direction of a94b100, a0b100, and aM94b100 are shown in Figure 14 at three different radii, namely 10r g , 15r g , and 20r g .For a nonspinning BH, depicted in the middle line of Figure 14, the results we obtained are similar to the results reported in Chatterjee & Narayan (2022).The radial component of the total angular momentum flux  J r total is negative in the disk region at all sampled radii, r = 10r g , 15r g , and 20r g .This indicates that as the matter accretes, it brings a net angular momentum to the BH.However, at θ ; 0.35π and 0.65π, the total angular momentum flux changes sign to become positive.This shows that angular momentum is being transported away from the accretion disk by the wind.The radial angular momentum flux converges to 0 as it approaches the poles, underlying the weakness of the jet for a = 0.In the disk, the contributions of  J r adv and  J r stress,R dominate and are always negative, suggesting that they are responsible for the inward  2 (first column), q b 2 (second column), and f b 2 (third column), and of the gas pressure (last column) for a = −0.94(top), a = 0 (middle), and a = 0.94 (bottom).The color scaling is the same for all spins and all components to display the differences better.The most striking difference is the small contribution of the azimuthal component b f to the total magnetic pressure for the nonspinning BH a = 0. We note that the white region in the second column for q b 2 corresponds to regions in which <  transport of angular momentum in the disk, as was previously demonstrated by Chatterjee & Narayan (2022).We note that the relative importance of the Reynolds components is high at small radii, but decreases toward large radii.This agrees with the results of Chatterjee & Narayan (2022) at r = 20r g .Furthermore, the Maxwell stress component  J r stress,M is always positive with two maxima, which are reached at θ ; 0.4π and θ ∼ 0.6π.At these high latitudes, the radial component of the total angular momentum flux is dominated by J r stress,M underlying the importance of magnetic fields in the production and dynamics of the wind above the disk in the transport of angular momentum.
The results of the BH with positive spin a = 0.94 are displayed in the top line of Figure 14.We find this case to be qualitatively similar to the nonspinning BH scenario: the radial component of the total angular momentum flux is negative in the disk and positive at high latitudes, with a change in sign at θ ∼ 0.4π and θ ∼ 0.6π.Outside the disk, we find that the Maxwell component strongly dominates and reaches its maxima symmetrically with respect to the equator at θ ∼ 0.25π and θ ∼ 0.75π.The radial components of the angular momentum flux remain large up to the poles.This indicates the existence of powerful magnetized winds and jets, through which the disk and the BH lose a substantial fraction of their angular momentum.We further see that the angular momentum flux contributed by advection is negative in the disk region.However, it is positive at high latitude, which is different from the case of the nonspinning BH.This may be due to the large contribution of the magnetic field energy density b 2 /2 in its expression.Finally, we find that the Reynolds stress component is much smaller than any of the other components and can be safely ignored.
For the negative-spin simulation a = −0.94displayed in the bottom line of Figure 14, we find that the radial angular momentum flux distribution with angle is significantly distinct from that of the nonrotating and positive-spin BHs.This is mostly due to the fact that the poloidal velocity changes sign: the plasma corotates with the BH close to the pole, while it corotates with the disk at the equator.First, we find that the radial component of the total angular momentum flux is always negative at all angles.As a result, the BH gains angular momentum and its spin decreases toward 0. This is consistent with the result presented in Figure 1.The contribution of the advection to the flux of angular momentum is similar to the cases of positive and null spins, namely, its contribution is maximum and negative at the equator.The Maxwell component contribution is positive at the equator, showing that magnetic fields are responsible for taking away radial angular momentum also in this case.However, this contribution is negative around θ ; 0.2π and 0.8π and dominates the radial total angular momentum flux at these angles.
The signs of angular momentum flux are mainly determined by the sign of the azimuthal velocity u f and of b r b f .In Figure 15 we show the fand time-averaged evolution of u f and of b r b f with the polar angle θ at the three radii selected for the angular momentum flux in Figure 14.For the a = 0 and a =0.94 cases, the toroidal velocity u f is positive for all angles θ as expected, and therefore, the sign of the advection contribution to the angular momentum flux  J adv depends on the sign of u r .On average, u r is negative in the disk region because the matter is accreted toward the BH.It is positive at high latitude, however, as matter is carried away in the disk wind and the jet.This change in sign explains the pattern of  J adv seen in Figure 14 for a = 0.94.For the nonspinning BH, the azimuthal velocity quickly goes to 0 at high latitude, which explains the null contrbution of the advection to the angular momentum flux in this region.On the other hand, the term b r b f is always negative at all angles θ for the BH with spin a = 0.94.In this case, the maxima (one in each hemisphere) of |b r b f | are reached closer to the poles than for the nonrotating BH, further underlying the different disk structure: rotating BHs have wider disks than nonrotating ones, as was already demonstrated in Section 3.7 and in Figures 10 and 11.The polar angles θ at which the maxima of b r b f are reached correspond to the angle at which the contribution of the Maxwell stress is maximum.
The actual sign of the component b r b f depends on our arbitrary choice of initial disk magnetization, and specifically, on the orientation of the initial magnetic field loop.A different orientation would lead to a different sign of the b r b f component.Indeed, the sign of b f is set by frame-dragging, and it is positive in the southern hemisphere and negative in the northern hemisphere for prograde disks.However, the sign of b r would be opposite if the initial magnetic field loop had a different orientation, as demonstrated in McKinney et al. (2012).We note here another difference between the a = 0 and a = 0.94 BHs: b r b f is negative close to the poles for the rotating BH, while it is very small for the nonrotating BH.
The situation is different for the retrograde disk, with a different evolution of the azimuthal velocity u f and of the term b r b f with polar angle θ, both shown in the third line of Figure 15 for a = − 0.94.The toroidal velocity u f is also positive in the disk region.However, it becomes negative close to the poles around θ ; 0.2π and θ ; 0.8π.The term b r b f is also different from that of the prograde disk.First, it has a different sign close to the pole: b r b f is positive for a retrograde disk since (i) b f is negative in the southern hemisphere and positive in the northern hemisphere because of frame-dragging, and (ii) b r is negative in the southern hemisphere and positive in the northern hemisphere.Here, the sign of b r is set by the orientation of the initial magnetic field loop.At the equator, however, b r b f has the same sign for both prograde and retrograde disks.This explains the sign of the Maxwell stress contribution to the radial angular momentum flux.
In order to understand how angular momentum is transported throughout the disk, we show in Figure 16 color maps of the angular momentum flux modulus, for the total, advection, and Maxwell stress components and the associated streamlines.The top line shows the total angular momentum flux, while the middle and bottom lines show the contribution of the advection and Maxwell stress, respectively.The left, middle, and right columns are for spins a= − 0.94, 0, and 0.94 respectively.We note that the scale of the colorcoding is the same in each line to facilitate comparison between the different BH spins.Furthermore, we point out that the equilibrium radius for these simulations is around r = 30r g , prompting caution for larger radii.First, it is clear that each of these figures is antisymmetric with respect to the equator, as expected.We note that these figures are obtained directly from the data from our simulation without imposing symmetry, as was done by Chatterjee & Narayan (2022).From the top line, we further see that for all three spins, angular momentum is lost by the disk through the disk wind.
There is one more difference between those figures: the angular momentum flux at the pole is negative for a = −0.94 and a = 0, while it is positive for a = 0.94.The contribution to the total angular momentum flux in this region also appears to be small compared to the flux at the equator, as seen in Figure 14.In addition, Figure 16 shows the presence of a transition layer in the simulation with a negative and a null spin.In this transition layer, the angular momentum flux changes direction, being oriented outward in the wind and inward in the jet.The position of this transition layer is at a larger angle from the pole for a = − 0.94 than for a = 0.This transition layer also appears in the Maxwell stress contribution.
The streamlines of the Maxwell stress contribution to the total momentum flux are shown in the last line of Figure 16.We find two interesting features.First, both the prograde and retrograde disks display a large Maxwell stress contribution to the angular momentum flux in the jet region, in contrast to the nonspinning BH.This underlines that the jets of a nonrotating BH in the MAD regime are weak and do not transport a substantial amount of angular momentum (and in fact, energy) to their surrounding environments.The largest outward contribution in the nonrotating BH comes from the magnetized wind.It is also clear that the jet of the prograde disk is stronger than that of the retrograde BH and faster deposits angular momentum into its environment.We further see a substantial Maxwell stress contribution to the wind, which shares the same characteristics as the total angular momentum flux, meaning that it is weaker for retrograde disks.Finally, we find that for the retrograde disk, a region of nearly zero Maxwell stress contribution separates the wind region from the disk region.This transition is associated with the change in sign of the toroidal velocity u f , shown by the red line in Figure 16 for a = −0.94.It is clear that the red line follows the transition region.

Conclusion
In this work, we performed several GRMHD simulations of thick accretion disks in the MAD regime around BHs characterized by different spin a with cuHARM.Our key results can be summarized as follows: 1. We studied the angular momentum flux for our simulations with a = −0.94,0, and 0.94 and underlined the differences.We found that (i) a substantial amount of angular momentum is transported away in the magnetized wind and that the Maxwell stresses are larger for a rotating than for a nonrotating BH, as expected because of frame-dragging.In fact, the amount of angular momentum transported by the jet for the nonrotating BH is small and negligible.These results were provided in Section 4 and are clearly displayed in Figure 16. 2. We did not find any correlation between the mass accretion rate and the MAD parameter.However, we did find an anticorrelation between their time derivatives, displayed for our simulations aM94b100 and a94b100 in Figure 7.We provided a heuristic explanation for this result in Section 3.4.3. We underlined the difference in the magnetic field component strengths for spinning and nonspinning BHs in Section 3.8.From Figures 12 and 13, it is seen that the θ component is always subdominant, while the relative importance of the f component depends on the spin of the BH.In the nonrotating case, the toroidal component is negligible compared to the radial component close to the horizon, while these components are comparable for the rotating BH.The toroidal component even dominates very close to the horizon for a = 0.94.We therefore recover the results from Begelman et al. (2022) and Chatterjee & Narayan (2022) on the importance of the toroidal component in regulating the accretion, and confirmed that the difference has its origin in the BH spin.4. We underlined the differences in the structure of the disks and the jets in the MAD state as a function of spin in Section 3.7.In particular, we found that retrograde disks are wider than the corresponding prograde disks, in agreement with the findings of Narayan et al. (2022).
Correspondingly, the jets of prograde disk are narrower than the jets of retrograde disks.These results are displayed in Figure 11. 5. We studied the MAD parameters and investigated their dependence on the spin and the initial magnetic field strength in Sections 3.2 and 3.6.We found that our numerical results agree very well with those of Narayan et al. (2022), namely, the MAD factor increases with spin until a = 0.5, after which it decreases.The fitting formula from Narayan et al. (2022), which we renormalized to account for the difference in definitions, accurately describes the distribution of Φ B with spin.Therefore, this result holds across resolution, simulation duration, and numerical method.The results of this analysis are summarized in the top row of Figure 4. We find that the MAD parameter does not depend strongly on the initial magnetic field strength parameter β 0 .It only takes longer for the simulation to reach the MAD state, namely, for the magnetic field to saturate to its final value.When the MAD state is reached, all simulations with different β 0 share the same characteristics, as shown in Figure 8.
6.We attempted to identify a characteristic variability time in the MAD regime by studying the temporal variation of Φ B via the Fourier transform in Section 3.3.In contrast to the 2D case (see, e.g., Chashkina et al. 2021), no clear period was unambiguously identified.This is because in a 3D simulation, accretion proceeds via nonaxisymmetric instabilities, such as the interchange instability (Spruit et al. 1995;McKinney et al. 2012;Begelman et al. 2022).
The typical variability time we find is around a few hundred to some thousand t g , which is comparable to the estimates from Lloyd-Ronning et al. (2016) and James et al. (2022).This timescale is also comparable to the timescale inferred by Wong et al. (2021), who studied the layer between the disk and the jet.
Our results shed light on the differences in the accretion dynamics of disks in the MAD state across spins, prompting a more detailed analysis of the transport of angular momentum by jets and winds, as well as of the the role of the toroidal magnetic field in (i) shaping the disk and jet, and (ii) providing the norm of the MAD parameter, and via the norm understanding the efficiency of conversion between accretion luminosity and BH spin energy deposited in the Poynting jets.

Figure 2 .
Figure 2. Time evolution of the mass accretion rate M .The left column pertains to retrograde disks, while the right columns corresponds to prograde disks.For convenience, we set the y-axes to identical scales to facilitate direct comparison.The horizontal lines are the time-averaged mass accretion rate  M , where the averaged is taken from 10 4 t g to 2 × 10 4 t g .The last row shows a discrepancy for a = 0 between our low-(a0b100) and high-resolution (a0b100h) simulations.It shows, in this case, that the resolution is too low to resolve the mass accretion rate accurately.However, we found that the MAD parameters are nearly identical, indicating that the properties of the MAD state should be similar.We provide more details on this issue at the end of Section 3.2.

Figure 4 .
Figure 4. Spin dependence of several disk and jet properties.Top left panel: mass accretion rate M at the horizon for different spins.The blue points are the results of low-resolution simulations, while the red points are the high-resolution simulations.Top middle panel: normalized MAD parameters for different spins.Plotted on top (dashed green line) are the fit results ofNarayan et al. (2022).Their formula is divided by p 4 2 to be consistent with our definition of the MAD parameter (see Section 3.2 for more details).Top right panel: jet efficiency as a function of spin.The green points are the predicted efficiencies, and the blue (red) points are the efficiencies derived from our simulation (red shows the higher resolution).The error bars represent the 1 sigma temporal variation of each quantity.η > 1 represents an efficiency higher than 100%.Bottom left: dependence of the slope( ) ( )  º F k d dt dM dt

Figure 5 .
Figure 5.The long-time evolution of the mass accretion rate and MAD parameters for simulation aM94b100h.The MAD parameter remains steady with time, although the mass accretion drops.

Figure 12 .
Figure12.The angular-and time-averaged gas pressure p (yellow) and magnetic pressure b 2 (blue) for aM94b100, a0b100, and a94b100.We further show each component of the magnetic pressure, b r 2 (red), q b 2 (green), and f b 2 (purple).These values are calculated by Equation (18) and are then time averaged.The dashed (light blue) curve represents the r −2 dependence.

Figure 13 .
Figure 13.Time and azimuthal averages of the magnetic pressure components, namely, b r2 (first column), q b 2 (second column), and f b 2 (third column), and of the gas pressure (last column) for a = −0.94(top), a = 0 (middle), and a = 0.94 (bottom).The color scaling is the same for all spins and all components to display the differences better.The most striking difference is the small contribution of the azimuthal component b f to the total magnetic pressure for the nonspinning BH a = 0. We note that the white region in the second column for q b 2 corresponds to regions in which < truncation errors in the calculation to produce this figure that arise because |b θ | = |b r |, |b f |.

Figure 14 .
Figure 14.Radial angular momentum flux components  J r as a function of the polar angle θ for a94b100 (top), a0b100 (middle), and aM94b100 (bottom) at three different radii, r = 10, 15, and 20r g .The solid blue line corresponds to the total radial flux  J r total , and the red and green dashed lines are for the advection  J r adv and for the Maxwell stress  J r stress,M contributions to the radial flux, respectively.Finally, the dotted purple line corresponds to the Reynolds stress contribution  J r stress,R .

Figure 15 .
Figure 15.Toroidal component of the four-velocity u f and b r b f as a function of θ for a94b100 (top), a0b100 (middle), and aM94b100 (bottom) at the radii at which the radial angular momentum flux is obtained.

Figure 16 .
Figure 16.Time-and faveraged maps of the angular momentum flux of the simulations with spin a = −0.94(left), a = 0 (middle), and a =0.94 (right).The color map represents the modulus of the angular momentum flux, ( )( )   + q J J r 2 2. The first row shows the total angular momentum flux, while the second and third rows show the contributions of the advection and Maxwell stress components, respectively.The color-coding for each row is the same to facilitate comparison.For a = − 0.94, we use the red lines to denote the transition radius where the toroidal velocity component u f changes sign.It is seen that its position corresponds to the region in which the radial component of the Maxwell stress changes sign as well.

Table 1
List of the Simulations Presented in This Paper with Their Initial Magnetization β 0 , Spin a, and Resolution