An Analytic Model For Magnetically-Dominated Accretion Disks

Recent numerical cosmological radiation-magnetohydrodynamic-thermochemical-star formation simulations have resolved the formation of quasar accretion disks with Eddington or super-Eddington accretion rates onto supermassive black holes (SMBHs) down to a few hundred gravitational radii. These 'flux-frozen' and hyper-magnetized disks appear to be qualitatively distinct from classical $\alpha$ disks and magnetically-arrested disks: the midplane pressure is dominated by toroidal magnetic fields with plasma $\beta \ll 1$ powered by advection of magnetic flux from the interstellar medium (ISM), and they are super-sonically and trans-Alfvenically turbulent with cooling times short compared to dynamical times yet remain gravitationally stable owing to magnetic support. In this paper, we present a simple analytic similarity model for such disks. For reasonable assumptions, the model is entirely specified by the boundary conditions (inflow rate at the BH radius of influence [BHROI]). We show that the scalings from this model are robust to various detailed assumptions, agree remarkably well with the simulations (given their simplicity), and demonstrate the self-consistency and gravitational stability of such disks even in the outer accretion disk (approaching the BHROI) at hyper-Eddington accretion rates.


INTRODUCTION
Accretion disks play a crucial role in many fields of astrophysics.In particular around supermassive BHs, these disks power active galactic nuclei (AGN) and quasars, and are believed to supply mass at rates exceeding ≳ 10 M ⊙ yr −1 to the BH, ultimately building up most of the mass in SMBHs (Schmidt 1963;Soltan 1982).The vast majority of the theoretical literature on quasar accretion disks has assumed as a starting point something like the Shakura & Sunyaev (1973, hereafter SS73) -disk model, which assumes that disks are geometrically thin (height  ≪ ), optically thick (black-body like), sub-sonically turbulent (sonic Mach number M  < 1), slowly cooling ( cool ≫  dyn = 1/Ω), gas and/or radiation-pressure-dominated (plasma1  ≡  2  / 2  > 1), radiatively efficient, well-ionized, and can be parameterized by an effectively constant- viscosity where  ∼ () 2 / 2  < 1 represents some Maxwell or Reynolds stresses and the kinematic viscosity scales as  visc ≡  2  /Ω.While many variations have been introduced (including radiatively inefficient, advection-dominated, magnetically-arrested, magnetically-elevated, "slim," and gravito-turbulent disks; for reviews see Frank et al. 2002;Abramowicz & Fragile 2013), for luminous quasars some form of SS73-like analytic model is still most often the "default" reference.
Recently, Hopkins et al. (2023, Paper I) and Hopkins (2023, Paper II) presented the first simulations to selfconsistently follow gas in a cosmological simulation from > Mpc to < 100 au scales (a few hundred gravitational radii) around an accreting SMBH, including the physics of magnetic fields (seeded from trace cosmological values), multiband radiation-hydrodynamics, non-equilibrium multi-phase radiative thermo-chemistry and cooling, self-gravity, star formation, and stellar evolution/feedback (jets, stellar mass-loss, radiation, supernovae).In these simulations, gas around the BHROI2 is tidally captured by the SMBH of mass  bh ∼ 10 7  ⊙ from larger-scale ISM gas complexes in the galaxy, and free-falls briefly before circularizing to form an accretion disk with Toomre  ∼  support Ω/  Σ gas ≫ 1 and little to no star formation or fragmentation on sub-pc scales.This disk evolves in quasi-steady-state and sustains super-Eddington accretion (up to  ∼ 20 − 30 M ⊙ yr −1 ) onto the SMBH for at least ∼ 10 5 dynamical times at the inner simulation boundary (the simulation duration).Crucially, in Paper II where the disk properties were studied in detail, it was shown that the disks are strongly magnetized, with  ∼ 10 −4 − 10 −2 in the midplane, in the form of primarily toroidal magnetic field (but with mean radial fields and quasi-isotropic turbulent fields only a factor of a few less strong) owing to amplification of magnetic flux accreted from the ISM.These stabilize the disk against catastrophic fragmentation and star formation: without magnetic fields, the disks were shown to be ordersof-magnitude less massive and support factor of ∼ 1000 lower accretion rates and higher star formation rates.The disks also have a flared structure (/ ∼ 0.1 − 0.3 at large radii) with weak vertical stratification owing to trans-Alfvénic, highly 2 Defined as the radius interior to which the BH dominates the potential over its host galaxy of characteristic velocity dispersion  gal , or  BHROI ∼   bh / 2 gal ∼ 5 pc in the reference simulations.
super-sonic turbulence, which is sustained by rapid cooling (M  ∼ /  ∼ 1, with M 2  ∼ 1/ ∼ 1/ cool Ω ≫ 1).These disks are therefore qualitatively distinct from SS73like -disks in most respects.As discussed in Paper II, they also appear to qualitatively differ in key respects from most or all historical models of "strongly magnetized" disks such as magnetically "elevated" or "arrested" disks (e.g.Bisnovatyi-Kogan & Ruzmaikin 1976;Miller & Stone 2000;Narayan et al. 2003), hence referring to them as "flux-frozen" and "hyper-magnetized."And their properties may resolve some crucial long-standing questions: most obviously, the wellknown problem that an SS73-like disk with quasar luminosities should be violently gravitationally unstable outside of just a couple hundred gravitational radii (Goodman 2003).
But, while simulations like those in Paper II are ultimately necessary to explore how, when, and "what kind" of accretion disks form (let alone to follow the highly non-linear, multiphysics evolution), they are complex and difficult to use for simple physical intuition and interpretation.Moreover, owing to computational expense, the simulations capture just one system at one moment in cosmic time, so it is not obvious what we can reliably extrapolate to other systems with different BH masses or inflow rates.Motivated by these considerations, in this manuscript ( § 2) we develop a simple self-similar analytic model for the super-sonically turbulent, rapidly-cooling, hyper-magnetized/flux-frozen disks seen in the simulations (building on previous work such as Pariev et al. 2003), modeling their density, magnetic, and angular momentum structure ( § 2.1), as well as thermal and radiation structure ( § 2.2), and validating this against the direct simulations.We discuss the sensitivity (or lack thereof) of the results to different assumptions ( § 3) and reasonable expectations for the model parameters, then compare to the scalings predicted by the classic SS73 disk theory ( § 4.1), as well as some previous, more similar literature models for "magnetically dominated" disks ( § 4.2) and magnetically elevated or arrested disks ( § 4.3).We discuss the upper limits to accretion and likely limitations of the model at small radii (approaching the horizon) in § 5, and summarize in § 6.Note that our primary focus here is on quasar accretion disks around SMBHs, as opposed to stellar-mass BH binaries (where it is less clear our assumptions apply).
2. THE MODEL With this motivation in mind, we develop here a simple approximate analytic similarity model for the "hypermagnetized" or "flux-frozen" accretion disks that form in the simulations.As we explain the model predictions below, we will compare them to the simulations from Paper II in Fig. 1,3 over the range where one can reasonably identify a rotationally-dominated "accretion disk" in such simulations ( ∼ (300 − 10 6 ) 2   BH / 2 ).4We stress that the model predictions are not fitted to the simulations -we simply adopt the boundary conditions (for values like  and  bh ) defined therein.

