# GRRMHD Simulations of MAD Accretion Disks Declining from Super-Eddington to
Sub-Eddington Accretion Rates

Brandon Curd,1,2 Ramesh Narayan1,2
1 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA
02138, USA
2 Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA
02138, USA

E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present two general relativistic radiation magnetohydrodynamics (GRRMHD)
simulations of magnetically arrested disks (MADs) around non-spinning
($a_{*}=0$) and spinning ($a_{*}=0.9$) supermassive black holes (BHs). In each
simulation, the mass accretion rate is decreased with time such that we sample
Eddington-scaled rates over the range
$3\gtrsim\dot{M}/\dot{M}_{\rm{Edd}}\gtrsim 0.3$. For the non-spinning BH
model, the total and radiative efficiencies increase as the accretion rate
decreases, varying over the range $\eta_{\rm{tot}}\sim 9-16\%$ and
$\eta_{\rm{rad}}\sim 6-12\%$, respectively. This model shows very little jet
activity. In contrast, the spinning BH model has a strong relativistic jet
powered by spin energy extracted from the BH. The jet power declines with
accretion rate such that $\eta_{\rm{jet}}\sim 18-39\%$ while the total and
radiative efficiencies are $\eta_{\rm{tot}}\sim 64-100\%$ and
$\eta_{\rm{rad}}\sim 45-79\%$, respectively. We confirm that mildly sub-
Eddington disks can extract substantial power from a spinning BH, provided
they are in the MAD state. The jet profile out to $100\,GM/c^{2}$ is roughly
parabolic with a power-law index of $k\approx 0.43-0.53$ during the sub-
Eddington evolution. Both models show significant variability in the outgoing
radiation which is likely associated with episodes of magnetic flux eruptions.
The $a_{*}=0.9$ model shows semi-regular variations with a period of $\sim
2000\,GM/c^{3}$ over the final $\sim 10,000\,GM/c^{3}$ of the simulation,
which suggests that magnetic flux eruptions may be an important source of
quasi-periodic variability. For the simulated accretion rates, the $a_{*}=0$
model is spinning up while the $a_{*}=0.9$ model is spinning down. Spinup-
spindown equilibrium of the BH will likely be achieved at
$0.5<a_{*,{\rm{eq}}}<0.6$, assuming continuous accretion in the MAD state.

###### keywords:

accretion, accretion discs - black hole physics - MHD - radiative transfer -
gamma-rays: galaxies - X-rays: galaxies

††pubyear: 2022††pagerange: GRRMHD Simulations of MAD Accretion Disks
Declining from Super-Eddington to Sub-Eddington Accretion Rates–B

##  1 Introduction

Accretion disks around black holes (BHs) provide powerful engines for
converting gravitational energy and BH spin energy into radiation emitted from
the viscously heated inner accretion flow or from a powerful jet. In essence,
BHs are engines that convert the accreting mass into other forms of energy
(radiation, kinetic energy, magnetic power). The efficiency of this energy
conversion depends intimately on the accretion disk properties, with one of
the key components determining the energetics being the mass accretion rate
itself (e.g see Abramowicz & Fragile 2013; Yuan & Narayan 2014). Various
models of accretion disks on the low and high end of accretion rates have
found success in applications to various astrophysical phenomena involving
BHs. When disks are accreting at super-Eddington rates, models such as
geometrically thick, radiatively inefficient slim disks may apply (Abramowicz
et al., 1988). For accretion rates in the range $\sim 10^{-2}-0.5$ times the
Eddington rate, a geometrically thin, radiatively efficient disk, formally
known as a Novikov-Thorne (NT hereafter) thin disk (Novikov & Thorne, 1973;
Shakura & Sunyaev, 1973) is likely present. Meanwhile, when the accretion rate
drops yet further, the disk evaporates and forms a two-temperature advection
dominated accretion flow, which is again radiatively inefficient and
geometrically thick (Narayan & Yi, 1995; Esin et al., 1997; Shapiro et al.,
1976; Ichimaru, 1977; Rees et al., 1982).

How exactly accretion disks evolve from one stage to another and how similarly
the transitions operate across different BH mass scales and astrophysical
settings is not well understood. Systems where the same BH experiences large
variations in the mass accretion rate $\dot{M}$ as a function of time are
particularly valuable in this regard. Accretion flows corresponding to
stellar-mass BH X-ray binaries (BHXBs) and tidal disruption events (TDEs)
around supermassive BHs (SMBHs), for example, show large changes in $\dot{M}$,
and likely undergo state transitions. In both cases, for some fraction of the
time the accretion rate is modestly sub-Eddington and the accretion disk is
likely in the NT thin disk state.

Many AGN launch powerful astrophysical jets. One of the most attractive models
for powering these jets is the Blandford-Znajek (BZ hereafter, Blandford &
Znajek 1977; Penrose 1969), which requires that the BH has a relatively high
spin and the BH is threaded by a large scale poloidal field. However, how
exactly BHs in AGN both acquire and maintain this magnetic flux is not clear.
In particular, poloidal fields in thin accretion disks are expected to rapidly
diffuse away radially, indicating that thin accretion disks may not power BZ
jets (Lubow et al., 1994; Lovelace et al., 1994; Esin et al., 1997; van
Ballegooijen, 1989). Additional support for this theoretical claim is provided
by the jetted TDEs Swift J1644+57 and Swift J2058+05, whose X-ray flux showed
a drastic decrease of more than two orders of magnitude several months after
the initial detection, falling below detection limits (Zauderer et al., 2013;
Pasham et al., 2015). Assuming both TDE events are the result of a SMBH
disrupting a main sequence star, this rapid decline is thought to have
occurred due to some critical accretion rate having been reached where the jet
shuts off. At the time of drastic X-ray decline, both systems were likely sub-
Eddington, which suggests that the disk may have indeed become geometrically
thin (Tchekhovskoy et al., 2014; Pasham et al., 2015).

How do jetted TDEs (and similar BH accretion systems) transition to a
quiescent state where no apparent jet is observed? In terms of geometrically
thick accretion disks, the magnetic flux and BH spin provide an explanation.
In simple terms, accretion disks threaded by weak magnetic fields accrete via
turbulence fueled by the magneto-rotational instability (MRI, Balbus & Hawley
1991) and are expected to produce only mildly relativistic outflows at best
(Sądowski & Narayan, 2015a; Curd & Narayan, 2019). Meanwhile, accretion disks
threaded by powerful poloidal fields are capable of launching powerful jets
(Curd & Narayan, 2019), presumably through the BZ process. It is common to
refer to disks with weak fields as “Standard and Normal Evolution” (or SANE,
Narayan et al. 2012) disks and those with powerful poloidal fields as
“Magnetically Arrested Disks” (or MAD, see Narayan et al. 2003; Igumenshchev
et al. 2003).

Although MADs provide an avenue for jet launching, theoretical studies, as
discussed earlier, suggests that geometrically thin accretion disks cannot
sustain a powerful jet. This fact led Tchekhovskoy et al. (2014) to conclude
that the disk in Swift J1644+57 must have become geometrically thin and
simultaneously shed its magnetic flux and become non-MAD leading up to the
X-ray shutoff. However, the assumption that the magnetic field diffuses away
from the BH horizon once the disk becomes geometrically thin requires some
scrutiny. Several works which simulated sub-Eddington accretion disks
(Sądowski, 2016; Avara et al., 2016; Morales Teixeira et al., 2018; Liska et
al., 2022) find that disks with $H/R\sim 0.1$ can be stabilized against
thermal collapse via powerful magnetic fields. Simulations of sub-Eddington
MADs (Avara et al., 2016; Morales Teixeira et al., 2018; Liska et al., 2022)
did not show significant loss of magnetic flux over the simulated time. If the
disk in Swift J1644+57 was indeed MAD, and was stabilized against collapse
when it went sub-Eddington by strong magnetic fields, then why did the X-rays
decline so suddenly, signifying that the jet shut off?

Models of AGN often assume a thin disk and the presence of some corona, which
produces significant X-ray emission through both direct emission and photons
which scatter off of the accretion disk. One possible source of the corona is
the jet/wind launched by the accretion flow. In particular, for spinning BHs,
the BZ process extracts power from the BH (Tchekhovskoy et al., 2010) at a
rate:

$P_{\rm{BZ}}=\dfrac{\kappa}{16\pi
c}\omega_{H}^{2}\Phi_{\rm{BH}}^{2}f(\omega_{H}).$ (1)

Here $\Phi_{\rm{BH}}$ is the magnetic flux threading the BH horizon, $\kappa$
is a constant related to the field geometry and is $0.044$ for a parabolic
field structure, $\omega_{H}=a_{*}/[1+(1-a_{*}^{2})^{1/2}]$ is the
dimensionless rotational frequency of the BH horizon, and
$f(\omega_{H})\approx 1$ for the spins $a_{*}=0,~{}0.9$ that we consider in
this work. Normalizing by the mass accretion rate $\dot{M}$, the efficiency
is:

$\eta_{\rm{BZ}}\equiv\dfrac{\langle P_{\rm{BZ}}\rangle}{\langle\dot{M}\rangle
c^{2}}=\dfrac{\kappa}{16\pi c}\omega_{H}^{2}\langle\phi_{\rm{BH}}^{2}\rangle
f(\omega_{H}),$ (2)

where the normalized flux $\phi_{\rm BH}$ is defined in equation (15). Thus
one expects a jet to be present so long as the BH has some rotation and the
disk has sufficient magnetic flux $\Phi_{\rm{BH}}$. This jet will likely
source X-ray emission.

Scaling relations from simulations suggests that thin disks maintain lower
magnetic flux. Geometrically thick disks reach a magnetically arrested state
once the normalized magnetic flux at the BH horizon (see Section 2.1 for
definitions) is $\phi_{\rm{BH}}\approx 70(1-0.38\omega_{H})[(H/R)/0.3]$
(Tchekhovskoy et al., 2011, 2012). Here $H/R$ is the disk scale height. If a
similar scaling applies as the accretion rate $\dot{M}$ and the disk scale
height $H/R$ decline, even sub-Eddington MAD disks may maintain substantial
magnetic flux provided $H/R$ does not become too small. The powerful magnetic
field in this scenario provides pressure support and also leads to a jet (and
thus X-ray emission) for spinning BHs. With these considerations in mind, the
rapid decline in the X-ray of Swift 1644+57 is quite puzzling when interpreted
as the result of a MAD accretion disk around a SMBH.

To further explore the viability of the MAD state in geometrically thin disks
and to understand its relation to relativistic jets, in this work we study the
properties of moderately thin, MAD accretion disks where the mass accretion
rate is steadily decreased from super-Eddington rates to sub-Eddington rates.
The systems described in this work are moderately thin ($H/R\approx 0.15$),
MAD accretion disks around a SMBH of mass $10^{6}\,M_{\odot}$ where the
accretion rate peaks at $3-6\,\dot{M}_{\rm{Edd}}$ and decreases to $\approx
0.3\,\dot{M}_{\rm{Edd}}$ by the end of the simulation. We perform one
simulation with BH spin $a_{*}=0$ and one with $a_{*}=0.9$ to study the
effects of BH spin. We study the efficiency, variability, jet structure, and
emission as the disk evolves.

We describe our numerical methods and choice of initial conditions in Section
2. We also describe various diagnostics used to study each model. In Section
3, we describe the time evolution of each model. We then analyze the outflow
properties and measure the jet profile. We also perform a cross-correlation
analysis to quantify the source of variability in the escaping radiation. In
Section 4, we quantify the spin evolution of each model as a function of
accretion rate. In Section 5, we compare our results with related work in the
literature, and also consider observational implications. We summarize our
findings in Section 6.

##  2 Numerical Methods

The simulations presented in this work were performed using the general
relativistic radiation magnetohydrodynamical (GRRMHD) code `KORAL` (Sądowski
et al., 2013, 2014; Sądowski et al., 2017; McKinney et al., 2014) which solves
the conservation equations in a fixed, arbitrary spacetime using finite-
difference methods. We solve the following conservation equations:

$\displaystyle(\rho u^{\mu})_{;\mu}$ $\displaystyle=0,$ (3)
$\displaystyle(T^{\mu}_{\nu})_{;\mu}$ $\displaystyle=G_{\nu},$ (4)
$\displaystyle(R^{\mu}_{\nu})_{;\mu}$ $\displaystyle=-G_{\nu},$ (5)

where $\rho$ is the gas density in the comoving fluid frame, $u^{\mu}$ are the
components of the gas four-velocity as measured in the “lab frame”,
$T^{\mu}_{\nu}$ is the MHD stress-energy tensor in the “lab frame”,

$T^{\mu}_{\nu}=(\rho+u_{g}+p_{g}+b^{2})u^{\mu}u_{\nu}+\left(p_{g}+\dfrac{1}{2}b^{2}\right)\delta^{\mu}_{\nu}-b^{\mu}b_{\nu},$
(6)

$R^{\mu}_{\nu}$ is the stress-energy tensor of radiation, and $G_{\nu}$ is the
radiative four-force which describes the interaction between gas and radiation
(Sądowski et al., 2013, 2014). Here $u_{g}$ and $p_{g}=(\gamma-1)u_{g}$ are
the internal energy and pressure of the gas in the comoving frame and
$b^{\mu}$ is the magnetic field four-vector which is evolved following the
ideal MHD induction equation (Gammie et al., 2003).

The radiative stress-energy tensor is obtained from the evolved radiative
primitives, i.e. the radiative rest-frame energy density and the four velocity
of this frame, using the M1 closure scheme modified by the addition of
radiative viscosity (Sądowski et al., 2013; Sądowski et al., 2015).

The interaction between gas and radiation is described by the radiation four-
force $G_{\nu}$. The opposite signs of this quantity in the conservation
equations for gas and radiation stress-energy (equations 4, 5) reflect the
fact that the gas-radiation interaction is conservative, i.e. energy and
momentum are transferred between gas and radiation. For a detailed description
of the four-force see Sądowski et al. (2017). We include the effects of
scattering via the electron scattering opacity ($\kappa_{\rm{es}}$). We also
account for absorption and emission via the free-free absorption opacity
($\kappa_{\rm{a}}$) and bound-free absorption opacity through the Sutherland
Dopita model (Sutherland & Dopita, 1993). We assume a Solar metal abundance
for the gas such that the Hydrogen mass-fraction $X=X_{\odot}=0.7381$, the
Helium mass-fraction $Y=Y_{\odot}=0.2485$, and the mass-fraction of all other
elements $Z=Z_{\odot}=0.0134$. We also include the effects of thermal
synchrotron and Comptonization (Sądowski & Narayan, 2015b; Sądowski et al.,
2017).

We quantify the resolution of the fastest growing mode of the MRI by computing
the quantities:

$\displaystyle Q_{\theta}=\dfrac{2\pi}{\Omega
dx^{\theta}}\dfrac{|b^{\theta}|}{\sqrt{4\pi\rho}},$ (7) $\displaystyle
Q_{\phi}=\dfrac{2\pi}{\Omega dx^{\phi}}\dfrac{|b^{\phi}|}{\sqrt{4\pi\rho}},$
(8)

where $dx^{i}$ (the size of the grid cell) and $b^{i}$ (the magnetic field)
are both evaluated in the orthonormal frame, $\Omega$ is the angular velocity,
and $\rho$ is the gas density. Numerical studies of the MRI have shown that
values of $Q_{\theta}$ and $Q_{\phi}$ in excess of at least 10 are needed to
resolve the fastest growing mode (Hawley et al., 2011).

