Rapid Black Hole Spin-down by Thick Magnetically Arrested Disks

Black hole (BH) spin can play an important role in galaxy evolution by controlling the amount of energy and momentum ejected from near the BH into the surroundings. We focus on radiatively inefficient and geometrically thick magnetically arrested disks (MADs) that can launch strong BH-powered jets. With an appropriately chosen adiabatic index, these systems can describe either the low-luminosity or highly super-Eddington BH accretion regimes. Using a suite of 3D general relativistic magnetohydrodynamic simulations, we find that for any initial spin, an MAD rapidly spins down the BH to the equilibrium spin of 0 < a eq ≲ 0.1, very low compared to a eq = 1 for the standard thin luminous (Novikov–Thorne) disks. This implies that rapidly accreting (super-Eddington) BHs fed by MADs tend to lose most of their rotational energy to magnetized relativistic outflows. In an MAD, a BH only needs to accrete 20% of its own mass to spin down from a = 1–0.2. We construct a semi-analytic model of BH spin evolution in MADs by taking into account the torques on the BH due to both the hydrodynamic disk and electromagnetic jet components, and find that the low value of a eq is due to both the jets slowing down the BH rotation and the disk losing a large fraction of its angular momentum to outflows. Our results have crucial implications for how BH spins evolve in active galaxies and other systems such as collapsars, where the BH spin-down timescale can be short enough to significantly affect the evolution of gamma-ray emitting BH-powered jets.


INTRODUCTION
It is well established that supermassive BHs (SMBHs) ranging in mass from 10 6 to 10 9 M ⊙ lie at the center of nearly every galaxy.Active galactic nuclei (AGN) are powered by the accretion of gas onto SMBHs.AGN release energy in the form of radiation, winds, and under the right circumstances, relativistic jets.The result is AGN feedback, and it is thought to play a key role in galaxy evolution (Ciotti & Ostriker 2007;Kormendy & Ho 2013).Scaling relations between SMBHs and properties of their host galaxies (e.g.stellar velocity dispersion, Ferrarese & Merritt 2000;King 2003;King & Nealon 2021;bulge mass, Marconi & Hunt 2003;Häring & Rix 2004; and total stellar mass, Reines & Volonteri 2015) Corresponding author: Beverly Lowell beverlylowell@u.northwestern.edu* Released on February, 1st, 2023 are powerful tools to understand the coevolution of BHs and galaxies.
Massive jets heat and ionize the gas in and around galaxies, preventing cooling and stifling star formation, and sweep gas out of the galactic bulge (Silk & Rees 1998;Di Matteo et al. 2005).Moreover, the energy deposited in the surrounding environment affects the fueling of the BH itself, consequently regulating its duty-cycle.
The outflows are likely powered in part by the accretion of gas onto the SMBH.The Eddington limit is an important reference point in understanding BH growth, and thus galaxy structure and evolution.Estimations of the Eddington ratio from observations are largely uncertain.Most AGN accrete below the Eddington luminosity (L Edd ), but there is evidence that some accrete at super-Eddington rates.In the super-Eddington regime, the flow becomes radiationdominated, and the accretion disk is geometrically-thick.While many details of this accretion regime remain unclear, super-Eddington accretion is expected, and even inevitable, in some circumstances.Super-Eddington accretion can be observed in relatively local AGN (Lanzuisi et al. 2016;Tsai et al. 2018) and tidal disruption events (TDEs, Rees 1988); it can be inferred in ultraluminous X-ray sources (ULXs, Sutton et al. 2013;Kitaki et al. 2017;Walton et al. 2018) and gamma-ray bursts (GRBs, Paczynski & Xu 1994).
Whereas in low redshift AGN super-Eddington accretion episodes are rare, the conditions for super-Eddington accretion were favorable at higher redshifts (Aird et al. 2018).The gas content of early galaxies correlates with the accretion rate of the SMBH, implying that a BH will accrete above the Eddington rate given sufficient gas supply (Izumi 2018).Furthermore, there is no theoretical evidence that super-Eddington accretion is not possible (Szuszkiewicz 2004;Fragile et al. 2012;S ądowski et al. 2014;Jiang et al. 2019).
Quasars at up to z ≈ 7 have been found to have BH masses of 10 9 M ⊙ (Mortlock et al. 2011;Yang et al. 2021).To grow to such a large mass at that redshift from a reasonable BH seed, a BH would need to spend at least some time accreting at or above the Eddington rate (Volonteri & Rees 2005;Volonteri et al. 2015;Begelman & Volonteri 2017).In fact, super-Eddington accretion episodes may be the most reasonable way to achieve 10 9 M ⊙ SMBHs by z ≈ 7, as maintaining accretion rates at or below the Eddington rate at the highest redshifts might require tuned conditions, although models that grow M BH = 10 9 M ⊙ SMBHs by z = 9 through mergers of 10 4 stellar-mass BHs, which formed at z ∼ 20 and had accreting most of their life at an Eddington rate since their formation, have been constructed (Yoo & Miralda-Escudé 2004).More recently, James Webb Space Telescope has revealed SMBHs already at redshift z = 10.1 (Bogdan et al. 2023), challenging the models of SMBH formation and growth (Goulding et al. 2023), and suggesting either the existence of massive BH seeds or highly super-Eddington accretion at such high redshifts.
For a given luminosity, the lower the BH mass, the less gas reservoir is needed to achieve the Eddington rate.It follows that high accretion rates will more likely be found in galaxies hosting low-mass SMBHs.Cosmological simulations from Anglés-Alcázar et al. (2017) showed that underweight BHs (those that lie below the M BH − M bulge relation) go through super-Eddington accretion episodes until the bulge reaches a critical mass, and the accretion rate then levels off and the BHs join the correlation.The short timescales involved with these episodes make them difficult to be observed.Similarly, the tightly-correlated M − σ relation may not represent all local AGN BH masses, but is instead an upper bound on BH mass and represents BHs accreting at or below the Eddington rate (King 2010).If that is the case, BHs that lie below the M − σ relation transiently accrete at super-Eddington rates: such outbursts may last a rather short time and thus are diffi-cult to catch in the act.Even if we were to catch one of these BHs feasting on gas at a super-Eddington rate, the observed luminosity would be likely limited to just a few Eddington luminosities due to the low radiative efficiency of super-critical accretion (Abramowicz et al. 1988).
An example of this may be the Narrow line Seyfert 1 (NLS1) galaxies, which are thought to host central BHs with masses lower than other Seyfert SMBHs.NLS1s are a subclass of broadline Seyfert galaxies.These are AGN detectable through both the active core and the radiation from the host galaxy.They are characterized by extreme multiwavelength properties, including narrow Balmer lines, strong Fe II emission, and super-soft X-ray emission with extreme variability.NLS1s are thought to accrete at close to Eddington and even super-Eddington rates in some cases (Pounds et al. 1995;Mathur 2000;Komossa et al. 2006;Jin et al. 2016Jin et al. , 2017)).In particular, the bolometric luminosity inferred for what is believed to be the most super-Eddington NLS1, RX J0134.2-4258, is firmly in the super-Eddington regime, L bol ∼ (6 − 16)L Edd (Jin et al. 2022(Jin et al. , 2023)).Furthermore, NLS1s feature steep X-ray spectra that are similar to highredshift quasars with super-Eddington luminosities.Indeed, NLS1s could be the early-stage of AGN evolution, specifically the low-redshift, low-luminosity versions of quasars that are actively undergoing super-Eddington accretion.
The centers of many galaxy clusters are observed to be hot, despite their central cooling timescales being much shorter than their lifetimes.BH feedback via jets could be the source of heating that is responsible for avoiding catastrophic cooling (see review by Gitti et al. 2012).Centered on the galaxy NGC 1275 is the Perseus cluster, the brightest known X-ray cluster.It is believed that jets inflate bubbles that displace the intracluster medium (ICM) to the North and South of the central BH.These jets, powered by the accretion and spin of the BH, may be the energy source for the observed luminosity (Sanders & Fabian 2007), by e.g.driving the turbulence that heats the gas (Zhuravleva et al. 2014).
The galaxy cluster MS 0735 contains an AGN with outflows that extend far beyond the size of the host galaxy, plowing through the ICM and forming large cavities.McNamara et al. (2009) argue that because the energy of the outflows is comparable to the rest mass energy of the SMBH, the outflows may be powered in part by the rotational energy of the BH.For example, the rotational energy of a maximallyspinning 10 9 solar mass BH is over 10 62 erg.This is large compared to the thermal energy of the surrounding X-ray atmosphere, enough to quench a cooling flow for several Gyrs.Accretion power alone is likely not enough to power the jets that heat the ICM, and therefore, they must be powered in part by BH spin energy.It is now generally accepted that spin power is coupled to the AGN feedback loop (McNamara et al. 2009).
It is widely thought that the rotational energy of the BH powers relativistic jets with the help of large-scale magnetic fields through the Blandford-Znajek (BZ, Blandford & Znajek 1977;Penrose 1969;Lasota et al. 2014) process.Naturally, the jet power depends on the value of BH spin, scaling approximately as P jet ∼ Φ 2 a 2 (MacDonald & Thorne 1982;Thorne et al. 1990;Narayan & McClintock 2012), where Φ is the magnetic flux on the BH.This scaling is consistent with observations of blazars (Chen et al. 2021) and X-ray binaries (Narayan & McClintock 2012;Steiner et al. 2013;McClintock et al. 2014;but see Fender et al. 2010).
Radio loudness, which is the total radio luminosity (core plus extended emission) normalized by the B-band nuclear luminosity, shows an apparent spread by factor of 10 3 between the radio-loud and radio-quiet AGN (Sikora et al. 2007).Differences in the accretion rate alone would not be sufficient to explain the dichotomy, because at the same optical luminosity, which is presumably the tracer of accretion rate, the radio-loudness appears to exhibit the dichotomy.One possibility to explain this dichotomy could be that the radio-loud AGN have much more rapidly-spinning BHs than radio-quiet AGN (Sikora et al. 2007;Tchekhovskoy et al. 2010;Martínez-Sansigre & Rawlings 2011;Wu et al. 2011).In this so-called "spin paradigm" the value of BH spin controls the jet power (Wilson & Colbert 1995;Blandford 1999), as opposed to the "accretion paradigm" (Blandford & Rees 1992) where BH accretion rate determines the spectral state of the accretion flow and the jet power.However, in order to achieve the full range of 10 3 in radio luminosity requires a range of a factor of 10 1.5 ≃ 30 in the spin, implying that radio-quiet AGN would have to have a spin ∼ 0.03, which is increasingly difficult to achieve unless the jet power dependence on spin can be steeper than a 2 (Tchekhovskoy et al. 2010).
A BH spin will evolve due to the torques acting on the BH, and over time, the spin will evolve toward an equilibrium value, a eq .BH spin evolution has been studied for various disk configurations (Bardeen 1970;Thorne 1974;Popham & Gammie 1998;Gammie et al. 2004;Narayan et al. 2022).Early 2D simulations by Gammie et al. (2004) found an equilibrium spin of a eq ≈ 0.93, lower than a disk with no large-scale magnetic field.Tchekhovskoy et al. (2012); Tchekhovskoy (2015); Narayan et al. (2022) found a much lower value of a eq than Gammie et al. (2004), and the effect is due to the presence of strong magnetic fields.
Thanks to recent advances in GRMHD global simulations of accretion disks, it is becoming clear that in the presence of a large amount large-scale vertical magnetic flux, accretion disks naturally end up in a magnetically arrested disk (MAD) state (Tchekhovskoy et al. 2011).Indeed, the largescale magnetic field threading the disk tends to be dragged inward by the accretion flow.Eventually, the magnetic flux saturates within the inner regions of the disk.This saturated state, the MAD state, is a robust result of GRMHD simulations of accretion disks.Furthermore, it can also occur when the initial magnetic field strength is weak (Jacquemin-Ide et al. 2021) or the magnetic flux is small (Christie et al. 2019).
Evidence is building that some AGN contain MADs with highly efficient jets (Zamaninasab et al. 2014;Zdziarski et al. 2015;Nemmen & Tchekhovskoy 2015;Tchekhovskoy 2015;Event Horizon Telescope Collaboration et al. 2021;Chen et al. 2021;Event Horizon Telescope Collaboration et al. 2022).GRMHD models of MADs with high BH spin produce strong jets (Tchekhovskoy et al. 2011;McKinney et al. 2012).Thus, modeling MAD accretion systems could allow us to study the effect the powerful jets have on the BH spin and, hence, jet power.In particular, maximizing jet power, as in MADs, extracts the maximum possible rotational energy from the BH.
In this paper we use the simulations of Tchekhovskoy et al. (2011); Tchekhovskoy et al. (2012) to investigate the dynamical differences between the MAD and the geometrically-thin Novikov-Thorne disk (Novikov & Thorne 1973), including the impact of jets and their underlying magnetic fields on the evolution of BH spin in time.Whereas the main application of this work is to highly super-Eddington MAD flows, we do not explicitly include radiation transport to model highly super-Eddington accretion disks.In such flows the radiation and gas cannot move relative to each other: effectively, they behave as a single fluid with a polytropic index γ = 4/3 characteristic of radiation-dominated ideal gas.Thus, we adopt an ideal gas equation of state with a polytropic index of Γ = 4/3.Future work we will include on-the-fly radiation transport to explicitly account for the effect of radiation feedback on the gas (Lowell et al. 2024).In Section 2 we give the equations of BH spin evolution and describe the numerical set-up for MAD simulations, where magnetic flux on the BH and jet power are maximized.In Section 3, we present our results and derive the semi-analytic model to reveal the physics behind the spin evolution in MADs.In Section 4, we discuss our results and conclude.