Solutions for the Accretion, Magnetic Field, and Disk
Structure Denote the thermal sound, Alfvén, circular (  ≡ √︁   enc (< )/ in terms of total enclosed mass  enc ), and effective "turbulent" speeds as   ,   ,   ,   , respectively.Motivated by the simulations, we will assume a disk geometry, and that the disk is magnetic-pressure-dominated, supersonically and trans-Alfvénically turbulent, efficiently cooling ( cool ≲  dyn ), and features angular momentum transfer dominated by some combination of Maxwell and Reynolds stresses (discussed below).We will neglect angular momentum transport via magnetically-driven outflows (these are negligible in the simulations).We validate the self-consistency of all these assumptions below.For simplicity, we focus on similarity solutions and denote the relevant O (1) dimensionless pre-factors with   ∼ 1 for different quantities .For generality, we will write our scalings in terms of some arbitrary orbital frequency Ω ≡  c /, but for the regime of greatest interest here we can safely assume the Keplerian limit Ω →  K / = √︁   bh /3 (i.e.  ≈  K = √︁   bh /).From these assumptions, we expect an effective viscosity5  visc ≈    2  /Ω giving rise to an inflow rate  = 3  visc Σ gas .Equivalently, the energy lost to turbulent dissipation, which is ultimately radiated away, is balanced by gravitational energy from accretion: so by definition 9     ≈ 2   (and recall   ∼   ∼   ∼ 1; for more detailed discussion see Pariev et al. 2003).
Since the disks are magnetically dominated, the disk scaleheight  ≈     /Ω ≡   /Ω √︁ 4  mid (  ∼ 1), in terms of the midplane density  mid ≡   Σ gas /2 , so (rearranging) we have  ≈ (  /Ω) 2 /(2   Σ gas ).From Eq. ( 1) we can re-arrange to write the gas surface density as Σ gas ≈  Ω/(3    2  ).In Paper II, it is shown that the magnetic-field strengths in the simulations arise from simple flux-freezing/flux-advection considerations, and indeed even the simplest-possible fluxfreezing argument, which predicts |B| ∝  2/3 in the midplane of the disk, works quite well at explaining the mean trend of |B()| and |B()|.We will therefore assume such a scaling here (though we discuss below the weak effects of adopting some different prescription for B()).Next, since we are interested in similarity solutions, we parameterize   ∝  − .Comparison with the simulations in Paper II and basic physical considerations suggest a plausible range of 0 ≲  ≲ 1/2, for which we will see that the solutions give a radial inflow velocity dominated "disk" are more like ∼ 1 pc (see Paper II), while ∼ 2 − 5 pc resembles more of a "free-fall zone." 5 Here, we introduce  2  which can have contributions O ( 2  ) from the Reynolds and (especially for primarily anti-correlated toroidal-radial fields as seen in the simulations) O ( 2  ) from the mean and/or fluctuating Maxwell stresses.These are comparable in the simulations and would appear with a fixed ratio here, so we can take  2  ∼  2  with their relative normalization and the appropriate O (1) prefactors subsumed into   .As noted below, this means that if the mean Maxwell   −   stress has the correct properties, the disks do not strictly need to be strongly turbulent, though they appear to be in practice., averaged over all simulation snapshots).We compare inflow rate  ( ) (here net inflow less outflow); magnetic field strength (with the simulation example showing the toroidal/poloidal/radial decomposition); gas density; scale height /; turbulent (1D direction-averaged rms ), thermal sound (  ) and Alfvén (  ) speeds; mean inflow velocity −⟨  ⟩; and effective Toomre  parameter (including thermal, magnetic, and turbulent support).Simulation quantities are volume-averaged in the midplane (|  | < ), and plotted versus cylindrical radius from the SMBH , from ∼ 80 au (the inner simulation accretion boundary) to ∼ 1 pc (exterior to which the physics become more ISM-like, with significant star formation in a GMC-like complex which is being tidally disrupted by the SMBH to form the accretion disk, biasing the statistics; see Paper I).For the models, we take the parameters  bh = 1.3 × 10 7  ⊙ ,  bh ≈ M ⊙ yr −1 (plotted; ≈ 100), and  ff ≈  BHROI ≈ 5 pc directly from the simulation, and assume trans-Alfvénic turbulence (  = 1,  = 1/3).We contrast a Shakura & Sunyaev (1973, SS73) model with  = 0.1 (thin dotted).The analytic model reproduces the key simulation properties reasonably well, especially compared to an SS73 model which predicts quantities like  differing by factors as large as ∼ 10 8 .
−  / K ≈   (  / K ) 2 which increases with .Therefore, it is convenient to normalize our solutions to an outer radius, which we call the "freefall radius"  ff , at which point   ( ff ) ≈ − K ; i.e. at  ∼  ff the gas is in near free-fall (at the Keplerian speed) onto the center.The solutions we derive therefore patch onto a solution of gas free-falling onto the nuclear region at  ≳  ff , as in the simulations.We can now write: (2)
Motivated by the results in Paper II, let us further assume that the turbulence/stresses are approximately trans-Alfvénic, being set by an initially turbulent field that is not order-ofmagnitude different from isotropic, so   ∼   .This is equivalent to taking   → (2/3  ) 1/2 ∼ 1 and  → 1/3, giving, Σ gas g cm −2 ∼ 0.01 where we define  ≡ /  edd ≡  S / bh as the accretion rate in Eddington units (with a Salpeter time  S ≡ 5 × 10 7 yr defined for a reference radiative efficiency ≈ 0.1),  7 ≡  bh /10 7 M ⊙ ,  ff,5 ≡  ff /(5 pc).The second scaling in each expression normalizes  to the gravitational radius of the BH,   ≡ /  with   ≡   bh / 2 , instead of normalizing  to  ff .6 We can immediately check the consistency of our assumptions:   ∼   ≲  K at all  <  ff (the regime of validity here), so the disk is thin/slim ( < ) and the accretion will not be magnetically "arrested"; the turbulence is trans-Alfvénic (by assumption); even at the ISCO (  ∼ 6),   / ∼ 0.04 ( 7 / ff,5 ) 1/6 ∼ 0.04  1/12 7 so for all interesting BH masses, it is reasonable to neglect relativistic corrections at this order; the total effective (magnetic + turbulent + thermal) Toomre  parameter  tot =  eff /  Σ gas (with  2 eff ≈  2  +  2  +  2  , and   small by assumption) is: so the disk is stable at all radii  <  ff (outside of  >  ff , ISM physics applies) so long as  < 3000  −1/4 7 -i.e. even for highly super-Eddington accretion rates.
6 If one wishes to scale to e.g.accreting binaries and stellar-mass BHs, it is helpful to note  7 / ff,5 ≈ 25 ( 1 / ff,day ) 2/3 where  1 = /10 M ⊙ and  ff, day is the orbital period at  ff in days.So Eqs. ( 5 We can also confirm that various more subtle properties observed in the simulations should hold for this model.Per the arguments in Paper II ( § 6.6 therein), if we assume the field |B| is initially dominated by some combination of radial/toroidal fields at ∼  ff (as seen in the simulations, and expected for free-falling/circularizing material tidally captured a turbulent ISM beyond the BHROI), then we would have fields with inverse coherence length   ,  ∼ 1/ at  ∼  ff .Given this, the model suggests that (1) the radial and toroidal fields will maintain   ,  ∼ 1/ as the flow migrates inwards (i.e. the field lines will be "stretched" inwards in the ,  plane, requiring  ln (−  /)/ln  > 0) so long as  ≳ −1/4 (true for all solutions we consider); ( 2 Comparing to the simulations in Fig. 1 (taking the simulation values of  bh , , and  ff ∼  BHROI directly, so there is no "adjustable parameter" to be fit here), the model Eqs.5-6 appears to describe the simulations remarkably well (at the order-of-magnitude level) and capture most of the important qualitative scalings (though of course there are obvious deviations from strictly self-similar, pure-power-law, steady-state behavior in the simulations).For comparison, the traditional Shakura & Sunyaev (1973) -disk model contrasted in Fig. 1 differs from the simulations by as much as factors of ∼ 10 8 in key properties such as the Toomre  or  or Σ gas .

Thermal & Radiation Structure
Note that because the disk is magnetically-dominated we did not need to invoke assumptions about its thermal properties to solve for its structure.However it is useful to do so for its own sake and to validate the consistency of our assumptions.Begin with the usual optically thick, thermal disk assumption:  ≈  gas, mid ≈  rad, mid with cooling flux  cool / ∼ 2      4 (where   is some function that accounts for optical depth effects, vertical stratification, etc.).Taking  cool / ∼  turb+mag /  we can solve for  or We can also solve for the ratio of cooling to dynamical or free-fall time  cool / dyn ∼ (3/4) (Σ gas  2  Ω)/(     4 ) ∼ (  /  ) 2 , sonic Mach number M  ∼   /  , or "equivalent Shakura-Sunyaev -parameter"  therm, equiv defined below, which for  ∼ 1/3 gives: We can immediately validate several assumptions: for any reasonable   ,  ≪ 1 is very small and depends weakly on radius so this is true at all radii from the freefall radius to the ISCO.7 Correspondingly, the turbulence is highly super-sonic (M  ≫ 1), and balanced by rapid cooling with  cool ≪  dyn .We can also calculate the "thermalonly" Toomre  therm ≡   /  Σ gas , which gives r −25/24 , so despite   ≪   the thermal  therm ≫ 1 at most radii (and even at extremely large radii of several pc, it only falls to unity or lower for  ≫ 10  −1/2 7 , which is consistent with the behaviors seen and discussed at radii > pc in Paper I).This is important in that it means the disk is only very weakly fragmenting even where e.g.collapse along field lines in magnetic field switches is possible, and becomes fully stable against e.g.gravito-turbulence at radii ≲ pc (for SMBH masses), as discussed and shown in Paper I (see also Pariev et al. 2003).
In the inner disk, to understand   and the vertical stratification (or lack thereof), first let us estimate .Scaling to the electron scattering opacity  es gives  ∼ 3000  ( ff,5 / 7 ) 1/3  −5/6  (/ es ), and estimating the absorption opacity as the Kramers opacity gives that the optical depth to thermalization is ≳ 1 where  ≫ 1.However, comparing the vertical mixing/transport timescale via turbulence,  mix ∼ /  ∼ Ω −1 from the above, to the radiation diffusion timescale required to establish vertical thermal stratification ( es /) ≪ 1.Thus we expect efficient vertical mixing to make the disk weakly thermally stratified (meaning that we do not expect a strong temperature gradient at || ≪ ), leading to   ∼ 1.We discuss this case further in § 5. Given these ratios, it also follows that the inner disk is stable against the usual viscous and thermal 7 If we used  −1  = 1 +  −1 + 3 /8 derived for a stratified LTE disk with  ≡   ( =  mid ) Σ gas ≲ 1, then this does not change our conclusions regarding  ≪ 1.For the outermost disk where  ≲ 1 is possible, this generically gives  decreasing towards smaller  from  ff (closer to what we actually see in the simulations at larger radii in Paper I).For  ≫ 1, if H − opacity dominates (as in the outer disk in Paper II), its strong  7.7 temperature dependence means again  decreases at smaller  from  ff ; at smaller , if Kramers and/or electron scattering opacity dominate and  ≫ 1 for this   , it only slightly changes our coefficients (giving e.g. ∝ r −1/9 instead of r −1/12 , not enough to change our conclusions).
instabilities, per the discussion in Begelman & Pringle (2007) § 4.1.3.In the outer disk, the surface densities are low and two things happen: first, from Eq. 6 we see as r → 1, Σ gas can be quite small and so even assuming full ionization with electron scattering opacities, one finds  ≲ 1 (let alone the absorption/Kramers opacity, which will be smaller).Moreover, the equilibrium temperature assuming blackbody-like cooling becomes quite low,  BB → 10 K r −3/4  1/4  1/8 7 , so the gas becomes primarily atomic at radii r ≳ (0.2−1) ×10 −4  1/6 7  1/3 -just as we see in our simulations (see Paper I).So the opacity drops significantly and the disk can even become multi-phase.At these radii the thermal physics can be influenced by optical depth effects, but are set by a more complicated balance of heating from external radiation (stars in the galaxy and the CMB, which at typical quasar redshifts may not be negligible), cosmic rays, thermal cooling, non-equilibrium chemistry, etc.These maintain a weakly-stratified, cool disk with  ≪ 1 per the discussion above (with warm atomic gas from a few hundred to ∼ 8000 K), and  rad ,  therm ≪  mag .
Finally, it is easy to verify from the temperatures and chemical structure above that the microphysical resistivity and nonideal MHD effects (e.g.ambipolar diffusion, the Hall effect, Ohmic resistivity) are negligible for all disk conditions of interest here, as expected (and verified in the simulations directly; see Paper II).

SENSITIVITY TO ASSUMPTIONS AND CAVEATS
It should be obvious that the scalings above represent a tremendously simplified model.We can see in Fig. 1 that in the "full physics" simulations, even over a modest dynamic range, are not precisely similarity solutions (the "effective power law index" of different quantities changes with scale); nor are they in exact steady state (power-law indices vary with time, and quantities like  are not exactly constant with  as required in a strictly steady-state model; see Paper II); nor, as discussed in Paper II, is the disk exactly axisymmetric (it has large coherent eccentricity and smaller-scale spiral modes).But in Fig. 1, we plot the analytic expectations from this model and show that they nonetheless do appear to capture the key qualitative behaviors and at least order-of-magnitude values of the most important quantities in the simulations.
We can also vary some of the detailed assumptions within our analytic model.From examination of the simulations we see  is not exactly constant in space or time and can vary between 0 ≲  ≲ 1/2.Varying over this range gives slightly different results but does not change any of our conclusions: for example, we still obtain  ≪ 1 with  ∝   with a very weak dependence on radius −5/28 <  < 3/28.We can vary   but so long as it is not extremely small (  ≪ 0.01), which would be inconsistent with the simulations motivating this model in the first place, we still obtain similar disk structure and  ≪ 1,  ≫ 1.We have also considered parameterizations of the form   ∝     1−  K for some , but this is just equivalent to a narrower range of 1/3 <  < 1/2 and   ∼ 1.We have also re-calculated the entire model assuming the alternative flux-freezing relation discussed in Paper II,   ∝  −1 (which is also broadly consistent with the simulations).For a given , assuming   ∝  −1 gives us identical power-law scaling for   / K and Σ gas ∝ r −3/2+2  , with  ∝ r −3/4+4  /3 , and for e.g.  / K ∼ / modifies the r dependence to r1/2−2  /3 .Comparing this and/or Σ gas to the simulations would, with this model, favor something more like  ∼ 3/8, modifying   / K to be ∝ r1/4 instead of ∝ r1/6 (and Σ gas ∝ r −3/4 ,   /  ∝ r −1/8 ), all fairly weak effects.8But even for this closure and more extreme assumptions like  ∼ 0 or   ∝  −1 , which would introduce a stronger (though still sub-linear) dependence of   / K on r (and deviate more notably from the simulation behaviors which motivate these models), this would not change our qualitative conclusions at large radii, but would extrapolate inwards differently and would lead to a thinner, more thermal-dominated disk at small  (approaching the inner disk and ISCO).This could produce a limit more like that discussed in Begelman & Pringle (2007) and § 4.2 below in the inner disk (with the parameter   / BP07 defined in § 4.2 falling to ≲ 1 as r → 0).
More generally, for power-law solutions like those here, we can consider the generalized "closure" relation   ∝  2 ∝   (our default model takes  = 4/3), and note this is equivalent to any arbitrary power-law function In Eq. ( 2), this modifies   / K →   r (7−5−4 +4 )/(2+2) (while   / K and Σ gas scale identically).Per the discussion above, if we still assume trans-Alfvénic turbulence (and/or a broadly similar range of ), the model predictions are quantitatively modified but remain qualitatively similar over a wide dynamic range in r for any closure in the plausible range 1 ≲  ≲ 2. For  ≳ 2 (depending on the exact value of  and the boundary condition   ),   / K and / will rise sufficiently steeply as  → 0 such that they could exceed unity before reaching the ISCO -this implies the disk would transition to a more MAD-like state at small radii (unless  ≳ 1/2, in which case they would instead transition back to free-fall solutions).For  ≲ 0.8 or so, such that the Alfvén speed   ∝  (−1)/2 decreases at higher densities and smaller radii within the disk,   / K will fall sufficiently rapidly as  → 0 so that  increases and exceeds unity before reaching the ISCO.This would imply a transition either to a SS73-like disk at small- or perhaps an intermediate disk more akin to the Begelman & Pringle (2007) model discussed in § 4.2.
Briefly, it is worth noting that a particularly simple and interesting model variant arises if we assume  = 5/3 and trans-Alfvénic turbulence, which gives: While not quite as good a fit in Fig. 1 (particularly at  ≳ 0.01 pc), in this case the flow is truly self-similar, is entirely parameterized by one dimensionless O (1) number   (note  ff entirely factors out of the solutions), and maps onto free-fall with   ∼  approaching the horizon.Given that the disk mass  gas (< ) <  bh within the BHROI (by definition), and that accretion cannot be faster than dynami-cal (  ≲    gas ( BHROI ) Ω( BHROI )), it becomes basically impossible to have  < 1 for any physical inflow rate at the outer boundary (the upper limit is set by supply from galactic scales).Also here  ∼  cool / dyn ∼ 1/M 2  ∼ 1/ therm, equiv ∼ 10 −5  −2  (   5 /    2 7 ) 1/4 ( 5 ≡ /5 pc) is similar to our default model in the outer disk and decreases weakly (for reasonable   ) as  → 0, so the disk becomes more magneticallydominated.All of our qualitative conclusions regarding consistency and relative importance of different terms for our "default" parameterization apply to this solution as well.
Per the discussion in § 2.2 and Fig. 1, the least robust part of our analytic model in general (comparing to the simulations) is probably the description of the thermal disk properties at the largest radii, although as noted there, this has no effect on the predicted accretion rates or mass profile since  ≪ 1.This owes primarily to the facts that LTE is a poor approximation, the disk is highly super-sonically turbulent, and (at large radii) the optical depths can become modest and the disk is no longer mostly ionized.More detailed thermochemical modeling is therefore warranted.
3.1.How Generic vs. Fine-Tuned Are the Parameters?Given the discussion above ( § 3), the qualitative behaviors in the model are quite robust to details of e.g. the detailed profile slopes (), scaling of   , etc., and as argued in Paper II the key dynamics and magnetic structure of the disk are set by simple flux-freezing/advection considerations and essentially independent of the detailed thermal structure.If those are general statements, then so long as the outer boundary conditions -encapsulated in our analytic model above by two parameters  ff and   -have "reasonable" values for a given  bh and , and are not especially "fine-tuned," it suggests this scenario may be robust.The "fine-tuned" question can be addressed immediately, as the dependence of key properties (/,   /  , etc.) on  ff and   is quite weak (see Eq. 2-6), so no fine-tuning is required.
Regarding "reasonable" values, recall  ff represents the radius at which the solution maps onto approximate free-fall onto the BH.It is notable, then, that our reference value of  ff , taken from the simulation, is very similar to the BHROI  BHROI ≡   bh / 2 gal where  gal is the velocity dispersion of the galaxy center/nucleus/bulge region.Indeed this is exactly what we expect (order-of-magnitude) for  ff if the accretion is ultimately driven by gravitational capture of gas from ambient complexes outside the BHROI moving at a characteristic velocity ∼  gal (precisely what is seen in the simulations, see Paper I).If we further assume that BHs lie on the observed galactic  bh −  gal relation (Kormendy & Ho 2013), we have 7 , which we can insert in the scalings above to eliminate  ff , making the dependence on  bh even weaker.
Next, consider   , which normalizes the magnetic field strength.For any system where the inflow velocities are trans or super-Alfvénic on large (galactic) scales (outside the BHROI),   ∼ 1 is a natural expectation.Moreover, recall we have only assumed flux-freezing, with  ∝  2/3 , so   is equivalent to the normalization of this relation.As demonstrated explicitly in Paper II, for the parameters in these simulations (these values of  bh , , and  ff ),   ∼ 1 is equivalent to  ∼ 8 G (/cm −3 ) 2/3 -i.e.completely typical ISM values (Crutcher et al. 2010;Ponnada et al. 2022).More generally, for arbitrary  and   if we assume  ∝  2/3 ≈  10 10 G (/cm −3 ) 2/3 strictly extrapolated from the ISM, we have   ∼ (2/3)  1/7  6/7 10  3/14 ff, 5  −5/14 7 (∼ 0.7  1/7  −1/4 7  6/7 10 if we assume  ff traces the BHROI, or ∼ 1.8  −2/7 7  6/7 10 if we assume  traces free-fall of an enclosed gas mass ∼  bh at  ff ).In other words, even if there were no amplification mechanism "driving" the system towards trans-Alfvénic turbulence (though in Paper I and Paper II we show there are multiple such mechanisms), we would still expect   ∼ 1 given the observed scalings of ISM magnetic field strengths.And to qualitatively change our conclusions regarding e.g. ≪ 1 or , we would need to lower the extrapolated  (e.g. 10 ) by at least two orders of magnitude -much lower than any reasonable models or expectations in the ISM.Further, even if this somehow happened, it does not necessarily mean the disks would not be magnetically-dominated, but instead, with weaker fields (such that e.g.|  | ≪ | BP07  | in Eq. 13) could lead to conditions more like the MRI-dominated Begelman & Pringle (2007) analytic disk models.
Formally speaking, there are of course physical criteria or conditions implicit in our assumption of how the magnetic fields are amplified (e.g. the models we consider with  2 ∝   where 1 ≲  ≲ 2 or  ∝ 1/).This includes, for example, the assumption that turbulent damping or buoyant escape of  will not strongly damp the midplane magnetic fields (compared to their rapid growth/supply via advection and flux-freezing).These and related qualifications are tested and discussed explicitly in the simulations at length in Paper II (see e.g. the discussion in § 6.6 as well as § 5.2 & 9.5 therein), where it is also shown that the analytic models here are internally self-consistent in these respects.Given that we simply adopt these relations here, we refer interested readers to the more detailed tests therein.
So given these considerations and the weak mass dependence, we expect our assumptions and model conditions to reasonably apply to essentially the entire SMBH mass range ∼ 10 4 − 1010 M ⊙ at high accretion rates (  ∼ 1).The upper limit (in accretion rate) to the validity of the assumptions above may be set by the stability condition,  < 3000  −1/4 7 -still allowing for extremely high accretion rates.The lower limit is less obvious but is likely set by similar physics as with a SS73-type disk; as  decreases, so does  mid ∝ Σ gas / ∝ , and the disk becomes more optically-thin and lower-density.Recall that even at  ∼ 1, Σ gas is significantly lower than a SS73 disk with the same , so this means the disks should approach the RIAF-type limit at some  ≪ 1 where the cooling becomes too inefficient to maintain the temperatures and dissipate the gravitational energy of the inflowing gas.This will lead to a completely different thermal structure of the disk.Naively comparing the optically-thin H+He cooling rates for a fully-ionized plasma to the values implied in the scalings above, one finds that this still requires  ≪ 0.1, but more detailed models would be needed in this limit.

What Drives the Turbulence?
In this paper, as in e.g.SS73, we are intentionally agnostic to the source of the turbulence, but simply assume it is trans-Alfvénic, motivated in part by the simulations in Paper II and by simple equipartition/non-linear saturation arguments.As discussed above, our conclusions are generally independent of this source, so long as it can sustain broadly order-ofmagnitude trans-Alfvénic turbulence with Alfvén Mach number M  not too strongly dependent on radius .In Paper II, we study the turbulence in much more detail, demonstrating that these statements are true in the simulations, and that there are many viable drivers of the turbulence including (but not limited to): inward propagation of large-scale/global gravitational  = 1 modes, gravito-turbulence in the outer disk, adiabatic compression/advection of turbulence with the accretion flow, and the Pessah & Psaltis (2005); Das et al. ( 2018) "Type II/III" or "SSMI/SHMI" radial magnetic buoyancy instabilities.Other buoyancy instabilities may be present, though the Parker-type modes discussed in Johansen & Levin (2008) are likely suppressed owing to the combination of strong turbulence, low-, and large / making their wavelengths ≳  (Paper II).Even the "traditional" linear MRI could operate in principle (despite the very low-) in regions where  ln   / ln  is sufficiently close to −1.9However, it is difficult to uniquely separate and identify a single "driver" in such highly non-linear, multi-physics simulations.
It is not even totally obvious that a turbulent driving mechanism is strictly necessary.For our purposes here, the angular momentum transport comes from some arbitrary combination of Maxwell+Reynolds stresses with  visc ∼  2  /Ω.This can come from a trans-Alfvénic turbulent Reynolds stress, or turbulent Maxwell stress, but in Paper II we show that in the inner disk (while there certainly is trans-Alfvénic turbulence present) this can actually be dominated by the mean-field −⟨    /4⟩ Maxwell stress as the toroidal field amplifies from the radial field producing a   −   anti-correlation giving rise to  visc ∼  2  /Ω.We also discuss therein that the turbulence itself may be largely advected from larger scales, rather than strictly locally amplified.Of course, one could have additional angular momentum transport via effects like outflows, but the simulations in Paper II find this is negligible on the scales we compare in Fig. 1.
Given this and our discussion of different variations for   and  above, it is worth discussing how different proposed saturation mechanisms for either the turbulence or toroidalradial mean field might scale.Since instabilities like the SSMI/SHMI are powered by the mean field magnetic tension, non-linear saturation around equipartition (M  ∼ 1) seems reasonable.But Begelman & Armitage (2023) propose some alternative analytic hypotheses for the saturated amplitudes of the MRI and/or Type III/SHMI instabilities in strongly-magnetized, toroidal-field dominated disks based on quasi-linear theory, or we could analogously speculate that saturation occurs when the linear growth time of these instabilities becomes comparable to other timescales such as the weak turbulence cascade time or the gas accretion timescale.We have re-derived our model scalings assuming each of these saturation scenarios in turn, for any given  = 4/3 − 5/3, assuming either the SHMI or MRI dominates the driving: these are all equivalent to very small changes in  (compared to our "default" trans-Alfvénic assumption) within the range −0.05 ≲  ≲ 0.07.Notably, this range of  is much smaller than the ad-hoc variations we already discussed above and showed do not change our conclusions.Thus, for the purposes of a global disk solution, most plausible saturation scenarios produce similar conclusions and, given the uncertainties, we do not feel it helpful to speculate further here.Future work using idealized numerical simulations will help to better constrain this physics.

COMPARISON TO OTHER ACCRETION DISK MODELS 4.1. Comparison to a Shakura & Sunyaev 𝛼 Disk
The above differs from the classic Shakura & Sunyaev (1973) (henceforth SS73) -model in many details, but the most important and fundamental is that SS73 assume a gaspressure dominated disk ( ≈   /Ω) with  ≈ (  /  ) 2 ∼ constant < 1.This in turn leads to many other differences: for example the disk structure and accretion depends on the opacity structure, the disks are strongly vertically stratified,  > 1 in the midplane,  cool ≫  dyn , and the surface densities Σ gas and Toomre  parameter can be very different.
To gain some insight, consider the "equivalent"  parameter  therm, equiv ≡    2  / 2  from our model above ( § 2.1) -defined as the  parameter in the SS73 model that would give rise to the same accretion rate as our fiducial model for a given   or Σ gas .Note that the thermal assumptions in § 2.2 are similar (where the disk is optically thick) so the zeroth-order temperature is similar.This means  therm, equiv ∼ 1/ as defined above -so we see  therm, equiv ≫ 1 -i.e. the accretion is much faster, for a given Σ gas , than the fastest-possible SS73 model (this is clear in ⟨  ⟩ in Fig. 1).This in turn means that for a given  in steady state, Σ gas at a given physical radius is lower in our fiducial model by a factor of ∼ ( SS / therm, equiv ) ∼  SS  relative to the Σ gas from an SS73 model with the same  and some  SS < 1.We see this in the density plotted in Fig. 1.Meanwhile /,   / K , and   / K are larger in our fiducial model relative to SS73 by factors of ≳  −1/2 ≫ 1.At a given , this means  (which is dominated by just the thermal contribution in the SS73 model) is larger in our fiducial model relative to SS73 by a factor ∝  eff /Σ gas , i.e. / SS ∼  −3/2  −1 SS , also shown directly in Fig. 1.For a canonical  SS ∼ 0.1, this can amount to a factor of ∼ 10 8 increase in !As a result, if one assumes a SS73  disk with  ∼ 0.1 − 1, then even for  below/around the Eddington limit one obtains the well-known result that the disk should fragment ( SS ≲ 1) outside of a few hundred gravitational radii.For the conditions in the simulations shown in Fig. 1 ( 7 ∼ 1 and super-Eddington  ∼ 20 M ⊙ yr −1 ) the prediction would be  SS ≲ 1 outside of a radius as small as ∼ 10 −5 pc.Instead, we see in the simulations and model that, at all radii out to a few pc,  ≫ 1 and even  thermal ≳ 1 (where we do not include magnetic or turbulent support but still see larger  because of the greatly reduced surface densities).
This also makes it clear that the simplified models we propose are not technically  models, since the effective  is not constant, though it does vary relatively weakly.
One non-intuitive result, comparing our default models to SS73 in Fig. 1, is that the rms magnetic field strength in absolute units is actually larger in an SS73 disk (assuming   ∼   )! Indeed, computing the ratio from the analytic expressions above and from SS73:10  SS73 / ∼ 12  9/80 7 r1/48  1/20 0.1  −3/40  −1/16 ff, 5 -i.e. an order of magnitude stronger fields appear in an SS73 disk, with extremely weak dependence on any parameter.But of course the models here feature much larger Alfvén speed (  ≡ |B|/ √︁ 4) and orders-of-magnitude lower  ∼  thermal / B ∼  2  / 2  com-10 We formally estimate the magnetic properties from SS73 by assuming  primarily arises from Maxwell stress, appropriate for e.g.MRI-driven transport.Adding e.g. a mean poloidal field to SS73 does not change our discussion.
pared to SS73:   ≫   in the models here, while   ≪   in SS73.But because of the enormous difference in density of the models (where  and  thermal in SS73 are larger than here by factors ∼ 10 6 ), even a relatively large |B| in SS73 (needed to explain some  ∼ 0.1) still gives very low   ≪   and therefore negligible magnetic disk support.This difference in |B|, however, has important consequences for formation of these disks, helping to emphasize the point made in Paper II that the magnetic flux from the ISM needed to supply these disks is not particularly large (corresponding to completely "mundane" ISM field strengths).This is also important to keep in mind when considering potential observable signatures.

Comparison to the Begelman & Pringle Model
To our knowledge the closest analytic model in the literature (of those which attempt to actually predict the scalings of () in the disk) to that proposed herein is the model discussed in Begelman & Pringle (2007) (hereafter BP07), which assumes an -like model for a disk with primarily toroidal magnetic fields and  ≪ 1.While there are a number of differences in detail regarding some of the assumptions above (e.g.regarding the thermal structure, accounting for turbulence and its effects on stratification, etc.), these mostly have weak effects on the qualitative behaviors and conclusions of the models.
The most important physical difference between the models proposed here and BP07 is the assumption of what sets the magnetic field strength.Paper II shows that in the simulations motivating our assumptions, this appears to be given by simple flux-freezing and advection of ISM magnetic flux, which gives rise to an analytic scaling broadly similar to |  | ∝  2/3 mid or ∝  −1 , and introduces the concept of the free-fall radius above because the magnetic flux, and hence the disk properties, explicitly depend on the outer boundary conditions (the magnetic flux being fed into the disk).In BP07, the authors essentially assume the opposite limit, where the "seed" or "initial" field is some small/trace field, which is then amplified by the MRI, until the MRI saturates and its linear growth ceases at approximately the threshold estimated in Terquem & Papaloizou (1996); Pessah & Psaltis (2005):11   ∼  sat, mri ∼ √    K .As shown explicitly in Paper II in multiple examples,   ∼  sat, mri is not, in fact, a good description of the magnetic fields in the simulations: their amplitude is everywhere orderof-magnitude larger than this threshold  sat, mri (  ≫  sat, mri or equivalently |B| simulation ≫  sat, mri ), they obey a different power-law scaling with BH-centric radius and density, and do not exhibit any clear dependence on sound speed at a given radius and density (unlike  sat, mri which depends directly on   ).But as shown in Paper II, all of these results are a natural expectation if the field strength determined by flux-freezing and the advected field is already above this nominal critical value (though it is also possible that MRI growth continues in more complicated, non-linear, radially-stratified systems like those here).We can validate this assumption in our analytic model above, calculating the ratio at each radius of the predicted field strength |  | to this threshold value  BP07 ≡  sat, mri (which should be ≳ 1): ∼ 4 ∼ 20  0.09 7  −0.06 r0.1 .
So indeed we predict order-of-magnitude larger fields, and this conclusion is relatively insensitive to the conditions and other detailed assumptions above.12Despite this, most of our qualitative conclusions are similar to BP07 -indeed, any magnetically-dominated model where   /  and   / K do not vary too strongly as a function of  would give many broadly similar conclusions to our proposed model regarding the scalings of , , Σ gas , etc.And while the detailed power-law dependencies are different, the "weak" dependencies remain "weak" (e.g.terms scaling here to the 1/16 power might scale as 1/36 or 1/12 in BP07).
Nonetheless, there is one important consequence of the order-of-magnitude difference as r → 1 in Eq. ( 13).Recall, the effective "" scales as ∝  2  ∝  −1 ∝  2 , and the effective Σ gas scales as 1/, and  ∝  eff /Σ gas ∝  3/2 .So at small radii (approaching the ISCO), the difference with BP07 is not so dramatic, but at large radii -noting the ratio in Eq. ( 13) increases with  -a factor of ∼ 10 in magnetic flux translates to factor ∼ 100 lower Σ gas and ∼ 1000 larger .Basically, in BP07, Σ gas decreases more slowly (relative to our model here) at large  because   ∝ √    K declines more rapidly.As BP07 note, this leads to a "bottleneck" where they predict their disk models would fragment before reaching the BHROI (roughly equal to our  ff here) for a critical  ≲ 1.8  3/4 70  −5/8 7 (where  70 scales the temperature to ∼ 70 K allowing even for the temperature at large radii to be set by complex ISM physics making it warmer than predicted by the blackbody scaling, as discussed above and in BP07, and we rescale from their definition of  to ours).Whereas we see from Eq. ( 7) (and directly in the simulations) that  > 1, and even  thermal > 1, out to the BHROI for  as large as ≳ 1000.
In short, while the BP07 model is already much less prone to fragmentation than a thermal-dominated model like SS73, by allowing for larger fields set by flux-freezing at larger radii (as opposed to only the field strengths at which the MRI saturates according to local linear theory in a non-stratified disk), our model allows the disk to support factor ∼ 1000 larger inflow rates at all radii from the BHROI inwards.

Comparison to Magnetically Elevated or Arrested Disks
In Paper II, we discuss and demonstrate extensively how the simulation disks are qualitatively and quantitatively different in their basic properties from historical models of magnetically "levitated" or "elevated" or "arrested" disks (see e.g.§ 9 therein).Rather than repeat all of the discussion and detailed comparisons therein, we briefly summarize how this also applies to the analytic models here.
First, consider levitated or elevated disks.The defining feature of these models (see e.g.Miller & Stone 2000) is that the magnetic fields are weak ( ≫ 1) in the disk midplane (where the structure is akin to an SS73-type disk), but rise in importance at larger vertical heights above the disk until  ≲ 1 at a few scale-heights (|| ≫ ).Clearly this is opposite our fundamental assumptions for  and the disk midplane structure, and produces predictions for disk properties like ,   ,   / K , , etc. akin to SS73 and completely different from our model in Fig. 1.Moreover in Paper II we show that the simulation disks are actually both weakly and inversely stratified ( is lowest, not highest, in the midplane) relative to the elevated assumption.
Second, consider magnetically arrested (MAD) disks, defined by dominant mean large-scale vertical/poloidal magnetic fields (B ≈ ⟨  ⟩ ẑ) whose magnetic tension balances gravity and prevents accretion (Bisnovatyi-Kogan & Ruzmaikin 1976;Narayan et al. 2003).It is easy to verify that our disk models never satisfy the standard MAD criterion of ⟨  ⟩ 2 mad /8 ≳   bh Σ gas /4  2 , even if we assumed the field was entirely in a coherent mean poloidal configuration.To see this, note the criterion can be re-arranged (using our implicit assumption for magnetically-dominated disks that / ≈   / K ) as ⟨  ⟩ 2 mad /|B| 2 ≳  K /  -which can only be realistically satisfied if both   ≳  K (the Alfvén speed rises to be super-Keplerian) and |⟨  ⟩| ≈ |B| (the field is primarily in a coherent mean vertical component).As noted above, even ignoring the field geometry considerations, achieving the former requires much more extreme amplification than seen in the simulations/assumed here ( ≳ 2, pushing towards the upper limit of what is possible for perfect flux-freezing in a laminar cylindrical vertical-field-only disk), together with the assumption that the turbulence becomes weak (  ≪   ) rather than remaining trans-Alfvénic (since in the trans-Alfvénic limit with such strong amplification, the solutions would extrapolate to free-fall, promoted by Maxwell+Reynolds stress, rather than arrested/static solutions supported by magnetic tension).Moreover, as discussed in Paper II at length, in the simulations motivating our models here, the fields are primarily toroidal, with secondary radial and quasi-isotropic turbulent components (so the mean vertical field |⟨  ⟩| ≪ |B| is generally the weakest component of B).This has important qualitative effects, changing the field amplification, role of magnetic pressure setting the vertical height and midplane structure of the disk, and the sign of the magnetic torques/Maxwell stresses (making the fields promote accretion, rather than suppress it) -all of which are implicit in our model assumptions.And unless one begins from boundary conditions with B ≈ ⟨  ⟩ ẑ, it is essentially impossible to generate this state via amplification within the disk (one cannot generically amplify the mean vertical field from trace to stronger than the turbulent field |⟨  ⟩| ≫ ⟨|  |⟩).
Finally, as discussed in detail in Paper II, although the criterion for a disk becoming MAD-like "eventually" is sometimes stated in terms of a critical magnetic flux at some arbitrarily large radius (e.g.some Φ crit ∼   2 |B| crit ∼ G pc 2 , as opposed to the local MAD criterion discussed above), this scaling makes a number of further assumptions -in addition to both the strong amplification ( ≳ 2) and mean-vertical field dominated (B ≈ ⟨  ⟩ ẑ) assumptions above -which are not valid in our model.Most notably, it also assumes that the field plays a negligible role in the disk structure (outside of the "arrested" radius) and so the disk can be treated as a  ≪ 1 SS73-like -disk, which we have shown is dramatically different from our model predictions.So while it is certainly possible that some accretion disks could, in general, approach a MAD-like state, it requires qualitatively different boundary/initial conditions, field amplification, and accretion disk structure from what we assume and model herein.Pariev et al. 2003 We note that the power-law scalings with radius (logarithmic slopes) that we obtain represent a generalization of those derived in Pariev et al. (2003) (P03): we obtain their scalings if we restrict to trans-Alfvénic turbulence, neglect advective radiation transport (assume   =  −1 ), and set  →  P03 /(3 P03 − 3) and  → 3 − 2 P03 in terms of their arbitrary parameter  P03 defined by their assumption that  turb () =  10, P03 (/20  ) −  P03 .This is despite the fact that P03 considered very different physical conditions (e.g.much lower accretion rates) and motivations, and assumed only "pure turbulent" magnetic fields (negligible flux freezing and mean fields).This similarity arises because the trans-Alfvénic assumption gives, by definition, |B 2 | 1/2 ∝ |⟨B⟩|, so should produce the same dimensional scalings for appropriate choice of  and .But in addition to extending what is presented in P03 by varying these assumptions and connecting the disks to their outer boundary at  ff or the BHROI, our model attempts to actually predict quantities like  P03 and  10, P03 (which are left arbitrary in P03) from the simulations and physical arguments.It is also worth noting that P03's discussion, which considered the absolute normalization of the magnetic field (in Gauss) around the ISCO required by the self-consistency of their assumptions, does not explain why  10, P03 would or should naturally lie in their "valid" range for a given  and  BH .By introducing the outer boundary conditions for the disk, we show that this is automatically ensured so long as   is not too small ( § 3).

LIMITING REGIMES OF THE MODELS AND
BEHAVIOR AT SMALL RADII The most basic requirement for accretion disks to be in the state predicted here is that their outer boundary conditions satisfy our assumptions: namely, gas is accreting which is relatively thermally "cold" (  ≪  K ) and well-magnetized ( ≲ 1) with (as detailed in Paper II) a steady supply of toroidal and radial magnetic flux.Beyond this, as noted above, at sufficiently large , the outer disks predicted here can still become Toomre unstable, but this will not occur until  ≳ 3000  −1/4 7 .At sufficiently low  (probably  ≪ 0.1), the thermal structure may be modified by inefficient cooling.Provided these conditions are met, the models here are designed to be extrapolated out to the BHROI, beyond which point, by definition, the physics must change since the potential changes, star formation becomes non-negligible, the thermal structure is set by diffuse ISM heating and cooling processes rather than accretion heating, etc.And indeed the simulations to which we compare are focused on this outer disk, from the BHROI at   ∼ 10 7 (and beyond) to   ∼ 600.But the behavior in the inner disk down to the ISCO/horizon (  ∼ a few) is less obvious.
Approaching the horizon, one has the boundary condition   ∼  K ∼ , while our "default" solution (Eq. 5) predicts   / K → (  / BHROI ) 1/3 ∼ ( gal /) 2/3 ≪ 1.Thus there would have to be some transition in physics between the radii of interest here and the near-horizon or "free-fall" region (i.e. one limitation of the analytic models here is that they must be modified near the inner boundary/horizon), but that is not particularly surprising (and similarly holds if one attempts to extrapolate most -disk models, including SS73 or BP07, from large radii to the horizon).Interestingly, the simpler solution for slightly "stiffer" magnetic amplification ( = 5/3; Eqs. ( 10)-( 12)) extrapolates to   ∼  K (for natural order-unity   ).Thus, the transition in the inner disk could simply amount to slightly more efficient field amplification.
At high accretion rates approaching or exceeding Eddington (  ≳ 1), it is also important to examine the role of radiation pressure in the disk.For the models in § 2, the ratio of midplane radiation pressure  rad ∼ (1/3)   4 to "other" (thermal+magnetic+turbulent) pressure  other is just 2).If we assume the standard optically thick diffusive-radiation-transport limit,  −1  ∼  ∼ Σ gas , then this is negligible in the outer disk:  rad / other ∼ 10 −6  ( 7 / ff, 5 ) (/ es ) r −7/6 (in addition, for the conditions in Fig. 1, the outer disk is primarily atomic, so  ≪  es ).But  rad / other will increase inwards until it crosses unity at some13   ≲ 150 ( ff, 5 / 7 ) 1/7 (/ es ) 6/7  6/7 .Note this is much is closer to the horizon than what one would obtain for a standard SS73 disk (where  rad ≳  other at  SS73  ≲ 9000  0.8 ( 7 ) 0.1 (/ es ) 0.9 ) owing to the much stronger magnetic pressure.But it is still well outside the horizon for large .Interior to this radius, naively the solutions would be modified (with radiation pressure "puffing up" the disk and setting ), and we could derive an appropriately modified solution which retains the key characteristics of our models herein (akin to the "region (a)" model in SS73).However, at even smaller radii, approximately (up to an O (1) constant) where the disk luminosity  disk ∼   2    4 eff exceeds  Edd , it becomes impossible to support the disk even with  ∼  and the radiation pressure forces exceed gravity, if the heat transport is diffusive with electron-scattering opacity (the minimum scale height effectively scales as (/) min ∼ (  −1  / es ) ( disk / Edd ), with  es ≡ Σ gas  es ).For  −1  ∼  es , this would occur at   ∼ 4 , well outside the horizon/ISCO for  ≫ 1.
However, the ratio  rad / other in the diffusive limit (  −1  ∼ ) in the models above is just ∼   /(/) -i.e. the same as the ratio of the advective (∼   ) to diffusive (/) heat transport rates.So in the models here, at the same radii where  rad would naively begin to exceed  other , the (vertical/local) heat transport becomes advective ( adv ∼   ≳  diff ∼ /), so the disk should become vertically mixed and weakly stratified with  −1  ∼ 1 ≪ .If we assume the maximum-possible stratification scales as  −1  ∼ /  , the disk would converge to  rad ∼  mag ∼  turb ≫  thermal (if instead  −1  → 1, then the disk remains with  mag ∼  turb >  rad ≫  thermal at all ) and the solutions we describe above would only be weakly modified by a systematic O (1) prefactor, and so  <  (from the original solutions) and continued accretion could be sustained down to near-horizon scales.Note that if we continue to assume diffusive radiation transport at  diff ∼ /, then the radial advective speed   ∼  2  / K becomes larger than the diffusive speed at   ≲ (10 − 20) (/ es ) , suggesting the system would become radiatively inefficient akin to the standard slim disk solution.However, because we just argued the radiation is not "locked" to the diffusive speed at these radii and  t >   always, the actual radiative efficiency of the inner disk when  ≫ 1 (and thus whether the quasars can appear super-Eddington or not, in the observable luminosity sense) will be set by some non-linear competition between accretion (radial advection) and convective/turbulent escape (vertical advection) of radiation.Given the somewhat speculative discussion above, it clearly would be valuable to directly simulate disks like those proposed here on scales smaller than those resolved in the Paper II simulations, resolving the radii where the key transitions above should take place with explicit radiationmagnetohydrodynamics, in order to properly model and understand whether the heat transport becomes advective (in both the vertical and/or radial directions) and how this modifies the disk structure, thermal properties, and radiative efficiencies.This could be done in the Newtonian or pseudo-Newtonian limits for sufficiently large , given the characteristic   above.But obviously, properly modeling scales of order the horizon (the inner boundary) would require GRMHD simulations.

SUMMARY & CONCLUSIONS
We present a simple analytic similarity model for stronglymagnetized flux-frozen accretion disks; specifically, disks with midplane  ≪ 1, super-sonic turbulence, and efficient cooling.We validate this model against recent numerical simulations that found such disks forming around SMBHs accreting at quasar-level accretion rates from cosmological initial conditions, and show that it can reasonably explain the qualitative behaviors and scalings in the simulations.
We show that in this limit, for reasonable assumptions (e.g.trans-Alfvénic turbulence), the model is internally consistent and entirely specified by the boundary conditions, specifically the accretion rate and magnetic flux around the "free fall" radius (the outer radius at which gas is first captured into the accretion disk, and the model maps to a solution of free-fall onto the SMBH, which should occur near the BHROI).In the outer accretion disk, which extends all the way to the BHROI (exterior to which the physics must become more ISM-like), we show these disks should remain consistent and maintain gravitational stability ( ≫ 1) up to extremely large Eddingtonscaled accretion rates  < 3000 ( bh /10 7 M ⊙ ) −1/4 .We show that most of the predicted structure and scalings are insensitive to the detailed scalings adopted for the magnetic field structure or other parameter choices, provided the limits above are satisfied (and a sufficient supply of coherent toroidal/radial magnetic flux exists in the first place).
These disks have qualitatively different structure from a Shakura & Sunyaev (1973)-like  disk, most notably with  larger by factors up to ∼ 10 8 for the simulations to which we compare.They are also qualitatively distinct from magnetically "elevated" or "arrested" disks.They share similarities with the model outlined in Begelman & Pringle (2007), but differ in some key respects, most notably in that the magnetic fields set by accretion of ISM magnetic flux (from typical ISM field strengths) can be much larger than the upper limit assumed in Begelman & Pringle (2007) (which takes the field to arise from the local unstratified, weak seed-field MRI only).This difference leads to different scalings and allows factors of ≳ 1000 larger accretion rates to be stably supported.Indeed, as discussed in Paper II extensively, a remarkable feature of the simulations (and implicit in these analytic models) is that such hyper-magnetized disks do not strictly need to be "powered" by the MRI, and can feature magnetic fields much stronger than the often-quoted upper limit for linear MRI growth from Pessah & Psaltis (2005).
Finally, we discuss implications for the inner disk and regimes where this model may require extensions at either much lower accretion rates or much smaller radii (approaching the ISCO).At small radii, our model, which has been primarily focused on the outer disk given the simulations and questions of gravitational stability, may require modification due to general relativistic effects and strong radiation pressure if  ≫ 1 (where we argue the heat transport may become advective, enabling sustained super-Eddington accretion).This will be an important subject for future work.It will also be important to explore the limits of the boundary conditions (in terms of accretion rate and magnetic flux) where such disks can be sustained.
That said, the robust nature of the analytic scalings we derive and their relative insensitivity to the details of the boundary conditions, BH mass, and accretion rates, together with the fact that the first (and thus far only) cosmological radiation-MHD simulations to resolve quasar accretion disks saw such hyper-magnetized disks form naturally, all suggests such disks may be common around quasars in nature.Given that they naturally appear to resolve the decades-old problem of gravitational instability in the outer accretion disk, it seems crucial to study their properties in more detail.Their spectral properties, detailed internal structure, and behaviors at small radii and at extreme (low and high ) regimes should all be explored more thoroughly in numerical accretion disk simulations.The key difference between such disks and the overwhelming majority of historical work on accretion disk simulations is the boundary/initial condition assumed for the magnetic flux (discussed in detail in Paper II) -we discuss how this is qualitatively distinct from e.g.models of magnetically arrested or elevated disks.Our hope is that the simple analytic models in this paper can inform and guide such studies, providing some guidance for extrapolation, but also be refined and improved by such studies in the future.
Fig. 1.-Comparison of the model predictions here (thick dashed lines, § 2) to the numerical simulations of quasar accretion disks forming from cosmological initial conditions in Paper II (solid lines show mean values at each  at one moment/snapshot, shaded range shows the 90% inclusion interval for all gas at that and Σ gas ∼ 3000 (  ff,day / 1 ) 2/9 ) the supply of   from advection of radial flux (d    ∼       + ...) will maintain the   −   anti-correlation needed to ensure the Maxwell stress transports angular momentum outwards (Paper II, § 5); and (3) the turbulent resistivity, d turb    ∼ − turb  2     ∼ − η (  / K ) 2 Ω   with η ≲ 1, will not damp the fields signficantly faster than they are advected inwards, since |d turb  | cannot be ≫ |d advection  | ∼      ∼   (  / K ) 2 Ω (though they can be of similar magnitudes).For a much more detailed discussion and justification of these points, see Paper II.