††thanks: ∗E-mail<EMAIL_ADDRESS>
,

# Multi-Phase Thermal Structure & The Origin of the Broad-Line Region, Torus,
and Corona in Magnetically-Dominated Accretion Disks

Philip F. Hopkins1,∗ 1TAPIR, Mailcode 350-17, California Institute of
Technology, Pasadena, CA 91125, USA

###### Abstract

Recent simulations have demonstrated the formation of “flux-frozen” and hyper-
magnetized disks, qualitatively distinct from both classical $\alpha$ disks
and magnetically-arrested disks, as a natural consequence of fueling gas to
supermassive black holes in galactic nuclei. We previously showed that the
dynamical structure of said disks can be approximated by simple analytic
similarity models. Here we study the thermal properties of these models over a
wide range of physical scales and accretion rates (from highly sub-critical to
super-critical). We show there are several characteristic zones: a dusty
“torus”-like region, a multi-phase neutral and then multi-phase ionized, broad
line-emitting region interior to the sublimation radius, before finally a
transition to a thermal accretion disk with a warm Comptonizing layer. The
disks are strongly-flared with large scale heights, and reprocess and/or
scatter an order-one fraction of the central disk emission. As a result, this
simple accretion disk model predicts phenomena including the existence of a
dusty torus and its covering factor, geometry, clumpiness, and dust
temperatures; a broad-line-region (BLR) with its characteristic sizes and
luminosities and ionization properties; extended scattering/reprocessing
surfaces producing cooler disk continuum and apparently large observed disk
sizes; and existence of warm Comptonizing layers and hard coronal gas.
Remarkably, these properties emerge without our having to introduce new
“components” or parameters: they are all part of the accretion flow if the
disks are in the hyper-magnetized limit.

###### keywords:

quasars: general — accretion, accretion disks — quasars: supermassive black
holes — galaxies: active — galaxies: evolution — galaxies: formation

##  1 Introduction

Table 1: Common variables in this manuscript (others are defined throughout where relevant). $M_{\rm BH}$, $m_{7}$ |  BH mass $M_{\rm BH}\equiv m_{7}10^{7}M_{\odot}$
---|---
$\dot{M}$, $\dot{m}$ |  Accretion rate $\dot{M}$ and accretion relative to critical $\dot{m}\equiv\dot{M}/\dot{M}_{\rm crit}$ with $\dot{M}_{\rm crit}\equiv L_{\rm Edd}/0.1c^{2}$
$\epsilon_{r}$, $\epsilon_{r,\,0.1}$ |  Integrated (bolometric) radiative efficiency $L\equiv\epsilon_{r}\dot{M}c^{2}$ with $\epsilon_{r}\equiv\epsilon_{r,\,0.1}\,0.1$
$R_{g}$, $x_{g}$, $R_{\rm ISCO}$ |  BH Schwarzschild/gravitational radius $R_{g}\equiv 2\,GM_{\rm BH}/c^{2}$, $x_{g}\equiv R/R_{g}$, and ISCO radius (for non-spinning BH $R_{\rm ISCO}=3R_{g}$)
$R_{\rm BHROI}$ |  BH radius of influence $R_{\rm BHROI}\equiv GM_{\rm BH}/\sigma_{\rm gal}^{2}$ (with $\sigma_{\rm gal}$ the galactic velocity dispersion)
$r_{\rm ff}$, $r_{{\rm ff},\,5}$ |  Outer boundary or “freefall” radius into Keplerian potential $r_{\rm ff}\approx R_{\rm BHROI}$, defined by $r_{\rm ff}\equiv r_{{\rm ff},\,5}5\,{\rm pc}$ (with $r_{{\rm ff},\,5}\sim m_{7}^{1/2}$ expected)
$R$, $r$, $\phi$, $\theta$, $z$ |  Cylindrical $R$ (spherical $r$) radii, azimuthal $\phi$ (polar $\theta$) angles, and vertical height $z$ (disk-aligned and centered on the BH)
${\bf B}$, $B_{i}$ |  Magnetic field ${\bf B}$ and components $B_{i}$ (e.g. radial, toroidal, poloidal components $B_{R}$, $B_{\phi}$, $B_{z}$)
${\bf v}$, $v_{i}$ |  Gas velocity ${\bf v}$ and components $v_{i}$ (e.g. radial, azimuthal, vertical components $v_{R}$, $v_{\phi}$, $v_{z}$)
$\rho$, $n$, $\Sigma_{\rm gas}$ |  Gas 3D density $\rho$ or number density $n_{\rm gas}$ (in ${\rm particles\,cm^{-3}}$), and projected surface density $\Sigma_{\rm gas}$
$T$, $c_{s}$ |  Gas temperature $T=T_{\rm gas}$ ($T_{\rm rad}$, $T_{\rm dust}$ denote radiation and dust temperatures) and thermal sound speed $c_{s}\equiv\sqrt{k_{B}\,T/\mu\,m_{p}}$
$v_{A}$, $v_{\rm turb}$ |  Alfvén speed $v_{A}\equiv|{\bf B}|/\sqrt{4\pi\rho}$, and typical turbulent velocity $v_{\rm turb}$
$\beta$, $\mathcal{M}_{s}$, $\mathcal{M}_{A}$ |  Plasma $\beta\equiv c_{s}^{2}/v_{A}^{2}$ parameter, sonic $\mathcal{M}_{s}\equiv|\delta{\bf v}|/c_{s}$ and Alfvén $\mathcal{M}_{A}\equiv|\delta{\bf v}|/v_{A}$ Mach numbers
$H$ |  Gas disk vertical scale-height $H$ (defined within a given annulus $R$)
$v_{\rm K}$, $\Omega$ |  Keplerian circular velocity $v_{\rm K}\equiv G\,M_{\rm BH}/r$ and frequency $\Omega\equiv v_{\rm K}/R$
$t_{\rm dyn}$, $t_{\rm cool}$ |  Disk dynamical time $t_{\rm dyn}\equiv\Omega^{-1}$, and gas cooling time $t_{\rm cool}$ (at a given radius and temperature, etc.)
$\kappa$, $\kappa_{s}$, $\kappa_{a}$, $\kappa^{\ast}$ |  Total, scattering, absorption and effective $\kappa^{\ast}\equiv\sqrt{\kappa_{a}(\kappa_{a}+\kappa_{s})}$ opacities, with corresponding optical depths $\tau$, $\tau_{s}$, $\tau_{a}$, $\tau^{\ast}$
$Z$, $x_{e}$, $f_{\rm mol}$, $\xi$ |  gas metallicity $Z\equiv\tilde{Z}Z_{\odot}$, free electron & molecular fractions $x_{e}\equiv n_{e}/n$, $f_{\rm mol}$, ionization parameter $\xi\equiv({\rm d}\dot{N}_{\rm ion}/{\rm dA})/(nc)$

Active galactic nuclei (AGN) and quasars are powered by accretion disks around
supermassive BHs (Schmidt, 1963; Soltan, 1982), with accretion rates exceeding
$\gtrsim 10\,{\rm M_{\odot}\,yr^{-1}}$ in the most luminous sources. It has
been known for decades that there must be multi-phase gas structure around
such BHs, with large covering factors or vertical extent above the midplane
(e.g. Davidson & Netzer, 1979; Peterson, 1997; Krolik, 1999, and references
therein), in order to explain many observational features of their spectra and
variability. This includes, for example (1) optical narrow line regions (NLR),
probably from more “typical” interstellar medium (ISM) gas in the AGN host
galaxy; (2) the infrared dusty “torii” of clumpy gas with $\mathcal{O}(1)$
covering factor, cool and well-shielded enough to host dust at a range of
temperatures in “clumpy” structures at $\sim 0.01-10\,$pc (Antonucci, 1993;
Urry & Padovani, 1995; Burtscher et al., 2013); (3) the optical/UV broad
emission line region (BLR), believed to come from partially-ionized gas at
$\sim 10^{4}$ K with a range of densities but relatively modest range of
ionization parameters at $\sim 1-100\,$ld distances, reprocessing $\sim
10-20\%$ of the light in a thick-disk type geometry (Kaspi et al., 2005;
Peterson, 2006; Gravity Collaboration et al., 2018); (4) an extended,
relatively cool central thermal optical/UV continuum region perhaps dominated
by scattering (potentially related to any extended scattering surface from a
thick disk, or outflows/“warm absorbers,” or patchy cloud-covering) to explain
both the weak dependence of AGN SEDs on BH mass and microlensing observations
indicating emission or reprocessing of continuum at large radii (Laor et al.,
1997; Dai et al., 2010; Giustini & Proga, 2019); (5) “warm” ($\sim 1\,$keV)
Comptonizing structures covering the thermal continuum source needed to
explain the EUV/soft X-ray excess (Kubota & Done, 2018; Liu & Qiao, 2022) and
a more diffuse/extended “hard” ($\sim 10-100\,$keV) X-ray corona, plus again
large covering-factor scattering structures needed to explain the X-ray
reflection spectra (George & Fabian, 1991; Haardt & Maraschi, 1991; Marinucci
et al., 2018).

Many theoretical papers have been written about the potential origins of these
various multi-phase structures. But notably, the vast majority assume as a
starting point that the accretion disk itself is thermal and/or radiation-
pressure dominated, and so is qualitatively something akin to a Shakura &
Sunyaev (1973)-like $\alpha$-disk (hereafter SS73), whether geometrically
“thin” or supercritical/“slim.” While there are many variant accretion disk
models in the literature, (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 the fundamental assumption of SS73, that magnetic
pressure is small compared to thermal pressure ($\beta\gg 1$), is still most
often the “baseline” for both analytic studies, observational forward-
modeling, or setting up initial conditions for idealized accretion-disk
simulations. Importantly, in this category of thermal-pressure dominated
$\alpha$-disk models, the multi-phase structure observed must arise “outside”
of the accretion disk: the disk itself cannot be multi-phase (for reasons we
review below), nor if it were, could it reproduce the observed properties of
these different multi-phase structures (for example, their large covering
factors to the central source). Moreover, it has been known for decades that
thermal-pressure dominated disks are violently gravitationally unstable at
most of these radii, and predict gas densities and temperatures many orders-
of-magnitude different (at the same distances from the BH) to those inferred
for the BLR/torus/Comptonizing and scattering surfaces/coronae (Goodman,
2003). This in turn has led to other many other models for the origins of this
multi-phase structure, including popular ideas such as the BLR being a part of
an outflow or “failed wind”/fountain-flow (Murray et al., 1995; Krolik &
Begelman, 1988; Elitzur & Shlosman, 2006; Naddaf et al., 2021). But these have
their own challenges, and some variants of such models may even be ruled out
by recent observations finding a rotating, thick-disk geometry (Gravity
Collaboration et al., 2018; GRAVITY Collaboration et al., 2020, 2021).

Recently, Hopkins et al. (2024b, Paper I) and Hopkins et al. (2024c, Paper II)
presented the first simulations to self-consistently 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), multi-band 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 black hole radius of influence (BHROI)111Defined as the radius interior to
which the BH dominates the potential over its host galaxy of characteristic
velocity dispersion $\sigma_{\rm gal}$, or $R_{\rm BHROI}\sim G\,M_{\rm
bh}/\sigma_{\rm gal}^{2}$ (about $\sim 5\,{\rm pc}$ in the reference
simulations). is tidally captured by the SMBH of mass $M_{\rm bh}\sim
10^{7}\,M_{\odot}$ from larger-scale ISM gas complexes in the galaxy, and
free-falls briefly before circularizing to form an accretion disk with $Q\gg
1$ and little to no star formation or fragmentation on sub-pc scales. This
disk evolves in quasi-steady-state and sustains super-critical accretion (up
to $\dot{M}\sim 20-30\,{\rm M_{\odot}\,yr^{-1}}$) onto the SMBH for at least
$\sim 10^{5}$ disk dynamical times (the simulation duration). Crucially, in
Paper II where the disk properties were studied in detail, it was shown that
these “hyper-magnetized” and “flux-frozen” disks have $\beta\sim
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 orders-of-
magnitude less massive and support factor of $\sim 1000$ lower accretion rates
and higher star formation rates. The disks also have a flared structure
($H/R\sim 0.1-1$ at large radii) with weak vertical stratification owing to
trans-Alfvénic, highly super-sonic turbulence, which is sustained by rapid
cooling ($\mathcal{M}_{A}\sim v_{\rm turb}/v_{A}\sim 1$, with
$\mathcal{M}_{s}^{2}\sim 1/\beta\sim 1/t_{\rm cool}\,\Omega\gg 1$). With it
now possible to capture these multi-scale ISM-to-disk conditions, since Paper
I these flux-frozen disks have been seen in a number of other simulations of
distinct parameter spaces including sub-Eddington accretion of galactic hot
gas onto much more massive BHs (Guo et al., 2024) and accretion onto
intermediate-mass to $\sim 10^{6}\,{\rm M}_{\odot}$ BHs in dense star clusters
(Shi et al., 2024a, b), as well as more idealized simulations with the
appropriate initial conditions (e.g. Gaburov et al. 2012, as well as Squire et
al., in prep., Guo et al., in prep., Tomar et al., in prep.).

Here, we show that such disks could present a natural solution to the puzzle
of the origin of the various multi-phase structures or components of the AGN
ecosystem reviewed above. We take the simple analytic flux-frozen disk model
in Hopkins et al. (2024d) (Paper III; which was shown therein to reasonably
reproduce the simulation properties from Paper I-Paper II), and explore the
opacity and thermal structure of the disk in more detail.222Note that Paper
II-Paper III briefly discussed an extremely simple estimate of the
characteristic disk temperatures, but this largely ignored the details of the
real opacity and thermo-chemistry in the disk, as it was only intended to
demonstrate that for even the most extreme plausible opacity structure
producing the highest possible disk temperatures, the disk models above would
still be self-consistent in that they would still maintain $\beta\ll 1$ in the
midplane down to horizon scales. Moreover their discussion was primarily
focused on the largest radii in the disk, and ignored much of the
phenomenology we discuss here. We find that – without introducing any new
“components” or parameters to the model – this naturally predicts both the
existence and properties (including densities, temperatures, scale
heights/covering factors, reprocessed light fractions, ionization parameters,
and clumping/density variations) of many different AGN ecosystem components
above including the dusty “torus”; the broad-line region; a cooler, more
spatially extended scattering effective central emitter (compared to SS73);
and the warm Comptonizing and hard coronae. The gravitational, thermal, and
buoyant stability of this gas, as well as its location/height above the
midplane, are naturally explained by the disk itself – in brief, they are
“held up” by magnetic fields. Crucially, these flux-frozen disks are (1)
vastly geometrically thicker and more strongly-flared at large radii, (2) much
lower density and lower mass/surface density, and (3) much more strongly
turbulent (in terms of the sonic Mach number $\mathcal{M}_{s}\equiv v_{\rm
turb}/c_{s}$) than thermal-pressure-dominated disks. These completely change
the predictions for their thermal structure.

In § 2, we outline some fundamental assumptions and scalings we will use
throughout (common variables are defined in Table 1 for reference). In
subsequent sections, we divide the disk into “zones” where different gas
phases or sources of illumination dominate, outlined in Fig. 1, which we
discuss in subsequent sections, including the galactic ISM (§ 3), dust torus
(§ 4), broad-line region (§ 5), neutral (§ 6) and multi-phase (§ 7) optically-
thin disks, central thermalized and blackbody-emitting disk (§ 8), corona (§
9), and extended scattering layers and surfaces (§ 10). We describe how
behavior should change at both super-critical and highly subcritical accretion
rates (§ 11), and briefly comment on outflows and jets there and in § 12. We
discuss the key differences from thermal or radiation-pressure dominated disks
in § 13, and conclude in § 14.