We use modified Kerr-Schild coordinates with the metric corresponding to a BH
with mass $10^{6}\,M_{\odot}$ and two values of spin, $a_{*}=0,\,0.9$. In each
simulation, we run the initial inflow of gas in 2D ($r,\theta$) coordinates
with a resolution $N_{r}\times N_{\theta}=256\times 128$. We then pause the
simulation and perform a regrid (which we later refer to as ‘regrid 1’) to 3D
($r,\theta,\phi$) with a resolution of $N_{r}\times N_{\theta}\times
N_{\phi}=256\times 128\times 128$ and evolve the simulation until the
accretion rate approaches sub-Eddington levels. At this time, we again pause
the simulation and perform an additional regrid (which we later refer to as
‘regrid 2’) in which we double the resolution in $\theta$ with a resolution of
$N_{r}\times N_{\theta}\times N_{\phi}=256\times 256\times 128$ and evolve the
simulation until the accretion rate approaches $~{}\sim 0.3$ times the
Eddington limit. We set the inner edge of the domain such that the inner six
domain cells are inside the BH horizon, which guarantees that the interior
grid boundary is causally disconnected from the exterior. The radial grid
cells are spaced logarithmically in radius and the cells in polar angle
$\theta$ have smaller sizes $\Delta\theta$ towards the equatorial plane (see
Tchekhovskoy et al. 2011). The cells are equally spaced in azimuth. At the
inner radial boundary ($R_{\rm{min}}$), we use an outflow condition while at
the outer boundary ($R_{\rm{max}}$) we use a similar boundary condition and in
addition prevent the inflow of gas and radiation. At the polar boundaries, we
use a reflective boundary. The polar boundaries are radius dependent such that
near the horizon the minimum/maximum polar coordinate is shifted slightly away
from the pole to reduce the minimum time step, while at larger radii the grid
extends closer to the poles (Tchekhovskoy et al., 2011). We employ a periodic
boundary condition in azimuth. For full 3D portions of each simulation, the
grid covers $-\pi\leq\phi\leq\pi$.

Model Name | $a_{*}$ | $R_{\rm{in}}$ | $R_{\rm{max}}$ | $R_{\rm{inj}}$ | Resolution | Duration
---|---|---|---|---|---|---
|  | $(r_{g})$ | $(r_{g})$ | $(r_{g})$ | $N_{r}\times N_{\theta}\times N_{\phi}$ | ($t_{g}$)
KIa00 | 0 | $1.5$ | $10^{5}$ | $40$ | $256\times 256\times 128$ | 29,000
KIa09 | 0.9 | $1.1$ | $10^{5}$ | $40$ | $256\times 256\times 128$ | 34,000
Table 1: Description of simulations presented in this work. Note that we give
the final resolution of each simulation. Both simulations were performed in
2D, and 3D with multiple resolution increases. We describe this in detail in
Section 2.2.

We performed two simulations, one with zero spin $a_{*}=0$ (KIa00) and one
with near maximal spin $a_{*}=0.9$ (KIa09). In this naming convention ‘KI’
identifies these models as Keplerian disk injection models (as described in
Section 2.2, the gas in these simulations is injected at a radius
$R_{\rm{inj}}$) and ‘a00’ or ‘a09’ identifies the BH spin. Details of each
simulation are given in Table 1.

###  2.1 Definitions

Throughout this work, we use gravitational units to describe physical
parameters. For distance we use the gravitational radius $r_{g}\equiv
GM_{\rm{BH}}/c^{2}$ and for time we use the gravitational time $t_{g}\equiv
GM_{\rm{BH}}/c^{3}$, where $M_{\rm{BH}}$ is the mass of the BH. Often, we set
$G=c=1$, so the above relations would be equivalent to $r_{g}=t_{g}=M_{\rm
BH}$111For a BH mass of $10^{6}\,M_{\odot}$, the gravitational radius and time
in cgs units are $r_{g}=1.48\times 10^{11}$ cm and $t_{g}=4.94$ s,
respectively.. Occasionally, we restore $G$ and $c$ when we feel it helps to
keep track of physical units.

We adopt the following definition for the Eddington mass accretion rate:

$\dot{M}_{\rm{Edd}}=\dfrac{L_{\rm{Edd}}}{\eta_{\rm NT}c^{2}},$ (9)

where $L_{\rm{Edd}}=1.25\times 10^{38}\,(M_{\rm{BH}}/M_{\odot})\,{\rm
erg\,s^{-1}}$ is the Eddington luminosity, $\eta_{\rm{NT}}$ is the radiative
efficiency of a NT thin disk around a BH with spin parameter $a_{*}$,

$\eta_{\rm{NT}}=1-\sqrt{1-\dfrac{2}{3r_{\rm{ISCO}}}}.$ (10)

In the above expression,
$r_{\rm{ISCO}}=3+Z_{2}-\sqrt{(3-Z_{1})(3+Z_{1}+2Z_{2})}$ is the radius of the
Innermost Stable Circular Orbit (ISCO, Novikov & Thorne 1973) in the Kerr
metric, where
$Z_{1}=1+(1-a_{*}^{2})^{1/3}\left((1+a_{*})^{1/3}+(1-a_{*})^{1/3}\right)$ and
$Z_{2}=\sqrt{3a_{*}^{2}+Z_{1}^{2}}$. For $a_{*}=0.0$, the efficiency is
$\eta_{\rm{NT}}=5.70\%$ while for $a_{*}=0.9$ the efficiency is
$\eta_{\rm{NT}}=15.58\%$.

In this section, we provide brief descriptions of quantities used to study the
simulation. We compute the net mass inflow rate as:

$\dot{M}(r)=-\int_{0}^{\pi}\int_{0}^{2\pi}\sqrt{-g}\,\rho\,u^{r}d\phi
d\theta.$ (11)

We treat the accretion of mass onto the BH as this integral taken at the
horizon $r_{H}$, and for simplicity we define $\dot{M}\equiv\dot{M}(r_{H})$
throughout the text. We also define the time averaged Eddington ratio (with
time averages over $3000\,t_{g}$):

$\lambda\equiv\dfrac{\langle\dot{M}\rangle}{\dot{M}_{\rm{Edd}}},$ (12)

which we use to differentiate regions of super-Eddington ($\lambda>1$) and
sub-Eddington ($\lambda<1$) accretion states. We commonly refer to $\lambda$
to describe stages in each simulation. We also consider the net mass outflow
rate at radius $r$:

$\dot{M}_{\rm{out}}(r)=\int_{0}^{\pi}\int_{0}^{2\pi}\sqrt{-g}\,\rho\,{\rm{max}}(u^{r},0)d\phi
d\theta.$ (13)

The magnetic flux is computed as:

$\Phi=\dfrac{1}{2}\int_{0}^{\pi}\int_{0}^{2\pi}|B^{r}(r)|\sqrt{-g}\,d\phi
d\theta,$ (14)

where $B^{r}(r)$ is the radial component of the magnetic field in Gaussian
units. We quantify the magnetic field strength at a given radius through the
normalized magnetic flux parameter (Tchekhovskoy et al., 2011):

$\phi(r)=\dfrac{\Phi(r)}{\sqrt{\dot{M}(r)}}.$ (15)

Throughout, we use $\phi_{\rm{BH}}\equiv\phi(r_{H})$ to describe the
normalized magnetic flux threading the BH.

To characterize the properties of the gas, we use the total pressure
$p_{\rm{tot}}\equiv p_{\rm{gas}}+p_{\rm{rad}}+p_{\rm{mag}}$, where the
individual terms are the gas pressure ($p_{\rm{gas}}$), radiation pressure
($p_{\rm{rad}}$), and magnetic pressure ($p_{\rm{mag}}$). It is also useful to
define the magnetic pressure ratio

$\beta\equiv\dfrac{(p_{\rm{gas}}+p_{\rm{rad}})}{p_{\rm{mag}}}.$ (16)

We quantify the disk thickness by estimating the density scale height:

$\dfrac{H}{R}=\sqrt{\dfrac{\int\int\rho\sqrt{-g}\,\,|\pi/2-\theta|^{2}\,d\theta
d\phi}{\int\int\rho\sqrt{-g}\,d\theta d\phi}.}$ (17)

The above integral is performed over all angles $\theta$ and $\phi$.

We estimate the electron scattering photosphere location for an observer at
infinity along the direction $(\theta,\phi)$ by integrating the optical depth
radially inward from the outer boundary of the grid. Far from the BH, the
curvature of spacetime is negligible, so we simply integrate at constant
$(\theta,\phi)$ in the “lab frame”:

$\tau_{\rm{es}}(r)=-\int_{R_{\rm{max}}}^{r}\dfrac{\rho\kappa_{\rm{es}}}{c}\left(u^{t}-u^{r}\right)\sqrt{g_{rr}}\,dr^{\prime},$
(18)

where $\kappa_{\rm{es}}=0.2(1+X)\kappa_{\rm{KN}}\,{\rm cm^{2}}$ is the
electron scattering opacity, $\kappa_{\rm{KN}}$ is the Klein-Nishina
correction factor for thermal electrons (Sądowski et al., 2017), and
$R_{\rm{max}}$ is the radius corresponding to the outer boundary of the grid.
For the gas and radiation temperatures in the simulations presented here, the
Klein-Nishina correction is negligible and the electron scattering opacity is
essentially equal to the Thomson opacity. In this work, we choose the location
of the photosphere as the $\tau_{\rm{es}}=1$ surface.

We also characterize the accretion flow using the magnetization parameter,
$\sigma_{M}\equiv b^{2}/\rho$, which quantifies the magnetic energy to rest
mass energy ratio. We define the magnetized jet as the region where
$\sigma_{M}>1$. To maintain numerical stability, we introduce mass in the
magnetized jet using a ceiling condition on the magnetization $\sigma_{M}\leq
100$ throughout each simulation. In measuring the jet profile, we make use of
the $\sigma_{M}=1$ contour. However, when computing the total energy flux in
the jet, we use a different definition of the jet boundary, based on the
Bernoulli parameter. We prefer this definition since gas which may be
accelerated to highly relativistic asymptotic speeds ($\gamma_{\infty}\gtrsim
2$) will lie in the $\sigma_{M}>1$ region. In addition, the Bernoulli
parameter profiles are less well defined at large radii and are not useful for
computing the jet shape. However, a significant fraction of gas which is too
fast ($\gamma_{\infty}>1.05$) to be called a ‘wind’ lies outside of the
$\sigma_{M}>1$ region so we make use of the Bernoulli definition to compute
the total energy of the jet.

We are particularly interested in the energy outflows. To that end, we define
the outgoing radiative bolometric luminosity by:

$L_{\rm{bol}}(r_{\rm{lum}})=\int_{0}^{\pi}\int_{0}^{2\pi}\sqrt{-g}\,\,{\rm{max}}(-R^{r}_{\
\,t},0)d\phi d\theta,$ (19)

which gives the flux of outgoing radiation energy through a surface at radius
$r_{\rm{lum}}=\max(30\,r_{g},r_{\tau_{\rm{es}}=1})$. By this definition, the
default $r_{\rm{lum}}=30\,r_{g}$, which lies inside of the injection radius
$R_{\rm{inj}}$, but it moves out to the photosphere if the latter lies outside
$30\,r_{g}$.

The MHD energy flux is computed as:

$L_{\rm{MHD}}(r)=-\int_{\theta_{1}}^{\theta_{2}}\int_{0}^{2\pi}\sqrt{-g}\,(T^{r}_{\
\,t}+\rho u^{r})d\phi d\theta,$ (20)

where we integrate the radial flux of energy carried by gas plus magnetic
field ($T^{r}_{\ \,t}$) and subtract out the rest-mass energy ($\rho u^{r}$,
since it does not lead to observational consequences for an observer at
infinity). Note that the integral is only over
$\theta_{1}\leq\theta\leq\theta_{2}$ since we consider the energy outflow in
jet and wind seperately. We define these regions below. When measured at large
radii, this quantity determines the total outflowing MHD energy extracted from
the accretion disk.

We define the accretion flow as three distinct regions (disk, wind, and jet)
based on the total energy via the Bernoulli parameter. We also make use of the
electron scattering opacity, $\tau_{\rm{es}}$, to determine where the region
of the fluid under consideration is optically thin or thick. In optically
thick regions ($\tau_{\rm{es}}\geq 1$), the radiation is advected with the
flow and can contribute to acceleration of the gas so we treat it as
contributing to the Bernoulli parameter. Meanwhile, in optically thin regions
($\tau_{\rm{es}}<1$) we assume that only the MHD components are relevant to
the total gas energy. Building on the definition given in Curd & Narayan
(2019), we define a Bernoulli parameter which smoothly connects between
optically thin and optically thick regions via:

$\rm{Be}=\begin{cases}-\dfrac{(T^{t}_{\ \,t}+\rho
u^{t})+(1-\tau_{\rm{es}}^{-1/2})R^{t}_{\ \,t}}{\rho u^{t}},\
\,\tau_{\rm{es}}\geq 1\\\ -\dfrac{(T^{t}_{\ \,t}+\rho u^{t})}{\rho u^{t}},\
\,\tau_{\rm{es}}<1\\\ \end{cases}$ (21)

We have adopted this definition to account for the fact that radiative
acceleration of the gas will be weak in the optically thin region. We define
the ‘disk’ region of the simulation as bound gas with $\rm{Be}<0$. Both the
‘wind’ and ‘jet’ are unbound and are generally radially out flowing. The wind
is defined as any fluid with $0<\rm{Be}\leq 0.05$. The jet is any fluid with
$\rm{Be}>0.05$. This choice of cutoff for wind vs. jet is based on the
velocity at infinity ($v_{\infty}$), with the wind having $v_{\infty}\lesssim
0.3c$ and the jet having $v_{\infty}\gtrsim 0.3c$.

We define the radiation, jet, and wind efficiencies as:

$\eta_{\rm{rad}}=\dfrac{1}{\dot{M}}L_{\rm{bol}}(r_{\rm{lum}}),$ (22)
$\eta_{\rm{jet}}=\dfrac{1}{\dot{M}}L_{\rm{MHD}}(r=30\,r_{g},{\rm{Be}}\geq
0.05),$ (23)

and

$\eta_{\rm{wind}}=\dfrac{1}{\dot{M}}L_{\rm{MHD}}(r=30\,r_{g},0.05>{\rm{Be}}>0).$
(24)

In each case, positive values correspond to energy being extracted from the
system. We define the total efficiency as:

$\eta_{\rm{tot}}=\eta_{\rm{jet}}+\eta_{\rm{wind}}+\eta_{\rm{rad}}.$ (25)

We are also interested in the total angular momentum flux flowing into the
system:

$\dot{J}(r)=-\int_{0}^{\pi}\int_{0}^{2\pi}\sqrt{-g}\,(T^{r}_{\
\,\phi}+R^{r}_{\ \,\phi})d\phi d\theta.$ (26)

The specific angular momentum flux into the BH is then

$j(r)=\dfrac{\dot{J}(r)}{\dot{M}}.$ (27)

In this work, we adopt $r=r_{H}$ to measure the specific angular momentum into
or out of the BH. If $j(r_{H})>0$, angular momentum is being gained by the BH,
and if $j(r_{H})<0$, angular momentum is being lost by the BH. We follow
Shapiro (2005) and define the BH’s dimensionless spin-up parameter:

$s\equiv\dfrac{d(a/M)}{dt}\dfrac{M}{\dot{M}}=j-2a_{*}(1-\eta_{\rm{tot}}),$
(28)

where positive values correspond to the BH spin increasing and negative values
correspond to the BH spin decreasing with time. More positive/negative values
imply more rapid changes in the BH spin.

###  2.2 Disk Injection

Figure 1: Example magnetic field at $t=15,000\,t_{g}$ in model KIa09. The gas
density of the disk (colors) is shown along with the magnetic field lines and
the $\sigma_{M}=1$ contour (blue line), which delineates the highly magnetized
region of the jet.

The simulation domain is initialized with a low density atmosphere with a
density profile that scales as $r^{-2}$ and has a maximum density of $6\times
10^{-17}{\rm{\,g\,cm^{-3}}}$ at the horizon. The atmosphere is initialized
with a constant radiation temperature of $T_{\rm{atm}}=10^{5}$ K.