Spin evolution
BHs are defined by three parameters: electric charge Q, angular momentum J, and mass M. We assume that BHs are formed from charge-neutral matter, so we work in the Kerr vacuum metric with two parameters, J and M.
We define dimensionless BH spin as where J is the BH angular momentum, and M is the BH mass.
Here and below, we use geometrical units such that c = G = 1, and assume that the spin vector is parallel or anti-parallel to the angular momentum vector of the accreting gas.
The change in BH mass-energy M and angular momentum J are given by dM = e in dm (2) and dJ = l in Mdm (3) (Bardeen 1970), where dm is the amount of accreted rest mass, and e in and l in are the the specific energy and angular momentum, respectively, of that mass.
Combining equations ( 2) and ( 3) gives an ODE that describes spin evolution as a function of BH mass: where log is a natural logarithm.Equations ( 2), (3), and (4) can be solved for the dependencies a(m) or a(t).Gammie et al. (2004) defined the dimensionless spin-up parameter, to describe BH spin evolution due to accretion: s is positive when the BH is spinning up, negative when the BH is spinning down, and zero when the BH has reached the equilibrium spin, a eq .Using equation (5) we recast equation (4) to find an expression for s as a function of the specific energy and angular momentum fluxes, Here, the angular momentum is positive when the disk is prograde and negative when the disk is retrograde relative to the BH spin.Thus, the accretion of angular momentum l in will lead to spin-up (spin-down) for prograde (retrograde) disks.Accreted energy is always positive, so the second term spins down (up) the BH for positive (negative) spin.The terms in equation ( 6) are intuitive in that the accretion of positive accreted angular momentum spins up the BH (increases a), but the accreted rest-mass spins down the BH.