Figure 1: Heuristic illustration of the accretion disk properties for a
magnetically-dominated, flux-frozen disk, with distinct “zones” with different
thermo-chemical properties. Black line shows the disk scale-height $z=H$
versus cylindrical radius $R$. We label critical radii including: the
gravitational radius/horizon $R_{g}$, ISCO $R_{\rm ISCO}$, radius interior to
which the midplane is thermalized and self-ionized $R_{\rm therm,\,ion}$,
radius at which the atomic disk is thermalized $R_{\rm therm,\,atomic}$,
radius of dust sublimation in the shielded midplane $R_{{\rm sub},\,T_{\rm
eff}}$, radius interior to which there is a Thompson-thick electron-scattering
intercepting surface illumination $R_{\rm es,\,ion}$, radius interior to which
the surface/illuminated layers have dust sublimated $R_{\rm sub,\,ext}$, and
BH radius of influence $R_{\rm BHROI}$. We divide the system into distinct
“zones” labeled (§ 3-10), described in Fig. 2.
Figure 2: Illustration as Fig. 1 of a flux-frozen disk, with descriptions of
each zone. Top right: Zone label (§ 3-10; as Fig. 1). Black line shows disk
scale-height $z=H$, versus cylindrical radius $R$. Bottom right: Phase of the
gas in each zone (e.g. “dusty” outside of the sublimation radii, or atomic, or
ionized). Bottom left: Whether each phase is effectively optically-thin or
thick to absorption, and multiphase or thermalized. Top left: Phenomenological
structures corresponding to each zone. Zones include the galactic ISM (1; §
3), dust torus (2; § 4), broad-line region (3; § 5), neutral (4; § 6) and
multi-phase (5; § 7) optically-thin disks, central thermalized and blackbody-
emitting disk (6; § 8), corona (7; § 9), and extended scattering layers and
surfaces (8; § 10).  Figure 3: Predicted disk surface densities $\Sigma_{\rm
gas}$ versus cylindrical radius $R$ (in units of gravitational radius $R_{g}$
on the bottom axis, or light-days on the top axis). Each panel shows a single
disk example (value of BH mass $M_{\rm BH}$ and accretion rate $\dot{m}$),
plotted from the ISCO ($R_{\rm ISCO}$) to the BHROI ($R_{\rm BHROI}$R), for
the flux-frozen magnetized ($\beta\ll 1$) models here (black solid line). We
contrast the prediction for the standard SS73 thermal-pressure-dominated
$\alpha$-disk model (pink dotted line). We compare different parameters: ${\rm
M}_{\rm BH}=10^{7}\,{\rm M}_{\odot}$, $\dot{m}=1$ (top left); ${\rm M}_{\rm
BH}=10^{9}\,{\rm M}_{\odot}$, $\dot{m}=1$ (top right); ${\rm M}_{\rm
BH}=10^{7}\,{\rm M}_{\odot}$, $\dot{m}=0.01$ (bottom left); ${\rm M}_{\rm
BH}=10^{7}\,{\rm M}_{\odot}$, $\dot{m}=100$ (bottom right). Shaded bars along
the horizontal axis represent the locations of different zones (colored as
Figs. 1-2) for each case, with their name (bottom right) and dividing radii
(bottom left) labeled. Their $y$-axis values are not meaningful (they only
highlight the range of radii for each zone). Bottom row shows the midplane
($|z|\ll H$) zones: thermalized disk (6a), ionized (5) and neutral (4)
optically-thin disks, and obscuring torus (2; with 2b the thermalized-dust
subregion). The row above corresponds to the disk illuminated surface
($|z|\sim H$) zones: warm comptonizing skin (6b), the scattering/reprocessing
(8) and optically-thin ionized illuminated disk (3) BLR-like region; and
illuminated warm dust-reprocessing torus (2c). Above this we show the range of
radii of the coronal gas ($|z|\gtrsim H$), zone (7). Galactic ISM (Zone (1))
resides to the right of the plot ($R>R_{\rm BHROI}$) and is not modeled here.
At most radii, masses, and accretion rates, the proportionally much stronger
Maxwell stresses in the flux-frozen disks translate to lower $\Sigma_{\rm
gas}$ compared to a thermal-pressure-dominated disk. Note the declining
central densities in SS73 at high $\dot{m}$ are a consequence of the inner
disks becoming radiation-pressure dominated. This only influences the models
here weakly, as there is a small change in slope in the $\dot{m}=100$ flux-
frozen case at $x_{g}\equiv R/R_{g}\sim 1500$ interior to which the disk
becomes radiatively inefficient (§ 11.1).  Figure 4: Predicted midplane 3D gas
densities $n\equiv\rho/m_{p}$, as Fig. 3. We compare the same BH mass and
accretion rates (panels), range of radii (ISCO to BHROI), zone locations
(shaded horizontal bars), and flux-frozen disk (thick black) versus thermal-
pressure-dominated (SS73; thin pink) prediction. SS73-like models assume disks
are weakly turbulent/laminar, so there is only one midplane density to plot.
Flux-frozen disks are supersonically turbulent and multi-phase, so we plot our
estimate of the volume-weighted mean midplane density $\langle\rho\rangle_{\rm
vol}\sim\Sigma_{\rm gas}/(2\,H)$,; mass-weighted mean midplane density
$\langle\rho\rangle_{\rm mass}\equiv M^{-1}\int\rho
dm\approx\langle\rho\rangle_{\rm vol}C^{1/2}$ where $C$ is the “clumping
factor”; and mass-weighted $\pm 2\sigma$ range of densities (grey shaded). We
assume a lognormal density distribution with the standard variance-Mach number
relation $S\approx\ln[1+(\mathcal{M}_{s}/3)^{2}]$ for supersonic turbulence to
compute these (see Konstandin et al., 2012). The discontinuities in
$\langle\rho\rangle_{\rm mass}$ appear at midplane zone transitions because we
use the analytic approximations from § 3-10 for the temperature/phase
structure and dominant opacities, hence sound speed and $\mathcal{M}_{s}$. The
mean profile crudely follows $\rho\propto R^{-2}$, with a broad range of
densities (i.e. “clumpy”/inhomogeneous structure) at all radii. In general the
densities are orders-of-magnitude lower than in an SS73 disk (owing both to
lower $\Sigma_{\rm gas}$, Fig. 3, and thicker disks with $H/R\sim 0.1-1$),
except for interior regions at large $\dot{m}$ where radiation pressure
modifies the magnetized solutions much more weakly. Note the similarity of the
densities for flux-frozen disks to observational inferences for the BLR (§ 5)
and dusty torus (§ 4) at corresponding radii (where the SS73 model is many
orders-of-magnitude more dense).  Figure 5: Total enclosed disk gas mass
$M_{\rm gas,\,disk}(<R)\equiv\int_{0}^{R}2\pi\Sigma_{\rm gas}\,R\,dR$ versus
radii $R$ for the flux-frozen models (thick black) or thermal-pressure-
dominated models (thin pink) with different BH masses and accretion rates
(linestyles labeled). Per Fig. 3, the flux-frozen disks predict much lower-
disk masses. The gas masses required for thermal pressure-dominated disks
become enormous at large radii, which in turn raises many theoretical
challenges (e.g. Goodman, 2003), and even comparable to the BH mass (so
Keplerian approximations break down, the disk cannot stabilize, and orbits
become strongly perturbed) at $\gtrsim 1000\,R_{g}$. In contrast the flux-
frozen disks require far less disk mass for the same BH mass and accretion
rate/luminosity.

##  2 Fundamental Model Assumptions & Scalings

###  2.1 Assumptions

In Paper III, we present a simple analytic model for the disk properties from
Paper II. This model makes two particularly important ansatz:

  1. 1.

Magnetic Pressure Dominates over Thermal Pressure ($\beta\sim P_{\rm
thermal}/P_{\rm mag}\ll 1$) in the Disk, with in-plane (toroidal, radial, and
turbulent) field amplification driven by flux-freezing. So $|{\bf B}|$ obeys a
relation like $|{\bf B}|^{2}\propto\rho^{\gamma}$ with $\gamma\sim 1-5/3$ or
$|{\bf B}|\propto H^{-1}$ or other similar scalings, depending on the assumed
geometry (but the detailed choice has little effect on our results). This is
the most important difference from classical $\alpha$-disk models which assume
$\beta\gg 1$ ($P_{\rm mag}\ll P_{\rm thermal}$).

  2. 2.

The turbulence is trans-Alfvénic ($\mathcal{M}_{A}\equiv v_{\rm
turb}/v_{A}\sim 1$), or more generally speaking the in-plane total stress
driving angular momentum transport is comparable to the Maxwell stress. We
note that this assumption (but not $\beta\ll 1$) is also made in many thermal-
pressure dominated $\alpha$ disk models, and is implicit in e.g. any model
where the MRI or the mean flux-frozen in-plane (toroidal-radial) field
regulates the stresses (e.g. it follows naturally from physics like those in
Balbus & Hawley 1998). The key difference is that, because $\beta\ll 1$
($v_{A}\gg c_{s}$), this means the turbulence should be highly supersonic
($\mathcal{M}_{s}\equiv v_{\rm turb}/c_{s}\sim\mathcal{M}_{A}\beta^{-1/2}\gg
1$), while for classical $\alpha$-disk models it would imply sub-sonic
turbulence.

###  2.2 Basic Disk Properties

Taking this together with the fact that the gravitational energy flux of
accreting material is $F_{\rm
grav}\approx\varpi_{g}\,(3/4\pi)\,\dot{M}\,\Omega^{2}$, Paper III showed that
this admits a simple similarity solution valid over all radii where the
potential is approximately Keplerian (from outside the ISCO out to the BHROI).
Taking the default ansatz (1) and (2) as therein, the resulting scaling for
the turbulent velocity $v_{t}$, Alfvén speed $v_{A}$, scale height $H$, and
gas surface density $\Sigma_{\rm gas}$ becomes:

$\displaystyle\frac{v_{t}}{v_{\rm K}}$ $\displaystyle\sim\frac{v_{A}}{v_{\rm
K}}\sim\frac{H}{R}\sim\left(\frac{R}{r_{\rm ff}}\right)^{1/6}\sim
0.07\,\left(\frac{m_{7}\,x_{g}}{r_{{\rm ff},5}}\right)^{1/6},$ (1)
$\displaystyle\frac{\Sigma_{\rm gas}}{{\rm g\,cm^{-2}}}$ $\displaystyle\sim
0.014\,\frac{\dot{m}\,m_{7}^{1/2}}{r_{{\rm
ff},5}^{1/2}}\,\left(\frac{R}{r_{\rm ff}}\right)^{-5/6}\sim
10^{4}\frac{\dot{m}\,r_{{\rm ff},5}^{1/3}}{m_{7}^{1/3}\,x_{g}^{5/6}}$ (2)

in terms of the cylindrical radius $R$ (also defined in terms of the usual
gravitational radius $x_{g}\equiv R/R_{g}$ with $R_{g}\equiv 2\,G\,M_{\rm
BH}/c^{2}$), Keplerian speed $v_{\rm K}^{2}\equiv G\,M_{\rm BH}/R$, Eddington-
scaled accretion rate $\dot{m}\equiv\dot{M}/\dot{M}_{\rm crit}$ (with
$\dot{M}_{\rm crit}\equiv M_{\rm BH}/5\times 10^{7}\,{\rm yr}=L_{\rm
Edd}/(0.1\,c^{2})$ as defined here, for a reference radiative efficiency
$\epsilon_{r,\,0.1}=0.1$), BH mass $m_{7}\equiv M_{\rm BH}/10^{7}\,{\rm
M_{\odot}}$, and “free-fall” radius $r_{\rm ff}\equiv r_{\rm ff,\,5}\,5\,{\rm
pc}$ defined as the radius where the solutions map onto free-fall into the BH
potential from larger radii outside of the Keplerian potential. For reasonable
models of accretion onto a SMBH from galactic scales via gravitational capture
(e.g. Bondi-Hoyle accretion, or tidal disruption of ISM clouds or star
clusters, or loss-cone capture, or gravitational torques from nested bars or
mergers; Shlosman et al. 1989; Jogee 2006; Hopkins & Hernquist 2006; Hopkins &
Quataert 2010b; Anglés-Alcázar et al. 2021) $r_{\rm ff}$ will correspond to
the BHROI,

$\displaystyle r_{\rm ff}\sim R_{\rm BHROI}\sim G\,M_{\rm BH}/\sigma_{\rm
galaxy}^{2}\approx 5\,{\rm pc}\,m_{7}^{1/2}$ (3)

in terms of the galactic nuclear velocity dispersion $\sigma_{\rm galaxy}$
(where the latter scaling inserts the observed $M_{\rm BH}$-$\sigma_{\rm
galaxy}$ relation; McConnell & Ma 2013). This scale, the gravitational
radius/ISCO, and other key scales, along with the (weak) scaling of $H/R$ with
radius, are illustrated in Figs. 1-2. The predicted surface density scalings
are shown in Fig. 3, with 3D midplane densities in Fig. 4 and disk mass in
Fig. 5.

From these, we can immediately derive other relevant scalings, e.g. the
volume-averaged midplane density $\rho\approx\Sigma_{\rm gas}/2\,H$ or
effective Toomre $Q$ (including magnetic+turbulent support):

$\displaystyle\frac{\rho}{\rm m_{p}\,cm^{-3}}$ $\displaystyle\sim
272\,\frac{m_{7}^{1/2}\,\dot{m}}{r_{\rm ff,\,5}^{3/2}}\,\left(\frac{R}{r_{\rm
ff}}\right)^{-2}\sim 7\times 10^{15}\frac{\dot{m}\,r_{\rm
ff,\,5}^{1/2}}{m_{7}^{3/2}\,x_{g}^{2}}$ (4) $\displaystyle Q_{\rm tot}$
$\displaystyle\sim 3000\,\frac{m_{7}^{1/2}}{r_{\rm
ff,\,5}^{3/2}\,\dot{m}}\,\left(\frac{R}{r_{\rm ff}}\right)^{-1}\gg 1$ (5)

These can be compared to the values in Fig. 4 (where we also show the mass-
weighted densities and range of densities, which depends on the turbulent Mach
numbers and therefore indirectly on the thermal state, calculated below). Both
$Q\gg 1$ as well as our assumption of an approximately Keplerian potential
follow from these densities and $M_{\rm disk}\ll M_{\rm BH}$ as shown in Fig.
5.

###  2.3 General Considerations

Particularly convenient for us, precisely because these disks are
magnetically-dominated, it means that unlike standard thermal-pressure
dominated disk models, changing the disk thermal or opacity structure has no
effect on the accretion rates or basic dynamical properties of the disk
derived in Paper III. So do not need to re-derive a self-consistent “global”
disk model, but instead solve for the thermal structure given a fixed
“background” disk dynamical structure. For example, the scale heights in Fig.
1, or surface and 3D densities in Figs. 3 & 4, are defined by Eqs. 1-2 above,
independent of the disk thermal/opacity properties.

For calculating irradiation effects, we will define the bolometric accretion
luminosity, dominated by the emission from the inner disk regions, as

$\displaystyle L_{\rm bol}\equiv\epsilon_{r}\,\dot{M}\,c^{2}\approx 1.1\times
10^{45}\,{\rm erg\,s^{-1}}\,m_{7}\,\dot{m}\,\epsilon_{r,\,0.1}$ (6)

for some arbitrary radiative efficiency.

Paper II and Paper III briefly considered upper limits to the disk temperature
to show that these assumptions are self-consistent, i.e. the disk maintains
$\beta\ll 1$ at all radii. This also means, with ansatz (2), that the
turbulence is supersonic at all radii

$\displaystyle\mathcal{M}_{s}\sim\frac{v_{\rm
turb}}{c_{s}}\sim\frac{v_{A}}{c_{s}}\sim\beta^{-1/2}\gg 1,$ (7)

and that the ratio of cooling time $t_{\rm cool}$ to dynamical time $t_{\rm
dyn}\sim 1/\Omega$ is $\ll 1$ (for cooling flux balancing $F_{\rm grav}$,
$t_{\rm cool}/t_{\rm dyn}\sim\beta\sim 1/\mathcal{M}_{s}^{2}\ll 1$). Other
implicit assumptions we will make here, like that ideal MHD and a collisional
fluid description is a reasonable approximation (the ionization fractions are
high, and the atomic/molecular/ion/electron mean free paths are all very small
compared to global length scales)333For example, the ion collision time in the
ionized disk is approximately $\sim 10^{-9}\Omega^{-1}(T/T_{\rm
eff})^{3/2}\,(n_{\rm midplane}/n)\,(\dot{m}x_{g})^{-5/8}\ll 1$. are easily
validated and checked explicitly in Paper I. We confirm this ourselves in e.g.
Fig. 6. In § 11.1, we also explicitly confirm that we can neglect radiative
viscosity and damping effects.

Implicit in this is that the gravitational flux/energy change (as gas flows
in) appears in kinetic (turbulent) and magnetic energy, via simple
infall+compression/flux-freezing and mixing via various instabilities and
nonlinear processes. Since this is supersonic and trans-Alfvénic the driving
scale is naturally $\mathcal{O}(H)$ and turnover time of the largest eddies
(containing most of the power) is $\sim\Omega^{-1}$. The heat transfer to gas
is therefore mediated by shocks and (turbulent) reconnection, before energy is
radiated away.

###  2.4 Self-Illumination

We will throughout consider cases where the illumination of the surface of the
outer disk by the central inner disk (which dominates the bolometric
luminosity) may be important. Assuming $L_{\rm bol}$ comes from the inner
disk, then the flux per unit cylindrical area for a flared disk scales as

$\displaystyle F_{\rm
illum}\sim\left(\frac{\partial\ln{H}}{\partial\ln{R}}-1\right)\left(\frac{H}{R}\right)\,F_{\rm
incident}\approx\frac{1}{6}\frac{H}{R}\,F_{\rm incident}$ (8)

(Eq. 1), where

$\displaystyle F_{\rm incident}\equiv\frac{L_{\rm
bol}\,f_{\theta}}{4\pi\,r^{2}}$ (9)

depends on $f_{\theta}$ which defines the illumination pattern as a function
of polar angle. For a flat, geometrically thin ($H/R\ll 1$), effectively
optically-thick ($\tau^{\ast}\gg 1$), and aligned (coplanar with the outer
disk) inner disk,

$\displaystyle f_{\theta}\approx
2\,\cos{\theta}=\frac{2\,(H/R)}{\sqrt{1+(H/R)^{2}}}.$ (10)

For an isotropic, or spherical/point-source-like, or optically thin central
source,

$\displaystyle f_{\theta}=1\ .$ (11)

In practice, we show below that the disk models here should be “in between”
these extremes: the inner disk is effectively optically-thick but not by a
large margin, and the extended inner-disk scattering layer is not; the inner
disk has $H/R<1$ but again not by a large margin, and is flared; and the inner
disk will often be mis-aligned with the outer disk (e.g. in the simulations,
the disk has warps of $\sim 40^{\circ}$ at multiple radii owing to accretion
of material with different angular momentum, and is typically mis-aligned with
the BH spin which will realign the inner disk; see Paper II). In these cases
an exact solution requires detailed radiative transfer calculations, so we
instead simply consider both extremes to bracket the range of possibilities
analytically.

Figs. 1-2 show, as we discuss below, the importance of self-illumination.
After calculating the relevant opacity effects in Fig. 7 (and pressure
profiles in Fig. 8), we illustrate this in more detail in Fig. 9.

###  2.5 Comparison to Thermal-Pressure Dominated (SS73) Disks

We will refer to SS73 disks throughout for comparison (e.g. Figs. 3, 4, 5, 6,
7, 8, 10), as a reference model for the category of thermal-pressure dominated
disks. We take their scalings directly (converting to our definitions of
$\dot{m}$, etc.), assuming the Maxwell stress is comparable to the total
stress (determined by the $\alpha$ parameter in their model which we will take
by default to be $\alpha\sim 0.1$). While details in some models differ, for a
radiatively efficient, thermal-pressure-dominated disk, the scalings in SS73
are all qualitatively similar (for our purposes) to other thermal-pressure-
dominated model variants in the literature (see Abramowicz & Fragile 2013),
and we review below how the most general conclusions from this comparison are
robust to other detailed assumptions once one assumes a midplane $\beta\gg 1$.

In this sense, one should note that in the scalings we presented above, there
are a couple of key differences in the structural properties of magnetically-
dominated disks compared to thermal-pressure-dominated disks, which will
appear repeatedly in our analysis. The magnetically-dominated disks are
flared, and extremely geometrically thick (with $H/R\sim 1$ in the outer
regions) compared to a thermal pressure-dominated disk (for which $H/R\ll
0.01$ in the outer regions). For the same accretion rate, magnetically-
dominated disks have both much lower 3D density $\rho$ (by factors up to $\sim
10^{7}$; Fig. 4) and much lower surface density $\Sigma_{\rm gas}$ (by a
factor $\sim 10^{5}$; Fig. 3) owing to the much stronger Maxwell stresses (for
a detailed comparison of these basic structural properties, see Paper III). As
a result, we will show that the thermal properties of the different disk
models are qualitatively distinct (e.g. Figs. 6, 7, 8).

##  3 Zone 1: The Galactic Nucleus/ISM/Starburst Environment

We now consider different regions or “zones” around the BH, roughly proceeding
“outside-in” with the accretion flow (Fig. 1). We divide the zones by their
qualitative thermochemical and radiative behaviors (Fig. 2), which are
determined by the dominant gas phases and heating/cooling processes (Figs. 3,
4, 6, 7, 8). Fig. 1 shows the definitions of these zones; Fig. 2 illustrates
the key properties of the dominant gas phases in each zone and their
connections to different classic structures known around AGN.