Using a similar scheme to Curd (2021), we introduce gas that is initially on a
circular orbit into the simulation domain via an injection boundary condition
in a small ring of cells at an injection radius $R_{\rm{inj}}=40\,r_{g}$ that
is interior to the maximum radius ($10^{5}\,r_{g}$) of the simulation domain.
The vertical extent of the injected disk is set to $(H/R)_{0}=0.05$, but as
the gas first flows into the simulation domain (prior to the formation of a
disk), it self-consistently puffs up due to the high magnetic and radiation
pressure. The gas density in the injection cells varies with vertical distance
from the mid plane ($\theta=\pi/2$) as
$\rho\propto\left[1-\left(\frac{|\theta-\pi/2|}{(H/R)_{0}}\right)^{2}\right]^{3}$.
The normalization is set to give the desired mass injection rate, which we
define later. The primary motivation for this approach is that it allows for
precise control of the accretion rate as the simulation evolves in time so
long as the timescale on which the accretion rate is varied is longer than the
viscous accretion time from $R_{\rm{inj}}$ to the horizon $r_{H}$. The
accretion time for a parcel of gas in a MAD, geometrically thick disk is
approximately given by (in gravitational units):

$\dfrac{T_{\rm{acc}}}{t_{g}}\approx\dfrac{r}{0.1v_{\rm{Kep}}}=10r^{3/2},$ (29)

where $r$ is the radius, and $v_{\rm{Kep}}=r^{-1/2}$ is the Keplerian orbital
velocity. Here we have assumed that the radial velocity of the inflowing gas
is 10% of the Keplerian velocity motivated by the findings of McKinney et al.
(2012).

For this numerical experiment, we chose $R_{\rm{inj}}=40\,r_{g}$, which
corresponds to an accretion time of $T_{\rm{acc}}\approx 2500\,t_{g}$. We
initialize the orbital properties of the gas by setting the radial velocity
$v_{r}=0$ and setting the angular momentum to that of a circular orbit at the
injection radius.

A nominal gas temperature of $T_{\rm{inj}}=10^{9}$ K is set at the injection
boundary. This is not the true temperature, but is a measure of the total
pressure. We use the initial gas pressure obtained from
$p^{\prime}=p_{\rm{gas}}(T_{\rm{inj}})$ to split the internal energy into gas
and radiation energy density by solving the condition
$p^{\prime}=p_{\rm{gas}}+p_{\rm{rad}}$ and finding a new gas and radiation
temperature which assumes thermal equilibrium of the gas. This gives
$T_{\rm{gas}}\approx 5\times 10^{5}$ K in the injection cells.

We set the magnetic field in the two radially interior injection boundary
cells to an initially nearly vertical field geometry
$B^{\theta}\propto\cos\left[\frac{\theta-\pi/2}{(H/R)_{0}}\right]$ and a
magnetic pressure ratio of $\beta=33$. This choice of $\beta$ leads to the
rapid formation of a MAD disk as the gas begins to accrete. We set
$B^{r}=B^{\phi}=0$ at the injection point. The divergence constraint
$\nabla\cdot\mathbb{B}=0$ is enforced through the flux constrained transport
algorithm introduced by Tóth (2000). We do not explicitly initialize the ghost
cells such that $\nabla\cdot\mathbb{B}=0$ is guaranteed, but the magnetic
field that ultimately gets injected with the gas has a divergence of the order
$|\nabla\cdot\mathbb{B}|/|\mathbb{B}|\sim 10^{-7}$ throughout the evolution.
We note that a similar magnetic field initialization was used in (Kaaz et al.,
2022) to simulate a magnetized Bondi-Hoyle-Lyttleton flow; however, they
inject gas from the outer boundary of the simulation. We show a snapshot of
the resulting disk and field loops in Figure 1.

The injection boundary is four cells wide in the radial direction. At the two
radially interior cells, we set the gas primitives using the quantities
described above. Meanwhile, for the two radially exterior cells we employ a
reflecting boundary condition in radius. This is to avoid any spurious inflows
or outflows behind the injection point. Gas moving with $v^{r}<0$ at
$r>40\,r_{g},\,|\theta-\pi/2|\leq(H/R)_{0}$ encounters a reflective boundary
and cannot flow further in.222We note that this does not affect the validity
of the numerical solution of the inflow/wind/jet since the gas behind the
injection point is outflowing after the disk forms.

One difficulty with treating the injection of gas inside of the simulation
domain is the choice of boundary conditions at the top and bottom ($\theta$)
edges and the radial ($r$) edges. Enforcing reflection at both the exterior
radius and the top/bottom edges in $\theta$ is non-trivial. For instance,
setting the primitives in the ghost cells at the radial face to prevent
inflow/outflow does not guarantee that no inflow/outflow will occur at the
top/bottom angular corner cells. This is particularly pernicious with MAD
disks as enforcing a reflecting boundary condition on the magnetic field can
lead to unwanted inflow of gas carrying the opposite magnetic field polarity
at the angular boundary. We sought to minimize mass from flowing out of the
simulation domain through the injection region while maintaining the desired
injected field topology at all times, so we opted to set the magnetic field
$B^{r}=B^{\theta}=B^{\phi}=0$ at the 2 radially exterior grid cells and
enforce a reflecting boundary condition only in radius on the other gas
primitives at the injection boundary. This guaranteed that no magnetic field
is injected if mass is allowed to flow in through a corner cell. We do not
believe this to be a major issue as the gas flowing around the boundary is
radially outflowing throughout the simulation. In addition, we compute the
simulation properties near the midplane at $r<R_{\rm{inj}}$ throughout this
work where the effect of these boundary conditions should be minimal. We
demonstrate in Appendix A that either reflecting $\mathbb{B}$ or setting
$\mathbb{B}=0$ at the 2 radially exterior cells does not impact the
conclusions drawn from the simulation.

Here we further detail the regrid process and how the simulation evolves
through each stage. As detailed in the previous subsection, we first inject
the gas in 2D with a resolution of 256 cells in $r$ and 128 cells in $\theta$.
We initially fix the injection rate to $\dot{M_{0}}=10\dot{M}_{\rm{Edd}}$
until the magnetic flux at the horizon reaches $\phi_{\rm{BH}}>40$.

After this stage, which takes $10,000\,t_{g}$ in each simulation, we perform
‘regrid 1’ and copy the 2D data onto a fully 3D grid with 128 cells in polar
coordinate $\phi$. We also introduce a 5% perturbation in the azimuthal
velocity before restarting the simulation to break the azimuthal symmetry. In
KIa00, we continue injecting mass at a fixed rate until $t=15,000\,t_{g}$ and
begin decreasing it after that time. For KIa09, we begin to decrease the mass
injection rate at $t=10,000\,t_{g}$, immediately after the regrid. The mass
inflow rate at the injection boundary decreases over time following a
$\dot{M}\propto t^{-5/3}$ profile. The exact description of the mass injection
is given by:

$\dot{M}_{\rm{inj}}(t)=\dot{M}_{0}\left[\dfrac{(t-t_{0})}{T_{\rm{acc}}(R_{\rm{inj}})}+1\right]^{-5/3}$
(30)

where $\dot{M}_{0}$ is the peak mass inflow rate, $t_{0}$ is the time at which
the decay starts, and $t$ is the time since the beginning of injection in
gravitational units. The power-law dependence that we have adopted in Equation
30 is arbitrary and is merely adopted for its similarity to the time evolution
in TDEs, whose mass fallback rate varies approximately as
$\dot{M}_{\rm{fb}}\propto t^{-5/3}$. A similar result could be obtained by
using a larger $T_{\rm{acc}}$ coupled with a shallower power-law as long as
the disk has time to adjust to changes at the injection boundary. It is worth
noting that $T_{\rm{acc}}\approx 2500\,t_{g}$ is substantially shorter than
the decay times of $\sim 28,000-180,000\,t_{g}$ adopted in Curd (2021). In
this work, our somewhat small choice of $T_{\rm{acc}}$ is motivated by the
computational requirements to probe a range of accretion states that differ by
roughly an order of magnitude over the duration of a single simulation.

After the disk becomes sub-Eddington, which happens around $21,000\,t_{g}$ for
KIa00 and $20,000\,t_{g}$ for KIa09, we perform ‘regrid 2’ and double the
resolution in $\theta$ from $N_{\theta}=128$ to $N_{\theta}=256$ to better
resolve the disk. In terms of resolution, the grid has $\Delta\theta\approx
0.005\,(0.01)$ near the mid-plane after (before) the regridding. As a result,
we resolve the disk with roughly 30 (15) cells per disk scale height (we
provide the precise disk scale height later in the text), even at the lowest
accretion rates we consider.

Figure 2: In the top row we show snapshots of the accretion disk in the
midplane spanning the evolution of KIa00 from $t=15,000-29,000\,t_{g}$ in a
$60r_{g}\times 60r_{g}$ region centered on the black hole. The gas density,
scaled by the maximum density in the disk ($\rho/\rho_{\rm{max}}$), is
indicated by colors in each panel. Gas primarily flows in along several spiral
arms and the low density regions (blue) of highly magnetized gas, sometimes
referred to as ‘magnetic islands,’ can sometimes make up a significant
fraction of the surface area in the disk. The full simulation history is shown
in the bottom panels. We show the time evolution of accretion rate $\dot{M}$,
normalized magnetic flux at the BH horizon $\phi_{\rm{BH}}$, and bolometric
luminosity $L_{\rm{bol}}$. The vertical gray lines indicate the time at which
we perform the regrid from 2D to 3D (regrid 1) and when we double the $\theta$
angular resolution (regrid 2). The accretion rate $\dot{M}$ and bolometric
luminosity $L_{\rm{bol}}$ roughly follow a $t^{-5/3}$ power law, indicated as
the dashed line in the top panel. The magnetic flux after $t>10,000\,t_{g}$
quickly decreases to $\phi_{\rm{BH}}\approx 30$ and is highly variable,
indicating that flux eruptions are frequent.

Figure 3: Here we present the time averaged radial profiles (top panel) for
the accretion rate (solid lines) and mass outflow rate (dashed lines), and the
density scale height as a function of time at selected radii (bottom panel),
for KIa00. The time averaging range for each profile is indicated at the top
of the panel. In the bottom panel, we show the density scale height throughout
the evolution at $5\,r_{g}$ (solid line), $10\,r_{g}$ (dashed line), and
$20\,r_{g}$ (dashed-dotted line).

##  3 Results

###  3.1 Time Evolution

####  3.1.1 Zero Spin Model: KIa00

We show the full time evolution of the accretion rate $\dot{M}$, normalized
magnetic flux at the BH horizon $\phi_{\rm{BH}}$, and the bolometric
luminosity $L_{\rm{bol}}$ for KIa00 in Figure 2. Also shown are selected
snapshots of the disk at the equatorial plane which span the full 3D
evolution.

The accretion rate very closely follows the mass injection rate. As detailed
in the previous section, we maintained a constant mass injection rate
$\dot{M}_{\rm{inj}}$ until $t\geq 15,000\,t_{g}$, and the accretion rate
$\dot{M}_{\rm{BH}}$ into the BH prior to this stage rises to a mean value of
about $\dot{M}/\dot{M}_{\rm{Edd}}\approx 3$ and then remains nearly constant
as a result. After we start to decrease the injection rate (see Equation 30),
the accretion rate similarly follows a $t^{-5/3}$ power law. By the end of the
simulation, we sample an accretion state with
$\dot{M}/\dot{M}_{\rm{Edd}}\approx 0.35$.

The magnetic flux grows rapidly in the initial 2D evolution from
$t=0-10,000\,t_{g}$ and reaches a maximum value of $\phi_{\rm{BH}}\approx 55$
where it ceases to grow. Before the disk becomes MAD, it is the MRI which
drives the angular momentum transport and hence the accretion. We find that
the density weighted MRI quality factor integrated over all angles is
$Q_{\theta}>100$ during the entire evolution so MRI is well resolved in our
simulations. After regridding to 3D, the magnetic flux escapes from the
horizon in frequent flux eruption events which flow out radially and partially
mix with the disk at larger radii. The most dramatic drop in flux occurs
around $11,000\,t_{g}$ when $\phi_{\rm{BH}}$ decreases from $\approx 55$ to
$\approx 20$. As a result of the flux eruption events, the disk becomes highly
non-axisymmetric throughout the entire evolution in 3D and the gas flows in
primarily along dense spiral arms with much lower density gas filling the area
in between, where the magnetic field is very strong. This behaviour is clearly
shown in the top panels in Figure 2. The net effect of the flux eruptions is
that the magnetic flux saturates at a lower value than its initial peak,
varying between $\phi_{\rm{BH}}\approx 20-40$ throughout. Flux eruptions are
well known phenomena in MAD accretion disks (for instance, see Igumenshchev
2008; Tchekhovskoy et al. 2012), but our disk appears to be more variable and
have more severe flux eruptions in comparison to other models in the
literature. This is likely due to the fact that the entire disk is MAD out to
$R_{\rm{inj}}$, and we employ a unique disk initialization which drives the
disk towards a large magnetic flux at all times. Gas flowing in from the
injection radius is evidently able to efficiently transport magnetic flux
inward and rapidly restores any lost magnetic flux at the horizon. The fact
that $\phi_{\rm{BH}}$ is lower than the typical $\phi_{\rm{BH}}\sim 50$ seen
in hot accretion flows in the MAD state (e.g. Narayan et al. 2022) is likely
because our systems are vertically thinner.

The luminosity behaves similarly to the accretion rate. As in the case of
SANE, super-Eddington disks (Sądowski & Narayan, 2016; Curd & Narayan, 2019),
the luminosity is marginally super-Eddington when $\lambda=2.8$ and maintains
a value of $L_{\rm{bol}}\approx 2.8L_{\rm{Edd}}$. At the end of the simulation
when $\lambda=0.35$, $L_{\rm{bol}}\approx 0.76L_{\rm{Edd}}$.

We also consider how the simulation evolves in a time averaged sense by
averaging the simulation data over $\phi$ and over chunks of $3000\,t_{g}$ in
time. This allows us to sample the accretion state of the disk at various time
intervals while ignoring the rapid variability; furthermore, the overall
change in the disk properties over a given $3000\,t_{g}$ chunk of time is not
substantial. As shown in Table 2, the magnetic flux is stable across accretion
rates with $\langle\phi_{\rm{BH}}\rangle\approx 32-33$. The luminosity does
not follow precisely the same decline as $\lambda$, and $\langle
L_{\rm{bol}}\rangle$ instead declines slightly more slowly. This has
implications for the efficiency of the system, which we discuss later.

The net accretion rate as a function of radius is shown in the left panel of
Figure 3. We find that the profiles are nearly constant out to $r\approx
30\,r_{g}$ with a net outflow rate that is generally much less than the net
inflow rate. Flat accretion rate profiles are typically seen in simulations
where the initial conditions are an equilibrium magnetized torus. In such
simulations, the gas flows in on the viscous accretion time and the disk is
said to have reached inflow equilibrium up to the maximum radius where the net
accretion rate is flat. The present simulations have a time dependent mass
injection rate $\dot{M}_{\rm{inj}}$. Nevertheless, the profiles that we
measure suggest that the disk is in a quasi inflow equilibrium state out to
$r=30\,r_{g}$.

We measure the density scale height of the disk following Equation 17 at radii
$r=5,10,20\,r_{g}$ (bottom panel in Figure 3, see also Table 2). The disk
exhibits substantial thinning at all radii as the accretion rate decreases.
The initial scale height is $H/R\approx 0.25$ during the super-Eddington phase
and increases slightly with radius. This is similar to the scale height of
previous super-Eddington disks (Sądowski et al., 2015; Sądowski & Narayan,
2016). When the accretion rate drops to $\lambda=0.35$, the scale height
reaches $H/R\approx 0.15$ but the behavior as a function of radius is
inverted, being smaller at $r=10,20\,r_{g}$ compared to $r=5\,r_{g}$. While
the disk appears to thin somewhat as the accretion rate drops, it is clearly
stable against runaway cooling and collapse. This is to be expected since the
magnetic field provides pressure support and accounts for nearly half of the
pressure in the disk. A similar result for thin MAD disks was also found by
Morales Teixeira et al. (2018); Liska et al. (2022).