Standard thin disk
In a "standard", razor-thin, radiatively-efficient disk model, or Novikov-Thorne (hereafter NT) model Novikov & Thorne (1973), gas moves in Keplerian orbits and experiences viscous forces that cause the gas to lose angular momentum and migrate inward to smaller radii until it reaches the inner-most stable circular orbit (ISCO), also known as the marginally-stable orbit.Inside of the orbit, of radius R ms , the gas plunges into the event horizon.In this case the inner radius of the disk is R in = R ms .The equations for e in , l in , and R ms describing the NT disk spin evolution are taken from Moderski & Sikora (1996).Analytic expressions for e in and l in can be found in Appendix A.
Physically, the primary origin of viscosity is thought to be the magnetorotational instability (MRI, Balbus & Hawley 1991;Balbus & Hawley 1998), driven by the magnetic field in the disk.The NT disk model describes the high/soft, or thermal, spectral state of accretion that occurs when the luminosity is 1 − 10% of L Edd .The model neglects the large-scale poloidal magnetic field: relativistic jets do not form and do not play any role in the spin evolution.The BH spin changes only due to mass accretion onto the BH.This can only increase the angular momentum and energy of the BH, resulting in the spin-up to the maximum spin, a eq = 1.In this work, we ignore the accretion of low angular momentum photons in the NT disk model: including them would limit the BH spin-up to a eq = 0.998 (Thorne 1974).

Simulations of Magnetically Arrested Disks
Unlike for the NT disk, there is no analytic model to describe MAD spin evolution.Therefore, numerical simulations of MADs must be used.We use non-radiative 3D GRMHD simulations run with the code HARM (Gammie et al. 2003;Noble et al. 2006) of BH accretion in the MAD regime with BH spins −0.9, −0.5, −0.2, 0.0, 0.1, 0.2, 0.5, 0.9, and 0.99 (Tchekhovskoy et al. 2011;Tchekhovskoy et al. 2012).In order to make inferences for super-Eddington accretion, we adopt an ideal gas law equation of state with a polytropic index of Γ = 4/3, characteristic of a radiationdominated gas.The initial conditions are such that an accretion disk will form with scale height h/r ≈ 0.3, and enough magnetic flux Φ BH will accumulate on the BH to obstruct the accreting gas and lead to a MAD.The simulation parameters are given in Table 1.Because BH mass and spin evolve over cosmological time, much longer than the duration of any GRMHD simulation, we "stitch" together simulations for different BH spin values to represent the BH spin evolution over cosmological time.Our simulations reach the quasi-steady MAD state by t ∼ 10, 000r g /c.From here we time-average quantities to the end of each simulation, as shown in Table 1, to reduce the effect of fluctuations due to magnetic flux eruptions.The simulations reach inflow equilibrium out to (20 − 50)r g within the time-average, where the distance depends on the simulation.
We calculate the specific angular momentum and energy fluxes from the simulations using the stress-energy tensor, defined as The hydrodynamic component is given by and is used to calculate specific fluxes e in and l in on the BH from the disk.Here, ρ is the gas density, u is the specific internal energy and p is thermal pressure of the gas, and u µ is the gas contravariant 4-velocity.Electromagnetic (EM) fluxes are calculated using the electromagnetic component of the stress-energy tensor, given by where is the magnetic pressure and b µ is the comoving magnetic field 4-vector.The mass, specific energy, and specific angular momentum fluxes are calculated as and where dA θϕ is the area element in the θ − ϕ plane.The torques on the BH are due to accreted fluxes, so in our simulations we calculate the fluxes at the BH event horizon, , where r g = GM/c 2 is BH gravitational radius.The total, electromagnetic, and hydrodynamic specific angular momentum fluxes are calculated at r H as l in = l(r = r H ), as is the electromagnetic specific energy flux.
However, when calculating hydrodynamic fluxes we must account for numerical floors.Realistically, there is no mass in the jet, but numerical MHD schemes do not work in a vacuum.To handle this, when mass or internal energy densities are very small, below a set floor value, internal energy or mass is added to the grid near the BH horizon to keep the numerical scheme stable.When calculating the mass, hydrodynamic energy and angular momentum fluxes, we eliminate contributions from the highly magnetized regions affected by the floors using a magnetization cutoff condition.This method gives the physical values at the BH horizon.Our treatment of the floors on the hydrodynamic fluxes is described further in Appendix B.

RESULTS
We compared the angular momentum supply at the BH event horizon radius for the thin disk model and the MAD simulations.We found that the hydrodynamic disk in the MAD simulations does not behave in the same way as the NT disk. Figure 1 shows that the infalling material in the inner accretion disk becomes strongly sub-Keplerian and sub-NT.Strikingly, specific angular momentum flux on the BH, computed using the hydrodynamic component of T r ϕ (Eq.12), makes up only 37% of the NT value, a consequence of strong magnetic fields in MADs.This value is strikingly low even for MADs, because usually MAD angular momentum is compared to the Keplerian value, with the typical degree of MAD sub-Keplerian rotation quoted at around 50% (Igumenshchev 2008;McKinney et al. 2012;Begelman et al. 2022).Hence, the disk has nearly three times less specific angular momentum than the NT expectation to offer the BH and might not be able to significantly contribute to spinning up the BH.

BH spin evolution for NT and MAD flows
To determine BH spin evolution, we solved the coupled ODEs, equations (2), (3), and (4), for BH spin a(m) and mass M(m), for both the standard NT disk and MAD flows.Figure 2 shows the BH spin evolution a(m) for the NT (Fig. 2a) and MAD (Fig. 2b) cases.The NT evolution is calculated using the equations described in Section 2.1.From each MAD simulation, we calculate the fluxes at the horizon, e in and l in .We then use spline interpolation to evaluate the simulated data points over the full range of spins, resulting in the interpolation functions, e in (a) and l in (M, a), that we plug into Non-relativistic Keplerian specific angular momentum (blue; lKep = r 1/2 ) and MAD (red) normalized angular momentum flux vs. radius for a = 0.9.The BH horizon radius rH is the lower x-axis limit, and the innermost stable circular orbit radius rISCO is shown by the vertical dashed line.The angular momentum in the hydrodynamic flow in the MAD, lHD,MAD, is sub-Keplerian as the gas falls onto the BH, indicating that some mechanism is robbing the disk of angular momentum thereby depriving the BH of the angular momentum supply.
the ODEs (see Fig. 4).We then solve the ODE equations to model the spin evolution.
For each disk type, we consider three values for the initial spin: a 0 = −1, a 0 = 0, and a 0 = 1.For any initial spin a 0 , the NT disk in Figure 2(a) will spin up the BH to equilibrium at the maximum a eq = 1 before the BH has accreted twice its initial mass, m/M 0 = 2.The infalling gas supplies angular momentum and energy to the BH, and in the NT model there is nothing to remove angular momentum or energy from the BH.Thus, a can only increase.For the MAD case shown in Figure 2(b), the BH spins down rapidly to the low value of a eq,MAD ≈ 0.07, and in doing so it has only accreted roughly 50% of its initial mass.
We translate this timescale into physical units at the upper x-axis in Fig. 2. The Eddington accretion rate is given by where σ T is the Thomson cross-section, m p is the proton mass, and ϵ is the radiative efficiency.We assume ϵ = 0.1.We show the cosmological spin-down timescale in units of Myrs over a fixed Eddington ratio, For a BH accreting at the Eddington rate (λ 0 = 1), the MAD BH will spin down from a = 1 to a = 0.2 in approximately 9 Myrs, by which time it will consume about 20% of its own mass.In order to explain how the equilibrium spin in the MAD case, a eq,MAD , can become so small, as well as why the spindown is so much faster compared to the NT case, we need to understand the physics of the torques acting on the BH.First, we know that MADs launch powerful electromagnetic jets, and from the BZ model we know that they are powered by BH spin.Additionally, the accretion disk will exert a torque on the BH.Each in part carry their own angular momentum and energy.We now identify the influence of each of these factors.