Beginning from the largest radii, therefore, at $R\gtrsim r_{\rm ff}\sim
R_{\rm BHROI}$, we are fully outside the regime where the BH dominates the
gravitational potential, so cannot in any meaningful way talk about being in
an “accretion disk.” Of course gas may lose angular momentum and flow into
this general region from the outer galaxy or intergalactic medium (IGM), but
this is the galactic ISM in the galaxy nucleus/bulge center, and such flows
are governed by galactic dynamics. This should be multi-phase, dusty, clumpy,
and can have a covering factor for optical obscuration of the AGN ranging from
nil (in e.g. “passive” elliptical host galaxies) to $\sim 100\%$ (in e.g.
starburst nuclei with isotropic, Compton-thick nuclear dust columns), as
discussed in e.g. Simcoe et al. (1997); Hopkins et al. (2005a, b, c, d);
Hopkins et al. (2006); Ghosh et al. (2007); Trump et al. (2009); Trump (2011);
Hopkins et al. (2012a); Gilli et al. (2022); Glikman et al. (2024). So
certainly part of the clumpy torus and some re-emission/reprocessing can arise
from these radii: e.g. dust directly exposed to an unobscured QSO sightline
can reach dust temperatures

$\displaystyle T_{\rm dust}\sim 200\,{\rm
K}\left(\frac{m_{7}\,\dot{m}\,\epsilon_{r,\,0.1}}{r_{\rm
ff,\,5}^{2}}\right)^{1/4}\left(\frac{r_{\rm ff}}{R}\right)^{1/2},$ (12)

modestly higher than the median $\sim 70-100\,$K expected from illumination
just by stars in a starburst environment (Sanders et al., 1988; Narayanan et
al., 2005; Veilleux et al., 2009; Younger et al., 2009b; Pope et al., 2008;
Greve et al., 2009; Younger et al., 2009a; Narayanan et al., 2011; Hayward et
al., 2011). And of course, the NLR will arise from gas on these scales (Stark
& Carlson, 1984; Bennert et al., 2002; Rice et al., 2006; Meena et al., 2022).
But modeling these scales in more detail depends on the galactic environment,
star formation properties, morphology, etc., and most of the gas properties
will depend on the detailed physics of star formation and stellar feedback
self-regulating the galactic ISM (see Paper I and Hopkins & Quataert 2010b;
Hopkins et al. 2016; Anglés-Alcázar et al. 2021; Hopkins et al. 2022a; Wellons
et al. 2023; Byrne et al. 2023a, b; Mercedes-Feliz et al. 2023; Cochrane et
al. 2023; Mercedes-Feliz et al. 2024), and so is outside the scope of our
focus here.

##  4 Zone 2: The Cold, Dusty Outer Disk/Torus

The first region where our model applies is just interior to $r_{\rm ff}\sim
R_{\rm BHROI}$ (Figs. 1-2). We expect, just like in the ISM, the medium to be
dusty (given the dust temperature scaling above, we are well exterior to the
sublimation radius $R_{\rm sub}$). For dust one can approximate the ratio of
absorption to total opacities as $\kappa_{a}/(\kappa_{a}+\kappa_{s})\sim
1-0.5\,(1+10^{5.72}\,(1+T_{\rm rad}^{2})^{-1})^{-1}$ (compare Weingartner &
Draine, 2001; Semenov et al., 2003; Draine, 2011) for a given radiation
temperature $T_{\rm rad}$, which is $\kappa_{a}/(\kappa_{a}+\kappa_{s})\approx
0.5$ for $T_{\rm rad}\gtrsim 1000\,$K and
$\kappa_{a}/(\kappa_{a}+\kappa_{s})\approx 1$ for $T_{\rm rad}\lesssim 1000$
K. So at the cold temperatures of interest in the outermost disk, the dust
opacity will be primarily absorption (and in the warmest parts the scattering
opacity is at most comparable to absorption). This means the “effective”
optical depth from the midplane to infinity,
$\tau^{\ast}\equiv\int_{0}^{\infty}\sqrt{\kappa_{a}\,(\kappa_{a}+\kappa_{s})}\,\rho(R,\,z)\,dz\approx\langle\sqrt{\kappa_{a}\,(\kappa_{a}+\kappa_{s})}\rangle\,\Sigma_{\rm
gas}/2$ (where the surface densities follow the fiducial scalings for these
disks from Eq. 2, illustrated in Fig. 3) and the “total” optical depth
$\tau\equiv\langle\kappa_{a}+\kappa_{s}\rangle\,\Sigma_{\rm gas}/2$ are
roughly equal to the absorption optical depth $\tau_{\rm
abs}=\kappa_{a}\,\Sigma_{\rm gas}/2$. The absorption opacity $\kappa_{a}$ can
be approximated (at the level of approximation here, ignoring detailed
substructure within the disk and chemical inhomogenity, etc.) as ${\rm
MIN}[5\,,3\,(T_{\rm rad}/100\,{\rm K})^{\beta}]\,\tilde{Z}\,{\rm
cm^{2}\,g^{-1}}$ with $\beta\approx 1.5$ and $\tilde{Z}\equiv
Z/Z_{\odot}$,444Our scalings can be extrapolated to any $\tilde{Z}$, but
motivated by modern numerical simulations and more detailed QSO full-spectrum
line modeling (e.g. Temple et al., 2021b), we will often refer to a reference
$\tilde{Z}\sim 1$, as opposed to the highly super-Solar $\tilde{Z}\sim 10-50$
sometimes invoked in older quasar models (these highly super-solar values are
based on more simplified older line models which assumed all observed lines
come from a single density and ionization parameter and source distance).
until $T_{\rm dust}\gtrsim 1500-2000\,$K where the dust sublimates
($\tilde{Z}\rightarrow 0$) (Semenov et al., 2003; Draine, 2011). These
tabulations for the opacities are used for the calculation of the opacities in
this zone shown in Fig. 7. Given that the free electron fractions in this
region are $\lesssim 0.01$ (Paper I), and metallicities are $\tilde{Z}\sim 1$,
this will dominate over other sources of opacity. So we have:

$\displaystyle\tau^{\ast}$ $\displaystyle\approx{\rm
MIN}[0.021T_{100}^{2},\,0.035]\frac{m_{7}^{1/2}\,\dot{m}\,\tilde{Z}}{r_{\rm
ff,\,5}^{1/2}}\left(\frac{R}{r_{\rm ff}}\right)^{-5/6}$ (13)

where $T_{100}\equiv T_{\rm rad}/100\,{\rm K}$. Fig. 1 shows the range of
scales (calculated below) over which this “zone” applies; Fig. 2 illustrates
the characteristic thermochemical properties qualitatively, and breaks the
zone up into sub-zones which we describe below. The quantitative values of the
temperature, opacity, and pressure properties, given the approximations in
this section, are shown for various values of the BH mass and accretion rate
in Figs. 6, 7, 8, respectively.

Note the covering factor reaches $\sim 0.7$ for the very outer edge of this
disk as $R\rightarrow r_{\rm ff}$, but only some of this is able to be
illuminated by the central disk (illustrated in Fig. 9). At the (order-unity)
values of $H/R$ predicted for the radii of interest here, it makes a
negligible difference for our conclusions whether we assume the central source
is an optically-thick flat aligned disk or assume quasi-isotropic illumination
(i.e. either form of the function $f_{\theta}$ in § 2.4 we assume gives nearly
identical results).

###  4.1 2a: The Outer, Optically-Thin Cold Dust Region

The outer regions will therefore be optically-thin to the cooling radiation of
the dust. We show its properties as Zone (2a) in Fig. 2, and quantitatively
see the effects of optically-thin cooling in e.g. the temperature calculation
in Fig. 6 and corresponding pressures (Fig. 8). The optically-thin cooling
flux $F_{\rm cool,\,thin}=4\pi\,\kappa_{\rm em}\,\Sigma_{\rm
gas}\,\sigma_{B}\,T_{\rm dust}^{4}$ exceeds $F_{\rm grav}$ so long as $T_{\rm
dust}\gtrsim 12\,{\rm K}\,m_{7}^{3/8}\,r_{\rm ff,\,5}^{-5/8}\,(\kappa_{\rm
em}/5\,{\rm cm^{2}\,g^{-1}})^{-1/4}\,(R/r_{\rm ff})^{-13/24}$, which is easily
maintained. In practice the midplane regions of this zone, while shielded from
the central source, will be externally heated by the massive star irradiation
and cosmic ray production in Zone (1), and maintain typical nuclear dust and
gas temperatures in equilibrium with one another at these densities (but not
necessarily with $T_{\rm rad}$) at

$\displaystyle T_{\rm dust}\sim T_{\rm gas}\sim 20-100\,K.$ (14)

The upper parts of the disk near $|z|\sim H$ will be directly illuminated by
the central disk, and maintain warmer dust temperatures as discussed below:

$\displaystyle T_{\rm gas}\sim T_{\rm dust}\sim 200\,{\rm
K}\left(\frac{m_{7}\,\dot{m}\,\epsilon_{r,\,0.1}}{r_{\rm
ff,\,5}^{2}}\right)^{1/4}\left(\frac{r_{\rm ff}}{R}\right)^{1/2},$ (15)

until sublimation (noted below).555If the irradiated layer can efficiently
transfer heat to the central layers below it and reach equilibrium, then the
optically-thin cooling balances heating for $T_{\rm dust}\sim 200\,{\rm
K}\,m_{7}^{1/8}\,r_{\rm ff,\,5}^{-3/8}\,\epsilon_{r,\,0.1}^{1/4}\,(\kappa_{\rm
em}/5\,{\rm cm^{2}\,g^{-1}})^{-1/4}\,(R/r_{\rm ff})^{-5/24}$, similar to the
directly heated value but much more weakly dependent on radius and other
parameters. It is not obvious how this would occur, but we note it for
completeness.

In most ways, this is akin to Zone (1) and “normal” GMCs at least
qualitatively (through quantitatively it is “more extreme” in most respects):
it is optically-thin to its own cooling radiation, the dust and gas are both
“cool” and (given the low densities) may not be in exact equilibrium. The
plasma $\beta\ll 1$ and the turbulence is super-sonic

$\displaystyle\mathcal{M}_{s}\equiv\frac{v_{\rm
turb}}{c_{s}}\sim\beta^{-1/2}\sim 100\,m_{7}^{1/2}r_{\rm
ff,\,5}^{-1/2}T_{100}^{-1/2}\left(\frac{r_{\rm ff}}{R}\right)^{1/3}$ (16)

(quantitative examples of $\mathcal{M}_{s}$ are shown in Fig. 10, for
different BH masses and accretion rates), with cooling times much shorter than
dynamical times $t_{\rm cool}\ll t_{\rm dyn}\sim 1/\Omega$. Just like in GMCs
the gas should therefore be multiphase (with the phases still neutral but akin
to a warm and cold atomic and cold molecular medium) with large density
contrasts owing to the supersonic motions, with the largest density contrasts
reaching $\sim\mathcal{M}_{s}^{2}\sim 10^{4}$ (for shocks along the field
lines, as given the rapid cooling these are radiative/isothermal shocks; see
Fig. 10). The characteristic size of such extreme overdensities will reach as
small as the sonic length or post-shock length in supersonic turbulence,

$\displaystyle R_{\rm sonic}\sim\frac{H}{\mathcal{M}_{s}^{2}}\sim
10^{15}\,{\rm cm}\,r_{\rm ff,\,5}^{2}\,m_{7}^{-1}\,T_{100}\,(R/r_{\rm
ff})^{11/6}.$ (17)

The variation of that size scale with radius (again for different accretion
rates and BH masses) is shown in Fig. 11. Note that (like in GMCs, but even
moreso here) it is not correct to assume that different phases or clumps are
in thermal pressure equilibrium, since thermal pressure is vastly subdominant
to both magnetic and turbulent ram pressure.

While some very small amount of star formation can (and does, in simulations;
see Hopkins et al. 2024a) persist at these radii, it is strongly suppressed
(compared to e.g. the expectations for classical $\alpha$-disks discussed in
Goodman 2003), as the Toomre $Q\gg 1$ (Eq. 5) and similarly the magnetic
critical mass is formally larger than the entire disk mass (so collapse can
only occur in very special regions, e.g. local magnetic field line polarity
switches with strong shocks along the field lines). But unlike in the ISM,
young/massive stars are not needed to power either the turbulence or
“maintain” the large scale heights here. The magnetic field (which can grow to
these values purely via flux-freezing) and turbulence powered by the fields
and gravitational flux, before being dissipated, supports the material
vertically and prevents denser gas from “sinking” or sedimenting to the
midplane (by definition in the models here, but easily verified from the
scalings of Eq. 1-5).

The covering factors of this zone

$\displaystyle f_{\rm
cover}\approx\cos{\left(\tan^{-1}\frac{R}{H}\right)}=\frac{H/R}{\sqrt{1+(H/R)^{2}}}\sim
0.5-0.7$ (18)

extend from $f_{\rm cover}\sim 0.5\,(m_{7}/r_{\rm
ff,\,5})^{1/10}\,(\dot{m}\,\tilde{Z})^{1/5}$ at its inner boundary to $f_{\rm
cover}\sim 0.7$ at its outer boundary (see Fig. 9). This is therefore
important for re-radiation/reprocessing (discussed below) and obscuration, but
the intrinsic gravitational cooling luminosity (i.e. distinct from the
reprocessed radiation) coming from this region is small, $\sim 10^{41}\,{\rm
erg\,s^{-1}}\,m_{7}^{7/5}\,r_{\rm
ff,\,5}^{-2/5}\,\tilde{Z}^{-6/5}\,\dot{m}^{-1/5}$ (note this can decrease with
$\dot{m}$ because of how the location of this region varies with $\dot{m}$), a
tiny fraction of the bolometric luminosity for the parameter space of
interest.

###  4.2 2b: The Inner, Thermalized Warm Dust Region

At radii interior to

$\displaystyle R\lesssim r_{\rm dust,\,therm}\approx 0.1\,{\rm
pc}\,m_{7}^{3/5}\,r_{\rm ff,\,5}^{2/5}\,(\dot{m}\,\tilde{Z})^{6/5},$ (19)

$\tau^{\ast}$ will become $\gtrsim 1$, and the dust will thermalize and
radiate as a blackbody (zone (2b) in Fig. 2). The gas densities here are

$\displaystyle n_{\rm gas}\sim 10^{6}\,\dot{m}^{-7/5}\,m_{7}^{-7/10}\,r_{\rm
ff,\,5}^{-3/10}\,\tilde{Z}^{-12/5}\,(R/r_{\rm dust,\,therm})^{-2},$ (20)

so dust and gas are efficiently collisionally coupled, and now $T_{\rm
gas}\approx T_{\rm dust}\approx T_{\rm rad}$. Equating the blackbody cooling
flux $F_{\rm cool,\,thick}=2\,\sigma_{B}\,T_{\rm eff}^{4}$ to $F_{\rm grav}$
gives an effective temperature $T_{\rm eff}\approx 200\,{\rm
K}\,m_{7}^{1/20}\,\dot{m}^{-13/20}\,r_{\rm
ff,\,5}^{-3/10}\,\tilde{Z}^{-9/10}\,(R/r_{\rm dust,\,therm})^{-3/4}$, so the
dust is at hundreds of ${\rm K}$ (justifying the opacity used to calculate
where $\tau^{\ast}>1$). But here, the re-radiation from the illuminated layer
above should indeed be reprocessed within the disk, so it is more appropriate
to equate to $F_{\rm illum}$ rather than $F_{\rm grav}$, giving

$\displaystyle T_{\rm eff}\approx 1000\,(\epsilon_{r,\,0.1}/\dot{m}\,r_{\rm
ff,\,5}\,\tilde{Z}^{2})^{1/4}\,(R/r_{\rm dust,\,therm})^{-5/12}$ (21)

(Fig. 6). Under these conditions the radiation transport is diffusive
($\tau^{\ast}\sim\tau\gtrsim 1$, and the vertical turbulent transport speed
$v_{\rm turb}$ as well as radial/inflow advection transport $v_{r}\sim v_{\rm
turb}^{2}/v_{\rm K}$ are much smaller than the diffusive speed $v_{\rm
diff}\sim c/\tau$ here), so the dust becomes somewhat more mono-phase at a
given position and can set up a stratified plane-parallel radiation
temperature gradient with midplane temperature $T_{\rm mid}\approx T_{\rm
eff}\,\tau^{1/4}\approx T_{\rm eff}\,(R/r_{\rm dust,\,therm})^{-5/24}$, which
is a modest enhancement to the midplane over effective temperature, over the
range of interest. For this reason, in Fig. 6, we show multiple temperatures
corresponding to different regions or viable phases of the gas. This plus the
calculated densities (Fig. 4) and radiation fluxes gives us the resulting
pressures in Fig. 8. The covering factor $f_{\rm cover}$ is

$\displaystyle f_{\rm cover}\sim 0.5\,(m_{7}/r_{\rm
ff,\,5})^{1/10}\,(\dot{m}\,\tilde{Z})^{1/5}\,(R/r_{\rm dust,\,therm})^{1/6}$
(22)

(Fig. 9). But the turbulence is still highly supersonic
$\mathcal{M}_{s}\gtrsim 100$ (Fig. 10) and $t_{\rm cool}\ll t_{\rm dyn}$ so
the medium can still be clumpy/inhomogeneous down to the sonic scale

$\displaystyle R_{\rm sonic}\sim$ $\displaystyle 7\times 10^{12}\,{\rm
cm}\,m_{7}^{1/10}\,r_{\rm
ff,\,5}^{13/20}\,\dot{m}^{39/20}\,\tilde{Z}^{17/10}\,\epsilon_{r,\,0.1}^{1/4}\,\times$
$\displaystyle\ \ \ \ (R/r_{\rm dust,\,therm})^{29/24}$ (23)

(Fig. 11).