Finally we come to the efficiencies $\eta$ of the various forms of energy
outflow. Namely, we measure the total efficiency ($\eta_{\rm{tot}}$), the
radiative efficiency ($\eta_{\rm{rad}}$), and the jet and wind efficiencies
($\eta_{\rm{jet}}$ and $\eta_{\rm{wind}}$, which both exclude radiation). We
find that $\eta_{\rm{tot}}$ increases with decreasing accretion rate. We also
observe the same trend in terms of $\eta_{\rm{rad}}$, which is initially only
$5.7\%$ when the disk is super-Eddington and increases to $12.2\%$ when
$\lambda=0.35$. We interpret this change as being related to the weaker
advection of radiation as the disk scale-height and optical depth changes. For
instance, the shrinking of the outer disk at $r\geq 10\,r_{g}$ to a smaller
scale height than the inner disk suggests that radiative cooling is acting to
remove thermal energy from the disk before it is advected across the BH
horizon. Ultimately, this leads to a larger escaping radiant luminosity,
relative to the accreted rest mass energy, and an increase in efficiency as
$\lambda$ decreases. At the same time, $\eta_{\rm{jet}}$ increases from
$1.9\%$ to $3.4\%$ while $\eta_{\rm{wind}}$ decreases from $1.1\%$ to $0.3\%$.
This indicates that the mechanical component of the outflow energy (jet and
wind) tilts from the wind to the jet as $\lambda$ decreases. Indeed, as we
discuss later, the outflow velocity near $30\,r_{g}$ increases as the
simulation evolves. See Table 2 for tabulated efficiencies as a function of
time and accretion rate.

Focusing on the sub-Eddington evolution, a total efficiency of
$\eta_{\rm{tot}}\approx 15.1-15.9\%$ is more than three times the NT thin disk
efficiency. Additionally, a radiative efficiency of $\eta_{\rm{rad}}\approx
11.6-12.2\%$ is more than two times the NT thin disk efficiency. Both metrics
for the extracted energy show substantial differences from the theoretical
value expected for a NT thin disk. We note that other authors have
demonstrated deviations from the NT thin disk model in GRRMHD for MAD disks
(Dexter et al., 2021; Liska et al., 2022). GRMHD simulations which utilized a
cooling function have also found significant deviations for MAD (Avara et al.,
2016) and non-MAD (Kinch et al., 2021) accretion disks. The additional
radiation that we find here may be provided by the binding energy of gas
plunging into the BH. In particular, low angular momentum gas will have more
binding energy available for dissipation. We find that the disk becomes
significantly sub-Keplerian at $r<R_{\rm{inj}}$ and $u_{\phi}$ is nearly half
the Keplerian angular momentum as it crosses the horizon, which suggests that
this could be a contributing factor to the deviations we measure. There is
also the matter of significant dissipation at radii inside the ISCO (which is
also present in our simulations), which is not predicted in the NT model. We
note that the inclusion of a two-temperature plasma at low accretion rates may
significantly reduce emission. For example, Kinch et al. (2020) find a coronal
luminosity reduction of $\sim 1/2$ when contrasting a single temperature
plasma with a two temperature plasma.

Model | Time | $\lambda$ | $\langle\phi_{\rm{BH}}\rangle$ | $\langle L_{\rm{bol}}\rangle$ | $(H/R)_{r=5\,r_{g}}$ | $\eta_{\rm{tot}}$ | $\eta_{\rm{rad}}$ | $\eta_{\rm{jet}}$ | $\eta_{\rm{wind}}$
---|---|---|---|---|---|---|---|---|---
| ($10^{3}\,t_{g}$) |  |  | $(L_{\rm{Edd}})$ |  |  |  |  |
KIa00 | 15-18 | 2.8 | 33 | 2.80 | 0.24 | 8.7% | 5.7% | 1.9% | 1.1%
| 18-21 | 1.2 | 33 | 1.88 | 0.19 | 11.5% | 8.7% | 2.5% | 0.3%
| 22-25 | 0.58 | 32 | 1.19 | 0.18 | 15.1% | 11.6% | 3.2% | 0.3%
| 26-29 | 0.35 | 32 | 0.76 | 0.16 | 15.9% | 12.2% | 3.4% | 0.3%
NT ($a_{*}=0$) | - | - | - | - | - | 5.70% | 5.70% | - | -
KIa09 | 14-17 | 2.2 | 36 | 6.70 | 0.18 | 88.0% | 48.2% | 38.7% | 0.6%
| 20-23 | 0.76 | 33 | 3.68 | 0.15 | 100.1% | 75.1% | 25.8% | 0.6%
| 23-26 | 0.44 | 29 | 2.21 | 0.15 | 98.3% | 78.8% | 18.9% | 0.6%
| 28-31 | 0.31 | 26 | 0.99 | 0.13 | 68.3% | 50.2% | 17.7% | 0.5%
| 31-34 | 0.29 | 27 | 0.85 | 0.12 | 64.3% | 45.4% | 18.4% | 0.4%
NT ($a_{*}=0.9$) | - | - | - | - | - | 15.58% | 15.58% | - | -
Table 2: Tabulated values of the time- and $\phi$-averaged magnetic flux at
the BH horizon ($\langle\phi_{\rm{BH}}\rangle)$, bolometric luminosity
($\langle L_{\rm{bol}}\rangle$), disk scale height ($H/R$), total efficiency
($\eta_{\rm{tot}}$), radiative efficiency ($\eta_{\rm{rad}}$), jet efficiency
($\eta_{\rm{jet}}$), and wind efficiency($\eta_{\rm{wind}}$), for each model
as a function of time and Eddington-scaled accretion rate. We also list the
corresponding NT thin disk efficiency for comparison. Note that the time range
for time-averaging is shown in units of $10^{3}\,t_{g}$ and $\lambda$ is the
average accretion rate in units of $\dot{M}_{\rm{Edd}}$ during the given time
range. The disk scale height presented here is computed at $r=5\,r_{g}$.

Figure 4: The same as Figure 2 but for KIa09. The simulation history in the
bottom panels is similar to that found for KIa00 but the luminosity is higher
because of the contribution from a relativistic jet. Also note that the
injection rate starts to decline at $t=10,000\,t_{g}$.

Figure 5: The same as Figure 3 but for KIa09.

####  3.1.2 High Spin Model: KIa09

As demonstrated in Figure 4, the evolution of KIa09 is quite similar to KIa00
in several respects. The accretion rate after we begin to decrease the
injection rate (at $t>10,000\,t_{g}$) closely tracks the $t^{-5/3}$ behaviour
of the injection rate, and the accretion rate decreases to $\approx
0.29\dot{M}_{\rm{Edd}}$ by the end of the simulation. The magnetic flux rises
to an initial saturation value of $\phi_{\rm{BH}}\approx 55$ but then declines
during the 3D evolution before leveling off at $\phi_{\rm{BH}}\approx 20-40$.
The disk has frequent flux eruptions, and the BH again primarily accretes gas
through dense spiral arms. We note that the initial flux eruption after
‘regrid 1’ occurs at $t\approx 10,500\,t_{g}$ and is both earlier and far less
dramatic than the first flux eruption in KIa00 after ‘regrid 1’. The time
averaged magnetic flux (Table 2) declines throughout the simulation starting
at $\langle\phi_{\rm{BH}}\rangle\approx 36$ when $\lambda=2.2$ but it declines
to $\langle\phi_{\rm{BH}}\rangle\approx 26-27$ by the time $\lambda\leq 0.31$.
This suggests that magnetic flux is advected less efficiently at lower
accretion rates in KIa09, perhaps due to the lower $H/R$ that we measure. This
is the opposite of KIa00, which showed essentially no change in magnetic flux
by the end of the simulation; however, we note that $H/R$ was also generally
larger in KIa00.

The accretion rate as a function of radius in KIa09 also behaves similarly to
that in KIa00 (compare the top panels of Figures 3 and 5), and we observe
profiles which are nearly constant out to $r\approx 30\,r_{g}$. The net
outflow rate of mass is again substantially lower than the mass inflow at
small radii, but at $r=40\,r_{g}$ it is nearly equal to the mass inflow. This
is likely due to the BZ process, which leads to a substantially higher jet
power than in the non-spinning case.

The two primary differences between KIa00 and KIa09 are that (i) the disk
starts out noticeably thinner in KIa09 (bottom panel in Figure 5), and (ii)
the total efficiency, jet efficiency, and radiative efficiency are all
substantially higher in KIa09 (Table 2). Both of these differences come about
because of the addition of a spinning BH and the fact that the disk is MAD.
For example, the theoretical energy extraction via the BZ power assuming
$a_{*}=0.9$ and $\phi_{\rm{BH}}=30$ in Equation 2 gives an efficiency due to
the BZ process of $\eta_{\rm{BZ}}\approx 31\%$. This extracted energy powers a
jet which supplies significant pressure and compresses the disk while
simultaneously producing radiation and driving a more energetic outflow in the
form of a relativistic outflow. This extracted energy is in addition to any
energy extracted from the disk and comes instead from the BH spin. We note
that our total efficiency is substantially higher than $\eta_{\rm{BZ}}$, but
efficiencies larger than $\eta_{\rm{BZ}}$ are not uncommon for MAD disks
(Tchekhovskoy et al., 2012; McKinney et al., 2012; Narayan et al., 2022; Liska
et al., 2022).

The inner disk gets compressed due to the jet, and we observe a peak scale
height of $H/R\approx 0.18$ (compared to $H/R\approx 0.25$ in KIa00, this is
similar to the results reported by Narayan et al. 2022). Unlike KIa00, the
disk does not have an inverted $H/R(r)$ profile as the accretion rate
approaches $\approx 0.3\dot{M}_{\rm{Edd}}$ and instead the outer disk reaches
$H/R(r=20\,r_{g})\approx 0.14$ and the inner disk reaches
$H/R(r=5\,r_{g})\approx 0.12$.

We find that $\eta_{\rm{tot}}$ and $\eta_{\rm{jet}}$ have a downward trend as
the accretion rate drops. Notably, $\eta_{\rm{tot}}\approx 88.0\%$ when
$\lambda=2.2$ and increases to $\eta_{\rm{tot}}\approx 100\%$ when
$\lambda=0.44-0.76$. However, as the accretion rate declines further to
$\lambda=0.29$, $\eta_{\rm{tot}}\approx 64.3\%$. On the other hand,
$\eta_{\rm{jet}}$ is highest ($\approx 38.7$%) when the accretion rate is
super-Eddington and declines with $\lambda$, reaching only $\approx 18.4\%$
when $\lambda=0.29$. The wind efficiency remains fairly stable at
$\eta_{\rm{wind}}\approx 0.4-0.6\%$ throughout the entire evolution. See Table
2 for tabulated efficiencies in KIa09. Here we again find a large deviation
from the NT thin disk, with the total efficiency exceeding NT by a factor of
nearly $4$ and the radiative efficiency exceeding NT by a factor of $3$ during
the lowest accretion states.

We note that the numerical ceiling implemented in our simulations which
requires $\sigma_{M}\leq 100$ likely impacts the radiative efficiency in
KIa09. Mass is injected into the jet to maintain $\sigma_{M}\leq 100$. This
ultimately leads to some dependence on the numerical floors in the radiation
field in the jet since the gas in this region is both hot and magnetized and
will thus produce radiation. This uncertainty is difficult to quantify even
where $1<\sigma_{M}<100$ since the jet is optically thin and the radiation
from regions effected by floors will spread throughout the jet (e.g. McKinney
et al. 2017; Chael et al. 2019). We performed a rerun of the 2D portion of
KIa09 to quantify the impact of the $\sigma_{M}$ ceiling on our results and
find that setting a ceiling of $\sigma_{M}\leq 20$ decreases the radiative
efficiency as $L_{\rm{bol}}$ is decreased by roughly a factor of 2. The
$\sigma_{M}$ ceiling is less of a concern in KIa00 as the region where the
numerical ceiling is applied is much smaller in extent (see Figure 18 in
Appendix B).

###  3.2 Outflow Properties and Jet Structure

Figure 6: Here we show snapshots of KIa00 for time-averages over $3000\,t_{g}$
and $\phi$-averaged. The time averaged accretion rate ($\lambda$) is indicated
above each panel. The gas density (colors), fluid velocity (streamlines),
jet/wind boundary (${\rm{Be}}=0.05$, dashed green line), wind/disk boundary
(${\rm{Be}}=0$, green line), electron scattering photosphere
($\tau_{\rm{es}}=1$, yellow line), and magnetized jet boundary
($\sigma_{M}=1$, blue line) are indicated in each panel. Figure 7: The same as
Figure 6 but for model KIa09. As the accretion rate declines, the jet power
declines substantially and the jet width declines.

Figure 8: Here we show angular profiles of the gas temperature (top), lorentz
factor (middle), and asymptotic lorentz factor (bottom) in the jet
(${\rm{Be}}\geq 0.05$) for KIa00 (left) and KIa09 (right) as a function of
polar angle $\theta$ at $r=30\,r_{g}$. Note that we have symmetrized the data
across $\theta=\pi/2$ such that quantities shown are an average of the top and
bottom jet. The solid line is for regions where $\sigma_{M}>1$ while the
dashed line signifies where $\sigma_{M}<1$ and ${\rm{Be}}\geq 0.05$. We mark
the location where $\sigma_{M}=1$ (solid circle) as well as where
${\rm{Be}}=0.05$ (open circle). The temperature of the outflow in KIa00
systematically increase with decreasing accretion rate. Both $\gamma$ and the
asymptotic lorentz factor $\gamma_{\infty}$ also generally increase with
decreasing $\lambda$ with the exception of the $\sigma_{M}>1$ region
(especially near the poles). KIa09 does not show significant differences in
the jet temperature, nor do $\gamma$ or $\gamma_{\infty}$ show clear trends in
their behaviour as $\lambda$ decreases. However, $\gamma$ during the
$3000\,t_{g}$ period that we average over where $\lambda=0.44$ is much smaller
than at other points in the simulation due to a prolonged drop in jet power,
as we discuss in the text. The peak asymptotic lorentz factor suggests the
outflow in KIa00 may accelerate to mildly relativistic outflow with
$\gamma_{\infty}\sim 1.6-2.2$ while the outflow in KIa09 may reach highly
relativistic speeds of $\gamma_{\infty}\gtrsim 30$ at all mass accretion rates
$\lambda$ that we have sampled. Figure 9: Here we show snapshots of the gas
density (colors) and magnetized jet boundary $\sigma_{M}=1$ (blue line) for
KIa09 between $t=23,000-26,000\,t_{g}$ when the jet power is high (top panel)
and when the jet power is low (bottom panel). At $t=24,000\,t_{g}$ (top
panel), the normalized magnetic flux at the BH horizon is
$\phi_{\rm{BH}}\approx 35$ and the instantaneous jet power is
$\eta_{\rm{jet}}\approx 22\%$. Here a powerful jet creates a low density
funnel where the magnetic pressure dominates. At $t=26,000\,t_{g}$ (bottom
panel), $\phi_{\rm{BH}}\approx 15$ and $\eta_{\rm{jet}}\approx 3\%$. The polar
region in this low jet power state is filled in with higher density gas with
$\sigma_{M}<1$ and the overall width of the jet narrows. The low jet power
state in the bottom panel is prolonged for $t=23,000-26,000\,t_{g}$
($\lambda=0.44$) and leads to a lower jet velocity in a time averaged sense.
Figure 10: Here we show fits to the jet radius as a function of jet height,
$w\propto z^{k}$, at different times during the evolution of KIa09 by fitting
the $\sigma_{M}=1$ profile over radii in the range $r=10-100\,r_{g}$ for both
the top and bottom jet. The power-law index $k$ appears to be slightly larger
during the super-Eddington phase of the simulation, which coincides with a
higher jet power.

Figures 6 and 7 show the disk, wind and jet structure as a function of
accretion rate for the two simulations.