BH power output
We showed that in the NT disk, the BH can only spin up with time.The BH eats energy and angular momentum, and there are no large-scale magnetic fields to efficiently carry away the angular momentum.In a MAD, energy and angular momentum do not travel just into the horizon, but also a > 0 : 106.3a 4 + 39.5a 2 a < 0 : − 19.8a 4 + 48.9a 2 η EM η NT Figure 3. Dimensionless efficiency η as a function of BH spin a for different models.The radiative efficiency ηNT for the Novikov-Thorne disk is shown in the blue dash-dotted line, and the electromagnetic outflow efficiency ηEM = PEM/ ṁc 2 at the event horizon is shown with dark blue dots.We provide separate fits to ηEM(a) dependence for positive (pink) and negative (orange) spins (see eq. 16).A rapidly-spinning BH fed by a MAD can achieve EM energy outflow efficiencies > 100%.These fits form the inputs to our semi-analytic MAD spin-down model.outward in the form of jets and winds.They derive their power in part from the EM BZ power of the rapidly spinning BH.We define the EM power efficiency as the ratio of the time-averaged EM energy flux leaving the BH to the timeaveraged rest mass flux entering the BH: We plot the EM efficiency η EM as a function of spin in Figure 3.The electromagnetic energy flux on the BH is used to calculate the jet efficiency for each simulation and is shown as black points.To produce the pink and orange curves, we use the least squares fits in the form C 1 a 4 + C 2 a 2 for each of the negative and positive values of a.This results in the fit for η EM : For a NT disk, η NT never exceeds 45%.In the work of Gammie et al. (2004), the torque on the central BH was provided by a disk model that reached a maximum jet efficiency of 16%.MADs are known to achieve enormous efficiencies of 150% (Tchekhovskoy et al. 2011) and even higher (McKinney et al. 2012).An efficiency of 100% means that all accretion power is processed into jet power.Thus an efficiency exceeding 100% implies that the jets are feeding on the angular momentum of the BH.
For the BH to launch jets, it needs to hold magnetic flux.Jet power is directly proportional to magnetic flux, Φ BH .BH generates jet power via the BZ mechanism at the efficiency ( Blandford & Znajek 1977;Tchekhovskoy 2015), where κ is a constant that depends on the magnetic field geometry (0.044 for a parabolic geometry), Ω H = ac/(2r H ) is the angular frequency of the BH horizon, r H is the radius of the BH horizon, and Φ BH is the magnetic flux on the BH.The approximation, f (Ω H ) = 1, can be used for a ≤ 0.95, but for rapidly-spinning BHs with a > 0.95, the following approximation maintains accuracy, f (Ω H ) ≈ 1 + 1.38(Ω H r g /c) 2 − 9.2(Ω H r g /c) 4 (Tchekhovskoy 2015).We use η BZ to guide our choice of fit for η EM in Figure 3 and equation ( 16).

Disk and jet contributions to spin-down
Both the disk and the jets exert a torque on the BH.To understand why the MAD spins down to such a low value of spin and why it happens so quickly, we now decompose the total torque into its EM and hydrodynamic constituents, which we compute via the corresponding fluxes at the event horizon.We distinguish the disk (from the jets) using a cutoff condition on the magnetization, σ ≡ 2p mag /ρc 2 < 30; in the simulations, we set the density floor, which limits the magnetization in the polar regions to max σ = 50, see Appendix B for more details.Here, p mag is the magnetic pressure.Using the above condition, we find that, unsurprisingly, most of the electromagnetic (EM) flux escapes through the jets, and most of the hydrodynamic flux enters through the disk.To keep things simple, from now on, we will assign the torque by the jets to be the entire electromagnetic component and the torque by the disk to be the entire hydrodynamic component.We note that despite this naming convention, some "jet" magnetic field lines may connect to the BH through the disk.
We plot specific angular momentum and energy fluxes vs. spin in Figure 4. Figure 4(a) shows that the NT disk (blue dashed line) supplies the BH the most angular momentum flux and energy flux at a = −1 and the least at a = 1, and always supplies positive energy and angular momentum.The total MAD energy flux computed from our simulations, e MAD , is shown in red, where each data point corresponds to a simulation.We find that its hydrodynamic constituent remains roughly constant, e in,HD ≈ 0.97, for all but the highest positive values of a.However, at such high values of a, it is the EM component that dominates e MAD .Figure 4(b) reveals that the hydrodynamic specific angular momentum (orange) supplied by the MAD is significantly less than that supplied by the NT disk.The difference in angular momentum supply might be related to the strong magnetic fields efficiently transporting the disk angular momentum outwards, either within the disk or in the form of large-scale outflows (e.g., Blandford & Payne 1982).Furthermore, Fig. 1 shows that the disk rotation becomes significantly sub-Keplerian, which can starve the BH of the hydrodynamic angular momentum.This sub-Keplerian rotation can be related to the vertical thickness of the disk but it is probably also related to the large scale field (Scepi et al. 2023).
It is striking that l in,HD remains roughly constant at ≈ 0.86 for all values of a.However, it can be understood based on the argument described above.The dominant contributions to l HD are the strong field and the sub-Keplerian motions.Both effects will reduce the importance of the ISCO, which is defined for Keplerian motions, leading to l HD being independent of spin and only dependent on the disk dynamics.
Interestingly, the MAD spin-down cannot be reproduced by simply taking the torque from the NT disk and adding the spin-down contribution from the jets.In fact, we cannot assume that the two disks supply the same angular momentum and energy to the BH.In fact, the MAD hydrodynamic torque is quite different than that of the NT disk.This is one reason for why the MAD spins the BHs down to a eq ≈ 0.07, and why it happens so quickly: at every stage of its evolution, a MAD supplies very little angular momentum to the BH.
We illustrate the MAD spin-down in Figure 5 through the dimensionless spin-up parameter (Gammie et al. 2004), which is defined in equation (5).When s = 0 the BH spin is no longer evolving and has reached equilibrium.The blue dashed line shows the spin-up parameter for the thin Novikov-Thorne disk, where the blue vertical band shows where s NT = 0 at a = 1, just as we showed in Figure 2. The red curve shows the total spin-up from the MAD simulations and was calculated using the form of equation ( 6) with equations ( 11) and ( 12).The red vertical band shows s MAD = 0 at a ≈ 0.1, just as we showed a eq,MAD ≈ 0.1 in Section 3.1.The spin-up parameter remains positive for a NT disk, since accretion supplies angular momentum and energy, increasing the BH spin.For a MAD with s < 0, the BH loses energy to relativistic jets and its spin decreases back down to the equilibrium value.
We note that several simulations, e.g.those with a = 0.2 and a = 0.5, have a relatively short duration, which limits the time interval over which we averaged them over to measure the spin-up parameter (Table 1).However, the simulations that bracket the equilibrium spin, those for a = 0.0 and a = 0.1, extend out to later times.To compute the equilibrium spin numerically, we use linear interpolation between the simulated values of s for these two longer simulations.
Orange points in Figure 5 show the hydrodynamic component of the spin-up parameter, s.Interestingly, s HD becomes negative at a ≳ 0.5.Although the disk supplies the BH with angular momentum, at high spin the magnetic field removes much of the angular momentum from the disk before it reaches the BH.From Eq. ( 6), we see that the value of s comes from the balance between the energy and angular momentum fluxes into the BH.Because in the definition of a (Eq. 1) the angular momentum enters in the numerator (J), the supply of positive angular momentum increases the BH spin.In contrast, BH mass (M) enters in the denominator, so its increase (through the supply of positive energy) reduces We show the NT model in blue.MAD hydrodynamic disk is shown in orange, and our model for the disk is shown by the dashed orange line.In a MAD, the disk component follows a similar trend as the NT disk, but its torque on the BH is much smaller except at the highest spin values.The EM component is shown in green, with our jet model shown by the green dashed line.The total MAD data is shown by red points, and our total MAD spin-down model is shown by the thick, red curve.The values of equilibrium spin aeq for MAD and NT disks are shown as vertical lines at a = 0.1 and a = 1, respectively: these are where the spin-up curves vanish, s = 0. We note that the thickness of the vertical lines does not correspond to uncertainty.Our semi-analytic sMAD model agrees well with the simulation data.Differences for positive spin are primarily caused by an imperfect fit for ηEM (as seen in Fig. 3), affecting sEM.Our semi-analytic model for MAD torques on the BH shows how and why MADs spins down their BHs to a very low value of aeq,MAD ≈ 0.07.MADs are significantly modified by the large-scale, dynamically-important magnetic fields.Black triangles are higher resolution simulations with Γ = 13/9, demonstrating that higher Γ leads to stronger spin-down.the BH spin.Because the energy supply term in Eq. ( 6) is proportional to BH spin, its effect increases with increasing spin and ends up dominating the angular momentum term at high spin: in other words, at high spin the disk angular momentum supply to the BH is too low to balance out the energy supply.This causes the hydrodynamic torque on the BH to become negative at high spin values.