When $T_{\rm dust}$ at the midplane $T_{\rm mid}$ exceeds the sublimation
temperature $T_{\rm sub}\sim 1500\,$K, it will begin to be destroyed: this
will begin at $R_{\rm sub,\,T_{\rm mid}}/r_{\rm dust,\,therm}\sim
0.12\,m_{7}^{6/115}\,r_{{\rm
ff},\,5}^{-36/115}\,\dot{m}^{-78/115}\,\tilde{Z}^{-108/115}\,T_{\rm
sub,\,1500}^{-24/23}$ ignoring external illumination, or $\sim
0.5\,\epsilon_{r,\,0.1}^{2/5}\,\dot{m}^{-2/5}\,r_{\rm
ff,\,5}^{-2/5}\,\tilde{Z}^{-4/5}\,T_{\rm sub,\,1500}^{-8/5}$ including it (the
illuminated surface layer is discussed more below). But this will lower the
midplane opacity, making the midplane cooler, and slowing sublimation.
Sublimation will therefore not be complete until it occurs at the effective
temperature $T_{\rm eff}\gtrsim T_{\rm sub}$, which occurs at $R_{\rm
sub,\,T_{\rm eff}}/r_{\rm dust,\,therm}\sim
0.07\,m_{7}^{1/15}\,\dot{m}^{-13/15}\,r_{\rm
ff,\,5}^{-2/5}\,\tilde{Z}^{-6/5}\,T_{\rm sub,\,1500}^{-4/3}$ or

$\displaystyle R_{\rm sub,\,T_{\rm eff}}\sim 2\times 10^{16}\,{\rm
cm}\,m_{7}^{2/3}\,\dot{m}^{1/3}\,T_{\rm sub,\,1500}^{-4/3}$ (24)

ignoring external illumination, or

$\displaystyle R_{\rm sub,\,T_{\rm eff}}\sim 10^{17}\,{\rm
cm}\,\epsilon_{r,\,0.1}^{3/5}\,m_{7}^{3/5}\,\dot{m}^{3/5}\,r_{\rm
ff,\,5}^{-1/5}\,T_{\rm sub,\,1500}^{-12/5}$ (25)

including it. This inner boundary appears in Figs. 1-2.

The covering factor at the inner radius $R_{\rm sub,\,T_{\rm eff}}$ of this
zone (Fig. 9) is $f_{\rm cover}^{\rm inner}\approx
0.33\,m_{7}^{1/9}\,\dot{m}^{1/18}\,r_{\rm ff,\,5}^{-1/6}\,T_{\rm
sub,\,1500}^{-2/9}$ (ignoring illumination) or $f_{\rm cover}^{\rm
inner}\approx 0.43\,(m_{7}\,\dot{m}\,\epsilon_{r,\,0.1}/r_{\rm
ff,\,5}^{2}\,T_{\rm sub,\,1500}^{4})^{1/10}$ (including it), remarkably
independent of the BH accretion properties (Fig. 12 illustrates how the
boundaries of these zones move in radius and $H/R$ with varying $\dot{m}$),
and again $\mathcal{O}(1)$ so clearly important for re-radiation and
obscuration, but the intrinsic (neglecting reprocessing) gravitational cooling
emission luminosity is $\sim 3\times 10^{41}\,{\rm
erg\,s^{-1}}\,m_{7}^{7/5}\,\dot{m}^{2/5}\,T_{\rm
sub,\,1500}^{12/5}\,\epsilon_{r,\,0.1}^{-3/5}$, still a small fraction ($\sim
0.1\%$) of the bolometric output. The maximum dust column density for a
sightline through the midplane (integrating from the sublimation radius to
infinity) is equivalent to a gas column of

$\displaystyle N_{\rm H}\sim 3\times 10^{23}\,{\rm
cm^{-2}}\,\dot{m}^{1/5}\,m_{7}^{-1/10}\,r_{\rm
ff,\,5}^{1/2}\,\tilde{Z}^{3/5}\,T_{\rm sub,\,1500}^{8/5}$ (26)

(the usual units quoted assuming a standard dust-to-gas ratio, or $A_{v}\sim
160\,\dot{m}^{1/5}\,m_{7}^{-1/10}\,r_{\rm
ff,\,5}^{1/2}\,\tilde{Z}^{3/5}\,T_{\rm sub,\,1500}^{8/5}$). This “midplane”
section of this zone will therefore not contribute much to the direct
emission/SED, but will be important for obscuration of the central source, as
viewed from angles intercepted by its covering factor. Notably these predicted
covering factors, clumping factors, and column densities are very similar to
the canonical “obscuring torus” properties inferred from observations
(Lawrence & Elvis, 1982; Antonucci, 1982; Krolik & Begelman, 1988; Mor et al.,
2009; Hatziminaoglou et al., 2009; Hönig & Kishimoto, 2010; Alonso-Herrero et
al., 2011; Koshida et al., 2014; Hönig, 2019; Cackett et al., 2021).

###  4.3 2c: The Directly-Illuminated Dusty Region

Because the disk is flared with a covering factor $\mathcal{O}(1)$ at these
radii, above the scale height of the inner sublimation zone $\sim 0.33$ (below
this is shielded by the entire disk with Compton-thick column densities),
there will a region near the top of the disk ($|z|\sim H$) “exposed” to the
direct radiation of the central source (zone (2c) in Figs. 1, 2, 9). This
region will be optically-thin or only modestly optically-thick to electron
scattering ($\tau_{\rm es}\sim 0.1-10$ even if fully-ionized, and much smaller
in the neutral part of this region where free electron fractions range from
$\sim 10^{-8}-10^{-2}$; Fig. 7), so this can be neglected here. Given the non-
linearly large $H/R$ and scattering regions (discussed below, see Fig. 9)
around the central source it makes little difference if we take the incident
flux $F_{\rm incident}=f_{\theta}\,L_{\rm bol}/(4\pi\,r^{2})$ to be isotropic
or follow $f_{\theta}=2\,\cos{\theta}$, along rays from the center (recall
though this is distinct from the heating rate per unit area, which scales as
$F_{\rm illum}$). In either case this will sublimate dust out to a radial
distance

$\displaystyle r_{\rm sub,\,ext}$ $\displaystyle\sim(L_{\rm
bol}\,\tilde{Q}/16\pi\,\sigma_{B}\,T_{\rm sub}^{4})^{1/2}$ (27)
$\displaystyle\sim 3\times 10^{17}\,{\rm
cm}\,(\epsilon_{r,\,0.1}\,m_{7}\,\dot{m}\,\tilde{Q}_{\rm abs})^{1/2}\,T_{\rm
sub,\,1500}^{-2}.$

Here $\tilde{Q}=Q_{\rm abs}/Q_{\rm em}$ depends on the ratio of the absorption
and emission coefficients which themselves depend on the incident spectral
shape, dust composition, and temperature, but in general this is not a large
correction. This is the more commonly-quoted “sublimation radius.” Generically
this is $\gg R_{\rm sub,\,T_{\rm eff}}$ (where $H/R=(H/R)_{\rm in}\sim 0.43$
above) but $\ll r_{\rm ff,\,5}$ (generally comparable to $r_{\rm
dust,\,therm}$), with $H/R$ at $r_{\rm sub,\,ext}$ giving a covering factor if
we include the entire disk up to this height of $f_{\rm cover}^{\rm
out,\,illum}\sim(H/R)_{\rm out}\sim
0.51\,(m_{7}\,\dot{m}\,\tilde{Q}\,\epsilon_{r,\,0.1}/r_{\rm
ff,\,5}^{2}\,T_{\rm sub,\,1500})^{1/6}$.

So the covering fraction specifically of directly illuminated “hot” dust,
where $T_{\rm dust}$ approaches an appreciable fraction of the sublimation
temperature is given by

$\displaystyle f_{\rm cover}^{\rm illum}\sim f_{\rm cover}^{\rm
out,\,illum}-f_{\rm cover}^{\rm inner}\sim 0.08$ (28)

at $H/R\sim 0.4-0.5$ (Fig. 9). So (again roughly independent of $f_{\theta}$)
it will re-radiate a similar fraction of $L_{\rm bol}$ at hot dust
temperatures, and a fraction $\sim 0.2-0.5$ at cooler dust temperatures
(corresponding to the illuminated surface reaching $H/R\sim 0.7$ at the outer
$r_{\rm ff,\,5}$, and as large as $\sim 1$ depending on if there is colder
dust in zone 1 from the ISM). Again, this is consistent with what is typically
inferred from observations of the mid-IR emission spectra of quasars (see
references above and Alonso-Herrero et al. 2021; Lyu & Rieke 2022) and/or more
recent spatially-resolved imaging studies (García-Burillo et al., 2019, 2021;
Cackett et al., 2021; GRAVITY Collaboration et al., 2021; Izumi et al., 2023)
and associated with the “warm/hot dusty torus.”

The non-sublimated but illuminated “wedge” at $R\gtrsim r_{\rm sub,\,ext}$
will otherwise have properties (besides dust temperature) similar to region
(2a), per Figs. 3-11. We discuss the sublimated, externally-illuminated wedge
at $R\lesssim r_{\rm sub,\,ext}$ next. Its properties, along with those of
potential line-driven outflows discussed in § 11.1.4, are predicted to scale
together with the inner boundary of the torus-like region (forming its outer
boundary), both moving outwards (therefore increasing $H/R$ and their
corresponding covering factors) at higher Eddington ratios $\dot{m}$, in a
manner that appears quite similar to recent observational inferences in e.g.
Temple et al. (2021a, 2023).

Figure 6: Predicted disk temperature profiles for the flux-frozen magnetized
($\beta\ll 1$) disks versus radius as Figs. 3 & 4. We compare different
temperatures at each radius calculated in § § 3-10: the mean midplane gas
temperature $T_{\rm midplane}$; effective temperature including heating by
gravitational dissipation and (approximate) self-illumination $T_{\rm eff}$;
characteristic maximum post-shock/reconnection temperature for the disk
internal turbulence $T_{\rm shock}\sim(3/16)\,m_{p}\,v_{\rm
turb}^{2}/k_{B}\sim(3/16)\,m_{p}\,v_{A}^{2}/k_{B}$; virial temperature $T_{\rm
vir}\sim(1/2)(G\,M_{\rm BH}/r)\,(\mu\,m_{p}/k_{B})$; temperature of optically-
thin, directly illuminated surface layers (at $|z|\approx H$; taking the
“warmer” solution when multiple quasi-stable solutions exist), $T^{\rm
warm}_{\rm illum,\,thin}$; and the estimated “hard” electron coronal
temperature $T^{\rm cor}_{e}$ where coronal solutions are supported (§ 9). We
see clear transitions corresponding to the predicted zone/phase boundaries.
Flux-frozen disks have cooler $T_{\rm midplane}$ than SS73 owing to their
lower optical depths, but more importantly exhibit obvious multi-phase
structure with gas at a range of temperatures supported. The similarity
between the predicted temperatures and those needed to explain e.g. the soft
excess in the warm skin, hard coronal X-rays, BLR emission lines, torus
reprocessed warm dust emission, and the thermal continuum/big-blue-bump are
discussed in § 3-10. Note that one cannot directly translate $T_{\rm eff}$
versus $R$ into a predicted emission region size or observed effective
temperature because most of the radiation from the inner region will be
absorbed or reprocessed by the large $H/R$, flared disk (Fig. 2 and below).
Figure 7: Opacity structure of the disks, as Fig. 3. We plot the vertically-
integrated scattering $\tau_{s}$, absorption $\tau_{a}$, and “effective”
absorption $\tau^{\ast}\equiv\sqrt{\tau_{a}\,(\tau_{a}+\tau_{s})}$ optical
depths to the midplane (linestyles labeled). We compare flux-frozen ($\beta\ll
1$; thick black) and thermal-pressure-dominated ($\beta\gg 1$; thin pink)
disks. For the sake of direct comparison note we have recalculated the SS73
models using our more detailed opacity models here (including e.g. dust,
atomic, and iron absorption line opacities), but this does not change the
qualitative results. The optical depths are generally much lower in the
magnetized disks as expected from Fig. 3, but the different temperature and
density structure leads to qualitatively different scalings of the opacities
in many regions. Discontinuities owing to phase changes (e.g. dust
sublimation) are pronounced here. The ratio of e.g. $\tau_{s}>\tau_{a}$ in the
thermalized disk is important for Comptonization (§ 8) while (given the large
$H/R$ and reprocessing; Fig. 2) the effective temperature of e.g. the big blue
bump is influenced by reprocessing at the radii around $R\sim 100\,R_{g}$
where the absorption optical depths first drop below unity, and low mean
absorption $\tau_{a}\lesssim 1$ depths in the BLR-like and scattering zones (3
& 8) are critical for observed line emission and scattering (§ 5 & 10). Below
$\dot{m}\lesssim 0.01$, optical depths become so low that some regions of the
disk cannot cool efficiently.  Figure 8: Midplane pressures of the disks, as
Fig. 3. We plot the compensated pressure, $P_{i}\times(R/R_{g})^{2}$ (for ease
of visualization), at the midplane ($|z|\ll H$), for thermal pressure ($P_{\rm
thermal}\equiv\rho_{\rm mid}\,c_{s,\,{\rm mid}}^{2}=n_{\rm mid}k_{B}T_{\rm
midplane}$; Figs. 4, 6), magnetic+thermal pressure (which scale together as
$P_{\rm mag}\sim|B_{\rm midplane}|^{2}/8\pi\sim(1/2)\rho_{\rm mid}v_{A,\,{\rm
mid}}^{2}\sim P_{\rm turb}\sim(1/2)\rho_{\rm mid}\,v_{\rm turb}^{2}$ in both
the fluz-frozen and SS73-like $\alpha$-disk models), and radiation pressure
($P_{\rm rad}$, accounting for finite optical depth effects by defining this
as the radiation pressure/force per unit area exerted on gas in the midplane,
$P_{\rm rad}\sim\langle\Sigma_{\rm gas}\kappa{F}_{{\rm rad,\,mid}}/c\rangle$
in terms of the midplane flux $F_{\rm rad,\,mid}$). Owing to the much lower
densities (Fig. 4), the pressures – even magnetic pressure/field strength
$|{\bf B}|^{2}$ – in the flux-frozen disks are almost everywhere much lower
than in thermal-pressure-dominated models (see Paper III). While SS73 models
become radiation-pressure-dominated at small radii even at modest accretion
rates $\dot{m}\sim 0.05-1$, the radiation pressure even at $\dot{m}\gg 1$
never becomes much larger in the flux-frozen models than the magnetic pressure
(§ 11.1).

##  5 Zone 3: The Directly-Illuminated, Ionized Upper Disk as the Broad Line
Region

As noted above and shown in Fig. 9, rays from the central disk with
$\cos{\theta}\gtrsim 0.7$ will not intersect the disk, but travel to Zone (1)
(the ISM). Rays with $0.5\lesssim\cos{\theta}\lesssim 0.7$ will intersect
between $r_{\rm sub,\,ext}\lesssim R\lesssim r_{\rm ff}$, where the dust will
not be sublimated, i.e. Zone (2c), where the dust will strongly shield the gas
against ionization by the central source. Rays with $\cos{\theta}\lesssim 0.5$
will intersect at $r\lesssim r_{\rm sub,\,ext}$, where the dust is efficiently
sublimated. Over an intermediate range of angle $\theta$ or disk height $H/R$,
this will create an ionized layer reaching out to $r_{\rm sub,\,ext}$.

###  5.1 3a: The “Fully” Ionized (Dust-Limited) Layers

For the directly illuminated, dust-sublimated (effectively dust-free) disk,
two effects will limit ionization. First, the ionization layer could be
photon-bounded. We can estimate this via a simple Stromgren type argument.
Consider all rays from the central source in some azimuthally-symmetric
$\Delta\cos{\theta}$ opening angle, with some number of ionizing photons per
unit bolometric luminosity $(Q/L)$ which depends on the spectrum of the inner
disk, so $\dot{N}_{\rm ion}/dA\sim
2\pi\,\Delta\cos{\theta}\,f_{\theta}\,(Q/L)\,L_{\rm bol}/(4\pi\,r^{2})$. If
this ionizes a “wedge”, there must be sufficient number of photons per unit
time per unit solid angle to offset recombination in the wedge assuming it is
fully-ionized, given by $\dot{N}_{\rm rec}\approx\int\alpha_{\rm
rec}\,n(r,\,\theta,\,\phi)^{2}\,r^{2}\,{\rm d}\cos{\theta}\,{\rm d}\phi\,{\rm
d}r$ where $\alpha_{\rm rec}\approx 4\times 10^{-13}\,{\rm cm^{3}}$ for the
conditions of interest. Comparing the two (noting $n\propto R^{-2}$ in the
disk model), we see that the “fully ionized” condition (requirement to ionize
to $r\rightarrow r_{\rm sub,\,ext}$) is completely dominated by the ability to
ionize the innermost radii where $n(R)$ is maximized (Eq. 4 and Figs. 3, 4,
7), and this can occur for rays intercepting the disk initially at $R\gtrsim
R_{\rm min,\,strom}$ (so at $\cos{\theta}\gtrsim H/R|_{R=R_{\rm
min,\,strom}}$) given by

$\displaystyle R_{\rm min,\,strom}$ $\displaystyle\approx 1.6\times
10^{16}\,{\rm cm}\,r_{\rm ff,\,5}\,\dot{m}^{6/7}\,\tilde{Q}^{-6/7}$ (29)
$\displaystyle\sim 6\,{\rm ld}\,(L_{\rm bol}/10^{45}\,{\rm
erg\,s^{-1}})^{1/2}\,\epsilon_{r,\,0.1}^{-1/2}\tilde{Q}^{-6/7}\dot{m}^{5/14}$

($\tilde{Q}\equiv\epsilon_{r,\,0.1}\,(Q/L)/(0.1/13.6\,{\rm eV})\sim 1$ for a
typical quasar spectrum; Shen et al. 2020), where

$\displaystyle H/R|_{R=R_{\rm min,\,strom}}\approx
0.3\,(\dot{m}/\tilde{Q})^{1/7}.$ (30)

At heights larger than this, but less than $\sim 0.5$ where the ray intersects
the disk outside $r_{\rm sub,\,ext}$ (Figs. 9 & 12), the wedge will be ionized
but only a fraction of the incident ionizing flux is “needed” and therefore
some of the ionizing radiation will reach the sublimation region and be
attenuated by dust, giving rise to intermediate zones with lower-excitation
ions.

For incident heights $\lesssim 0.3\sim H/R|_{R=R_{\rm min,\,strom}}$, the
ionizing photons will be fully-consumed in a surface layer of radial size $\ll
R$, beyond which (noting that $R_{\rm therm,\,atomic}<R_{\rm
min,\,strom}<R_{\rm sub,\,T_{\rm eff}}$) the layer will resemble the
underlying atomic Zone (4) described below, until reaching the dusty Zone
(2c). But note that the Stromgren approximation assumes negligible electron
scattering or continuum (Kramers) absorption optical depth. The latter
assumption is quite good over the specific range of radii here, but the former
assumption becomes invalid at small $R$, as electron scattering can further
limit ionization by scattering photons out of the narrow “wedge” and into
polar angles where they can escape. Assuming a fully-ionized surface “layer”
at some $R$ and $H$, the electron scattering optical depth to some depth
$d\ell\ll R$ is $\sim 0.35\,\rho(R)\,d\ell$, if we solve for $\ell_{\rm es}$
where this exceeds a factor of $\tau_{\rm crit}\sim$ a couple, we have
$\ell_{\rm es}/R\sim 0.26\,\tau_{\rm crit}\,(R/10^{16}\,{\rm
cm})\,m_{7}^{-1/2}\,r_{\rm ff,\,5}^{-1/3}\,\dot{m}^{-1}$. When $\ell_{\rm
es}\ll R$ ($R\ll R_{\rm es,\,ion}$) where