For KIa00, we find that at high accretion rates ($\lambda\approx 2.8$), the
photosphere is extended in the polar region, extending to nearly $80\,r_{g}$.
As the accretion rate decreases to $\lambda\sim 1$, the polar region becomes
optically thin which allows for radiative acceleration close to the BH
horizon. The overall scale of the wind region shrinks considerably as the
accretion rate drops, which agrees with the fact that the measured wind
efficiency declines from 1.1% to 0.3% by the time $\lambda=0.35$ while the jet
efficiency increases from 1.9% to 3.4%. The enhanced energy deposited in the
jet is due to the shrinking of the photosphere, which is nearly coincident
with the disk boundary for sub-Eddington values of $\lambda$. Radiative
coupling allows greater acceleration as the optical depth drops, which leads
to larger total energy in the gas streaming away from the disk surface. The
${\rm{Be}}=0$ surface, which is the disk-wind boundary, also reflects our
results for the density scale height (Figure 3), and the disk shows
considerable changes in height as the accretion rate declines.

We present angular profiles in $\theta$ of the gas temperature, lorentz
factor, and asymptotic lorentz factor $\gamma_{\infty}$ as measured at
$r=30\,r_{g}$ in Figure 8. We compute the asymptotic lorentz factor from the
Bernoulli (Equation 21) by assuming all other forms of energy are converted to
kinetic energy far from the BH. Since $u_{t}\approx-\gamma$ at $r\gg r_{g}$,
we compute the asymptotic lorentz factor using:

$\gamma_{\infty}={\rm{Be}}+1.$ (31)

The initial gas temperature at $\lambda=2.8$ reaches $T_{\rm{gas}}\approx
10^{7}-10^{8}$ K, but when $\lambda=0.35$, the temperature has increased to
nearly $T_{\rm{gas}}\approx 10^{8}-10^{9}$ K. The outflowing gas in KIa00
achieves speeds in excess of $0.23c$ by $r=30\,r_{g}$, and the gas velocity
increases monotonically as the accretion rate declines, with the exception of
the $\sigma_{M}>1$ region which shows no clear trend across $\lambda$. We find
that when $\lambda=2.8$, the peak outflow velocity is $\approx 0.23c$, but
when $\lambda=0.35$ it is $\approx 0.33c$. The gas temperature in the jet
increases as the accretion rate shrinks. The asymptotic velocity also shows a
monotonic increase in the outflow velocity as $\lambda$ decreases aside from
near the poles in the $\sigma_{M}>1$ region. Despite the lack of BZ power in
KIa00, $\gamma_{\infty}$ implies the outflow may achieve a mildly relativistic
asymptotic velocity of $\sim 0.78-0.86c$.

In the case of KIa09, the pole is always optically thin down to the BH horizon
for the accretion rates we study in this work. In addition, the jet features
aren’t strongly affected by the accretion rate. The gas temperature in the jet
is generally $T_{\rm{gas}}>10^{8}$ K, with the hottest gas residing in the
$\sigma_{M}>1$ region where $T_{\rm{gas}}\sim 10^{9}-10^{10}$ K. The jet
achieves more relativistic velocities in this model compared to KIa00 (see
Figure 8) with a peak of $v\approx 0.65c$ when $\lambda=0.76-2.2$ (where the
jet power was observed to be highest). For $\lambda=0.29-0.44$ the peak
velocity reaches $v\approx 0.6c$; however, there is a noticeable decline to
$v\approx 0.4c$ at $\lambda=0.44$ which we discuss in the following paragraph.
The $\sigma_{M}=1$ boundary becomes visibly less wide over time. This agrees
well with the shrinking jet power. The profiles of $\gamma_{\infty}$
illustrate definitively that the relativistic jet does not shut off in KIa09.
Across all values of $\lambda$ that we sample, $\gamma_{\infty}$ reaches
$\gtrsim 30$. Unlike in KIa00, there is no clear trend in $\gamma_{\infty}$,
measured at $r=30\,r_{g}$, as $\lambda$ decreases for model KIa09. It is
possible that interaction with the wind as the jet expands outward could
modify the angular profiles of $\gamma_{\infty}$ as a function of radius and
slow some of the outflow into a slower moving jet sheath (as in Chatterjee et
al. 2019), but we are restricted to measuring the outflow properties at small
radii due to our implementation of a mass injection scheme.

When $\lambda=0.44$, the disk goes through two of its longest periods (at
$t\approx 23,500\,t_{g}$ and $t\approx 25,500\,t_{g}$) where $\phi_{\rm{BH}}$
reaches a minimum of $14-19$, substantially lower than the mean of $\approx
29$ (Table 2). Here $\phi_{\rm{BH}}$ takes longer to recover ($\sim
400\,t_{g}$ and $\sim 700\,t_{g}$, respectively) than during other magnetic
flux eruptions, where the magnetic flux only takes roughly $100-200\,t_{g}$ to
increase back to its mean value. The rather extreme nature of the magnetic
flux eruptions during this stage is well captured at $t=25,000\,t_{g}$ (3rd
panel from the left in Figure 4), where the accretion flow is essentially an
inner disk surrounded by low density gas where magnetic islands are orbiting
at $r\gtrsim 20\,r_{g}$. We show snapshots of the simulation when the jet
power is high versus low in Figure 9. The $\sigma_{M}=1$ boundary is closer to
the pole and the jet velocity is substantially lower when the jet power is
weak. Indeed, the instantaneous jet efficiency at these points is
$\eta_{\rm{jet}}\sim 3-7\%$. In addition, when the jet power temporarily drops
for extended periods of time, higher density gas launched by the disk wind
fills in the polar regions and mass loads the jet. As a result, the velocity
(in a time averaged sense) over the time range $t=23,000-26,000\,t_{g}$ is
noticeably smaller with a maximum of $v\approx 0.4c$ at $30\,r_{g}$. The
profiles of $\gamma_{\infty}$ at $\lambda=0.44$ suggest that at low jet power
states the highly relativistic component is not entirely shut off. Rather, a
smaller fraction of the jet gets accelerated to $\gamma_{\infty}\sim 30$
(between $\theta\sim 0-0.25$) while the mass loaded component at larger angles
reaches only $\gamma_{\infty}\lesssim 6$.

Note that the $\sigma_{M}=1$ boundary does not extend to large radii in KIa00
but closes off at $r\sim 20-80\,r_{g}$ (Fig. 6), whereas it extends to
infinity in KIa09 (Fig. 7). This is indicative of the lack of a magnetized
relativistic jet in KIa00. Hence, we only measure the jet profile in KIa09. To
measure the jet profile of KIa09, we treat the relativistic jet as the region
with $\sigma_{M}>1$ and take the polar boundary of the top/bottom jet to be
where $\sigma_{M}=1$. We assume that the profile takes the form $w\propto
z^{k}$, where $w=r\sin\theta_{\rm{jet}}$ is the jet radius,
$z=r\cos\theta_{\rm{jet}}$ is the jet height, and $k$ is a constant. We
measure the jet profile between $r=10-100\,r_{g}$ using time averages of
$\Delta t=3,000\,t_{g}$. This choice of radial limits is well motivated as the
jet profile within this range is approximately linear in $\log z-\log w$
space. We perform a linear fit in $\log z-\log w$ space to obtain the constant
$k$ as a function of time. As shown in Figure 10, the jet profile appears to
be initially slightly less parabolic for super-Eddington accretion rates with
$k\approx 0.58$ for both the top and bottom jet. As the accretion rate
declines and the jet power decreases, the jet settles into a very nearly
parabolic profile with $k\approx 0.43-0.49$ for the top jet and $k\approx
0.49-0.53$ for the bottom jet. The asymmetry between the top and bottom jet is
not unusual, as other GRRMHD simulations have observed similar asymmetrical
jets in BH accretion disks (e.g. see Morales Teixeira et al. 2018).

These results reflect a similar trend to the recently reported results of
Narayan et al. (2022) who also found that the jet parameter $k$ decreased as
the jet power decreased. However, our jet profile constants are slightly
higher in general, e.g. their $a_{*}=0.9$ model had $k\approx 0.43$. It is
interesting that our measured jet profiles are similar to those of large scale
AGN jets. For example, Kovalev et al. (2020) report jets with $k\approx$ 0.39
- 0.56 while Boccardi et al. (2021) find $k\approx$ 0.39-1.86. The jet profile
of KIa09 is comparable to 3C 120, which is a high state AGN whose Eddington
ratio may be as high as $\lambda\approx 0.3$ (Chatterjee et al., 2009) with a
profile at small radii ($r<1$ parsec) with $k\approx 0.586\pm 0.047$ (Kovalev
et al., 2020). As was noted by Narayan et al. (2022), the observed jet
profiles of AGN are derived from emission produced by the jet sheath, which
likely does not coincide with $\sigma_{M}=1$. Thus, our jet profiles may be
narrower than those inferred directly from observations. Ray traced radio
images may better quantify jet profiles from GRRMHD simulations of thin disks,
but this is beyond the scope of the present work.

Our results suggest that jets from a thin, highly magnetized accretion disk
can explain the jet profile at small radii seen in some AGN. Our model with
$\langle\phi_{\rm{BH}}\rangle\approx 30$ tends towards a slightly less
magnetized disk than geometrically thick, MAD ADAF models, which have
$\langle\phi_{\rm{BH}}\rangle\approx 50$ (e.g. Narayan et al. 2022). In
addition, our simulation differs from other works which have applied accretion
disk simulations to study the structure of AGN jets (e.g. Chatterjee et al.
2019) in that we do not have a large scale torus which may provide confinement
at larger radii. Instead, all of the confinement happens due to the wind
launched from the disk.

An important point to note is that our simulation setup forces the system to
be in the MAD state. It is suspected that thin disks in AGNs will allow
magnetic flux to diffuse away from the BH, thereby preventing the high
magnetic flux that we report in this work (Lubow et al., 1994; Heyvaerts et
al., 1996). This is in conflict with the presence of jets in AGN with
luminosities corresponding to the thin disk regime, which suggests that some
of these AGN accretion disks are capable of advecting magnetic field towards
the BH. For AGN that do maintain a significant magnetic flux, the total jet
and wind power will be non-negligible and confinement similar to what we find
in this work may be expected. Further simulations of steady state, thin, MAD
accretion disks are warranted.

###  3.3 Variability

MAD accretion disks are by nature highly variable and have been cited as
potential sources of quasi-periodic oscillations in AGN and other
astrophysical sources with BH accretion as the central engine (Tchekhovskoy et
al., 2012, 2014). In general, it would be ideal to study variability from
physically motivated simulations such as presented in this work. As was noted
by Sądowski & Narayan (2016), while numerical simulations can provide detailed
information on the radiation field, the accuracy of the radiation field is
affected by the numerical schemes used and the treatment of radiation. In
particular, the M1 closure scheme which is applied by KORAL assumes a gray
approximation in the radiative transfer and a simple angular distribution of
the radiation field (Sądowski et al., 2013).

With those complications in mind, we attempt to measure the variability
directly from the KORAL simulation. One would ideally measure the radiative
flux at a large distance (i.e. at the outer boundary of the simulation domain)
as these rays would be expected to reach a distant observer. This is not
accurate though if the system is not in equilibrium out to the radius at which
the flux is measured. If the radiation is decoupled from the gas the radiation
field should show similar variability at each radius with a time shift due to
the difference in the light crossing time. The present simulations are
particularly complicated because although the jet and magnetized wind are both
optically thin, the disk properties are changing on roughly the viscous time
so, by construction, the system is never in equilibrium.

Figure 11: Here we arbitrarily scale the accretion rate $\dot{M}$, bolometric
luminosity $L_{\rm{bol}}$, and magnetic flux at the BH horizon
$\Phi_{\rm{BH}}$ to compare their variability in KIa00 (top) and KIa09
(bottom). The accretion rate, bolometric luminosity, and magnetic flux appear
to show similar variability (e.g. we observe large dips at similar times);
however, visually, the bolometric luminosity and magnetic flux show the
strongest correlation. Figure 12: Here we show the fraction of the disk
surface area within $r\leq 10r_{g}$ that is optically thick $a_{\rm{eff}}$
over time in KIa00 and KIa09. We also scale the normalized magnetic flux at
the BH horizon $\phi_{\rm{BH}}$ for ease of comparison between the two
parameters. In both simulations, we find that the area of the optically thick
gas shows similar variability to the magnetic flux. The drop in $a_{\rm{eff}}$
appears to roughly coincide with drops in the magnetic flux due to the low
density regions caused by flux eruptions.

As shown in Figure 11, the variability of the bolometric luminosity in both
KIa00 and KIa09 is clearly also matched by variability in the magnetic flux.
We also note that KIa09 appears to show slightly more coherent quasi-periods
after $t>25,000\,t_{g}$. The largest dips appear in both the radiative flux
and magnetic flux, and are separated by quite long times ($\delta t\approx
2000\,t_{g}$, or nearly $10^{4}$ s for our $10^{6}M_{\odot}$ BH). The
variability in the luminosity results in a decay to nearly 25% of the peak
luminosity at the minima. This suggests that MAD accretion flows may also be
sources of longer period QPOs, in addition to the often cited short period
QPOs associated with the spin frequency at the BH horizon (e.g. Tchekhovskoy
et al. 2012, 2014).

As the apparent correlation between $\Phi_{\rm{BH}}$ and $L_{\rm{bol}}$ is
present in both KIa00 and KIa09, it is possible that the reduction in emitting
area of the disk (where the disk is optically thick at the midplane) is what
causes drops in the bolometric luminosity. To test this hypothesis, we compute
the fraction of the disk surface area that is optically thick between $r_{H}$
and $10\,r_{g}$ as

$a_{\rm{eff}}=\dfrac{\int_{r_{H}}^{10r_{g}}\int_{0}^{2\pi}f(\tau_{\theta})r\,dr\,d\phi}{\int_{r_{H}}^{10r_{g}}\int_{0}^{2\pi}r\,dr\,d\phi},$
(32)

where

$\tau_{\theta}=\int_{0}^{\pi/2}\dfrac{\rho\kappa_{\rm{es}}}{c}u^{t}\sqrt{g_{\theta\theta}}\,d\theta$
(33)

is the angle integrated electron scattering optical depth at the midplane and
$f(\tau_{\theta})=1$ if $\tau_{\theta}>1$ and $f(\tau_{\theta})=0$ otherwise.
We only integrate up to $r=10\,r_{g}$ since the disk tends to become optically
thick at larger radii, which would tend to damp oscillations in $a_{\rm{eff}}$
near the horizon if included in the integral. As we show in Figure 12,
$a_{\rm{eff}}$ generally decreases around the same time that $\phi_{\rm{BH}}$
decreases due to magnetic islands disrupting the inner disk.

Model | Variable 1 | Variable 2 | Correlation
---|---|---|---
KIa00 | $L_{\rm{bol}}$ | $\Phi_{\rm{BH}}$ | $0.59\pm 0.14$
| $L_{\rm{bol}}$ | $\dot{M}$ | $0.24\pm 0.18$
| $\Phi_{\rm{BH}}$ | $\dot{M}$ | $0.56\pm 0.22$
| $L_{\rm{bol}}$ | $a_{\rm{eff}}$ | $0.25\pm 0.15$
| $\Phi_{\rm{BH}}$ | $a_{\rm{eff}}$ | $0.38\pm 0.13$
KIa09 | $L_{\rm{bol}}$ | $\Phi_{\rm{BH}}$ | $0.75\pm 0.13$
| $L_{\rm{bol}}$ | $\dot{M}$ | $0.32\pm 0.17$
| $\Phi_{\rm{BH}}$ | $\dot{M}$ | $0.37\pm 0.18$
| $L_{\rm{bol}}$ | $a_{\rm{eff}}$ | $0.25\pm 0.12$
| $\Phi_{\rm{BH}}$ | $a_{\rm{eff}}$ | $0.38\pm 0.12$
Table 3: Cross-correlation coefficients (CCCs) of various quantities in which
we identify significant variability.