Model for MAD spin-down
Here we construct a model that reproduces the MAD simulation data.Jet power must play an important role in BH spin-down because the jets extract BH rotational energy at high spin.Moderski & Sikora (1996) described the spin evolution, including the effect of jet-like outflows, analytically by taking the BH spin evolution equations (Eqns. 2 and 4) and including the terms for jet power.We write them as and where k = Ω F /Ω H , the ratio of the angular frequency of the magnetic field lines to the angular frequency of the BH event horizon, Ω H .The first terms now only describe the disk (hydrodynamic) contributions.
As we did for the NT disk and MAD, we similarly recast equation ( 18) in the form of equation ( 5) to now determine how the disk and jets contribute to BH spin-up.The new expression for s is then We now use approximate fits for k, Ω H , η EM , l in,HD , and e in,HD (see Figs. 3, 4, and 6) to construct our semi-analytic model.Namely, we approximate l in,HD with the constant value, l HD,avg = 0.86, and e in,HD with the average value of 0.97, both shown by the dashed orange lines in Figure 4. We use the fit, η EM , in the form η ∝ C 1 a 4 + C 2 a 2 as we show in Figure 3 and equation ( 17).
The standard monopole BZ model has k = 1/2, which gives the maximum EM efficiency.We calculate k at the horizon and then average in t, θ and ϕ. Figure 6 shows that for our semi-analytic model, we approximate this dependence with the average value, k = 0.23, for a < 0. For a > 0, we find a good fit for k, k = 0.23, a < 0, min(0.1 + 0.5a, 0.35), a > 0. (21) Similar to Penna et al. (2013), k is smaller than what is predicted by the standard BZ model.We could have modeled the angular momentum extraction by the jet EM torques by fitting l EM instead of k.However, k is a useful quantity in BH jet theory because it connects the angular momentum extraction efficiency to the better-known and better-calibrated EM energy extraction efficiency, η EM .Moreover, from the consideration of the standard BZ effect, we expect k ∼ 0.25 − 0.5, for a wide range of jet geometries ranging from monopolar to parabolic (Tchekhovskoy et al. 2010).Therefore we expect k to show less variation with spin than η EM or l EM .As an additional perk, fitting k can inform us of the deviations from the standard BZ effect.
For a = 0, there is no energy extraction along the field lines because Ω H = 0, so we set the second term in equations ( 18), (19), and (20) to 0. This might not be obvious in equations ( 18) and ( 20).However, because Ω H = ac/2r H and η EM ∝ a 2 , the Ω H ∝ a term in the denominator cancels out: thus, for zero spin, this results in zero contribution from the jet terms (which contain P EM and η EM ).
With these constraints in our semi-analytic model we produce the thick red curve in Figure 5 .We plot k = ΩF/ΩH and its dependence on spin, which is crucial for converting jet power into angular momentum in our semi-analytic model.The red dashed line shows the fit that we used to construct our semi-analytic model for a > 0. We did not find a strong dependence on spin for a < 0 and so used a constant value, k = 0.23, in our semi-analytic model for a < 0.
MAD spin-up s MAD curve.This fits the data well, and we expect this semi-analytic model to be useful for modeling BH spin evolution.See Section 4 for more details.The dashed lines in Fig. 2(b) show the spin evolution, due to the integration of Eqs. ( 2) and ( 4), computed using a spline interpolation of the simulated data points, e in (a) and l in (M, a), shown in Fig. 4. In Fig. 2(b) we also show, in thin solid lines, the spin evolution computed using our semi-analytic model developed above.We notice a remarkable agreement between the spin-evolution computed using the interpolated data and the semi-analytic model.There are some slight difference for a 0 = 1 that could be related to our model being less accurate around a ∼ 0.5.We stress that even with those slight differences our semianalytic model reproduces the spin evolution with good accuracy and is adaptable for modeling BH spin in different accretion disk regimes and astrophysical contexts (e.g., Sec 4.2).
In Figure 7 we show the spin, mass and jet efficiency evolution in MAD flows vs. time over a characteristic BH massdoubling timescale τ .The left hand column shows the prograde BH case with initial spin a 0 = 1, and the right hand column shows the retrograde case with a 0 = −1.The timescale τ is defined such that t/τ = m/M 0 , where m is the accreted mass and M 0 is the initial BH mass.In the top panels of Fig. 7 we have fit the spin evolution to, (black dotted line), where τ MAD = τ /9.5 and a 0 is the initial spin.
The second row of panels from the top shows the dependence of BH mass M vs. t.The total BH mass, M is shown by the red lines, and the irreducible BH mass, M irr = M( 1 2 r H (a)) 1/2 , is shown by the dashed blue lines.The irreducible mass increases with time.The total BH mass in-creases with time for the retrograde BH, but for the prograde BH, we find that the total mass find that both the total mass and irreducible mass of the BH overall increase with time.While the BH is rapidly spinning, the total and irreducible BH masses are distinct because the BH also possesses rotational mass.As the BH spins down, M irr and M converge, after t/τ ≈ 0.1.
The bottom panel shows the efficiency η vs. massdoubling time.In terms of the specific energy flux e in onto the BH, we can write down the total efficiency with which the BH converts rest-mass energy into power: In terms of the hydrodynamic and electromagnetic specific energy fluxes, we can write it as The energy fluxes e HD and e EM are computed using the fits derived above.As the BH spins down, the EM efficiency drops.In comparison, the hydrodynamic efficiency, |1 − e HD |, remains roughly constant, consistent with the fact that in MADs, the disk torques do not significantly affect BH spin evolution with time.
Interestingly, the inset in the left middle panel of Fig. 7 shows that for the prograde case, the total BH mass slightly decreases before increasing.We see that during the short time the BH mass decreases, the rapidly rotating BH is spinning down so fast that the total efficiency η, the red curve in the bottom panel of Fig. 7, exceeds 100%.This is the first demonstration that, as a function of time, the BH actually loses mass due to the jets outshining accretion.The cause of this phenomenon is that the BH rapidly loses rotation energy at high spin electromagnetically (through the jets), faster than the irreducible mass increases due to the energy supplied by accretion (through the disk).

Analytic model to calculate a eq in MADs
Because the BH reaches a eq when s = 0, we can work backward and use our semi-analytic model to find the value of a eq analytically.We set equation ( 20) to 0 and solve for a eq .As before, we use the average values in our model for the disk fluxes: ⟨l in,HD ⟩ ≈ 0.87 and ⟨e in,HD ⟩ ≈ 0.97.We use the positive fit η EM ≈ 0.4a 2 .For k we use the fit for low positive spin, k = 0.1 + 0.5a.For small a, Ω H ≈ a/4.Using these approximations, dropping any terms that are higher order than ∝ a 2 (they are small for low spin), we get a quadratic equation for a eq that gives us the positive solution, a analytic eq,MAD ≈ 0.06.This is remarkably close to the equilibrium spin value we have obtained directly from the simulation results, a simulated eq,MAD ≈ 0.07 (see Sec 3.1).A MAD BH will spin down to aeq ≈ 0.07 by the time it has accreted a small fraction of its initial mass.For instance, a maximally spinning prograde BH spins down to a = 0.2 by t/τ = 0.2, i.e., after consuming 20% of its own mass, and to a = 0.1 by t/τ = 0.35, i.e., after consuming ≃ 35% of its own mass.[Middle panel]: Total BH mass (red solid) and irreducible BH mass (blue dashed) vs. time.At t/τ ≲ 0.04, the mass of the prograde BH briefly decreases before increasing again, as seen in the inset (left middle panel).This shows that BH loses rotational energy faster via the BZ mechanism than the accretion can replenish it (i.e., η > 100%).At later times, BH mass increases for both prograde and retrograde BHs.[Bottom panel]: Electromagnetic (EM) energy efficiency.In the prograde case, the EM efficiency drops by more than three orders of magnitude from a = 1 to late times, when BH spin reaches its equilibrium value, aeq ≈ 0.07.Whereas for the integration in time to obtain the a(t) dependence we use the interpolated values of s measured from the simulations, we evaluate the EM efficiency, shown in the bottom panels, using the fits given in Fig. 3 and computed in Sec.3.4.