$\displaystyle R_{\rm es,\,ion}\approx 3.8\times 10^{16}\,{\rm
cm}\,\dot{m}\,(m_{7}\,r_{\rm ff,\,5})^{1/2}\,\tau_{\rm crit}^{-1},$ (31)

this effectively means that the wedge further exterior is shielded by a
relatively narrow Thompson-thick scattering layer. At $R<R_{\rm es,\,ion}$,
the maximum possible reprocessed ionizing luminosity will scale $\propto
R^{4/3}$ or steeper. Thus the reprocessed ionizing luminosity will come
predominantly from the radii outside the maximum of either $R_{\rm es,\,ion}$
or $R_{\rm min,\,strom}$. But it is notable that these are very similar to one
another, so basically $R_{\rm es,\,ion}$ serves to further “cut off” the
fraction reprocessed at $R\lesssim R_{\rm min,\,strom}$ more rapidly. This
defines the boundaries in Figs. 1, 2, 3 and others.

Note that under some conditions – mostly notably when $\dot{m}\gg 1$
(especially if the accretion becomes radiatively inefficient as expected in
this regime) – $R_{\rm min,\,strom}$ and $R_{\rm es,\,ion}$ can expand
outwards to become larger than the outer radius $r_{\rm sub,\,ext}$ of this
zone. We stress that this does not mean the BLR-like region would vanish.
Rather it simply means that, as above, the fractional contribution from each
radius $R<r_{\rm sub,\,ext}$ will fall off as $\propto(R/r_{\rm
sub,\,ext})^{4/3}$ or so as the gas able to freely absorb and emit such lines
will be confined to a narrower surface layer. This simply means that the BLR-
like emission will be more strongly dominated by the gas just inside the dust
sublimation radius $r_{\rm sub,\,ext}$ when $\dot{m}\gg 1$. Qualitatively (as
discussed in § 11) as these move to larger radii and therefore larger $H/R$,
this implies a larger covering/reprocessing factor for hot dust and the BLR,
and stronger outflows, associated with weaker X-rays, at higher $\dot{m}$.

Depending on the photon energies of interest, there will be a region from
$\sim R_{\rm min,\,strom}$ to $r_{\rm sub,\,ext}$, corresponding to a covering
factor of

$\displaystyle f_{\rm cover}^{\rm 3a}\sim f_{\rm cover}^{\rm r_{\rm
sub,\,ext}}-f_{\rm cover}^{R_{\rm min,\,strom}}\sim 0.2$ (32)

(from $H/R\sim 0.3-0.5$; Figs. 9 & 12) which can reprocess a fraction up to
$f_{\rm reprocessed}\sim f_{\theta}f_{\rm cover}^{\rm 3a}\sim 0.12-0.16$ of
the bolometric luminosity from the central source (depending on how isotropic
or not its emission is). The total mass or volume subtended by this will be
dominated by the largest radii, $R_{\rm max,\,3a}\sim 100\,{\rm ld}\,(L_{\rm
bol}/10^{45}\,{\rm erg\,s^{-1}})^{1/2}$ (ld = light-days), while the smaller
radii where the efficiency of reprocessing peaks correspond to $R_{\rm
min,\,3a}\sim 6\,{\rm ld}\,(L_{\rm bol}/10^{45}\,{\rm erg\,s^{-1}})^{1/2}$.
Taking a typical linewidth to be twice the Keplerian velocity (for disky
material), this range of radii corresponds to $\Delta v\sim 2\,v_{K}$ in the
range

$\displaystyle 1400\,m_{7}^{1/4}\,\dot{m}^{-1/4}\,T_{\rm
sub,\,1500}\lesssim\Delta v\lesssim
6200\,m_{7}^{1/4}\,\dot{m}^{-3/7}\,\tilde{Q}^{3/7}.$ (33)

The temperatures in this zone will be set by photonionization equilibrium to
$\sim 10^{4}\,$K (Fig. 6), so the sonic Mach numbers will be $\sim 100$ (for
more details, see Fig. 10) and there should be inhomogeneous/clumpy/turbulent
structure down to scales of order the sonic length $R_{\rm sonic}\sim
H/\mathcal{M}_{s}^{2}$ as noted above (see Fig. 4), which also corresponds (by
definition) roughly to the Sobolev length in a supersonically turbulent
medium, relevant for the coherence length of resonant line emission/absorption
even in the smoother or lower-density gas (Fig. 11). Calculating this properly
with our model scalings,

$\displaystyle R_{\rm sonic}\sim 0.5\times 10^{12}\,{\rm
cm}\,\dot{m}^{11/7}\,\tilde{Q}^{-11/7}\,m_{7}^{1/12}\,r_{\rm ff,\,5}^{1/6}$
(34)

at the inner radii and $\sim 10^{14}\,{\rm cm}\,\dot{m}^{11/12}\,T_{\rm
sub,\,1500}^{-11/3}$ at the outer radii of this zone. The volume-weighted mean
disk densities in the zone range from $\sim 3\times 10^{8}\,{\rm
cm^{-3}}\,\tilde{Q}^{12/7}\,r_{\rm
ff,\,5}^{1/2}\,m_{7}^{-1/2}\,\dot{m}^{-5/7}$ (near the inner region of the
zone where most of the emission occurs) to $\sim 10^{6}\,{\rm
cm^{-3}}\,(Q_{\rm em}/Q_{\rm abs})\,T_{\rm
sub,\,1500}^{4}\,\epsilon_{r,0.1}^{-1}\,(m_{7}\,r_{\rm ff,\,5})^{-1/2}$ (at
the outer zone boundary), but the high sonic Mach numbers mean that the local
gas density fluctuations should follow something like the usual lognormal for
these mach numbers (with $\mathcal{M}_{A}\sim 1$; see Ostriker et al. 2001;
Hopkins 2013b; Beattie et al. 2021), so $\sim 90\%$ of the gas mass will be
between $\sim 0.2-6000\times\langle n_{\rm gas}(R)\rangle_{\rm vol}$. The gas-
mass weighted mean or median will be around

$\displaystyle\langle n\rangle_{\rm mass}\sim 50\,\langle n_{\rm
gas}(R)\rangle_{\rm vol}\sim 10^{10}\,{\rm cm^{-3}}\,\tilde{Q}^{12/7}\,r_{\rm
ff,\,5}^{1/2}\,m_{7}^{-1/2}\,\dot{m}^{-5/7}$ (35)

giving $\sim 10\%$ of the gas in densities more like

$\displaystyle n_{10\%}^{\rm inner}\gtrsim 2\times 10^{12}\,{\rm
cm^{-3}}\,m_{7}^{1/4}\,\tilde{Q}^{16/7}\,\dot{m}^{-9/7}$ (36)

(at the radii which dominate the total reprocessed emission of the zone) or

$\displaystyle n_{10\%}^{\rm outer}\gtrsim 10^{9}\,{\rm
cm^{-3}}\,m_{7}^{1/4}\,T_{\rm sub,\,1500}^{16/3}\,\dot{m}^{-1/3}$ (37)

(at the outer zone boundary), as shown more explicitly in Fig. 4. The
dimensionless version of this in terms of $\rho/\rho_{\rm midplane}$ or
corresponding “clumping factors” of the gas is shown in Fig. 10 in more
detail. Given this, the gas will be multi-phase for the same reasons and with
the same character as in Zone (5) below.

Because of how the various disk properties scale, the volume-weighted mean
ionization parameter for directly-illuminated gas depends quite weakly on
radius or on any other parameter except $\tilde{Q}$ (the shape of the actual
ionizing spectrum) and the local density relative to mean: $\xi\equiv({\rm
d}Q/{\rm d}A)/(c\,n_{e})\sim 3\,(f_{\theta}/0.5)\tilde{Q}\,(n_{\rm
gas}/\langle n_{\rm gas}(R)\rangle_{\rm mass})^{-1}\,(m_{7}/r_{\rm
ff,\,5})^{1/2}\sim\tilde{Q}\,(m_{7}\,R/10^{16}\,{\rm cm})^{1/6}\,(n_{\rm
gas}/\langle n_{\rm gas}(R)\rangle_{\rm mass})^{-1}$ (where in the latter
expression we take the expected $r_{\rm ff,\,5}\sim m_{7}^{1/2}$ and the
geometric mean of the two limiting cases for $f_{\theta}$ given in § 2.4, as
the difference between the two is small at these radii of interest). So $\xi$
will vary primarily with these same density fluctuations, giving $\xi\sim
0.01-1\tilde{Q}$ for the denser gas in the $1-2\sigma$ overdensity clumps
(i.e. the $\sim 10-30\%$ mass range, at densities $n\sim 10^{12}\,{\rm
cm^{-3}}$ at the smallest BLR radii and $n\sim 10^{9}\,{\rm cm^{-3}}$ at the
largest BLR radii) quoted above. We show this computed more explicitly, for
directly ionized gas at the surface layer, in Fig. 11.

As discussed below, these scalings are largely robust (modulo order-unity
coefficients) to whether this upper layer is in inflow with the disk or
outflow (if the incident luminosity is high enough to drive e.g. a line-driven
wind off the surface; discussed quantitatively below).

All of these properties are remarkably similar to the canonical
observationally-inferred properties of the BLR. Specifically, this includes
the emitting gas densities, temperatures, and ionization parameters (see
typical observed values in Peterson 2006 and references therein); strength of
the density fluctuations or range of densities at a given radius (compare
Dexter & Agol, 2011); observed effective “emitter” or “cloud” size which
should be similar to the predicted Sobolev length in both the discrete-cloud
scenario (per Krolik et al. 1981; Krolik 1999) and smooth emitting-structures
case (per Arav et al. 1998; Laor et al. 2006); linewidths (compare Laor,
1991); BLR covering factors and/or scale-heights constrained by direct imaging
and reverberation mapping in e.g. Gravity Collaboration et al. (2018);
luminosities/reprocessing fraction (compare Kaspi et al., 2005; GRAVITY
Collaboration et al., 2024); inner and outer characteristic BLR radii (Du et
al., 2015; GRAVITY Collaboration et al., 2020, 2021; Woo et al., 2023);
emitting gas mass (compare e.g. the BLR emitting gas mass of $\sim 40\,{\rm
M_{\odot}}$ for NGC 3227 in Devereux 2021 to the predicted $\sim 30\,{\rm
M}_{\odot}$ here from Fig. 5 for similar luminosity and BH mass); and thick-
disk like geometry and global kinematics (as observed in Gravity Collaboration
et al., 2018; GRAVITY Collaboration et al., 2021, 2024).

All of these will be the subject of more detailed study in future papers in
preparation, in which we will make more detailed quantitative comparisons with
specific well-observed AGN. A rigorous comparison, however, requires detailed
radiation transport calculations performed on the full numerical simulations,
in order to properly compare observables, and so is outside the scope of our
predictions here, as we discuss further in § 13. But compared this to, for
example, a thermal-pressure-dominated SS73-like $\alpha$-disk. In that case,
at the radii of observed BLRs, every disk property is many orders of magnitude
different from the properties observed in the BLR (the thermal-pressure-
dominated disk prediction is that the disk would be vastly more dense, hotter,
geometrically thinner, more massive, less turbulent, etc.). We note in § 13
that the same is true of e.g. marginally-self-gravitating disk models. So in
this model classes, it is simply not possible that the BLR comes from the disk
itself. Of course, in those models, as noted in § 1, one could posit that the
BLR is something different entirely which sits “above” the disk, held up by
different physics. Of course, if one invokes an arbitrary magnetic field or
other physics to hold up BLR “clouds” and fits or adjusts the model to
reproduce all the same observables as above, then there will be no obvious way
to observationally discriminate between these models (though there may still
be measurable differences in the kinematics owing to the effects of the thin
disk in the midplane). The key difference between those models from § 1 and
the model here is that the BLR properties we discuss above are predictions of
the disk model itself – they arise inevitably from an accretion disk that
obeys our two simple ansatz in § 2.2.

###  5.2 3b: The Shielded, Mostly Neutral Layer

The inner boundary of the fully ionized zone above corresponds roughly with
the radius $R_{\rm sub,\,T_{\rm eff}}$ – i.e. the radii where at the midplane,
the gas transitions to mostly neutral but dust-free. At $|z|\sim H$ at these
radii, the gas will be partially ionized but the electron-scattering effect
above shields most of the layer which would naively be directly illuminated
(Fig. 7). Thus the neutral zone can extend through nearly to the disk surface,
and this zone is basically identical to Zone (4) we calculate next, except for
the shielding layer (Fig. 2).

##  6 Zone 4: The Cool, Neutral Dust-Free Disk Midplane

Inside of the midplane (shielded/un-illuminated) sublimation radius, the disk
is dust free but still relatively cool with temperatures $T\sim
10^{3}-10^{4}\,$K (Fig. 6) so the gas should be mostly neutral, with free
electron fractions $x_{e}\sim 0.01\,x_{e,\,0.01}$ typical of gas in this
temperature range (see Paper I, Wolfire et al. 1995; Tielens 2005; Draine 2011
for more detailed discussion). The optical depth to starlight or ionizing
photons through the midplane from the central source (directly through the
midplane of the disk) is enormous (Compton-thick) and so this remains shielded
to those sources of irradiation. This zone will be less interesting
observationally, though we calculate its properties for the sake of
completeness and defining its boundaries.

Here the scattering opacity is now dominated by electron scattering with
$\kappa_{\rm es}\sim 0.35\,x_{e}$, but the absorption opacity is largely
dominated by H- and warm molecular opacities. In this temperature range we can
adopt the usual approximations $\kappa_{H^{-}}\approx 1.1\times 10^{-25}\,{\rm
cm^{2}\,g^{-1}}\,Z^{1/2}\,\rho^{1/2}\,T^{7.7}\sim 10^{-7}\,(n_{\rm
gas}/10^{8}\,{\rm cm^{-3}})^{1/2}\,\tilde{Z}^{1/2}\,T_{3000}^{77/10}$, and
$\kappa_{\rm mol}\sim 0.0014\,\tilde{Z}\,f_{\rm mol}$ where typical $f_{\rm
mol}\sim 0.001\,f_{\rm mol,\,0.001}$ at the midpoint of the temperature range
(Paper I and Hollenbach & McKee 1979; Tielens 2005). So we have
$\kappa_{s}\gg\kappa_{a}$ and the effective opacity for thermalization is
$\kappa^{\ast}=\sqrt{\kappa_{a}\,(\kappa_{s}+\kappa_{s})}\approx\sqrt{\kappa_{a}\,\kappa_{s}}$.
This gives a vertically-integrated scattering optical depth

$\displaystyle\tau_{s}\approx\tau\sim 0.0017\,x_{e,\,0.01}\,(\dot{m}\,r_{\rm
ff,\,5}/\epsilon_{r,\,0.1})^{1/2}\,T_{\rm sub,\,1500}^{2}\,\hat{R}^{-5/6},$
(38)

scaling to the outer radius $R_{\rm sub,\,T_{\rm eff}}$ (defined including
illumination, though this makes no important differences here) with
$\hat{R}\equiv R/R_{\rm sub,\,T_{\rm eff}}$, and effective opacity

$\displaystyle\tau^{\ast}\approx$
$\displaystyle\,\sqrt{\tau_{s}\,\tau_{s}}\sim 3.3\times
10^{-5}\,(x_{e,\,0.01}\,f_{\rm mol,\,0.001}\,\dot{m}\,r_{\rm
ff,\,5}\,\tilde{Z})^{1/2}$
$\displaystyle\times\epsilon_{r,\,0.1}^{-1/2}\,T_{\rm
sub,\,1500}^{2}\,\hat{R}^{-5/6}$ (39)

or (if H- dominates)

$\displaystyle\tau^{\ast}\sim$ $\displaystyle\,6\times
10^{-8}\,T_{1000}^{3.9}\,x_{e,\,0.01}^{0.5}\,\tilde{Z}^{0.25}\,\dot{m}^{0.45}\,m_{7}^{-0.18}\,r_{\rm
ff,\,5}^{0.73}$ $\displaystyle\times T_{\rm
sub,\,1500}^{3.2}\,\epsilon_{r,\,0.1}^{-0.8}\,\hat{R}^{-4/3},$ (40)

all shown in Fig. 7. So the disk is optically-thin, so long as temperatures
remain below $\sim 10^{4}$ K (given the steep H- opacity dependence).

Initially, given the $T_{\rm eff}\sim 1500\,$K at the sublimation radius, this
would give the result that the absorption, and therefore continuum emission
opacity is dominated by optically-thin molecular transitions. The optically-
thin cooling rate from said transitions, $F_{\rm
cool,\,thin}=4\pi\,\kappa_{\rm em,\,mol}\,\Sigma_{\rm gas}\,\sigma_{B}\,T_{\rm
mol}^{4}\sim 5\times 10^{6}\,f_{\rm mol}\,\dot{m}^{1/2}\,r_{\rm
ff,\,5}^{1/2}\,T_{\rm
sub,\,1500}^{2}\,\epsilon_{r,\,0.1}^{-1/2}\,\tilde{Z}\,\hat{R}^{-5/6}\,(T_{\rm
mol}/1500\,{\rm K})^{4}$ (in cgs) is much less than the gravitational heating
flux $F_{\rm grav}\sim 4\times 10^{6}\,m_{7}^{1/5}\,r_{\rm
ff,\,5}^{3/5}\,\dot{m}^{-4/5}\,T_{\rm
sub,\,1500}^{36/5}\,\epsilon_{r,\,0.1}^{-9/5}\,\hat{R}^{-3}$ at $R<R_{\rm
sub,\,T_{\rm eff}}$, even for $f_{\rm mol}\rightarrow 1$. For more realistic
lower $f_{\rm mol}$ the cooling luminosity would be dominated by CII at these
temperatures,666We adopt the standard expressions for effectively optically-
thin cooling: $\Lambda_{\rm CII}\sim 9.6\times
10^{-24}\,\tilde{Z}\,x_{e}\,T^{-1/2}\,e^{-92\,{\rm K}/T}$ and $\Lambda_{\rm
Ly\alpha}\sim 1.8\times 10^{-11}\,x_{e}^{2}\,T^{-1/2}\,e^{-120000\,{\rm K}/T}$
in ${\rm erg\,s^{-1}\,cm^{3}}$ (see e.g. Field et al., 1969; Dalgarno &
McCray, 1972). but that also gives too-low a cooling luminosity $F_{\rm
cool,\,thin}^{\rm CII}\sim 7500\,x_{e,\,0.01}\,m_{7}^{-7/10}\,r_{\rm
ff,\,5}^{7/5}\,\dot{m}^{-3/10}\,\tilde{Z}\,\epsilon_{r,\,0.1}^{-17/10}\,T_{\rm
sub,\,1500}^{34/5}\,\hat{R}^{-17/6}\,T_{1500}^{-1/2}$, so the molecular and
atomic gas would quickly heat to much larger equilibrium temperatures $\sim
2\times 10^{4}\,{\rm K}\,m_{7}^{1/72}\,\dot{m}^{-13/72}\,r_{\rm
ff,\,5}^{-1/12}\,\tilde{Z}^{-1/4}\,T_{\rm sub,\,1500}^{13/18}\,f_{\rm
mol,\,0.001}^{-1/4}\,\hat{R}^{-13/24}$ to balance heating, but as the
temperature rises the molecules dissociate further and the CII cooling rate
actually decreases, while H- opacity further increases. The medium therefore
becomes primarily atomic and the temperature quickly rises until it reaches
the equilibrium where optically-thin H cooling balances the gravitational heat
flux. The cooling at these temperatures will be dominated by a combination of
H- and Ly-$\alpha$, with either or both in combination giving an extremely
similar, nearly-constant temperature (owing to the extremely strong
temperature dependence of both) around