To quantify the correlation between variability in the above quantities, we
compute the cross-correlation between the variables $L_{\rm{bol}}$,
$\Phi_{\rm{BH}}$, $\dot{M}$, and $a_{\rm{eff}}$. For this analysis, we use
data from $t=15,000-29,000\,t_{g}$ for KIa00 and $t=15,000-34,000\,t_{g}$ for
KIa09. We detrend the data using a moving average so as to leave only the
short time scale variability. We then perform a stationary block bootsrap
cross-correlation (e.g. see Effron 1979) using $13-14$ blocks333Our chosen
number of blocks is on the low end since we also had the additional
requirement of preserving the auto-correlation. One would typically desire
$\geq 30$ blocks, where the Central Limit Theorem is satisfied for both
Gaussian and non-Gaussian sample distributions; however, we found that 13
blocks was adequate to obtain approximately Gaussian distributions of the
statistic. of equal duration in time for each model. We performed 5000 re-
samples with replacement.

We report the cross-correlation coefficient (CCC) and errors in Table 3. We
saved snapshots of the simulation every $\delta t=100\,t_{g}$, so we do not
report lags since the sampling rate is longer than the minimum lag we expect
between the mass accretion rate and luminosity, which is roughly the light
crossing time $\sim 30\,t_{g}$ for emission near the BH horizon. Although it
is commonplace in statistical analysis to define a range for the CCC as
strong, moderate, or weak, we do not adopt such a definition here as it is
somewhat arbitrary. We find it more useful to consider the relative value of
the CCC.

For KIa00, we find a significant correlation between $L_{\rm{bol}}$ and
$\Phi_{\rm{BH}}$ with a CCC of $0.59\pm 0.14$ and between $\Phi_{\rm{BH}}$ and
$\dot{M}$ with a CCC of $0.56\pm 0.22$. However, $L_{\rm{bol}}$ and $\dot{M}$
show a weak correlation with a CCC of $0.24\pm 0.18$. For KIa09, the cross-
correlation indeed finds that $L_{\rm{bol}}$ and $\Phi_{\rm{BH}}$ are strongly
correlated with a peak CCC of $0.75\pm 0.13$. Meanwhile, there is a much
weaker correlation between $\dot{M}$ and $L_{\rm{bol}}$ (and
$\Phi_{\rm{BH}}$).

For both KIa00 and KIa09, we find that $\Phi_{\rm{BH}}$ and $a_{\rm{eff}}$
have a CCC of $0.38$; however, $L_{\rm{bol}}$ appears to be very weakly
correlated with $a_{\rm{eff}}$ with a CCC of $0.25$. This suggests that the
variability is not caused by thermal emission from the disk. It is possible
that the variability in $L_{\rm{bol}}$ may instead be due to the flux
eruptions modifying the emission properties in the jet/wind, but a more in
depth investigation of the emission properties of the simulation are needed in
order to test this.

Since increasing the $\sigma_{M}$ ceiling also increases $L_{\rm{bol}}$ in the
jet region, the larger correlation between $L_{\rm{bol}}$ and $\Phi_{\rm{BH}}$
in KIa09 when compared to KIa00 may be due to the ceiling model. However, the
fact that KIa00, which did not have a jet, is less effected by the ceiling on
$\sigma_{M}$ (Figure 18) suggests that our finding of a high CCC for
$L_{\rm{bol}}$ and $\Phi_{\rm{BH}}$ in both models, when contrasted with the
CCC between other variables, is a robust result.

Figure 13: Here we show time-averaged specific angular momentum flux into the BH $j$ (top panel, see Equations 26, 27) and the dimensionless spin-up parameter $s$ (bottom panel, see Equation 28) as a function of the accretion rate $\lambda$ for spin $a_{*}=0$ (KIa00) and $a_{*}=0.9$ (KIa09). Model | Time | $\lambda$ | $j$ | $s$
---|---|---|---|---
| $(10^{3}\,t_{g})$ |  |  |
KIa00 | 15-18 | 2.8 | 1.70 | 1.70
| 18-21 | 1.2 | 1.91 | 1.91
| 22-25 | 0.58 | 1.79 | 1.79
| 26-29 | 0.35 | 1.84 | 1.84
NT $(a_{*}=0)$ | - | - | 3.46 | 3.46
KIa09 | 14-17 | 2.2 | -2.80 | -3.02
| 20-23 | 0.76 | -2.14 | -2.11
| 23-26 | 0.44 | -1.52 | -1.55
| 28-31 | 0.31 | -0.66 | -1.23
| 31-34 | 0.29 | -0.64 | -1.29
NT $(a_{*}=0.9)$ | - | - | 2.10 | 0.58
Table 4: Tabulated values of the specific angular momentum flux into the BH
($j$), and spin-up parameter ($s$) for each model as a function of time and
accretion rate ($\lambda$).

##  4 Spin Evolution

We show the specific angular momentum $j$ of the accreted gas and the
dimensionless BH spin-up parameter $s$ of each model in Figure 13. We also
present tabulated values of $j$ and $s$ as a function of time and $\lambda$ in
Table 4. For KIa00, we find that both $j$ and $s$ are positive, meaning that
angular momentum and net spin are being gained by the BH. Furthermore, both
quantities are nearly independent of the accretion rate. In the case of KIa09,
both $j$ and $s$ have the opposite sign, indicating that angular momentum is
being lost by the BH. Thus, rapidly spinning BHs, even with thin disks, spin
down if accretion occurs in the MAD state. For comparison, the traditional NT
disk spins up the BH for all $a_{*}$ values up to $0.998$ (Thorne, 1974). Note
that $s$ is correlated with $\lambda$ in KIa09 and becomes less negative with
decreasing $\lambda$. This suggests that, if $\lambda$ is sufficiently small,
the BH may switch from spin-down to spin-up.

Figure 14: Here we estimate the time for the spin $a_{*}$ of the BH to
increase/decrease by $0.1$ as a function of accretion rate $\lambda$ for spin
$a_{*}=0$ (KIa00) and $a_{*}=0.9$ (KIa09).

An interesting question to consider is the time evolution of the BH spin. We
follow the analysis in Narayan et al. (2022), which is based on Shapiro
(2005), and calculate the spin evolution via:

$\dfrac{\Delta a_{*}}{\Delta
t}=\lambda\dfrac{\dot{M}_{\rm{Edd}}}{M}s(\lambda).$ (34)

To calculate the change in spin over a time interval $\Delta t$, we assume
$\lambda$ and $s$ are constant. Figure 14 shows the time $\Delta t_{0.1}$ for
the spin to change by $|\Delta a_{*}|=0.1$ as a function of $\lambda$ for the
two models. The time required decreases as the accretion rate increases, as
expected. For KIa00, this change is approximately inversely proportional to
$\lambda$ since $s(\lambda)$ is nearly constant. For KIa09, the dependence is
steeper; $\Delta t_{0.1}$ is nearly the same as for KIa00 for large $\lambda$,
but 3 times longer for $\lambda\sim 0.3$ because $s$ decreases by a factor of
3. For $\lambda>1$, the time for the BH to spin up/down for both models is
roughly 1 Myr. When $\lambda\sim 0.3$, $\Delta t_{0.1}$ is $\approx 4$ Myr for
KIa00 and $\approx 15$ Myr for KIa09. These times are much longer than the
duration of a TDE, and hence not relevant for these events. The results are,
however, of great interest for AGN evolution in general, and we turn next to
this topic.

Since KIa00 has $s>0$ and KIa09 has $s<0$, there must be a spin value
$a_{*,\rm eq}$ in between $a_{*}=0$ and 0.9 where a thin MAD disk will reach
spin equilibrium with neither spin-up nor spin-down. We obtain a crude
estimate of the equilibrium spin $a_{*,{\rm{eq}}}$ using the following
approximations: (i) we assume that, for a given $\lambda$, $j$ varies with BH
spin as $j=c_{1}+c_{2}a_{*}^{2}$, where $c_{1}$ and $c_{2}$ are constants;
this approximation is motivated by the fact that angular momentum loss via the
BZ mechanism should scale as $a_{*}^{2}$, given that the BZ power varies as
$\omega_{H}^{2}\sim a_{*}^{2}$ (see Eq 1), (ii) we assume that
$\eta_{\rm{tot}}=c_{3}+c_{4}a_{*}^{2}$ for similar reasons.

We consider three accretion rates ranging from super-Eddington to sub-
Eddington: $\lambda=1.2,0.58,0.35$. For these values of $\lambda$, we have
estimates of $j$ and $\eta_{\rm tot}$ directly from Table 2 for KIa00. In the
case of KIa09, we interpolate $j$ and $\eta_{\rm{tot}}$ from the tabulated
values. We then compute where $s(\lambda,a_{*})=0$ using fits to
$j(\lambda,a_{*})$ and $\eta_{\rm{tot}}(\lambda,a_{*})$ vs $a_{*}$, using the
assumed $a_{*}$ dependence described in the previous paragraph. We thereby
estimate that the equilibrium spin lies roughly in the range
$0.5<a_{*,{\rm{eq}}}<0.6$ for $\lambda=0.35-1.2$. The fact that the spin-up
parameter of Avara et al. (2016), whose model was also thin and MAD with
$a_{*}=0.5$, was only $s\approx-0.4$ suggests that our estimate is not too far
from what a full parameter space study may find. Our estimated equilibrium
spin is approximately half the spin value expected with a traditional NT thin
disk model, which is expected to drive a BH to an equilibrium spin of
$a_{*,{\rm{eq}}}=0.998$ (Thorne, 1974).

Eddington-level or mildly sub-Eddington accretion during a SMBHs accretion
history is often cited to explain the occurrence of very massive SMBHs at high
redshifts $z\sim 6$ (Haiman, 2004; Yoo & Miralda-Escudé, 2004; Shapiro, 2005;
Volonteri & Rees, 2005, 2006; Pelupessy et al., 2007; Tanaka & Haiman, 2009;
Madau et al., 2014; Pacucci et al., 2015; Lupi et al., 2016). Our models
suggest that MAD accretion disks, even for thin disks, may play an important
role in the spin evolution of these BHs. The time that we estimate for a BH to
spin up or down is not substantial even for the sub-Eddington stages. For
example, even assuming a modest duty cycle $\mathcal{D}\approx 0.01$ for the
BH’s accretion history and Eddington-scaled accretion rates $\lambda\approx
0.1-1$, the time for the BH to spin up or down to its equilibrium spin is only
of the order of a few Gyr (ignoring other contributions to spin evolution such
as mergers).

##  5 Discussion

###  5.1 Comparison with other thin disk simulations

In recent years, studies of radiative GRMHD simulations have begun to explore
$\dot{M}/\dot{M}_{\rm{Edd}}<1$ accretion disks. In this regime, it was
initially unclear if thermally stable solutions with moderate vertical extent
existed, but several studies have shown that stability is indeed possible at
such accretion rates (Sądowski, 2016; Avara et al., 2016; Morales Teixeira et
al., 2018; Lančová et al., 2019; Mishra et al., 2020; Liska et al., 2022).
Sądowski (2016) was the first to demonstrate that stable disks can occur at
$\lambda\approx 0.8$ if a strong magnetic field threads the disk. This initial
solution was turbulent with a density scale height of $0.15$ while the
photosphere height was substantially larger due to magnetic pressure support.
Accretion occurred both through the much more dense disk (within a density
scale height) and the optically thick, but less dense gas inside the
photosphere.

Lančová et al. (2019) extended the range of simulations of magnetically
supported disks to $\lambda\approx 0.6$. They again found that a stable,
turbulent disk formed with accretion occuring both through the denser disk at
the midplane and above/below the disk in magnetically supported gas. The
density scale height was $\approx 0.1$ while the photosphere height was $H\sim
r$. They named this mode of accretion a ‘Puffy Disk’. Both Sądowski (2016) and
Lančová et al. (2019) found that the magnetic pressure was nearly 50-60% of
the total pressure and radiation pressure was far greater than the gas
pressure. Additionally, the radiative efficiency of these solutions is similar
to that of a Novikov & Thorne (1973) thin disk solution. We note that both
authors employed non-rotating BHs with SANE accretion disks.

Figure 15: Here we show time-$\phi$ averaged radiation energy density (colors)
and radiative flux (arrows) for both KIa00 (top) and KIa09 (bottom) at
Eddington ratios similar to the ‘Puffy Disk’ models presented in Sądowski
(2016); Lančová et al. (2019). The electron scattering photosphere (yellow
line), and magnetized jet boundary ($\sigma_{M}=1$, blue line) are indicated
as well.

In comparing our $a_{*}=0$ model KIa00 with the models presented in Sądowski
(2016); Lančová et al. (2019), we find that our thin MAD solution is
significantly more efficient, with $\eta_{\rm{rad}}\approx 11.6-12.2\%$ for
sub-Eddington phases of the evolution. The density scale height of KIa00
during the sub-Eddington stages ($H/R\approx 0.15-0.17$) is also slightly
larger in general. We also find a considerable difference in the behaviour of
the radiative flux (see Figure 15). Both KIa00 and KIa09 exhibit radiative
flux which escapes nearly vertically from the mid-plane even below the
electron scattering photosphere. This is the opposite of the behaviour in
Sądowski (2016); Lančová et al. (2019), who both found that the radiative flux
was advected towards the BH beneath the photosphere. This difference is
possibly due to the different angular momentum transport mechanisms, as
accretion in MAD disks is driven by magnetic Rayleigh-Taylor instabilities
(but see Begelman et al. 2022 who argue instead that radial
convective/interchange instabilities dictate the flow structure in MADs)
whereas SANE models such as the ‘Puffy Disk’ model accrete through MRI.

It is interesting that our radiative efficiency is higher, but there are
important differences between our models and the ‘Puffy Disk’. For instance,
they found that additional radiation from heightened viscous heating inside
the ISCO cannot escape due to radial advection inward, which does not appear
to be the case in our models. Also, the photosphere scale height in our models
is significantly smaller than $H\sim r$, which implies that the conversion of
radiative energy into kinetic energy will be less efficient than in a ’Puffy
Disk’, so more radiative energy will ultimately escape. The magnetized wind
carries gas which holds a significant amount of radiation from the disk with
it and it seems that unlike in the super-Eddington case (see McKinney et al.
2015), this radiation escapes freely with very little absorption due to the
low density of the wind. We note that for our model KIa09 the spin is slightly
higher ($a_{*}=0.9$ versus $a_{*}=0.8$) than the simulation of McKinney et al.
(2015). Furthermore, the fact that we must measure the energy outflow at small
radii ($r\geq 30\,r_{g}$) may impact our results. For instance, McKinney et
al. (2015) cite an order of two decrease in their radiative efficiency when
they account for the energy conversion out to large radii.

Morales Teixeira et al. (2018) presented radiative GRMHD simulations of thin,
MAD accretion disks around rotating ($a_{*}=0.5$) BHs. They simulated various
disks spanning accretion rates from $\lambda\approx 0.2-0.5$ and confirmed
that their disks were stable during the simulated time. They find total
efficiencies of $\sim 15-20$% with the lowest efficiency corresponding to the
lowest accretion rate model and vice versa, which is substantially larger than
the Novikov & Thorne (1973) thin disk efficiency of 8.2%. They also cite
radiative efficiencies of $\sim 2.4-2.9\%$, making their models rather
radiatively inefficient (that is, most of their large total efficiency is in
mechanical energy outflow).

Our results for KIa09 are similar to the findings of Morales Teixeira et al.
(2018) in the sense that the extraction of spin energy leads to a total
efficiency that is significantly larger than that of the NT thin disk
efficiency even at highly sub-Eddington rates. However, we observe a large
difference with regards to the photsophere scale height and the overall energy
output. Morales Teixeira et al. (2018) found a photosphere height that was
significantly larger than the density scale height, which suggests that
radiation energy would have been efficiently converted into kinetic energy.
This appears to be true of their simulations as the sum of their reported jet
and wind efficiency for their model RADHR ($14.7\%$) is nearly five times
larger than the radiative efficiency ($2.9\%$). It is worth noting that the
disk at $r\gtrsim 16\,r_{g}$ exhibited turbulence in all of their models since
the disk was only MAD up to $\sim 16\,r_{g}$. The application of disk
injection at $R_{\rm{inj}}=40\,r_{g}$ in our simulations, and the fact that
our entire disk below the injection radius becomes MAD, may be a possible
explanation for why our simulation behaves differently. Perhaps the extent to
which the disk is MAD plays a role in the photosphere scale height and the
total radiative efficiency in thin MAD disks.