Semi-analytic Model of MAD Spin-down
We have used a suite of 3D GRMHD simulations of MADs from Tchekhovskoy et al. (2011);Tchekhovskoy et al. (2012); Tchekhovskoy (2015) to study the BH spin evolution in MADs.Tchekhovskoy et al. (2012) were the first to find that a MAD will spin down a BH to an equilibrium spin of a eq ≈ 0.07.Here, using the equations describing the BH spin evolution in time under the torques from the accretion disk and the jets, as measured in a suite of 3D GRMHD simulations, we showed that MADs will cause the BH to rapidly spin down, such that the BH reaches its equilibrium spin after accreting only 50% of its initial mass (see Fig. 2).Indeed, BH spin-down by a MAD disk is much more efficient than the BH spin-up by a NT disk (see Fig. 2).
Such a rapid spin-down might come as a surprise, because intuitively we think of accretion as spinning up the BH, by supplying it with the angular momentum.To understand the physics behind the MAD spin-down, we split up the torques acting on the BH into their hydrodynamic and magnetic components, as seen in Fig. 5. Figure 1 shows that some mechanism, likely magnetic fields threading the inner disk, robs the disk gas of its angular momentum before it can reach the BH: the MAD specific angular momentum lies below the Keplerian one and, in fact, it decreases towards the BH within the innermost stable circular orbit (ISCO), inside of which the specific angular momentum is assumed to be constant in the NT model.We thus end up with a strongly sub-Keplerian disk whose specific angular momentum decreases even further inside the ISCO and hence carries into the BH very little angular momentum, which is available for spinning up the BH.
In addition to the magnetic fields extracting angular momentum from the BH gas supply, strong and dynamicallyimportant magnetic fields, which thread the BH, extract the angular momentum from the BH and send it out to large distances in the form of jets.Thus, the low equilibrium BH spin is a consequence of both of these effects: the lower-thanexpected angular momentum supplied by the disk gas pales in comparison with the large magnetically-mediated angular momentum extraction out of the BH by the large-scale magnetic flux (see Figs. 4 and 5).
To quantify this, we have developed the first semi-analytic model that quantitatively describes what causes the BHs in the radiatively-inefficient MAD state to spin down so efficiently.We find that the hydrodynamic angular momentum and energy fluxes are essentially independent of the BH spin, as seen in Fig. 4. Thus, to the accuracy we are interested in, the hydrodynamic angular momentum (l HD ) and energy (e HD ) flux values -which we take to be constants in our semi-analytic model (see Fig. 4) -are unique properties of the MAD state, at least for the typical disk thickness that we have studied here (h/r ≈ 0.3).Using equation (20), we now can compute the hydrodynamic contribution to the spin-up parameter, s HD , which we plot with the orange dashed line in Fig. 5.It is remarkable that the gas retains so little angular momentum by the time it reaches the BH: in fact, s HD vanishes at a = 0.5, which means that at a > 0.5, the hydrodynamic torques spin down the BH.This means that if we were to consider a hypothetical hydrodynamic disk model with the structure identical to MAD (but without the magnetic fields), this model would result in a eq ≃ 0.5.
The final step in building our semi-analytic model is to account for the EM torques, which we do via including the energy and angular momentum EM fluxes: we compute the values of these fluxes using the fitting formulas for the spindependence of the jet efficiency (η EM , Fig. 3) and the angular velocity of the magnetic field lines (k = Ω F /Ω H , Fig. 6).

Whodunnit?
Equation (20) gives the sum of the hydrodynamic and EM contributions, thus completing our semi-analytic model and giving us the spin-up parameter, s, which is shown in Fig. 5 with the thick light red curve.The model is in good quantitative agreement with the numerical values for s, shown with dark red squares.This agreement is even more remarkable given that in order to keep the model as simple as possible, we purposefully used quite rough approximations for the free parameters in the model (e.g., l HD , e HD , and k), as we discussed in Sec.4.1.
Figure 5 shows that MADs have |s EM | ≫ |s HD |, i.e., the EM torques dominate the hydrodynamic ones, over a wide range of spin, a ≳ 0.2.This is made uniquely possible by the strongly sub-NT angular momentum of MADs: the specific angular momentum of the MAD (orange dashed line) lies much lower than that of the NT disk (blue dash-dotted line).
To gauge the relative importance of the various processes for the rapid BH spin-down in MADs, we consider hypothetical disk models that differ in crucial aspects from our MAD simulations.For instance, in order to quantify the importance of sub-Keplerian nature of MADs on the spin-down, let us consider a disk model whose EM structure is that of the MAD, but whose hydrodynamic structure is that of the NT model.This combination might provide a rough description of luminous MADs, whose disks are thinner and hence closer to NT.For such disks, the equilibrium spin would shift to higher a values to satisfy s NT + s EM = 0. From Fig. 5, we can approximately see that this happens at a eq ≈ 0.3.This means that even if the accretion flow were Keplerian (instead of significantly sub-Keplerian), the jets would still be strong enough to spin down the BH to a rather low, albeit higher, spin.If, additionally, the EM outflow efficiency drops substantially in luminous MADs compared to non-radiative ones (e.g., η EM in Eq. 20 is lower by a factor of 2 than in our simulations, similar to what is seen in radiation-transport simulations of luminous MADs at L/L Edd = 0.35, Liska et al. 2022a), then the spin would evolve to s NT + 0.5s EM = 0, resulting in the equilibrium spin of a eq ≲ 0.5.
The rather low equilibrium spins in the above hypothetical disk models can have interesting implications for luminous MADs, ones in which the large-scale magnetic flux is saturated on the BH, but those that can efficiently cool down to small values of scale height, h/r ≪ 1.Such disks are particularly important because they can power luminous jetted quasars.The rotation profile of such thinner disks can be much closer to the Keplerian one than of our thick MADs.However, from the above arguments, it follows that even in this case the equilibrium spin can remain significantly below unity (e.g., a eq ∼ 0.2 − 0.5), defying the standard expectation of the "standard" NT disk model.