$\displaystyle T\sim 6000\,{\rm K}\,T_{\rm
sub,\,1500}^{0.24}\,(\dot{m}\,\hat{R})^{-0.1}$ (41)

or

$\displaystyle T\sim 6100\,{\rm
K}\,(1-0.04\,\ln{[x_{e,\,0.01}^{2}\,\dot{m}]}),$ (42)

respectively. Thus the temperature rapidly jumps within this zone to $\sim
6000\,$K, then rises more slowly with decreasing $R$ (Fig. 6), with $f_{\rm
ion}\sim 0.01-0.1$, akin to the warm neutral medium (WNM) in the ISM (Field et
al., 1969; Wolfire et al., 1995, 2003). Note that in the above, we have used
the gravitational flux – again since the region is extinction and effectively
optically-thin, the upper layer of the disk should be directly illuminated
(discussed below) but the midplane should see little reprocessed indirect
illumination flux, and this will not be effective at heating it. Even if we
did assume that $F_{\rm illum}$ were thermalized, the results for the
temperature and chemical structure would be the nearly identical, given the
strong temperature dependence of the cooling/opacity physics. Moreover at the
inner regions of this zone, the effective illumination flux $F_{\rm illum}$
becomes comparable to and smaller than $F_{\rm grav}$, with $F_{\rm
grav}\gtrsim F_{\rm illum}$ for $\hat{R}\lesssim
0.03\,\dot{m}^{-0.6}\,m_{7}^{0.15}\,r_{\rm ff,\,5}^{0.45}$.

In any case, the warmer temperatures (which mean much larger H- opacities)
plus rising surface densities as $R$ decreases means that $\tau^{\ast}$
increases, and crosses unity at $R\lesssim R_{\rm therm,\,atomic}$ where
$R_{\rm therm,\,atomic}\sim 0.014\,R_{\rm sub,\,T_{\rm
eff}}\,\tilde{Z}^{0.05}\,m_{7}^{0.003}\,\dot{m}^{0.03}\,r_{\rm
ff,\,5}^{0.36}\,\epsilon_{r,\,0.1}^{-0.6}\,T_{\rm sub,\,1500}^{2.4}$ or

$\displaystyle R_{\rm therm,\,atomic}\sim 1.3\times 10^{15}\,{\rm
cm}\,m_{7}^{0.6}\,\dot{m}^{0.63}\,r_{\rm ff,\,5}^{0.16}\,\tilde{Z}^{0.05}$
(43)

(Fig. 7). The disk will then thermalize and radiate with an effective
temperature balancing the gravitational accretion flux, giving

$\displaystyle T_{\rm eff}(R<R_{\rm therm,\,atomic})\sim$
$\displaystyle\,11000\,{\rm K}\,m_{7}^{0.05}\,r_{\rm
ff,\,5}^{-0.12}\,\dot{m}^{-0.22}\,\tilde{Z}^{-0.04}$
$\displaystyle\times(R/R_{\rm therm,\,atomic})^{-3/4}$ (44)

(Figs. 6, 8). The H- rapidly dissociates and the disk begins to efficiently
“self-ionize” from its own thermal emission at these temperatures. Thus the
midplane radial dynamic range of this zone in radius is modest but significant
($\sim 1\,$dex; Figs. 1-2). At this self-ionization radius, the covering
factor of the disk is $H/R\sim 0.21\,m_{7}^{0.1}\,r_{\rm
ff,\,5}^{-0.14}\,\dot{m}^{0.1}\,\tilde{Z}^{0.008}$ (Fig. 9). For the reasons
above, these scalings are nearly identical with or without including external
illumination effects.

The intrinsic gravitational luminosity from this region is $\sim 2\times
10^{43}\,{\rm erg\,s^{-1}}\,m_{7}^{1.4}\,\dot{m}^{0.37}\,r_{\rm
ff,\,5}^{-0.16}\,\tilde{Z}^{-0.05}\sim 0.017\,L_{\rm
bol}\,m_{7}^{0.32}\,\dot{m}^{-0.63}$, so still a small fraction of $L_{\rm
bol}$. Since most of this comes from the inner (and therefore warmer) regions,
it will primarily emerge in a mix of Ly$\alpha$ and H- continuum. The former
will be very strongly scattered in the external medium around this zone and is
much smaller than the reprocessed Ly$\alpha$ luminosity of the directly-
illuminated regions, while the latter will be much smaller than the thermal
continuum at the same wavelengths from the warm disk at smaller radii, so both
will be largely undetectable. Given the cooling rates and temperatures, a
smaller fraction $\sim 6\times
10^{-5}\,(x_{e}/0.1)\,\dot{m}^{1/2}\,\tilde{Z}\,m_{7}^{-0.16}$ of $L_{\rm
bol}$ will emerge in CII lines.

Figure 9: Illustration (as Figs. 1-2) of covering factors for the different
zones. We plot the relative height $z/R$ (top and middle) and polar angle from
the assumed disk axis $\cos{\theta}=z/r=z/\sqrt{R^{2}+z^{2}}$ (bottom) of the
zones ($z/R=H/R$ as solid black line). Zones are labeled (interior to horizon
in grey shaded), versus radius with their numbers (top/bottom) and
observational descriptors (middle). On the linear scale here, the predicted
flaring of the disk and large scale-heights ($H/R\approx(R/R_{\rm
BHROI})^{1/6}$) are more obvious. As discussed in § 10 and seen here, this
leads to large $\sim 10-20\%$ covering factors for specific regions like the
BLR and directly illuminated hot torus, as well as an extended scattering
zones (6-8). Together we see that most of the light from the innermost disk
should be scattered or reprocessed in some form (§ 10).

##  7 Zone 5: The Warm Ionized Multi-Phase Disk

Figure 10: Radial profile of the predicted midplane turbulent sonic Mach
number $\mathcal{M}_{s}\equiv\langle v_{\rm turb}\rangle/\langle c_{s,\,{\rm
mid}}\rangle$ (as Fig. 3) and corresponding expected gas clumping factor
$\langle\rho^{2}\rangle\sim 1+(b\mathcal{M}_{s})^{2}$ (where $b\approx 1/3$
reflects the compressive component expected; Fig. 4) given the predicted
trans-Alfvénic, supersonic turbulence in the flux-frozen $\beta\ll 1$ disks.
Despite various changes in temperature and structure, typical sonic Mach
numbers for a range of conditions tend to be around $\sim 100$, corresponding
to $\sim 5-10\%$ of the gas mass being in multi-phase denser structures with
densities $\sim 10^{3}-10^{4}$ times the volume-averaged mean density within
the disk at a given radius.

Figure 11: Top: Radial profile of the Sobolev velocity coherence or sonic
length (the two are similar in supersonic turbulence) of the disk, for
different parameters (labeled). Note this is defined in the midplane, so is
slightly different in the surface/illuminated layers, but the order-of-
magnitude values are similar. At BLR-like radii there is a plateau around
$\sim 10^{10}-10^{13}\,{\rm cm}$, while in the inner disk there is another
extended region where this is around $\sim 10^{10}\,{\rm cm}$. These will
correspond to the characteristic sizes, in optically-thin regions like the
BLR, of individual velocity-and-density coherent gas emission structures (e.g.
BLR “clouds,” loosely defined). Bottom: Same for the mass-weighted median
surface ionization parameter $\xi\equiv\langle(Q/L)\,F_{\rm
illum}/(n\,c)\rangle_{\rm mass}$, where $Q/L$ gives the ratio of ionizing
photon production rate to total luminosity, $F_{\rm illum}$ is the mean
directly-illuminated flux, and $n$ the gas density at $|z|\sim H$, all
evaluated assuming illumination with a flat-disk pattern
($f_{\theta}\propto\cos{\theta}$) and ignoring scattering/reprocessing and
optical depth effects interior to the disk. The shaded range shows the gas-
mass-weighted $90\%$ inclusion interval owing to the range of densities
predicted (Fig. 4). The range of $\xi$ is remarkably stable, only dropping
weakly for $\dot{m}\gg 1$ because the disks become radiatively inefficient
while being denser, while the scatter is large and similar to that observed in
e.g. the BLR (§ 5).

Now consider the midplane zone interior to the thermalization radius of Zone
(4), which occurs (as we noted before) very close to where the effective
temperature exceeds $\sim 10^{4}$ K (compare Figs. 2 & 6). Interior to this
radius $R_{\rm therm,\,atomic}$, the increasing temperatures rapidly
dissociate any atomic hydrogen. For convenience we will in this section scale
the radii to the outer boundary of the zone: $\tilde{R}\equiv R/R_{\rm
therm,\,atomic}$. Again this midplane zone is of less direct interest
observationally, as we note below it does not much influence the observed
spectrum.

The scattering opacity in this (now mostly-ionized) regime is dominated by
electron scattering with $\kappa_{s}\approx\kappa_{\rm es}\approx 0.35$.
Absorption opacity comes primarily from a combination of the usual Kramers
continuum opacity $\kappa_{\rm K}\approx 7\times
10^{25}\,(0.014\,\tilde{Z}+0.001)\,\rho\,T^{-7/2}\sim
0.0006\,m_{7}^{-0.71}\,r_{\rm
ff,\,5}^{0.19}\,\dot{m}^{-0.26}\,\tilde{Z}^{0.90}\,T_{4}^{-7/2}\,\tilde{R}^{-2}$
(where $T_{4}\equiv T/10^{4}\,{\rm K}$; Mihalas & Mihalas 1984; Rybicki &
Lightman 1986; Krolik 1999), plus iron line opacities (which effectively cover
the continuum of interest here), which over the temperature range
$10^{4}\,{\rm K}\lesssim T\lesssim 10^{6}\,{\rm K}$ can be reasonably
approximated by $\kappa_{\rm Fe}\sim 0.00036\,(\rho/10^{10}\,{m_{p}\,{\rm
cm^{-3}}})\,T_{4}^{-2}\,\tilde{Z}$ (Jiang et al., 2015) – larger by a factor
$\sim 2.2\,T_{4}^{3/2}$ compared to $\kappa_{\rm K}$. We are again in the
regime with $\kappa_{s}\gg\kappa_{a}$, so the effective opacity is
$\kappa^{\ast}\approx\sqrt{\kappa_{a}\,\kappa_{s}}$. This gives a scattering
optical depth

$\displaystyle\tau_{s}\sim 6\,\dot{m}^{0.47}\,r_{\rm
ff,\,5}^{0.2}\,m_{7}^{-0.003}\,\tilde{Z}^{-0.04}\,\tilde{R}^{-5/6},$ (45)

and effective optical depth (assuming Fe dominates since we are near
$10^{4}\,$K)

$\displaystyle\tau^{\ast}\sim 0.36\,r_{\rm
ff,\,5}^{0.3}\,m_{7}^{-0.36}\,\dot{m}^{0.34}\,\tilde{Z}^{0.41}\,T_{4}^{-1}\,\tilde{R}^{-11/6}$
(46)

(Fig. 7).

From this, we see that there can be a relatively small continuation of the
multi-phase structure from Zone (4), as the ionization causes a drop in the
absorption opacity to reduce $\tau^{\ast}<1$. In this region, the optically-
thin line cooling around $T\sim 10^{4}-10^{5}\,$K ($\Lambda\sim 10^{-22}\,{\rm
erg\,s^{-1}\,cm^{-3}}$) can easily offset the gravitational heating $F_{\rm
grav}$,777Even if we ignored the dominant H and He recombination line cooling
mechanisms and just used the effective iron continuum opacity, this would
balance $F_{\rm grav}$ for a slightly higher temperature $\sim 2\times
10^{4}\,$K. and we noted above for Zone (4) that even if external illumination
were possible it is now at these radii smaller than $F_{\rm grav}$ (but here
we are shielded by the scattering layer since $\tau_{s}>1$, and effectively
optically-thin, so this radiation should not effectively heat this region
anyways). The regime will therefore be akin to the usual optically-thin,
supersonically turbulent ISM at these temperatures, and unstable to multi-
phase structure (the sonic scale can reach as small as

$\displaystyle R_{\rm sonic}\sim H/\mathcal{M}_{s}^{2}\sim 10^{11}\,{\rm
cm}\,m_{7}^{0.11}\,r_{\rm
ff,\,5}^{0.45}\,T_{4}\,\dot{m}^{1.16}\,\tilde{Z}^{0.09}\,\tilde{R}^{11/6};$
(47)

Fig. 10) with shocks forming hot gas at temperatures

$\displaystyle T_{\rm shock}\sim\mathcal{M}_{s}^{2}\,10^{4}\,{\rm K}$ (48)

(Fig. 6) which can cool but slowly enough (given the decreasing cooling time
for hotter gas) that a non-negligible fraction of such gas persists while the
rest cools to $\sim 10^{4}\,$K below which temperature the cooling rate drops
rapidly.888We can estimate the steady-state fraction $f_{\rm hot}^{\rm
midplane}$ of still-hot post-shock gas by computing the ratio of optically-
thin cooling time to inter-shock timescale $\Delta t_{\rm shock}$. Given the
turbulence is supersonic with turnover scale $\mathcal{O}(H)$, the latter is
necessarily about equal to the dynamical time $\Delta t_{\rm shock}\sim t_{\rm
dyn}\sim 1/\Omega$ (Hopkins, 2013a). So computing $t_{\rm cool}^{\rm
shock}/\Delta t_{\rm shock}=\Omega\,(3/2)\,k_{B}\,T_{\rm
shock}/(4\pi\,\kappa_{\rm K}\,\sigma_{B}\,T_{\rm shock}^{4}\,m_{p})$ for
typical post-shock gas, assuming an initial strong shock $T_{\rm
shock}\approx(3/16)\,(m_{p}/k_{B})\,v_{\rm turb}^{2}$, $\rho_{\rm
shock}\approx 4\,\rho_{\rm gas}$ (before cooling), gives $f_{\rm hot}^{\rm
midplane}\sim t_{\rm cool}^{\rm shock}/t_{\rm dyn}\sim
0.005\,\tilde{R}^{1/6}\,m_{7}^{0.6}\,r_{\rm
ff,\,5}^{-0.64}\,\dot{m}^{-0.89}\,\tilde{Z}^{-1}$.

Given the $\tau_{\ast}$ calculated above, as we move inwards in radius $R$,
the disk will (for high accretion rates at least) soon exceed $\tau^{\ast}\sim
1$ and thermalize (with $T\rightarrow T_{\rm eff}$) at $\tilde{R}\lesssim
0.36\,\dot{m}^{0.52}\,m_{7}^{-0.37}\,r_{\rm ff,\,5}^{0.38}\,\tilde{Z}^{0.41}$,
or $R\lesssim R_{\rm therm,\,ion}$ where

$\displaystyle R_{\rm therm,\,atomic}\sim$ $\displaystyle 4.8\times
10^{14}\,{\rm cm}\,m_{7}^{3/13}\,\dot{m}^{15/13}\,r_{\rm
ff,\,5}^{7/13}\,\tilde{Z}^{6/13}$ $\displaystyle\sim
160\,R_{g}\,\dot{m}^{15/13}\,r_{\rm
ff,\,5}^{7/13}\,m_{7}^{-10/13}\,\tilde{Z}^{6/13}$ (49)

(Fig. 2). Note that at sufficiently large $\dot{m}$ or $m_{7}$ this expression
will over-estimate $R_{\rm therm,\,ion}$ because the scalings we used will
extrapolate beyond the radii where $T_{\rm eff}$ drops too low to maintain the
disk as fully-ionized and we need to change our opacity expressions. In that
limit we would have instead

$\displaystyle R_{\rm therm,\,ion}\sim$ $\displaystyle 7.7\times 10^{14}\,{\rm
cm}\,m_{7}^{9/22}\,r_{\rm ff,\,5}^{7/22}\,\dot{m}^{9/11}\,\tilde{Z}^{3/11}$
$\displaystyle\sim 260\,R_{g}\,m_{7}^{-13/22}\,r_{\rm
ff,\,5}^{7/22}\,\dot{m}^{9/11}\,\tilde{Z}^{3/11}$ (50)

(Fig. 12). These expressions are not strongly modified by illumination: for
the most extreme “total” illumination case anticipated at very large $\dot{m}$
(e.g. physically interior to an extended scattering layer), where we neglect
gravitational flux but assume all of the interior light is reprocessed, we
obtain $R_{\rm therm,\,ion}^{\rm illum}\rightarrow 230R_{g}m_{7}^{-5/8}r_{\rm
ff,\,5}^{7/16}\dot{m}^{15/16}\tilde{Z}^{3/8}\epsilon_{r,\,0.1}^{-3/16}$, very
similar for most practical purposes. So in any case zone (5a) is small at high
$\dot{m}$.

At this thermalization radius – the transition radius between zones (5) and
(6) – the effective temperature will be

$\displaystyle T_{\rm eff}(R\sim R_{\rm therm,\,ion})\sim 2\times 10^{4}\,{\rm
K}\,m_{7}^{0.33}\,r_{\rm ff,\,5}^{-0.4}\,\dot{m}^{-0.62}\,\tilde{Z}^{-0.35}$
(51)