A recent study of truncated disks by Liska et al. (2022) provides a similar
model to KIa09. Their numerical model was run with $a_{*}=0.9375$,
$\lambda\approx 0.35$, and the magnetic flux at the BH horizon saturated at
$\phi_{\rm{BH}}\approx 30$, as in the case of our simulations. However, the
magnetic pressure supported inner disk is somewhat thinner in their
simulations with $H/R\sim 0.05-0.2$. They also evolved simulations using a two
temperature electron model, which leads to cooler electrons in the jet. They
confirm the possibility of large total efficiency in a thin MAD. They find
$\eta_{\rm{tot}}\sim 90\%$, which is not very different from the values
$\eta_{\rm{tot}}\approx 70-100\%$ that we find in model KIa09. However, they
find a smaller radiative efficiency $\eta_{\rm{rad}}\sim 20\%$ compared to our
$\eta_{\rm{rad}}\sim 50\%$ at $\lambda\approx 0.3$. The results of Kinch et
al. (2020) suggest that the rather large difference in $\eta_{\rm{rad}}$
between our work and Liska et al. (2022) is possibly due to the fact that they
modeled a two temperature plasma, whereas our simulations assume a single
temperature. However, we also demonstrated that the numerical ceiling on
$\sigma_{M}$ impacts the radiative efficiency in KIa09. In addition, our
choice of $\sigma_{M}\leq 100$ is larger than the maximum values of
$\sigma_{M}\leq 25$ and $50$ implemented by Liska et al. (2022). An increase
in radiative efficiency with larger $\sigma_{M}$ was also shown in their work
(see Figure 8 in Liska et al. 2022). Therefore, we believe the different
numerical floors and the inclusion of a two-temperature plasma can account for
the factor of $\gtrsim 2$ difference in our radiative efficiency when compared
with Liska et al. (2022).

###  5.2 Variability

Simulations of MADs (e.g. Igumenshchev 2008) show significant variability in
the mass accretion rate. Because of magnetic flux eruptions, a quasi-periodic
accretion rate (and thus light curve) may be expected. Although the gas can
slip past the magnetic field through 3D RT instabilities and accrete, there
may be a connection between variability and the radius to which the disk is
MAD since the radius at which magnetic buoyancy can support the gas grows with
magnetization.

Our simulations are unique in comparison to other MAD disk simulations in the
sense that we inject a Keplerian disk which is forced to be MAD. Due to the
boundary conditions implemented in this work, the disk can only be MAD up to
$R_{\rm{inj}}$ by construction. The fact that we found structured variability
in the bolometric luminosity for KIa09 that closely matched that of the
magnetic flux suggests that quasi-periods linked to flux eruptions may appear
in thin MAD accretion disks.

If MAD disks can maintain large scale poloidal fields up to radii beyond the
horizon, our model KIa09 suggests quasi-periods of the order $2\times
10^{3}\,t_{g}$, accompanied by a luminosity variation with minima of $\lesssim
25\%$ of the peak value, may be possible. This could have implications for
well-known BH X-ray binaries (BHXBs) that show a variety of time-variable
states, e.g., GRS 1915+105 (Belloni et al., 2000) as well as some TDEs.

Type-B QPOs are long frequency QPOs in BHXBs that coincide with the soft
intermediate state (SIMS) and appear with a typical frequency ranging from
$4-6$ Hz and fractional root mean squared (rms) amplitude of less than $10\%$
(Belloni 2010). The SIMS and the occurrence of type-B QPOs is possibly
associated with the presence of transient jets (Fender et al., 2004; Corbel et
al., 2005), which suggests our model KIa09 may be analogous to the SIMS.
Scaling the BH mass in KIa09 to $10\,M_{\odot}$ suggests a QPO frequency of
$\sim 10$ Hz for a stellar mass BH analog, which is remarkably similar to the
frequency range of type-B QPOs. We note, however, that the fractional rms
amplitude in $L_{\rm{bol}}$ is $\sim 50\%$ in KIa09, which suggests that our
model may be too violently variable to explain type-B QPOs.

The fact that the variability period is $\delta t\approx 2\times
10^{3}\,t_{g}$ and is similar to $T_{\rm{acc}}(R_{\rm{inj}})\approx
2500\,t_{g}$ suggests a connection between the period and the truncation
radius of the MAD disk. However, we did not explore the effects of different
choices of $R_{\rm{inj}}$ on the variability in this work. Assuming
$M_{\rm{BH}}=10\,M_{\odot}$, a model where $\delta t\sim T_{\rm{acc}}$ would
require the radial edge of the MAD disk to extend to larger radii in order to
explain lower frequency QPOs.

The most well studied TDE which is thought to have gone MAD is Swift J1644+57
(Tchekhovskoy et al., 2014), which showed variability on a range of timescales
with both short period QPOs at $\sim 200$ s (Reis et al., 2012) and long
period dips in the light curve on time scales of $\sim 10^{6}$ s (Saxton et
al., 2012). The QPO period in KIa09 is $\approx 10^{4}$ s, but scaling $\delta
t$ with $M_{\rm{BH}}$ and $R_{\rm{inj}}$ allows for longer periods capable of
explaining the $10^{6}$ s variability in Swift J1644+57. In addition, the flux
dips of the long period variability had minima of $\lesssim 50\%$ of the peak
value, which is similar in magnitude to the dips that we find in KIa09.

One complication with the hypothesis that a MAD disk may drive longer period
QPOs in the context of TDEs is that the magnetic flux required for the disk to
be MAD beyond the horizon is quite large. This is a substantial astrophysical
hurdle that makes the MAD scenario for TDEs difficult to begin with since the
flux provided by the star is too small (Kelley et al., 2014). Nevertheless,
given the success of MAD disks in explaining events such as Swift J1644+57,
one can assume substantial flux was indeed brought in, whether from a fossil
disk or from a dynamo process.

###  5.3 Implications for Swift J1644+57

Our model KIa09 is the first GRRMHD simulation of a sub-Eddington MAD disk
around a SMBH. Our results demonstrate that a substantial jet power efficiency
of $\eta_{\rm{jet}}\sim 18\%$ can be expected even when $\lambda\sim 0.3$,
provided the disk is MAD and the BH spins rapidly. At such low accretion
rates, model KIa09 has a lower scaled magnetic flux of $\phi_{\rm BH}\sim
25-30$ compared to $\phi_{\rm BH}\sim 50$ for highly super-Eddington accretion
rates (Curd & Narayan, 2019). This does cause a modest decrease in jet
efficiency (see the BZ power scaling in Eq 1) at the lowest sub-Eddington
rates that we have explored. However, the reduction in $\eta_{\rm{jet}}$ is
not sufficient to explain the apparent shut-off of the jet in Swift J1644+57
assuming the fallback rate after $\sim 500$ days was similar to the minimum
mass accretion rate we probe in our simulations.

One possibility is that the disk in Swift J1644+57 did maintain enough
magnetic flux to remain MAD, but the disk became much thinner than in our
simulation. The fitting formula provided by Avara et al. (2016) suggests an
association between $H/R$ and the extracted jet power:
$\eta_{\rm{jet}}\propto(H/R)^{2}$. This is also suggested by the BZ power in
Equation (2) since various authors demonstrated through numerical simulations
that $\phi_{\rm{BH}}\propto H/R$ for thick accretion disks (Penna et al.,
2010; Tchekhovskoy et al., 2011, 2012).

Based on the disk height evolution in KIa09 (see the bottom panel in Figure
5), the above scenario would only be realized at even lower accretion rates
than we have considered. In this context, the accretion rate at the time that
the X-ray flux shut off in Swift J1644+57 may have been significantly smaller
than the cutoff mass accretion rate of $\lambda\sim 0.3$ assumed by
Tchekhovskoy et al. (2014).

For instance, if we consider $H/R\sim 0.03$, which is nearly an order of
magnitude smaller than the disk height at $\lambda=0.29$ in KIa09, the jet
efficiency assuming $\eta_{\rm{jet}}\propto(H/R)^{2}$ is only $\sim 2\%$. The
trend that we find in $H/R$ with decreasing $\lambda$ suggests
$(H/R)\propto\lambda^{0.5}$ for $\lambda=0.29-0.44$. Assuming this behavior
holds for yet smaller accretion rates, the jet power may decrease to
$\eta_{\rm{jet}}\lesssim 2\%$ if $\lambda\lesssim 0.02$. A potential
difficulty with this scenario is that, while a lower accretion rate and
correspondingly smaller $H/R$ predicts a steep reduction in the jet power, we
expect the mass accretion rate itself to decline smoothly with time in a TDE.
The tidal stream is assumed to feed mass into an accretion disk at a rate
$\dot{M}_{\rm{fb}}\propto(t/t_{\rm{fb}})^{-5/3}$, where $\dot{M}_{\rm{fb}}$ is
the fallback rate and $t_{\rm{fb}}$ is the fallback time. Assuming
$\lambda\propto(t/t_{\rm{fb}})^{-5/3}$, it would take nearly
$4t_{\rm{fb}}\approx 124$ days for $\lambda$ to decrease from $0.3$ to $0.02$
if $M_{\rm{BH}}=10^{6}\,M_{\odot}$. It is hard to see how this can lead to as
abrupt a decline in the X-ray flux as was seen in Swift J1644+57 since the
X-ray flux declined by a factor of 15 in only 25 days (Sbarufatti et al.,
2012), which is substantially shorter than $4t_{\rm{fb}}\approx 124$ days.

We caution that the $\eta_{\rm{jet}}$ fits in Avara et al. (2016) that we used
in the above calculation were only provided over data sampling $H/R\gtrsim
0.05$ so we have assumed that $\eta_{\rm{jet}}\propto(H/R)^{2}$ at lower disk
scale heights. Also of note is the fact that Liska et al. (2022) report
$\eta_{\rm{jet}}\sim 50\%$ for a disk with $H/R\sim 0.05$, which suggests that
the correlation of $\eta_{\rm{jet}}$ with $H/R$ may not hold when MAD disks
are simulated with extreme resolution. Liska et al. (2022) implement a grid
with substantially higher resolution in the disk midplane (more than a factor
of 10 in $\theta$) compared to simulations presented in both Avara et al.
(2016) and this work.

The other possibility is that the magnetic flux diffuses away from the BH
horizon to larger radii once $\lambda$ falls below some threshold value. That
is, the MAD state is no longer permitted. Such diffusion has not been seen in
any simulation so far. However, it is possible that diffusion occurs on a
longer time scale than the simulations can explore. For instance, KIa09 was
run for a mere $\sim 2$ days (and spent $\sim 0.3$ days at $\lambda\lesssim
0.3$) while the X-ray flux in Swift J1644+57 declined by a factor of 15 in
$\sim 25$ days (Sbarufatti et al., 2012). In this scenario, the disk would
switch from MAD to SANE below a critical $\lambda$ and the jet and its
associated X-ray emission would be lost. Considering that all stable sub-
Eddington $\lambda\lesssim 1$ disk simulations so far have included strong
magnetic fields (even the ‘Puffy Disk’ model which is technically not MAD),
one wonders whether the loss of magnetic field might cause the disk to become
thermally unstable. If it does, the disk will likely collapse quickly to a
very small scale height since magnetic pressure support would be lost. Since
the time scale for the disk to collapse due to becoming thermally unstable is
short, this scenario is currently the most favorable for explaining the sudden
decline in the X-ray in Swift J1644+57.

In summation, we believe that diffusion of the magnetic flux away from the
horizon and the disk no longer being MAD must coincide with the disk becoming
modestly sub-Eddington to explain the X-ray shut off in Swift J1644+57 (and
possibly Swift J2058+05). It is conceivable that the required magnetic flux
diffusion would become evident with longer simulations, and we plan to
consider this in a future study.

##  6 Conclusions

We have simulated the evolution of a small scale, thin, MAD accretion disk
around both non-spinning (model KIa00) and spinning (model KIa09) SMBHs and
tracked the evolution from mildly super-Eddington to sub-Eddington mass
accretion rates. Our conclusions are as follows:

  1. 1.

Throughout the evolution, the disk is stable against thermal collapse because
of the large magnetic pressure in the disk.

  2. 2.

Both of our models are highly radiatively efficient during the sub-Eddington
phase. The spin $a_{*}=0$ model KIa00 reaches $\eta_{\rm{rad}}\approx
11.6-12.2\%$ and the spin $a_{*}=0.9$ model KIa09 reaches
$\eta_{\rm{rad}}\approx 45.4-50.2\%$. These are significantly larger than the
NT thin disk efficincies $\eta_{\rm{rad}}=5.72\%$ and
$\eta_{\rm{rad}}=15.58\%$, respectively, for the two spins. The total
efficiency also shows significant deviations from NT during the sub-Eddington
phase, with $\eta_{\rm{tot}}\approx 15.1-15.9\%$ for KIa00 and
$\eta_{\rm{tot}}\approx 64.3-68.3\%$ for KIa09.

  3. 3.

The jet efficiency in KIa09 decreases from $\eta_{\rm{jet}}\approx 38.7\%$ to
$\eta_{\rm{jet}}\approx 17.7-18.4\%$ as the disk becomes sub-Eddington,
meaning that, even at low accretion rates, there is substantial jet power in a
mildly relativistic jet driven by the extraction of spin energy. We only
measure the lorentz fact of the jet at $r=30\,r_{g}$ due to the our
application of the mass injection scheme used to vary the mass accretion rate
over time. However, the asymptotic lorentz factor suggests that the outflow
may reach highly relativistic speeds of $\gamma_{\infty}\gtrsim 30$ at $r\gg
r_{g}$ across all $\lambda$. This jet is powered by the relatively strong
magnetic flux $\phi_{\rm{BH}}\approx 30$ in spite of the decreasing accretion
rate and disk scale height.

  4. 4.

For KIa09, we find a best fit jet profile power-law $k\approx 0.43-0.53$,
which is approximately parabolic, within $r<100\,r_{g}$ when $\lambda<1$. This
is consistent with measurements of high state AGN jets on the sub-parsec
scale.

  5. 5.

Both KIa00 and KIa09 are highly variable. We find significant correlation
between the escaping luminosity and the magnetic flux in both models, which
suggests that the flux eruptions characteristic of the MAD state may be an
important source of QPOs. We found that the fraction of the disk surface area
that is optically thick within $10\,r_{g}$, $a_{\rm{eff}}$, is very weakly
correlated with $L_{\rm{bol}}$ which suggests the cause for the correlation
between $L_{\rm{bol}}$ and $\Phi_{\rm{BH}}$ is not blackbody disk emission.

  6. 6.

KIa09 exhibits coherent quasi-periodic variability with a period of $\sim
2000\,t_{g}$ in the final $\sim 10,000\,t_{g}$ of the simulation. The
luminosity decreases to $\lesssim 25\%$ of its peak during each dip. Such
variability may be useful in interpreting long period variability in BHXBs and
jetted TDEs such as Swift J1644+57.

  7. 7.

For all mass accretion rates $\lambda$ that we have studied, the BH in KIa00
($a_{*}=0$) spins up with a positive spin-up parameter $s\sim 2$. In contrast,
the BH in KIa09 ($a_{*}=0.9$) spins down with $-3\lesssim s\lesssim-1$. The
spin-down rate becomes larger with increasing $\lambda$.

  8. 8.

For accretion disks in the MAD state with mass accretion rates
$\lambda=0.35-1.2$, we estimate that the BHs will, with continued accretion,
reach an equilibrium spin in the range $0.5\lesssim a_{*,{\rm{eq}}}\lesssim
0.6$. This value of $a_{*}$ is only about half the equilibrium spin,
$a_{*,{\rm{eq}}}=0.998$, expected with the standard NT thin disk model.

## Acknowledgements