Comparison to other works
There has been plentiful other theoretical work studying BH spin evolution.As we mentioned above, Bardeen (1970) was the first to compute an equilibrium spin for a BH surrounded by a turbulent disk, a eq = 1.Fully turbulent rel-ativistic disk solutions were also calculated by Popham & Gammie (1998).They obtained smaller equilibrium spins, as low as a eq = 0.8.They achieved this by modifying the inner boundary condition for the angular momentum.They found that deviations from the Keplerian angular momentum profile controlled the equilibrium spin.Finally, the turbulent viscosity acting on the disk controls the angular momentum profile.We find striking similarities with our results as deviations from a Keplerian angular momentum profile are also important in our model.Gammie et al. (2004) were the first to use simulations of magnetized disks to measure the spin-up parameter and try to determine the equilibrium spin.They computed a large equilibrium spin when compared to our simulations, a eq ≃ 0.93.A possible explanation for this discrepancy is that their simulations are only 2D.Furthermore, compared to MADs, the jet power they measure in their simulations is weak (McKinney & Gammie 2004).In this work, we find that the magnetized jet is quintessential for the rapid BH spin-down.Thus, we expect the difference in the equilibrium spin to be primarily the consequence of the difference in the jet power.Narayan et al. (2022) ran long-duration GRMHD simulations of MADs.They calculate the BH spin-up parameter as a function of BH spin.Using a fifth-degree polynomial fit for a, they measured an equilibrium spin of 0.04.This value is a factor of 2 smaller than ours.However, they only ran one simulation with a BH spin close to their a eq , which may limit the accuracy to which a eq can be determined.The simulations of Narayan et al. (2022) use an adiabatic index of Γ = 13/9 and tested the effect of changing Γ in their simulations and reported no significant difference in their results.We have run new simulations with the GPUaccelerated GRMHD code H-AMR (Liska et al. 2022b) with Γ = 13/9 to study the effect of changing the adiabatic index.The spin-up parameter values for these simulations are shown by the black points in Figure 5.The spin-up parameter values for all spins are shifted down from the Γ = 4/3 simulations in red.These Γ = 13/9 simulations are consistent with Narayan et al. (2022).Therefore, the difference in equilibrium spin may also be due to a different choice of Γ.
It is also possible that the spin-up parameter and the equilibrium spin weakly depend on time.Indeed, Narayan et al. (2022) ran their simulations for a much longer time than ours.We investigated this possibility by computing the spin-up parameter as a function of time in our simulations.We found that the spin-up parameter does not show obvious trends in time, at least for the duration of our simulations.
In comparison to our work, Narayan et al. (2022) focus on low-luminosity AGN and do not account for the change of BH mass when computing the time evolution of BH spin.This limits their ability to evolve the system to late times, when the BH spin approaches its equilibrium value.

Observational implications
Our 3D GRMHD simulations of radiatively-inefficient accretion flows can be applied to the super-Eddington accretion regime: in this regime, the accretion flow is so dense and optically-thick that the radiation does not have the time to escape out of the accretion flow and is instead advected into the BH or ejected as part of the outflow.Although radiativelyinefficient GRMHD simulations are usually interpreted in the context of very sub-Eddington accretion ( Ṁ ≪ 0.01 ṀEdd ), we can make inferences for super-Eddington accretion by treating radiation-dominated optically-thick highly super-Eddington ( Ṁ ≫ ṀEdd ) flows as an ideal gas with the polytropic index, Γ = 4/3, that corresponds to radiationdominated gas.In an upcoming work (Lowell et al. 2024), we will show that the spin-down torques in this work are consistent with those in radiation transport simulations of super-Eddington BH flows.
Figure 7 shows the time evolution of the BH spin and the associated quantities over cosmological timescales.The bottom panels reveal that the spin-down leads to a considerable decrease in the jet efficiency, η EM (a = 1)/η EM (a = a eq ) ≃ 10 3 .This can have significant observational consequences.
In the context of AGN, this could naturally explain the reported 10 3 spread in the radio loudness of AGN (Sikora et al. 2007).Indeed, radio-quiet quasars can have low spin, whereas radio-loud quasars can have high spin, as expected in the spin paradigm (Sec.1).The typical challenge facing this paradigm is that in order to explain such a large spread requires radio-quiet AGN to have very small spin, much lower than expected in the chaotic accretion model (Volonteri et al. 2007;Tchekhovskoy et al. 2010).If AGN undergo spectral state transitions similar to microquasars, their central BHs would undergo periods of spin-down during the jetted MAD stages and spin-up during the jet-less stages.Thus, the radio-quiet AGN can form as the direct consequence of BH spin-down due to MADs.We note that other possibilities, such as the changes in the strength of the magnetic flux on the BH (relative to the MAD flux), changes in the ambient medium, or the presence of selection effects could contribute to the radio-loud/quiet dichotomy.
Gravitational wave detections (Abbott et al. 2016(Abbott et al. , 2019a(Abbott et al. , 2021) ) have led to new ways of constraining BH spin.The gravitational wave signals contain information on the angular momentum of the BHs and their orbits.However, disentangling the BH spin distribution from this signal requires careful Bayesian analysis and is subject to multiple biases.Many teams have performed such type of analysis (Abbott et al. 2019b;Wysocki et al. 2019;García-Bellido et al. 2021;Edelman et al. 2022;Callister et al. 2022;Tong et al. 2022;Abbott et al. 2023).They consistently find that the BH spin distribution is centered around small spins, a peak < 0.3 (see however Galaudage et al. 2021).This distribution is consistent with our measurement of the equilibrium spin.Whereas it is unknown whether super-Eddington accretion or the MAD state are present in the formation and evolution LIGO/VIRGO/KARGA BHs, in certain BH formation scenarios it is possible to infer the presence of a MAD state and BH spindown.For instance, Jacquemin-Ide et al. ( 2023); Gottlieb et al. (2023) argue that BH spindown is inevitable in the context of collapsars, dying stars whose core collapses into a BH.Conversely, other evolutionary pathways, where the MAD state and super Eddington accretion are absent, could lead to a BH spin distribution centered around small spins.In this context, it is intriguing that BH spins measured in X-ray binaries and AGN have populations of both slowly and rapidly spinning BHs (Reynolds 2021;Fishbach & Kalogera 2022).
Gamma-ray bursts are some of the most powerful astrophysical phenomena in the universe.The collapsar model describes GRBs as the consequence of a failed supernova where all the energy of the gravitational collapse is funneled into a jet (Woosley 1993).The jet is fed by the rotation of a central engine, a BH, through the BZ mechanism.Even though they are short events, lasting around 100 seconds, they are subject to very high accretion rates and are highly super-Eddington.
Recent large-scale GRMHD simulations of collapsars bridging the gap between the BH scales and the photosphere were run by Gottlieb et al. (2022a,b).They showed that their simulations achieve the MAD state and launch powerful jets.They also found that, for realistic conditions, the accretion rate can reach 0.01 solar masses per second.This accretion rate would lead to around 1 solar mass accreted for a GRB lasting around 100 seconds.For a BH of 2 solar masses, this would lead to a considerable spin-down (see Fig. 2).We should expect the remnant BHs to have a spin distribution centered around the equilibrium spin.The spin evolution during GRBs will be the subject of future work.
Furthermore, such a powerful spin-down could shut off the jet produced by the BH engine (see jet efficiency in Fig. 7, and recall that m/M 0 = t/τ ).If the jet power is reduced considerably, the weaker jet could have trouble escaping the stellar envelope.However, the engine activity may be controlled by the magnetic reservoir in the stellar envelope instead of the angular momentum reservoir in the central BH (Gottlieb et al. 2022a).Thus, our results imply that collapsar mod-els need to take into account the effects of BH spin-down to model the evolution of the jet engine.
We do not expect the accretion disk to accrete at super-Eddington rates most of the time.Therefore, it is important to determine how radiative effects influence the spin evolution of BHs.Cooler disks will be thinner and will have different properties when compared to thicker super-Eddington disks: they tend to more closely resemble a Keplerian disk as the scale height decreases.Cooling will likely modify the angular momentum profile and possibly the jet structure.Therefore, a thinner MAD will have a different equilibrium spin.However, as described in Sec.4.2, the equilibrium spin is still expected to be much lower than unity, even in the presence of radiative effects.Future work will focus on the effects of cooling on the spin evolution of MADs.This will provide insight into the physics and observational signatures of luminous quasars accreting at sub-Eddington rates, 1 − 10% of L Edd .where and Within the marginally-stable orbit, the energy and angular momentum fluxes are conserved.Using these, we calculate the specific energy flux on the BH horizon to be The specific angular momentum flux on the BH is L in = 2M 3 3/2 (1 + 2(3R ms − 2) 1/2 ).(A5)

B. HANDLING NUMERICAL FLOORS
This appendix describes how we avoid numerical floor contamination of the hydrodynamic fluxes computed in Section 3. In the integrals, Eqs. ( 11) and ( 12), we only consider regions where σ < 30.This procedure excludes cells in the jet that are subject to numerical gas floors.
Figure 8 shows the σ profile as a function of theta measured at r = r H and azimuthally averaged.The gas is the dominant component in a thin sheet around θ = π/2.Hence, we choose a cutoff of σ = 30 that removes the floor contamination, visible in the polar regions, at |θ/π − 0.5| ≳ 0.1.Different snapshots in time show similar trends.
In Fig. 9 we show the accretion rate (Eq.10) computed with and without the σ cutoff.The accretion rate, corrected using our σ cutoff, has a flatter profile at the inner radii when compared to the uncorrected one.We can remove the accretion boost due to the numerical floor contamination of density.We can then use this σ cutoff to compute the fluxes for quantities that steeply depend on radii like l HD .We apply the sigma cutoff to all hydrodynamic fluxes ( ṁ, l HD , e HD ).However, in practice ṁ and e HD are constant in radius after applying the σ cutoff and could be safely measured at 5r g .The electromagnetic fluxes, l EM and e EM , require no correction as they do not include quantities contaminated through the numerical floors (enthalpy, density, pressure).
We have run an extra simulation with a different density floor value, given by a smaller σ max = 5 instead of σ max = 50 that is used throughout the manuscript.This factor of 10 increase in the density floors, coupled with a lower cutoff value, σ < 5, did not lead to any significant differences in the measured values of the fluxes.10) for a simulation with a = 0.9 with (green line) and without (purple line) using the σ < 30 cutoff.Notice that once the cutoff is applied the accretion rate profile becomes flatter at the inner radii.This way, we can remove the density floor contamination from mass and energy accretion rates.