(Fig. 6), and the gravitational cooling luminosity is $\sim 3\times
10^{43}\,m_{7}^{23/13}\,r_{\rm
ff,\,5}^{-7/13}\,\dot{m}^{-2/13}\,\tilde{Z}^{-6/13}\sim 0.03\,L_{\rm
bol}\,m_{7}^{10/13}\,\dot{m}^{-15/13}\,r_{\rm
ff,\,5}^{-7/13}\,\tilde{Z}^{-6/13}\,\epsilon_{r,\,0.1}^{-1}$. So the direct
emission is small but not completely negligible. However, most of the line
emission will be strongly scattered and sub-dominant to the line emission from
the directly illuminated region, and potentially Comptonized, as we discuss
below. But the continuum radiation can contribute non-negligibly to the total
NIR emission, potentially. For these reasons it is also worth noting that at
the inner boundary of this zone the scale height will be

$\displaystyle H/R|_{R=R_{\rm therm,\,ion}}\sim 0.18\,m_{7}^{1/26}\,r_{\rm
ff,\,5}^{-1/13}\,\dot{m}^{5/26}\,\tilde{Z}^{1/13}$ (52)

(Fig. 9).

##  8 Zone 6: The Central Disk

###  8.1 6a: The Thermalized, Modified Black-Body, UV-Emitting Disk

Interior to $R\lesssim R_{\rm therm,\,ion}$, in the midplane, the disk is
effectively optically-thick (Fig. 7) and ionized (Fig. 6). At these radii we
have $F_{\rm grav}\gg F_{\rm illum}$ (Figs. 2 & 9), so direct illumination
owing purely to the disk flaring is negligible as a heating mechanism for the
disk at a given radius (though reprocessing in general, from scattering by the
gas at $|z|>H$, may not be negligible, as we discuss below).

Now scaling by $x_{g}\equiv R/R_{g}$ for these smaller radii, it is useful to
note that the scattering opacity is still primarily from Thompson
$\kappa_{s}\approx\kappa_{\rm es}\approx 0.35$, with optical depth

$\displaystyle\tau_{s}\approx\tau\approx 1000\,\dot{m}\,m_{7}^{-1/3}\,r_{\rm
ff,\,5}^{1/3}\,x_{g}^{-5/6}$ (53)

($\gg 1$ at all radii of interest out to $R_{\rm therm,\,ion}$ for large
$\dot{m}$, but not always so if $\dot{m}$ is very small, as we discuss below;
see Fig. 7). The effective temperature assuming (for now) local balance of
cooling with $F_{\rm grav}$ and no illumination (we will free this assumption
below) is

$\displaystyle T_{\rm eff}^{\rm local,\,grav-only}\sim 10^{6}\,{\rm
K}\,m_{7}^{-1/4}\,\dot{m}^{1/4}\,x_{g}^{-3/4}.$ (54)

So except for fairly small BHs (especially if we include reprocessing effects
as noted below), we should be in the regime where we can reasonably
approximate the opacity at $T_{\rm eff}$ from the iron opacities plus Kramers,
giving $\kappa_{a}\sim\kappa_{\rm Fe}\sim 0.02\,(\dot{m}\,r_{\rm
ff,\,5}/x_{g})^{1/2}\,\tilde{Z}\,m_{7}^{-1}\,(T/T_{\rm eff}^{\rm local})^{-2}$
which is smaller than $\kappa_{s}$.999Note if we do have small BHs where this
would naively appear to exceed $\kappa_{s}$, this means the central $T_{\rm
eff}^{\rm local}$ becomes large, and at sufficiently small $x_{g}$ we should
use the usual Kramers opacity expression, giving $\kappa_{a}\rightarrow
10^{-5}\,r_{\rm
ff,\,5}^{1/2}\,\dot{m}^{1/8}\,m_{7}^{-5/8}\,x_{g}^{5/8}\,(T/T_{\rm eff}^{\rm
local})^{7/2}$, which changes some of our expressions below but not the
qualitative behaviors of interest. So the effective opacity is
$\kappa^{\ast}\approx\sqrt{\kappa_{a}\,\kappa_{s}}$. We therefore have for the
massive BH regime

$\displaystyle\tau^{\ast}\approx 250\,\dot{m}^{5/4}\,m_{7}^{-5/6}\,r_{\rm
ff,\,5}^{7/12}\,\tilde{Z}^{1/2}\,(T_{\rm eff}^{\rm local}/T)\,x_{g}^{-13/12}$
(55)

(Fig. 7).

The disk is geometrically thinner than at much larger radii but still flaring
and quite geometrically-thick compared to e.g. a thermal-pressure supported
disk, with

$\displaystyle H/R\sim 0.1\,(m_{7}/r_{\rm ff,\,5})^{1/6}\,(x_{g}/3)^{1/6}$
(56)

(Fig. 9; assuming no additional radiation pressure support, which may be
important in the supercritical limit discussed below).

If the disk is radiatively efficient (which should apply at intermediate
$\dot{m}$; see discussion in § 11), then (like any $\alpha$-disk model) we
predict a spin-dependent $\epsilon_{r}\sim 0.1$ as the total luminosity
increases integrating from large radii towards the ISCO, with the majority of
the intrinsic luminosity emitted coming initially from gas within $R\lesssim
R_{\rm therm,\,ion}$. The emission from this zone will give rise to the usual
NIR-optical-UV continuum and big blue bump, with an intrinsic spectrum in good
agreement with that observed in quasars (see e.g. Kishimoto et al. 2008) given
the standard $T_{\rm eff}\propto R^{-3/4}$ relation. However, just like in
observed quasars, a majority of this light is reprocessed and/or scattered by
material at larger radii (Kishimoto et al., 2008; Lawrence, 2012; Chelouche et
al., 2019), which we discuss further below. Given the outermost radii, we
expect the thermalized continuum to transition to the Rayleigh-Jeans tail (of
the temperature at $R_{\rm therm,\,ion}$) at wavelengths

$\displaystyle\lambda\gtrsim\lambda^{\rm char,\,R_{\rm
therm,\,ion}}\sim{\rm\mu m}\,m_{7}^{-0.1}\dot{m}^{0.6}\tilde{Z}^{0.3},$ (57)

quite weakly dependent on AGN properties, except that the spectra are
predicted to be slightly flatter/cooler at higher $\dot{m}$. However at these
wavelengths the emission will be filled in by re-emission from the illuminated
torus (zone (2c), in particular), also as observed (Kishimoto et al., 2008).
In very low accretion rate systems $\dot{m}\lesssim 0.01$, this would predict
the thermal optical/NIR continuum begins to vanish, as the thermalized disk
radius moves inwards; this involves a fundamental state change in the
accretion disk which we discuss below (§ 11).

####  8.1.1 Illumination/Reprocessing and the Disk Effective Temperature

At high $\dot{m}$, as noted above, we expect the thermalized disk to be highly
self-illuminated and interior to an extended scattering atmosphere/corona.
This therefore represents a thermalized “re-processing zone,” with size $\sim
R_{\rm therm,\,ion}^{\rm illum}$ (i.e. size given by the radius interior to
which the scattering optical depth above the disk and effective absorption
optical depth of the disk are both $>1$), as illustrated quantitatively in
Fig. 9. Thus the emergent thermal emission from this region will have an
effective temperature given not by $T_{\rm eff}$ at each interior radius $R$,
but $T_{\rm eff}(R_{\rm phot})$ at the effective photospheric radius of this
region $R\rightarrow R_{\rm phot}\sim R_{\rm therm,\,ion}^{\rm illum}$.
Calculating this for the illuminated limit gives:

$\displaystyle T_{\rm eff}(R_{\rm therm,\,ion}^{\rm illum})\sim 4\times
10^{4}\,{\rm
K}\,\epsilon_{r,\,0.1}^{11/32}\,m_{7}^{-3/64}\,\dot{m}^{-7/32}\,\tilde{Z}^{-3/16}$
(58)

(Fig. 6) – very weakly dependent on any of the properties of the disk, and
much cooler than classic SS73-like $\alpha$-disk prediction which is $T_{\rm
eff}^{\rm SS73}\sim 4\times 10^{5}\,{\rm
K}\,(\dot{m}\,\epsilon_{r,\,0.1}/m_{7})^{1/4}$ (which we also compare in Fig.
6). Crucially, in the latter, the extremely geometrically-thin scale-heights
(which we also show in Fig. 14) mean that there is negligible reprocessing at
larger radii (the radiation flux from the innermost disk directly escapes to
infinity). Note that even if we neglected the details of illumination as
calculated above in shifting this boundary, and simply assumed most of the
emission came from the largest thermalized radii within the disk, we would
obtain similar expressions: $T_{\rm eff}\sim 4\times 10^{4}\,{\rm
K}\,\epsilon_{r,\,0.1}^{1/4}\,m_{7}^{-3/88}\,\dot{m}^{-7/44}\,\tilde{Z}^{-3/22}$
or $T_{\rm eff}^{\rm alt}\sim 5\times 10^{4}\,{\rm
K}\,\epsilon_{r,\,0.1}^{1/4}\,\dot{m}^{-17/52}\,\tilde{Z}^{-3/13}$, given the
scalings in § 7.

This predicted effective temperature agrees quite well with the canonical
values for bright quasars (Peterson, 1997; Vanden Berk et al., 2001; Richards
et al., 2006; Shen et al., 2020; Cai & Wang, 2023) especially after correcting
for dust reddening (Hopkins et al., 2004; Gaskell & Benker, 2007). But it is
not just the value of the effective temperature itself, but the scaling which
is important. It is well-established that the simple prediction of an
SS73-like razor-thin disk model (in which any such reprocessing is impossible)
extending to the ISCO, namely $T_{\rm eff}\propto L^{1/4}\,M_{\rm BH}^{-1/2}$
(as above) is not observed in AGN, with most exhibiting a much weaker, and
oppositely-signed mass dependence, see e.g. Bonning et al. 2007, 2013; Davis
et al. 2007; Hall et al. 2018; Antonucci 2023; Cai & Wang 2023.101010Or,
conversely, fitting AGN to model SS73-like $\alpha$-disks which assume $T_{\rm
eff}\propto L^{1/4}\,M_{\rm BH}^{-1/2}$, the fit appears to prefer
systematically higher $M_{\rm BH}$ and spin values than inferred from any
other methods, see Capellupo et al. 2015. For the flux-frozen disks, because
of the reprocessing effects above, we predict (depending weakly on exactly
which of the above scalings we adopt, and how $\epsilon_{r}$ and $\tilde{Z}$
scale with mass or accretion rate) $T_{\rm eff}\propto L^{-0.2}\,M_{\rm
BH}^{0.1-0.2}$, a much weaker and oppositely-signed scaling (much closer to
the behavior observed).

####  8.1.2 Vertical Structure within the Thermalized Zone

For a laminar ($\mathcal{M}_{s}\ll 1$), stable, hydrostatic, optically-thick,
geometrically thin disk, we would naively expect the vertical heat flux to be
governed by radiative diffusion with diffusive speed

$\displaystyle v_{\rm diff}\approx c/\tau,$ (59)

which establishes a radiation energy density gradient with midplane $T_{\rm
mid}\sim T_{\rm eff}\,\tau^{1/4}$.

However, as discussed in Paper II and Paper III, turbulence also transports
radiation advectively in the vertical direction with speed $v_{\rm
adv,\,z}\approx v_{\rm turb}$, and by definition the supersonic, trans-
Alfvénic turbulence means the eddies reach heights $\approx H$ on turnover
times $\sim 1/\Omega$, which means they carry past the effectively optically-
thick photosphere of the disk (i.e. they can carry trapped radiation to
heights where $\tau^{\ast}\ll 1$), to heights where it can become effectively
optically-thin radiation transport in the outer disk. The criterion for
advection to be faster than diffusion is

$\displaystyle 1<v_{\rm turb}/v_{\rm diff}\sim 52\,\dot{m}\,(r_{\rm
ff,\,5}/m_{7})^{1/6}\,x_{g}^{-7/6},$ (60)

i.e.

$\displaystyle x_{g}\lesssim x_{g}^{\rm turb,z}\approx 30\,(r_{\rm
ff,\,5}/m_{7})^{1/7}\,\dot{m}^{6/7}$ (61)

(noting $r_{\rm ff,\,5}\approx m_{7}^{1/2}$, this is basically $\sim
30\,\dot{m}^{6/7}$). Interior to these radii, the system will become
vertically well-mixed by turbulence. If the medium is sufficiently clumpy or
inhomogeneous, or the turbulent transport behaves advectively, turbulent
mixing will basically set $T_{\rm mid}\sim T_{\rm eff}$ throughout the disk.
If we more conservatively allow for some vertical stratification and assume
the disk remains smooth on average with vertical turbulent transport being
quasi-diffusive with “turbulent diffusivity” $\sim\ell_{\rm turb}\,v_{\rm
turb}$, then solving the vertical steady-state flux equations gives a local
maximum possible midplane temperature relative to effective temperature of
$T_{\rm mid}^{\rm max}/T_{\rm eff}\approx{\rm MIN}[(c/v_{\rm
turb}),\,\tau]^{1/4}\sim{\rm MIN}[2\,(r_{\rm
ff,\,5}/m_{7})^{1/24}\,x_{g}^{1/12},\,\tau^{1/4}]$ with this ratio reaching
its global maximum at $x_{g}\approx x_{g}^{\rm turb,\,z}$ of

$\displaystyle T_{\rm mid}^{\rm max}/T_{\rm eff}\lesssim
2.7\,\dot{m}^{1/14}\,(r_{\rm ff,\,5}/m_{7})^{3/56}.$ (62)

We adopt the latter for calculating the midplane temperatures in this limit in
Fig. 6.

As a result, basically independent of accretion and BH properties, $T_{\rm
mid}$ never exceeds a couple times $T_{\rm eff}$ at the surface (Fig. 6). As
discussed below, these physics mean that even in the highly super-critical
limit, radiation pressure should never inflate the disk by more than a factor
$\sim 1.7$ or so in $H$ (shown quantitatively in the dependence of $H/R$ on
scale and $\dot{m}$ in Figs. 12-13 as well as directly in pressure in Fig. 8),
which does not fundamentally alter any of our conclusions, but could (since
the inflation would be highest at small radii) make the radiation from the
central disk somewhat more isotropic ($f_{\theta}\sim 1$) instead of
vertically collimated ($f_{\theta}\sim 2\,\cos{\theta}$). It also means that
$\beta$ remains $\ll 1$ and the turbulence remains highly supersonic and
compressible (sonic Mach numbers $\sim 50-90$, weakly dependent on disk
parameters; Fig. 10) and likewise the sonic scale $\sim H/\mathcal{M}_{s}^{2}$
remains very small (Fig. 11).

###  8.2 6b: The Warm, Comptonizing Skin

At radii interior to $R\lesssim{\rm a\ few}\times R_{\rm es,\,ion}$ and
$R_{\rm therm,\,ion}$, the disk necessarily will have an electron-scattering
optically-thick but absorption-optically-thin layer or “skin” on top of the
thermalized effectively-optically thick disk (Fig. 2). The depth from the
surface of the disk to thermalization, between which this layer exists, is
approximately (for a constant density disk midplane as we have assumed) $\sim
H/\tau^{\ast}$, which means that the scattering optical depth through the
layer is (Fig. 7)

$\displaystyle\tau_{s}^{\rm skin}$ $\displaystyle\approx\kappa_{s}\,\rho_{\rm
gas}\,H/\tau^{\ast}=\sqrt{\kappa_{s}/\kappa_{a}}$ (63) $\displaystyle\approx
8\,(x_{g}/\dot{m}\,r_{\rm ff,\,5})^{1/4}\,(m_{7}/\tilde{Z})^{1/2}\,(T_{\rm
mid}/2\,T_{\rm eff}^{\rm local})\sim 1-20$

for a wide range of parameters in this zone (recalling from above that within
this zone of the disk $T_{\rm mid}\sim 2\,T_{\rm eff}^{\rm local}$ is a
typical value; Fig. 6 & § 8.1.2).

The gas in this skin is, by definition, effectively optically-thin, so should
cool according to the optically-thin rates. Those cooling rates still give
$t_{\rm cool}\ll t_{\rm dyn}$, however, so the gas temperature should quickly
reach some equilibrium even if (again as noted above) it is being advected up
from the midplane. The heating rate per particle, whether mediated by
viscosity or turbulent dissipation or shocks or magnetic reconnection, is the
necessarily the same in the disk models here up to an $\mathcal{O}(1)$
constant, $\sim v_{\rm turb}^{2}/t_{\rm eddy}\sim v_{\rm turb}^{2}\,\Omega\sim
F_{\rm grav}\,m_{p}/\Sigma_{\rm gas}\sim((B^{2}/8\pi)/\rho)/t_{\rm
reconnect,\,turb}$, and recall from above that these eddies transport said
energy fully into $|z|\gtrsim H$. The gas can cool according to either the
optically-thin Kramers+Fe rates, or via Compton cooling111111Note that despite
the relative importance of magnetic pressure, we can largely neglect
synchrotron/cyclotron cooling. As calculated in detail in Appendix A of Pariev
et al. (2003), the electrons are highly non-relativistic here (except perhaps
in the hottest regions of the “hard” Corona) and so cooling via cyclotron
emission is both inefficient and strongly self-absorbed at the resonant
frequencies where most of the radiation would be emitted. But even if we
assumed ultra-relativistic electrons so that synchroton and Compton cooling
scaled identically with the magnetic and radiation energy densities, in this
particular environment, close to the surface of the disk, if $u_{\rm rad}\sim
a\,T_{\rm eff}^{4}$ (a plausible lower limit) and $u_{\rm B}\sim u_{\rm
B}^{\rm mid}\sim B_{\rm mid}^{2}/8\pi$ (its midplane value) is uniform up to
the upper disk layer, than synchrotron cooling (already much larger than
cyclotron) could only be more efficient than Compton cooling by at most a
factor of a few, which translates to an equilibrium skin temperature
difference of a factor $\lesssim 2$, insufficient to change our conclusions.
from the photons in the disk below. Technically if we include all of these and
simply equate heating and cooling rates, there are two solution branches, the
first branch being the one relevant here121212The second solution branch
involves much lower temperatures $\sim 3\times 10^{4}\,{\rm K}$ dominated by
optically-thin Fe line and H recombination cooling, which is relevant for
cooler gas much further out in the disk or raining back into the disk near the
photosphere, but not of interest here since we are interested in the warm gas
skin heated by shocks or reconnection, where the characteristic post-shock or
reconnection temperature is $\sim(3/16)\,m_{p}\,v_{\rm
turb}^{2}/k_{B}\sim(3/16)\,m_{p}\,v_{A}^{2}/k_{B}\sim 6\times 10^{9}\,{\rm
K}\,(m_{7}/r_{\rm ff,\,5})^{1/3}\,x_{g}^{2/3}$, so the skin will reach the
hotter solution first and remain there in a quasi-stable state. is dominated
by free-free cooling at larger radii ($x_{g}\gtrsim 8\,(m_{7}/r_{\rm
ff,\,5})^{8/7}\,\dot{m}^{-6/7}\,\tilde{Z}^{-12/7}\,(u_{\rm rad}^{\rm
skin}/a\,T_{\rm eff}^{4})$) and Compton cooling at smaller radii. This gives