We are grateful to Koushik Chatterjee, Josh Speagle, and Matthew Liska for
useful comments and suggestions. We also thank the anonymous reviewer for many
useful suggestions which improved the quality of this work. This work was
supported in part by NSF grant AST-1816420, and made use of computational
support from NSF via XSEDE resources (grant TG-AST080026N). This work was
carried out at the Black Hole Initiative at Harvard University, which is
supported by grants from the John Templeton Foundation and the Gordon and
Betty Moore Foundation.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abramowicz & Fragile (2013) Abramowicz M. A., Fragile P. C., 2013, Living Reviews in Relativity, 16, 1
  * Abramowicz et al. (1988) Abramowicz M. A., Czerny B., Lasota J. P., Szuszkiewicz E., 1988, ApJ, 332, 646
  * Avara et al. (2016) Avara M. J., McKinney J. C., Reynolds C. S., 2016, MNRAS, 462, 636
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Begelman et al. (2022) Begelman M. C., Scepi N., Dexter J., 2022, MNRAS, 511, 2040
  * Belloni (2010) Belloni T. M., 2010, in Belloni T., ed., , Vol. 794, Lecture Notes in Physics, Berlin Springer Verlag. p. 53, doi:10.1007/978-3-540-76937-8_3
  * Belloni et al. (2000) Belloni T., Klein-Wolt M., Méndez M., van der Klis M., van Paradijs J., 2000, A&A, 355, 271
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Boccardi et al. (2021) Boccardi B., et al., 2021, A&A, 647, A67
  * Chael et al. (2019) Chael A., Narayan R., Johnson M. D., 2019, MNRAS, 486, 2873
  * Chatterjee et al. (2009) Chatterjee R., et al., 2009, ApJ, 704, 1689
  * Chatterjee et al. (2019) Chatterjee K., Liska M., Tchekhovskoy A., Markoff S. B., 2019, MNRAS, 490, 2200
  * Corbel et al. (2005) Corbel S., Kaaret P., Fender R. P., Tzioumis A. K., Tomsick J. A., Orosz J. A., 2005, ApJ, 632, 504
  * Curd (2021) Curd B., 2021, MNRAS, 507, 3207
  * Curd & Narayan (2019) Curd B., Narayan R., 2019, MNRAS, 483, 565
  * Dexter et al. (2021) Dexter J., Scepi N., Begelman M. C., 2021, ApJ, 919, L20
  * Effron (1979) Effron B., 1979, Ann. Stat., 7, 1
  * Esin et al. (1997) Esin A. A., McClintock J. E., Narayan R., 1997, ApJ, 489, 865
  * Fender et al. (2004) Fender R. P., Belloni T. M., Gallo E., 2004, MNRAS, 355, 1105
  * Gammie et al. (2003) Gammie C. F., McKinney J. C., Tóth G., 2003, ApJ, 589, 444
  * Haiman (2004) Haiman Z., 2004, ApJ, 613, 36
  * Hawley et al. (2011) Hawley J. F., Guan X., Krolik J. H., 2011, ApJ, 738, 84
  * Heyvaerts et al. (1996) Heyvaerts J., Priest E. R., Bardou A., 1996, ApJ, 473, 403
  * Ichimaru (1977) Ichimaru S., 1977, ApJ, 214, 840
  * Igumenshchev (2008) Igumenshchev I. V., 2008, ApJ, 677, 317
  * Igumenshchev et al. (2003) Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ, 592, 1042
  * Kaaz et al. (2022) Kaaz N., Murguia-Berthier A., Chatterjee K., Liska M., Tchekhovskoy A., 2022, arXiv e-prints, p. arXiv:2201.11753
  * Kelley et al. (2014) Kelley L. Z., Tchekhovskoy A., Narayan R., 2014, MNRAS, 445, 3919
  * Kinch et al. (2020) Kinch B. E., Noble S. C., Schnittman J. D., Krolik J. H., 2020, ApJ, 904, 117
  * Kinch et al. (2021) Kinch B. E., Schnittman J. D., Noble S. C., Kallman T. R., Krolik J. H., 2021, ApJ, 922, 270
  * Kovalev et al. (2020) Kovalev Y. Y., Pushkarev A. B., Nokhrina E. E., Plavin A. V., Beskin V. S., Chernoglazov A. V., Lister M. L., Savolainen T., 2020, MNRAS, 495, 3576
  * Lančová et al. (2019) Lančová D., et al., 2019, ApJ, 884, L37
  * Liska et al. (2022) Liska M. T. P., Musoke G., Tchekhovskoy A., Porth O., Beloborodov A. M., 2022, arXiv e-prints, p. arXiv:2201.03526
  * Lovelace et al. (1994) Lovelace R. V. E., Romanova M. M., Newman W. I., 1994, ApJ, 437, 136
  * Lubow et al. (1994) Lubow S. H., Papaloizou J. C. B., Pringle J. E., 1994, MNRAS, 267, 235
  * Lupi et al. (2016) Lupi A., Haardt F., Dotti M., Fiacconi D., Mayer L., Madau P., 2016, MNRAS, 456, 2993
  * Madau et al. (2014) Madau P., Haardt F., Dotti M., 2014, ApJ, 784, L38
  * McKinney et al. (2012) McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, MNRAS, 423, 3083
  * McKinney et al. (2014) McKinney J. C., Tchekhovskoy A., Sadowski A., Narayan R., 2014, MNRAS, 441, 3177
  * McKinney et al. (2015) McKinney J. C., Dai L., Avara M. J., 2015, MNRAS, 454, L6
  * McKinney et al. (2017) McKinney J. C., Chluba J., Wielgus M., Narayan R., Sadowski A., 2017, MNRAS, 467, 2241
  * Mishra et al. (2020) Mishra B., Begelman M. C., Armitage P. J., Simon J. B., 2020, MNRAS, 492, 1855
  * Morales Teixeira et al. (2018) Morales Teixeira D., Avara M. J., McKinney J. C., 2018, MNRAS, 480, 3547
  * Narayan & Yi (1995) Narayan R., Yi I., 1995, ApJ, 452, 710
  * Narayan et al. (2003) Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ, 55, L69
  * Narayan et al. (2012) Narayan R., SÄ dowski A., Penna R. F., Kulkarni A. K., 2012, MNRAS, 426, 3241
  * Narayan et al. (2022) Narayan R., Chael A., Chatterjee K., Ricarte A., Curd B., 2022, MNRAS, 511, 3795
  * Novikov & Thorne (1973) Novikov I. D., Thorne K. S., 1973, in Black Holes (Les Astres Occlus). pp 343–450
  * Pacucci et al. (2015) Pacucci F., Volonteri M., Ferrara A., 2015, MNRAS, 452, 1922
  * Pasham et al. (2015) Pasham D. R., et al., 2015, ApJ, 805, 68
  * Pelupessy et al. (2007) Pelupessy F. I., Di Matteo T., Ciardi B., 2007, ApJ, 665, 107
  * Penna et al. (2010) Penna R. F., McKinney J. C., Narayan R., Tchekhovskoy A., Shafee R., McClintock J. E., 2010, MNRAS, 408, 752
  * Penrose (1969) Penrose R., 1969, Nuovo Cimento Rivista Serie, 1, 252
  * Rees et al. (1982) Rees M. J., Begelman M. C., Blandford R. D., Phinney E. S., 1982, Nature, 295, 17
  * Reis et al. (2012) Reis R. C., Miller J. M., Reynolds M. T., Gültekin K., Maitra D., King A. L., Strohmayer T. E., 2012, Science, 337, 949
  * Saxton et al. (2012) Saxton C. J., Soria R., Wu K., Kuin N. P. M., 2012, MNRAS, 422, 1625
  * Sbarufatti et al. (2012) Sbarufatti B., Burrows D. N., Gehrels N., Kennea J. A., 2012, The Astronomer’s Telegram, 4398, 1
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shapiro (2005) Shapiro S. L., 2005, ApJ, 620, 59
  * Shapiro et al. (1976) Shapiro S. L., Lightman A. P., Eardley D. M., 1976, ApJ, 204, 187
  * Sądowski (2016) Sądowski A., 2016, MNRAS, 459, 4397
  * Sądowski & Narayan (2015a) Sądowski A., Narayan R., 2015a, MNRAS, 453, 3213
  * Sądowski & Narayan (2015b) Sądowski A., Narayan R., 2015b, MNRAS, 454, 2372
  * Sądowski & Narayan (2016) Sądowski A., Narayan R., 2016, MNRAS, 456, 3929
  * Sądowski et al. (2013) Sądowski A., Narayan R., Tchekhovskoy A., Zhu Y., 2013, MNRAS, 429, 3533
  * Sądowski et al. (2014) Sądowski A., Narayan R., McKinney J. C., Tchekhovskoy A., 2014, MNRAS, 439, 503
  * Sądowski et al. (2015) Sądowski A., Narayan R., Tchekhovskoy A., Abarca D., Zhu Y., McKinney J. C., 2015, MNRAS, 447, 49
  * Sądowski et al. (2017) Sądowski A., Wielgus M., Narayan R., Abarca D., McKinney J. C., Chael A., 2017, MNRAS, 466, 705
  * Sutherland & Dopita (1993) Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253
  * Tanaka & Haiman (2009) Tanaka T., Haiman Z., 2009, ApJ, 696, 1798
  * Tchekhovskoy et al. (2010) Tchekhovskoy A., Narayan R., McKinney J. C., 2010, ApJ, 711, 50
  * Tchekhovskoy et al. (2011) Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS, 418, L79
  * Tchekhovskoy et al. (2012) Tchekhovskoy A., McKinney J. C., Narayan R., 2012, in Journal of Physics Conference Series. p. 012040 (arXiv:1202.2864), doi:10.1088/1742-6596/372/1/012040
  * Tchekhovskoy et al. (2014) Tchekhovskoy A., Metzger B. D., Giannios D., Kelley L. Z., 2014, MNRAS, 437, 2744
  * Thorne (1974) Thorne K. S., 1974, ApJ, 191, 507
  * Tóth (2000) Tóth G., 2000, Journal of Computational Physics, 161, 605
  * Volonteri & Rees (2005) Volonteri M., Rees M. J., 2005, ApJ, 633, 624
  * Volonteri & Rees (2006) Volonteri M., Rees M. J., 2006, ApJ, 650, 669
  * Yoo & Miralda-Escudé (2004) Yoo J., Miralda-Escudé J., 2004, ApJ, 614, L25
  * Yuan & Narayan (2014) Yuan F., Narayan R., 2014, ARA&A, 52, 529
  * Zauderer et al. (2013) Zauderer B. A., Berger E., Margutti R., Pooley G. G., Sari R., Soderberg A. M., Brunthaler A., Bietenholz M. F., 2013, ApJ, 767, 152
  * van Ballegooijen (1989) van Ballegooijen A. A., 1989, in Belvedere G., ed., Astrophysics and Space Science Library Vol. 156, Accretion Disks and Magnetic Fields in Astrophysics. p. 99, doi:10.1007/978-94-009-2401-7_10

##  Appendix A Effects of Boundary Condition Choice Behind Injection Region

In this section, we quantify whether or not the choice of boundary conditions
has a significant impact on the solution of the simulation in the regions
where we analyze the output. We perform a rerun of KIa09 in 3D from
$t=30,000-34,000\,t_{g}$ with a reflecting boundary condition on $\mathbb{B}$
in the 2 radially exterior ghost cells of the injection region
($\mathbb{B}_{\rm{GC}}$ henceforth). We then compare diagnostics with the same
time period in KIa09 with $\mathbb{B}_{\rm{GC}}=0$ to evaluate whether or not
key results presented in this work would be significantly impacted by this
choice.

In Figure 16, we show the evolution of $\dot{M}$, $\phi_{\rm{BH}}$, and
$L_{\rm{bol}}$ for both boundary condition choices using $\phi$-averaged data.
As expected, the evolution is similar although there are minor deviations in
the rerun as the system is turbulent. The most significant deviation is in
$\phi_{\rm{BH}}$ and $L_{\rm{bol}}$ around around $33,000\,t_{g}$ due to a
more significant flux eruption in the $\mathbb{B}_{\rm{GC}}=0$ simulation
data. This is not unusual as the system is turbulent and reruns of the same
state will yield slightly different behavior on short time scales.

We also compare the $\sigma_{M}=1$ profile of the jet using data averaged in
$\phi$ and time averaged over $t=30,000-34,000\,t_{g}$ (Figure 17). The
overall jet profile is similar, though we note that the sharp drop in
$\phi_{\rm{BH}}$ around $33,000\,t_{g}$ leads to a slightly less wide jet in a
time-averaged sense due to the jet narrowing during low power states similar
to Figure 9. With that minor deviation aside, the jet behaves in a similar
manner during the rerun. Indeed, we find that a measure of the jet profile for
the top (bottom) jet over the time period $t=30,000-34,000\,t_{g}$ yields
$k\approx 0.45$ (0.47), which is in the range of values shown in Figure 10 for
the sub-Eddington phase.

Lastly, we measured the disk scale height at $r=5\,r_{g}$ and find
$(H/R)_{r=5r_{g}}\approx 0.12$, which is similar to what we find over the same
time period in Table 2.

As expected, the effect of our boundary condition choice for
$\mathbb{B}_{\rm{GC}}$ on the evolution of the disk and jet appears to be
minimal. As long as a large flux of poloidal magnetic field is maintained (and
a MAD state is maintained), the behaviour is generally the same. As noted in
the text, this choice of boundary conditions $\mathbb{B}_{\rm{GC}}=0$ was made
to avoid the injection of magnetic field of the opposite sign through the
corner ghost cells.

Figure 16: Here we compare $\dot{M}$, $\phi_{\rm{BH}}$ and $L_{\rm{bol}}$ in
KIa09 for boundary condition choices $\mathbb{B}_{\rm{GC}}=0$ (solid line) and
$\mathbb{B}_{\rm{GC}}$ reflected (dashed line). As described in the text, we
find that the simulations have similar behaviour, but there are some
deviations due to the turbulent nature of the system. Figure 17: Here we
compare the $\sigma_{M}=1$ contour for the top jet in KIa09 for boundary
condition choices $\mathbb{B}_{\rm{GC}}=0$ (solid line) and
$\mathbb{B}_{\rm{GC}}$ reflected (dashed line). Note that we have time-
averaged over $t=30,000-34,000\,t_{g}$ and $\phi$-averaged over $2\pi$. We
find that the jet shape is roughly the same, although the
$\mathbb{B}_{\rm{GC}}=0$ jet is slightly less wide in a time averaged sense
due to the more extreme flux eruption event around $t=33,000\,t_{g}$.

##  Appendix B Comparison of $\sigma_{M}$ in Models

We compare the behaviour of $\sigma_{M}$ as a function of polar angle $\theta$
at $r=5,15,30\,r_{g}$ between models KIa00 and KIa09 in Figure 18. Model KIa00
does not have a relativistic jet and is not significantly impacted by the
addition of mass in the funnel since $\sigma_{M}<100$ up to $r=5\,r_{g}$.
Model KIa09 on the other hand is subject to density floors at radii as large
as $r<15\,r_{g}$. The gas density $\rho$ and magnetization $\sigma_{M}$
profiles are flat in $\theta$ where the ceiling is active (i.e. see the
profile at $r=5\,r_{g}$). The radiation properties ($L_{\rm{bol}}$ and
$\eta_{\rm{rad}}$) in these regions becomes dependent on the choice of ceiling
model, as we have illustrated in Section 3.1.2.

Figure 18: We show angular profiles of the gas density ($\rho$, top) and
magnetization ($\sigma_{M}$, bottom) for model KIa00 (left) and KIa09 (right).
We show the maximum $\sigma_{M}=100$ as a dashed horizontal line. In KIa09,
the implementation of the numerical ceiling $\sigma_{M}\leq 100$, especially
at $r=5\,r_{g}$, leads to artificially high density gas in the funnel and a
flat density profile in $\theta$. KIa00 generally does not artificially add
mass in the jet since $\sigma_{M}$ is below 100 throughout most of the funnel
even at radii as small as $r=5\,r_{g}$.