Figure 2 .
Figure2.Dimensionless spin a(m) and a(t) for different values of initial spin: a0 = −1, 0, 1.The lower x-axis gives the mass accreted over the initial BH mass: when m/M0 = 1, the BH has accreted its entire initial mass.The upper x-axis gives spin evolution in Myrs over an initial Eddington rate λ0.[panel (a):] The BH with a Novikov-Thorne (NT) disk gains angular momentum and energy, and will spin up to aeq = 1.[panel (b):] A MAD will spin down the BH to aeq ≈ 0.07.The MAD spin-down timescale is very short (much shorter than NT spin-up timescale): the BH reaches its equilibrium spin after accreting only half of its initial mass.For a fixed λ0, this spin-down timescale is independent of initial BH mass.Dashed lines show BH spin evolution based on our GRMHD simulation results (computed using a spline interpolation of simulated data points, ein(a) and lin(M, a), shown in Fig.4) and solid thin lines show the spin evolution computed using the semi-analytic model described in Sec.3.4: our model provides a good description of the simulation results.

Figure 4 .
Figure 4. Specific energy and angular momentum fluxes in the Novikov-Thorne and MAD models.[panel (a):] Specific energy flux: The Novikov-Thorne energy flux, eNT, on the BH remains positive and decreases rapidly at high prograde spin (blue dash-dotted curve).In a MAD, the hydrodynamic portion, eHD, stays close to unity (orange dashed line), although it shows a slight decrease near a = 1 (orange squares).Additionally accounting for the EM contribution leads to the total specific MAD energy flux, eMAD (red dots connected by light red line).[panel (b):] In the NT model, the specific angular momentum flux, lin,NT, decreases with increasing spin (blue dash-dotted line).Surprisingly, in a MAD, lin,HD is much smaller than the NT value (orange squares) and roughly constant at all spins: the MAD specific angular momentum flux is significantly sub-NT.In fact, the hydrodynamic component of the angular momentum flux is essentially constant across all values of spin (orange dashed line).In our semi-analytic model of MAD spin-down, we use the two (constant) fits shown with the dashed orange lines, lHD and eHD.The thickness of the red lines following lMAD and eMAD data points does not represent uncertainty.The red lines show the spline interpolation through the simulated data points, and horizontal orange lines in panels (a) and (b) represent the best-fits for eHD and lHD, respectively, that we use in our semi-analytic model shown in Fig. (5).

Figure 5 .
Figure5.Dimensionless spin-up parameter s as a function of BH spin a for MAD and NT disk models.We show the NT model in blue.MAD hydrodynamic disk is shown in orange, and our model for the disk is shown by the dashed orange line.In a MAD, the disk component follows a similar trend as the NT disk, but its torque on the BH is much smaller except at the highest spin values.The EM component is shown in green, with our jet model shown by the green dashed line.The total MAD data is shown by red points, and our total MAD spin-down model is shown by the thick, red curve.The values of equilibrium spin aeq for MAD and NT disks are shown as vertical lines at a = 0.1 and a = 1, respectively: these are where the spin-up curves vanish, s = 0. We note that the thickness of the vertical lines does not correspond to uncertainty.Our semi-analytic sMAD model agrees well with the simulation data.Differences for positive spin are primarily caused by an imperfect fit for ηEM (as seen in Fig.3), affecting sEM.Our semi-analytic model for MAD torques on the BH shows how and why MADs spins down their BHs to a very low value of aeq,MAD ≈ 0.07.MADs are significantly modified by the large-scale, dynamically-important magnetic fields.Black triangles are higher resolution simulations with Γ = 13/9, demonstrating that higher Γ leads to stronger spin-down.
Figure6.We plot k = ΩF/ΩH and its dependence on spin, which is crucial for converting jet power into angular momentum in our semi-analytic model.The red dashed line shows the fit that we used to construct our semi-analytic model for a > 0. We did not find a strong dependence on spin for a < 0 and so used a constant value, k = 0.23, in our semi-analytic model for a < 0.

Figure 7 .
Figure7.Spindown of a maximally-spinning prograde BH down to the equilibrium spin results in the reduction of jet power by a factor of ∼ 10 3 .BH spin, mass, and efficiency vs. time for prograde (left column) and retrograde (right column) MADs with initial spins of a0 = ±1, respectively.The time is measured in units of τ , the timescale for the BH to accrete its own mass, t/τ = m/M0.[Top panel]: BH spin exponentially decays in time, closely following the fitting formula (22).A MAD BH will spin down to aeq ≈ 0.07 by the time it has accreted a small fraction of its initial mass.For instance, a maximally spinning prograde BH spins down to a = 0.2 by t/τ = 0.2, i.e., after consuming 20% of its own mass, and to a = 0.1 by t/τ = 0.35, i.e., after consuming ≃ 35% of its own mass.[Middle panel]: Total BH mass (red solid) and irreducible BH mass (blue dashed) vs. time.At t/τ ≲ 0.04, the mass of the prograde BH briefly decreases before increasing again, as seen in the inset (left middle panel).This shows that BH loses rotational energy faster via the BZ mechanism than the accretion can replenish it (i.e., η > 100%).At later times, BH mass increases for both prograde and retrograde BHs.[Bottom panel]: Electromagnetic (EM) energy efficiency.In the prograde case, the EM efficiency drops by more than three orders of magnitude from a = 1 to late times, when BH spin reaches its equilibrium value, aeq ≈ 0.07.Whereas for the integration in time to obtain the a(t) dependence we use the interpolated values of s measured from the simulations, we evaluate the EM efficiency, shown in the bottom panels, using the fits given in Fig.3and computed in Sec.3.4.

Figure 8 .
Figure8.Example of an instantaneous polar profile of an azimuthal average of b 2 /ρ, as evaluated at the BH event horizon, for σmax = 50.At θ/π = 0.5 we see the density dominated disk, whereas in the polar regions, |θ|/π ≳ 0.1, we notice the contamination due to density floors.Due to the steepness of the profile a σ cut, σ < 30, allows to select most of the disk while removing the floor contamination in highly magnetized polar regions.

Figure 9 .
Figure9.Time-average accretion rate as a function of radius, as computed using Eq.(10) for a simulation with a = 0.9 with (green line) and without (purple line) using the σ < 30 cutoff.Notice that once the cutoff is applied the accretion rate profile becomes flatter at the inner radii.This way, we can remove the density floor contamination from mass and energy accretion rates.

Table 1 .
Simulation details. lMAD and eMAD are the t, θ, ϕ -averaged angular momentum and energy fluxes, respectively, on the horizon.