$\displaystyle T_{\rm skin}\approx 10^{7}\,{\rm K}\,(m_{7}/r_{\rm
ff,\,5})^{5/3}\,\dot{m}^{-2}\,\tilde{Z}^{-2}\,x_{g}^{-1/3}$ (64)

and

$\displaystyle T_{\rm skin}\approx 10^{6}\,{\rm K}\,(m_{7}/r_{\rm
ff,\,5})^{1/3}\,\dot{m}^{-1}\,x_{g}^{5/6}\,(u_{\rm rad}^{\rm skin}/a\,T_{\rm
eff}^{4}),$ (65)

respectively, with a maximum temperature of

$\displaystyle T_{\rm skin}^{\rm max}\sim 5\times 10^{6}\,{\rm
K}\,\left(\frac{m_{7}}{r_{\rm
ff,\,5}}\right)^{9/7}\tilde{Z}^{-10/7}\dot{m}^{-12/7}\left(\frac{a\,T_{\rm
eff}^{4}}{u_{\rm rad}^{\rm skin}}\right)^{2/7}$ (66)

and minimum (at the outermost radii of the thermalized zone) of

$\displaystyle T_{\rm skin}^{\rm min}\sim{\rm MIN}[$ $\displaystyle 2\times
10^{6}\,{\rm K}\,(m_{7}/r_{\rm
ff,\,5})^{1.9}\,\tilde{Z}^{-2.2}\,\dot{m}^{-2.4}\ ,$ $\displaystyle\ \ \
7\times 10^{7}\,{\rm K}\,m_{7}^{-1/4}\,\tilde{Z}^{5/13}]$ (67)

(shown in Fig. 6).

Given this characteristic temperature, we can estimate the Compton $y$
parameter of the skin as $y=(k_{B}\,T_{e}^{\rm
skin}/m_{e}\,c^{2})\,(\tau_{s}^{\rm skin})^{2}$ giving

$\displaystyle y\sim$
$\displaystyle\,\frac{0.21}{\tilde{Z}}\,\left(\frac{T_{e,\,{\rm
skin}}}{2\times 10^{6}\,{\rm K}}\right)\left(\frac{T_{\rm mid}}{2\,T_{\rm
eff}^{\rm local}}\right)^{2}\left(\frac{x_{g}}{100\,\dot{m}\,r_{\rm
ff,\,5}}\right)^{1/2}$

At small radii approaching the ISCO, the skin will not Comptonize the gas, but
at most radii covering the thermal gas disk out to $\sim 200\,R_{g}$, this is
appreciable.

Comparing to the electron temperatures ($k_{B}\,T_{e,\,{\rm skin}}$; Hubeny et
al. 2001; Czerny et al. 2003); optical depths $\tau_{s}^{\rm skin}$ and
Compton $y$ parameter (Czerny et al., 2003; Kubota & Done, 2018; Petrucci et
al., 2018); covering factors and geometry (Wilkins & Gallo, 2015; Petrucci et
al., 2018) observationally inferred from AGN spectra – most notably required
to explain the soft excess – the predicted properties all appear consistent.
Note that while the radial extent of the Comptonizing region is predicted to
expand with $L$ or $\dot{m}$ as does the disk itself (Fig. 12), as implied by
observations (Palit et al., 2024), the temperatures above become lower and
therefore the excess becomes softer at higher $\dot{m}$, again consistent with
observed behaviors (Done et al., 2012; Krawczyk et al., 2013; Stevans et al.,
2014; Boissay et al., 2016; Ballantyne et al., 2024). Thus a Comptonized soft
excess similar to that observed seems inevitable in these flux-frozen disks.

####  8.2.1 6c: Towards the ISCO and Near-Horizon Region

At sufficiently small radii $x_{g}$ approaching and interior to the ISCO (Zone
(6c) in Fig. 1), the solutions above will need to be modified to map to the
different orbit structure and “plunge” or rapid infall, as well as including
the appropriate non-Newtonian terms in the potential and angular momentum.
However, at such radii, the detailed results will depend on properties like
the BH spin, and almost certainly on components like a jet and the detailed
behavior of the magnetic fields in the near-horizon regime which we are not
modeling here (though we discuss jets briefly in § 12). Moreover as shown in
Paper I, the accretion disks are often misaligned with the BH spin: this can
introduce more radical changes like strong precession, warps, or even disk
“tearing,” none of which we model here (see e.g. Kaaz et al., 2022). Clearly,
detailed modeling of this regime requires GRMHD simulations to extend the
Newtonian simulations in Paper I and Paper II.

However, to gain some insight, we can attempt to simply estimate how our
solutions should be modified assuming an inner free-fall boundary at the ISCO
for a simple quasi-Newtonian approximation around a non-rotating BH (assuming
there is no warp or precession in the disk), as done in most $\alpha$-disk
models. There is no conceptual difficulty here and it does not contradict any
of our important assumptions from Paper III. For example, retaining the
foundational two assumptions from § 2.1, we can replace the potential with
Paczynski-Wiita (PW) which multiplies the orbital frequency $\Omega$ by
$\tilde{\Omega}\equiv{x_{g}/(x_{g}-1)}$, and multiply the angular momentum
term by $\tilde{g}(x)\equiv 1-(x_{0}/x_{g})^{1/2}$ (where $x_{0}=3$ for the
ISCO of a non-rotating BH, as in Shakura & Sunyaev 1973). In principle, we can
represent any quasi-Newtonian potential (e.g. the galactic potential with
contribution from stars and dark matter in Zone (1)) and zero angular momentum
disk “cutoff” via appropriate $\tilde{\Omega}$ and $\tilde{g}$, respectively,
with such an approach (though of course this cannot capture various
relativistic effects). In fact, we use this in all of our more detailed
calculations in Figs. 3-11.

Within such a simple approximation, this gives us the exact same global disk
solutions from § 2 except we modify:

$\displaystyle v_{t}$ $\displaystyle\sim v_{A}\rightarrow v_{A}^{\rm
Keplerian}\,\tilde{\Omega}^{2/9}\,\tilde{g}^{1/9},$ (68)
$\displaystyle\Sigma_{\rm gas}$ $\displaystyle\rightarrow\Sigma_{\rm gas}^{\rm
Keplerian}\,\tilde{\Omega}^{5/9}\,\tilde{g}^{7/9}.$ (69)

Or, for the modified (advection-dominated) solutions in § 11.1.3, $v_{t}\sim
v_{A}\propto v_{\rm c}\propto\tilde{\Omega}$ and $\Sigma_{\rm
gas}\propto\tilde{\Omega}^{-1}\,\tilde{g}$ At large $x_{g}$ these reduce to
the expressions we have used throughout. At small $x_{g}\lesssim 10-20$, these
introduce modest quantitative changes to the solutions, most notably,

$\displaystyle T_{\rm eff}^{\rm local}\rightarrow T_{\rm eff}^{\rm
local,\,Keplerian}\,\tilde{\Omega}^{1/2}\,\tilde{g}^{1/4}.$ (70)

(again, included in Fig. 6 and those like it for both our models and SS73
comparisons), and the optical depths change by a small amount. None of this
changes our conclusions above about the disk being optically-thick, having
$\kappa_{s}>\kappa_{a}$, being strongly turbulent, etc.131313If anything they
tend to slightly strengthen some of the above conclusions. For example, the
upper limit – assuming optically-thick diffusive transport only – of the ratio
of midplane radiation pressure to magnetic+turbulent pressure scales as
$\propto\tilde{\Omega}^{7/9}\,\tilde{g}^{8/9}<1$ (a factor of $\sim 0.5$ at
$x_{g}\approx 10$, decreasing to $0$ as $x_{g}\rightarrow 3$). Similarly, the
disk flaring and vertical structure is nearly identical, except that it
becomes slightly more strongly flared ($H/R$ decreases with smaller $R$ more
rapidly) in the very narrow range $3\lesssim x_{g}\lesssim 5$. Of course, very
close to $x_{g}\rightarrow 3$, these correction terms strongly modify the
solution (formally extrapolating to a vacuum; Figs. 3-4), but this occurs so
close to $x_{g}\approx 3$ that it is, for most purposes, akin to simply
truncating the disk. This is precisely the regime where our SS73-like
assumption of zero angular momentum at the ISCO is suspect and requires more
detailed fully-GRMHD models to calibrate (see discussion in Abramowicz &
Fragile 2013).

In practice, the most important consequence of these correction terms for our
results here, like with a standard thermal-pressure dominated $\alpha$-disk,
is that they significantly suppress both the relative amount of direct
radiation flux coming from small $x_{g}$ and its effective temperatures, so
the thermal emitted continuum will be softer than if we ignored the term
$\tilde{g}(x)$.

##  9 Zone 7: The “Hard” Corona

The Comptonizing skin (zone (6b); § 8.2) is made up of gas within one disk
scale height $H\,(1-1/\tau^{\ast})<|z|<H$ above the midplane of the inner,
thermalized disk (see Figs. 1, 2, 9). Now consider the gas at somewhat larger
$|z|\gtrsim H$, at similar and smaller radii (Figs. 12 & 13). The vertical
profiles of these disks at $|z|>H$ are discussed in detail in Paper I: the gas
density (and rms vertical turbulent velocity) follows an approximate
$\rho\propto{\rm sech}^{2}(z/H)$ or Gaussian “core” profile (at $|z|\lesssim
H$) with an extended power-law falloff at larger radii, while the magnetic
field strength $|{\bf B}|$ falls of less rapidly than $\rho$ (especially at
$|z|\gtrsim H$), and the temperatures rise in the more rarified optically-thin
coronal gas above the disk. The behavior of the vertical profiles – most
notably, that $|{\bf B}|$ falls off slowly compared to $\rho$ – is expected
from basic vertical equilibrium considerations given our ansatz. Consider:
equating the (primarily toroidal) field support for the gas at $|z|\gtrsim H$,
$(\nabla\times{\bf B})\times{\bf B}/4\pi\sim B_{\phi}^{2}/z$, and gravity,
$-\rho\nabla\Phi_{\rm grav}=\rho\Omega^{2}z$, gives
$B_{\phi}(r,z)\propto\rho(r,z)^{1/2}z$ at a given $r$. Provided this
qualitative behavior holds, we find the details of the vertical profile are
not so important, so consider a generalized $\rho(z)=\tilde{\rho}(|z|\gtrsim
H)\,\rho_{\rm mid}$, $|{\bf B}(z)|=\tilde{B}(|z|\gtrsim H)\,|{\bf B}|_{\rm
mid}$ (with $\tilde{B}\sim\tilde{\rho}^{1/2}\,(|z|/H)$). Assume (since we are
not very far above the disk) the gas is illuminated primarily by the flux from
the disk surface, and define its electron temperature $T_{e}^{\rm cor}$ and
ion temperature $T_{i}^{\rm cor}$ (shown in Fig. 6). In steady-state, this gas
should be heated by magnetic reconnection (or shocks or turbulent dissipation,
which give a similar rate since $v_{\rm turb}\sim v_{\rm shock}\sim
v_{A}>c_{s}$) with time-averaged heating flux $F_{\rm heat}^{\rm cor}\sim|{\bf
B}|^{2}\,v_{A}(z)/(8\pi)$ and characteristic post-shock/reconnection
temperatures

$\displaystyle T_{\rm shock}^{\rm max}\sim T_{\rm recon}^{\rm
max}\sim\frac{3m_{p}}{16k_{B}}v_{\rm
shock}^{2}\sim\frac{3m_{p}}{16k_{B}}(2\,\alpha_{2}\,v_{A})^{2},$ (71)

of which a fraction $\delta\approx 0.5$ and $(1-\delta)$ go into the electrons
and ions respectively (Bisnovatyi-Kogan & Lovelace, 1997, 2000). Gas at these
temperatures is cooled by a combination of Compton, cyclotron/synchrotron (if
the electrons become highly relativistic), and Kramers (Bremsstrahlung, at the
temperatures of interest) cooling of the electrons $F_{\rm cool}^{\rm cor}\sim
F_{\rm Compton}^{\rm cor}+F_{\rm synch}^{\rm cor}+F_{\rm Brems}^{\rm cor}$. In
the inner regions, since the radiation flux is strong and falls off less-
rapidly than other quantities above with height, Compton cooling will be most
important.

As we move above the disk, since $\tilde{\rho}$ decreases more rapidly than
$\tilde{B}$, “hard” X-ray coronal solutions can appear at small enough $R$ and
$H\lesssim|z|\lesssim r$. Specifically, at a given sufficiently low
$\tilde{\rho}$, the timescale for full ion-electron temperature equilibration
via Coulomb interactions $t_{\rm ie}^{\rm full}$ (Spitzer, 1962; Narayan et
al., 1998; Cao, 2009; Faucher-Giguère & Quataert, 2012) at the characteristic
post-shock/reconnection temperature $T_{\rm shock}^{\rm max}$ becomes longer
than the characteristic shock crossing/reconnection timescale $t_{\rm
cross}\sim|z|/v_{\rm char}(z)\sim|z|/v_{A}(z)\sim\Omega^{-1}\sim t_{\rm dyn}$.
If that occurs, the ions will be heated faster than they can cool and will
therefore virialize at

$\displaystyle T_{i}^{\rm cor}\sim T_{\rm vir}\sim\frac{1}{2}\frac{GM_{\rm
BH}}{r}\,\frac{\mu\,m_{p}}{k_{B}},$ (72)

at which point the scale-height of the hard corona is $\sim r$, with $|{\bf
B}|$ and $T_{e}$ sub-virial (so the gas remains bound with $c_{s}$ and $v_{A}$
modestly sub-virial). If the gas density in this region is initially more
tenuous, it will build up rapidly with inflow (since this virialized component
is not rapidly accreting) until reaching a quasi-steady-state density where
continued heating is balanced by cooling for electrons and ions (at $T_{\rm
mid}\ll T_{e}^{\rm cor}\leq T_{i}^{\rm cor}$).

We can then calculate the approximate coronal properties following Cao (2009)
and equating the Coulomb, Compton, and heating rates (expressions given
therein) for ions and electrons in steady-state, assuming a virialized ion
system with $H/R\sim 1$ and electron-scattering optical depths in the corona
of order unity. In the inner regions at sufficiently low $\dot{m}$, self-
consistent coronal solutions can extend to $\tilde{B}\sim 1$, i.e. to the
scale-height $|z|\sim H$, with ion

$\displaystyle T_{i}^{\rm cor}\sim T_{\rm vir}\sim 10^{12}\,{\rm K}/x_{g},$
(73)

electron

$\displaystyle T_{e}^{\rm cor}\sim 10^{9}\,{\rm K}\,(r_{\rm
ff,\,5}/m_{7})^{1/19}\,(x_{g}/\dot{m})^{4/19}\sim 50-100\,{\rm keV},$ (74)

saturated density

$\displaystyle\rho_{\rm cor}^{\rm sat}/\rho_{\rm mid}\sim 0.001\,(m_{7}/r_{\rm
ff,\,5})^{7/19}\,(x_{g}/\dot{m})^{10/19},$ (75)

and approximate Thompson optical depth

$\displaystyle\tau_{\rm es}^{\rm cor}\sim 13\,(r_{\rm
ff,\,5}/m_{7})^{5/38}\,(\dot{m}/x_{g})^{9/19}\sim 4\,(x_{g}/10\dot{m})^{-0.5}$
(76)

and Compton

$\displaystyle y\sim 10\,m_{7}^{-0.2}\,(x_{g}/10\dot{m})^{-0.8}.$ (77)

This regime is valid approximately between radii $10\,\dot{m}\lesssim
x_{g}\lesssim x_{g}^{\rm cor,\,H}$. At smaller radii $x_{g}\lesssim
10\,\dot{m}$ (only relevant at higher $\dot{m}$), radiation pressure and
optical depth effects become important as expected (see discussion below, §
11). At larger radii and/or lower heights above the midplane at large $x_{g}$
and $\dot{m}$, more specifically where

$\displaystyle x_{g}\gtrsim x_{g}^{\rm cor,\,H}\sim 24\,(m_{7}/r_{\rm
ff,\,5})^{30/33}\,\alpha_{2}^{130/33}\,\dot{m}^{-6/11}\,(|z|/H)^{20/33},$ (78)

coronal solutions cannot extend to the disk scale height $|z|\sim H$, or
equivalently must have $\tilde{B}\ll 1$ (or else e.g. the implied Alfvén
speeds become super-virial and no self-consistent solution is possible). In
this regime the coronal gas will reside at larger vertical heights $|z|$ above
the disk. In this regime, from $x_{g}^{\rm cor,\,H}\lesssim x_{g}\lesssim
x_{g}^{\rm cor,\,z}$, using the virial limit for $\tilde{B}$ and
$\tilde{\rho}$, we obtain the modified relations for the corona at
sufficiently large height:

$\displaystyle T_{e}^{\rm cor}$ $\displaystyle\sim 0.4\times 10^{9}\,{\rm
K}\,(x_{g}/10)^{1/2}\,\dot{m}^{-1},$ (79) $\displaystyle\tilde{\rho}$
$\displaystyle\sim 0.0002\,(m_{7}/r_{\rm
ff,\,5})^{1/2}\,(x_{g}/10)^{5/4}\,\dot{m}^{-5/2},$ (80)
$\displaystyle\tau_{\rm es}^{\rm cor}$ $\displaystyle\sim
0.3\,(x_{g}/10)^{1/4}\,\dot{m}^{-3/2},$ (81) $\displaystyle y$
$\displaystyle=4\,k_{B}\,T_{e}\,\tau_{\rm es}/(m_{e}\,c^{2})\sim
0.07\,(x_{g}/10)^{3/4}\,\dot{m}^{-5/2}.$ (82)

In Fig. 6 we compare these coronal electron and ion temperature profiles to
the midplane temperatures at different radii, BH mass, and accretion rates.
Note that both the electron temperature (hardness of the corona) and Compton
$y$ depend inversely on $\dot{m}$.