# Magnetic field evolution in high and low $\beta$ disks

Payton E. Rodman Institute of Astronomy, University of Cambridge, Madingley
Rd, Cambridge CB3 0HA, UK  Christopher S. Reynolds Institute of Astronomy,
University of Cambridge, Madingley Rd, Cambridge CB3 0HA, UK  Dept. of
Astronomy & Joint Space Science Institute (JSI), University of Maryland,
College Park, MD 20742

(Received September 14, 2023)

###### Abstract

We present results from a pair of high resolution, long timescale ($\sim
10^{5}~{}GM/c^{2}$), global, three dimensional magnetohydrodynamical accretion
disk simulations with differing initial magnetic plasma $\beta$ in order to
study the effects of initial toroidal field strength on production of large-
scale poloidal field. We initialize our disks in approximate equilibrium with
purely toroidal magnetic fields of strength $\beta_{0}=5$ and $\beta_{0}=200$.
We also perform a limited resolution study. We find that simulations of
differing field strength diverge early in their evolution and remain distinct
over the time studied, indicating that initial magnetic conditions leave a
persistent imprint in our simulations. Neither simulation enters the
Magnetically Arrested Disk (MAD) regime. Both simulations are able to produce
poloidal fields from initially-toroidal fields, with the $\beta_{0}=5$
simulation evolving clear signs of a large-scale poloidal field. We make a
cautionary note that computational artifacts in the form of large-scale
vortices may be introduced in the combination of initially-weak field and
disk-internal mesh refinement boundaries, as evidenced by the production of an
$m=1$ mode overdensity in the weak field simulation. Our results demonstrate
that the initial toroidal field strength plays a vital role in simulated disk
evolution for the models studied.

accretion — black hole physics — magnetic fields — magnetohydrodynamics (MHD)

††software: Athena++ (Stone et al., 2020), numpy (Harris et al., 2020), pandas
(Wes McKinney, 2010) scipy (Virtanen et al., 2020), VisIt (Childs et al.,
2012)

##  1 Introduction

Accretion disks are highly turbulent and complex astrophysical systems formed
by infalling gas circularizing around a massive central body such as a young
star, a white dwarf, a neutron star, or a black hole. In the case of accretion
onto a supermassive black hole, the disk physics governs both the growth of
the black hole and the form of the resulting energy output and thus has wide-
reaching implications not just for the immediate environment but on galaxy
scales (e.g. the $M$-$\sigma$ relation; King, 2003). The strong gravitational
forces, relativistic velocities, and magnetic fields present make these
accreting black hole systems an invaluable laboratory for studying a wide
range of extreme physical processes, including the properties of strong
gravitational fields, particle acceleration, and the interaction of magnetic
fields and fluid turbulence, the last of which forms the primary focus of this
paper.

For an accretion disk to actually accrete it must find a way to remove angular
momentum from orbiting fluid elements. This could in theory occur through
local exchanges between gas parcels within the disk and/or through a wholesale
loss of angular momentum from the system via, for example, magnetic braking
from a disk wind. Given the extremely high Reynolds numbers of most
astrophysical disks, intrinsic microscopic viscosity can be immediately ruled
out as an angular momentum transport mechanism as it would drive accretion
over timescales that are far too long to explain the observed luminosity of
disks (see Balbus, 2003, for derivation). The seminal work of Shakura &
Sunyaev (1973) acknowledged that a form of anomalous viscosity could exist as
a result of turbulence, potentially related to the presence of magnetic
fields, but the mechanism for generating this turbulence was then unable to be
specified. Simply invoking a very high Reynolds number as a driver of
turbulence is not sufficient as Keplerian disks are linearly and non-linearly
stable to shearing instabilities even at very high Reynolds numbers (Hawley et
al., 1995, 1996; Balbus et al., 1996). Later work by Balbus & Hawley (1991,
1998) identified a powerful magnetohydrodynamic (MHD) instability which can
drive the necessary turbulence. This magneto-rotational instability (MRI) has
since become the paradigm for accretion disk theory.

Magnetohydrodynamic (MHD) and (more recently) General Relativistic
Magnetohydrodynamic (GRMHD) simulations are now common tools for studying the
internal dynamics of the disk due to the highly non-linear nature of the MRI-
driven turbulence (see e.g. Balbus, 2003; Davis & Tchekhovskoy, 2020, for
review). In spite of much attention in recent years, however, many questions
still remain about the evolution of magnetic fields within accretion disks and
the nature of the turbulence that is formed. Equally importantly,
uncertainties remain about the influence of computational realities (finite
spatial resolution, memory of initial conditions in a limited temporal span,
and non-helicity preserving MHD algorithms) on the outcome of accretion disk
simulations.

Choices in the initial vertical magnetic field were shown to have important
implications for turbulence and overall disk evolution even in early
simulations (Hawley et al., 1995), with turbulent stresses increasing with
increasing net vertical field. Later studies of both shearing box (e.g.
Salvesen et al., 2016; Bai & Stone, 2013) and global (e.g. Mishra et al.,
2020) simulations showed that the presence of an initially vertical field
affects the structure of the simulated accretion flow. If the vertical field
in the inner regions of the disk is sufficiently strong, it may cause the disk
to enter a Magnetically Arrested Disk (MAD; Narayan et al., 2003) state in
which mass accretion is suppressed or even temporarily halted due to magnetic
pressure and, if the central black hole is spinning, jets are launched via the
Blandford-Znajek mechanism (Blandford & Znajek, 1977; Tchekhovskoy et al.,
2011). MADs have been invoked to explain the low-luminosity and powerful jet
in M87 and GRMHD models of MADs can quantitatively match the horizon-scale mm-
band image obtained by the Event Horizon Telescope (Event Horizon Telescope
Collaboration et al., 2019). The counterpart to MAD is the Standard And Normal
Evolution (SANE) state, which is effectively characterized by the absence of a
MAD state. Recent work by Begelman et al. (2022), however, suggests that while
the presence of a net vertical field is important, it plays a secondary role
to that of the toroidal field, with the toroidal component being critically
important in the saturation of magnetic flux. To add further complexity to
this discussion, the net vertical magnetic flux may be zero on a global scale
but non-zero locally (Sorathia et al., 2010). A zero-net-vertical-flux (ZNVF)
global disk may well be composed of many net-vertical-flux (NVF) local
“shearing box” patches, magnetically connected at high-altitudes through the
low density corona of the disk.

In addition to questions on the role of vertical fields and fluxes, the final
state of a simulation appears to depend significantly on the particular
initial field morphology used, even where the initial net vertical field is
kept constant. Early work by Machida et al. (2001) suggested that the initial
magnetic field configuration could affect the mass outflow rate in
Radiatively-Inefficient Accretion Flows (RIAFs), disks in which the radiative
cooling is negligible and local cooling driven by advection — such disks are
generally geometrically-thick and have low-accretion rate. Such flows lack a
large separation of scales between their inflow/viscous times, thermal times,
and dynamical times and so seem to retain memory of their initial conditions
during the inflow (White et al., 2020). Further, since they lack radiative
cooling, the sign of the Bernoulli parameter (and hence the degree to which
regions of the disk are unbound) is sensitive to both the feeding prescription
and the detailed internal energetics associated with the initial field. Later
simulations by Narayan et al. (2012) found that a MAD state could be produced
by varying the number and direction of poloidal loops in the initial magnetic
field, with a single-loop configuration being most conducive to building up
vertical field necessary for the MAD state. White et al. (2020) then built
upon this work, showing that within SANE models variations in initial magnetic
field configuration have strong effects on the energetics and structure of
RIAFs, with smaller initial poloidal loops generating a less coherent radial
field (leading to less horizon-penetrating magnetic flux) and producing polar
outflows while larger poloidal loops showed the opposite. Similar work by
Jiang et al. (2019) on thin, radiation-dominated disks also found factors of
three differences in accretion rate and vertical disk structure depending on
whether the initial poloidal field was dipolar or quadrupolar.

These combined results clearly demonstrate that the properties of simulated
accretion disks are dependent on the artificial and somewhat arbitrary initial
magnetic conditions used in their construction, and such a dependence raises
serious questions as to the validity of simulations as emulations of real
astrophysical disks.

Magnetic fields with net vertical flux or with significant large-scale
vertical/polar components play an important role not only for the emergence of
MAD states, but also for the production of astrophysical jets through the
Blandford-Znajek (Blandford & Znajek, 1977) and Blandford-Payne (Blandford &
Payne, 1982) processes. An open question remains whether, and under what
conditions, a disk can eventually forget the details of the initial field
configuration and grow a self-consistent large-scale field via dynamo action
given sufficient time. Fragile & Sądowski (2017) found that initially strongly
magnetized disks (where the ratio of thermal-to-magnetic pressure is only
$\beta=P_{\rm gas}/P_{\rm mag}\sim 10$) with purely toroidal fields were not
locally self-sustaining, albeit the simulations were run over a short
timescale ($700~{}GM/c^{3}$) that may not have allowed enough time for
poloidal growth. Global simulations by Liska et al. (2020) run for much longer
times ($1.3\times 10^{5}~{}GM/c^{3}$) successfully generated a strong enough
poloidal field in situ to enter the MAD regime, but only from an initially
strong field ($\beta=5$), which is likely not reflective of initial conditions
in nature.

Given the long times needed for a disk with an initially weak toroidal
magnetic field to reach equilibrium out to reasonable radii, most simulations
with initially toroidal fields are initialized with a low magnetic plasma
$\beta$ (i.e. a ‘strong’ magnetic field), meaning little is known of the
amplification process for $\beta>100$ and few simulations cover the
intersection of weak and toroidal initial field configurations.

In this work, we study the ability of simulated disks initialized with a
purely toroidal field to grow a large-scale poloidal field. We agree with and
extend the previous findings of Liska et al. (2020), finding large-scale field
growth for disks initialized with either a strong ($\beta=5$) or a weak
($\beta=200$) magnetic field. The core of this work is based on two long-
duration simulations with initially-toroidal fields, a strong-field
($\beta=5$) and a weak-field ($\beta=200$) case. Our longest simulation (the
weak-field case) is run out to $t\sim 4\times 10^{5}~{}GM/c^{3}$ in order to
capture the long-term evolution of the magnetic field. Given our focus on
understanding the influence of numerics on results, we also run two additional
simulations that have identical set-up except for higher resolution. We see
evidence of poloidal flux generation in all cases, with the weaker initial
field case saturating with weaker and more disorganized poloidal field. All of
our model disks display dynamo-cycles in the toroidal field even though the
disks are geometrically-thick. We also identify the subtle ability of internal
mesh-refinement boundaries to drive vorticity-induced density asymmetries.

The basic equations of the problem and computational setup are given in
Section 2. The analysis of our four simulations, with special emphasis on the
growth of large-scale poloidal field, are presented in Section 3. We discuss
and summarize our findings with our conclusions in Section 4.

##  2 Basic Model and Computational Setup

In this paper, we model the evolution of two global accretion disks with
initially toroidal magnetic fields, one weak-field case and one strong-field
case (see Section 2.2). We additionally produce two more simulations with the
same initial magnetic field strengths but with double resolution in the
$\theta$\- and $\phi$-directions. As the goal of this work is to study the
evolution of the magnetic field, we keep all other aspects fixed between
simulations. To reduce the computational resources needed, we neglect any
physics which does not contribute directly to the evolution of the magnetic
field; specifically this work does not include radiative physics (formally
making our disks RIAFs) nor full general relativity (instead employing a
pseudo-Newtonian potential). Instead, we focus all resources into evolving the
equations of non-relativistic ideal-MHD at high resolution for as long as
possible.

All lengths, times, and velocities are given in natural units scaled to the
mass of the black hole, formally setting $G=M=c=1$.

###  2.1 Governing Equations and Simulation Code

Our accretion disk models are evolved using the MHD code Athena++ v21.0 (Stone
et al., 2020) on a spherical polar grid $(r,\theta,\phi)$, with the piece-wise
linear method (PLM) for spatial reconstruction, second-order (VL2) time-
stepping, a HLLD Riemann solver, and a Courant, Friedrichs, & Lewy (CFL)
number of $0.3$. The code is parallelized over 1008 Cascade Lake cores on the
Cambridge Service for Data Driven Discovery (CSD3) system Peta4, achieving a
full-node performance of $\sim$ 21-22 MZone-cycles/core-second.

Athena++ solves the equations of ideal MHD, given in conservative form as
(Stone et al., 2008)

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot\left(\rho\boldsymbol{v}\right)$ $\displaystyle=0,$ (1)
$\displaystyle\frac{\partial\rho\boldsymbol{v}}{\partial
t}+\nabla\cdot\left(\rho\boldsymbol{vv}-\boldsymbol{BB}+{\boldsymbol{P}^{*}}\right)$
$\displaystyle=-\rho\nabla\Phi_{\rm PW},$ (2) $\displaystyle\frac{\partial
E}{\partial
t}+\nabla\cdot\left[\left(E+P^{*}\right){\boldsymbol{v}}-{\boldsymbol{B}}\left(\boldsymbol{B\cdot
v}\right)\right]$ $\displaystyle=0,$ (3)
$\displaystyle\frac{\partial{\boldsymbol{B}}}{\partial
t}-\nabla\times\left({\boldsymbol{v}}\times{\boldsymbol{B}}\right)$
$\displaystyle=0,$ (4)

where $\boldsymbol{P^{*}}=(P+B^{2}/2)\boldsymbol{I}$, $\boldsymbol{I}$ is the
diagonal unit tensor, $B^{2}=\boldsymbol{B}\cdot\boldsymbol{B}$, and all other
variables having their usual meaning. The total energy density is given by

$E=\frac{P}{\gamma-1}+\frac{1}{2}\rho v^{2}+\frac{B^{2}}{2}+\rho\Phi_{\rm
PW},$ (5)

where we have imposed a $\gamma$-law equation of state relating the internal
energy density $\mathcal{E}$ to the pressure, $P=(\gamma-1)\mathcal{E}$. We
use the pseudo-Newtonian Paczyński–Wiita gravitational potential, $\Phi_{\rm
PW}$,

$\Phi_{\rm PW}(r)=-\frac{GM}{r-2r_{g}},\qquad r_{g}\equiv GM/c^{2},$ (6)

which features an innermost stable circular orbit (ISCO) at $r_{\rm
ISCO}=6r_{g}$ and is an adequate approximation (within Newtonian physics) of
the gravitational field around a non-rotating black hole. The solenoidal
condition on the magnetic field,

$\nabla\cdot\boldsymbol{B}=0,$ (7)

is enforced by upwind constrained transport (CT; Evans & Hawley, 1988).

The system of equations 1–7 is closed and may be solved numerically. During
the numerical solution, we enforce density and pressure floors at $\rho_{\rm
floor}=3\times 10^{-5}$ and $P_{\rm floor}=1\times 10^{-10}$. Due to the high
Alfvén speeds that arise in the low-density regions of our b5_hi run, we
increased the density floor in this simulation to $\rho_{\rm floor}=3\times
10^{-4}$ at $t\sim 2\times 10^{4}~{}GM/c^{3}$, representing a $0.16$% increase
in the total mass at that moment which does not appreciably affect our
results.

###  2.2 Simulation Setup

Our computational domain spans $r\in\left[5r_{g},290r_{g}\right]$,
$\theta\in[0,\pi]$, and $\phi\in[0,2\pi)$. Our fiducial pair of simulations
have a base resolution of $(N_{r},N_{\theta},N_{\phi})=(112\times 32\times
32)$ and 3 additional levels of static mesh refinement (SMR), with each level
double the resolution of the one below (Figure 1). This gives us a maximum
effective resolution of $(N_{r},N_{\theta},N_{\phi})=(896\times 256\times
256)$ for $r<100r_{g}$ and height $\lesssim 2H$, resulting in 20 cells per
scale height $H$. Using SMR in this way, our resulting mesh (Figure 1) is de-
resolved near the poles, increasing cell sizes in this region and avoiding
untenably small timesteps that would arise from the CFL condition. As a result
of this mesh construction we are able to simulate the full spherical domain
including the poles. Our two high-resolution runs use a base resolution of
$(N_{r},N_{\theta},N_{\phi})=(112\times 64\times 64)$, i.e. double the
resolution in $\theta$ and $\phi$, but are otherwise identical. This gives a
maximum effective resolution of $(N_{r},N_{\theta},N_{\phi})=(896\times
512\times 512)$ and 40 cells per scale height $H$.

Figure 1: Visualization of the fiducial-resolution simulation mesh in 3D, with
a cutout to show internal zones.

The base grid is logarithmically-spaced in $r$ and linearly spaced in $\theta$
and $\phi$. We use “polar” boundary conditions (BCs) at the poles (see Stone
et al., 2020, for details), periodic BCs in the $\phi$-direction, and outflow
BCs on the inner and outer $r$ boundaries. Our polar boundary condition allows
magnetic field lines to correctly pass through the coordinate singularity, and
allows a polar flux bundle to stay in place without (unphysically) being
destroyed by the boundary.

###  2.3 Initial Conditions

To contrast with previous work and further explore robustness to initial
conditions, we initialize our disks with an approximate disk-like
configuration following that of Hogg & Reynolds (2018a), rather than the
commonly-used Fishbone-Moncrief torus (Fishbone & Moncrief, 1976).

We initialize our disks with a power-law density profile

$\rho(R)=\rho_{0}R^{-3/2}e^{-(z^{2})/(2c_{s}^{2}R^{3})}\,F_{\rm s}(R),$ (8)

where $R$ is the cylindrical radius, $F_{s}(R)$ is a function motivated by the
zero-torque inner boundary condition of Shakura & Sunyaev (1973),

$F_{s}(R)=1-\left(\frac{r_{\rm ISCO}}{R}\right)^{1/2}.$ (9)

$\rho_{0}=100$ is a scaling factor set to give ${\rm max}\\{\rho(r)\\}\sim 1$,
and $c_{s}$ is the isothermal sound speed in a thin disk. The resulting
density scale-height of the initial disk $H_{i}$ is given by,

$c_{s}=v_{\rm K}H_{i}$ (10)

for Keplerian rotational velocity $v_{K}$.

To this density profile, we add a non-axisymmetric perturbation of the form
$\sin(10\phi)$ at the 1% level. Breaking the formal axisymmetry of the initial
conditions lessens the impact of (transitory) channel modes and hastens the
development of saturated turbulence.

The initial velocity profile is that of a test-particle on a circular orbit
(i.e. “Keplerian”) in the Pseudo-Newtonian potential,

$v_{\phi}=\frac{\sqrt{R}}{R-2},$ (11)

with the other components of velocity set initially to zero
($v_{r}=v_{\theta}=0$). This initializes the disk close to, but not exactly
matching, its natural equilibrium state, reducing the time needed for the disk
to settle into a quasi-steady state and the simulation data to become useful.

The initial magnetic field is purely toroidal field, set by

$B_{\phi}=\sqrt{\frac{2P}{\beta_{0}}},$ (12)

with all other components equal to zero, i.e. the magnetic pressure is set to
be a constant fraction $1/\beta_{0}$ of the gas pressure $P$. We calculate $P$
from the density before the introduction of the non-axisymmetric
perturbations. Thus the corresponding magnetic field is axisymmetric and,
being toroidal, is hence divergence free everywhere. We explore two cases of
field strength: a stronger field scenario ($\beta_{0}=5$) and a weaker field
scenario ($\beta_{0}=200$). Combined with the two resolutions that we explore
(2.2), this defines the four simulations presented in this paper and
summarized in Table 1.

Table 1: Simulation parameters

Simulation ID | Max. Resolution | Duration | ISCO orbits | $\beta_{0}$ | $\langle\langle\beta\rangle\rangle$ | $\langle\langle Q_{\theta}\rangle\rangle$ | $\langle\langle Q_{\phi}\rangle\rangle$ | $\langle\langle\theta_{B}\rangle\rangle$
---|---|---|---|---|---|---|---|---
| ($r\times\theta\times\phi$) | ($GM/c^{3}$) |  |  |  |  |  |
b5 | $896\times 256\times 256$ | $1.67\times 10^{5}$ | 2704 | 5 | $9.39$ | $20.56$ | $63.22$ | $11.85$
b200 | $896\times 256\times 256$ | $4.02\times 10^{5}$ | 6535 | 200 | $23.22$ | $9.22$ | $27.38$ | $10.62$
b5_hi | $896\times 512\times 512$ | $3.15\times 10^{4}$ | 511 | 5 | $7.87$ | $39.89$ | $106.56$ | $11.14$
b200_hi | $896\times 512\times 512$ | $1.54\times 10^{5}$ | 2502 | 200 | $20.90$ | $10.75$ | $47.82$ | $10.75$

Note. — Averages denoted by $\langle\langle X\rangle\rangle$ are taken over
$r\in[5,100]$, $\phi$, and three scale heights in $\theta$, at late times in
the disk as defined in Section 2.4.

###  2.4 Diagnostics and Convergence

MHD accretion disks are highly turbulent systems and thus, even when they have
reached a “quasi steady state”, will be characterized by strong spatial and
temporal fluctuations in all quantities. Useful diagnostics of disk structure
necessarily require some averaging process. In this paper, spatial averages
are denoted by $\langle X\rangle$; by default, this means averaging over
$r<100r_{g}$, all $\phi$, and across $\theta={\rm constant}$ shells within
$\pm 3$ scale heights of the mid-plane. Quantities denoted by $\langle\langle
X\rangle\rangle$ are additionally averaged over time; by default, we average
our weak field simulations (b200 and b200_hi) for $t\gtrsim 2000T_{\rm ISCO}$
($t>1.2\times 10^{5}~{}GM/c^{3}$) and our strong-field simulations (b5 and
b5_hi) for $t\gtrsim 500T_{\rm ISCO}$ ($t>3.0\times 10^{4}~{}GM/c^{3}$). Any
deviations from these defaults are indicated by $\langle X\rangle^{\star}$ or
$\langle\langle X\rangle\rangle^{\star}$ and details specified within the
text.

In addition to looking at the evolution of quantities already described
($\beta,\rho,v_{r}$), we calculate and track several additional standard
quantities. The mass accretion rate across some closed surface $S$ is

$\dot{M}=\iint_{S}\rho\boldsymbol{v}\cdot d\boldsymbol{S},$ (13)

which we calculate across the sphere $r=6$, namely the innermost circular
stable orbit (ISCO) of this potential.

The geometrical scale height $h(r)$ is calculated as in Hogg & Reynolds
(2018b), with the scale height at each radius given by

$\frac{h(r)}{r}=\Biggl{\langle}\frac{\int[\theta(r)-\bar{\theta}(r)]^{2}\rho
d\Omega}{\int\rho d\Omega}\Biggr{\rangle}$ (14)

for solid angle $d\Omega=\sin\theta\,{\rm d}\theta\,{\rm d}\phi$ and

$\bar{\theta}(r)=\frac{\int\theta(r)\rho d\Omega}{\int\rho d\Omega}.$ (15)

The radial bins from equation 14 are then averaged by weighting them according
to the radial width at each shell (as cells are logarithmically spaced in
$r$), giving a disk-averaged scale height $H$. To prevent confusion with other
parameters within the text, we use $H$ to refer to the disk-averaged
geometrical scale height rather than $\langle\frac{h(r)}{r}\rangle$.

With the focus of this paper on the magnetic evolution of the disk, a
parameter of key interest is the horizon penetrating magnetic flux,

$\Phi_{\rm BH}=\frac{\sqrt{4\pi}}{2}\iint_{S_{H}}\boldsymbol{B}\cdot
d\boldsymbol{S},$ (16)

where $S_{H}$ is the upper half-sphere near the inner boundary of our domain
$(r=6,\theta>0)$. This is then re-normalized by the late-time-averaged
accretion rate to $\phi=\Phi_{\rm BH}/(\langle\langle\dot{M}\rangle\rangle
r_{g}c^{2})^{-1/2}$ (Tchekhovskoy et al., 2011). For a disk scale height
$H\sim 0.3$ and BH spin $a=0$, the transition from SANE to MAD states is
expected to occur for $\phi\sim 50$ (Tchekhovskoy et al., 2012)111We note that
$\phi\sim 50$ is the saturation value in code or Gaussian units. When using
Lorentz-Heaviside units the factor of $\sqrt{4\pi}$ in equation 16 disappears,
giving $\phi\sim 15$ instead..

We also calculate the Shakura-Sunyaev alpha (Shakura & Sunyaev, 1973),

$\alpha_{\rm SS}=\frac{T_{\hat{r}\hat{\phi}}}{P}=\frac{T_{{\rm
Max},\hat{r}\hat{\phi}}+T_{{\rm Rey},\hat{r}\hat{\phi}}}{P}$ (17)

where $T_{\hat{r}\hat{\phi}}$ is the fluid frame stress, and $T_{{\rm
Max},\hat{r}\hat{\phi}}=-B_{r}B_{\phi}$ and $T_{{\rm
Rey},\hat{r}\hat{\phi}}=\rho u_{r}\delta u_{\phi}$ are the Maxwell and
Reynolds components of that stress, where $\delta u_{\phi}\equiv
u_{\phi}-r\Omega_{\rm Kep}$ is the fluctuating component of the
$\phi$-velocity. For thin disks, $\alpha_{\rm SS}\sim 0.1$, however this value
is known to vary with disk thickness, initial magnetic field
morphology/strength, and radius within the disk (e.g. Penna et al., 2013).

A key question is the degree to which a numerical model of an accretion disk
is converged. There are two distinct senses in which we must consider
convergence. Firstly, since our disks are initialized in a laminar state that
is close to hydrodynamic equilibrium but strongly unstable to the MRI, they
must first undergo evolution before settling into a (turbulent) statistical
steady state. Only then can the simulation data be considered physically
meaningful (out to the radius where we achieve inflow equilibrium). Diagnosing
this state informs the times on which we focus our analysis. Secondly, as we
are interested in studying the behavior of the magnetic field we need to be
sure that the underlying physical mechanisms that govern its evolution are
spatially well-resolved within the disk. To assess these criteria, we turn to
two measures of convergence: the quality factors $Q_{\theta}$ and $Q_{\phi}$,
and the magnetic tilt angle $\theta_{B}$. Additionally, we are able to do a
very limited resolution study by comparing our low and high resolution
simulations.

The quality factors, $Q$, defined by Noble et al. (2010) as

$\displaystyle Q_{\theta}$ $\displaystyle=\frac{\lambda_{{\rm
MRI},\theta}}{r\Delta\theta}$ (18) $\displaystyle Q_{\phi}$
$\displaystyle=\frac{\lambda_{{\rm MRI},\phi}}{R\Delta\phi},$ (19)

are a direct measure of the resolvability of the linear MRI modes via a
comparison of the characteristic wavelength of the MRI modes $\lambda_{\rm
MRI}=2\pi v_{A}/\Omega$, where $v_{A}$ is the Alfvén velocity, with the size
of the basic simulation voxel. The $Q$’s are simply a measure of the number of
resolution elements that span one wavelength of the fastest growing MRI mode.

Values as low as $Q_{\theta}=6-8$ are sufficient to capture the linear MRI
(Flock et al., 2010) in low-spatial order methods, but must rise to
$Q_{\theta}>10$ with $Q_{\phi}>20$ to capture the non-linear interactions
between MRI modes (Hawley et al., 2011, 2013; Sorathia et al., 2012), although
a higher $Q_{\phi}$ or $Q_{z}$ (when in cylindrical coordinates) may offset a
slightly lower $Q_{\theta}$.

Our second diagnostic tool is the magnetic tilt angle, $\theta_{B}$. First
used by Guan et al. (2009), $\theta_{B}$ is a measure of the correlations
between the radial and toroidal fields,

$\theta_{B}=-\arctan\left<\frac{B_{r}}{B_{\phi}}\right>.$ (20)

Guan et al. (2009) originally found this to have a characteristic value of
$\theta_{B}\sim 15^{\circ}$ for a well-resolved saturated MRI-driven
turbulence with later works on local (Hawley et al., 2011) and global
(Sorathia et al., 2012; Hawley et al., 2013; Hogg & Reynolds, 2016, 2018a)
simulations revising this number down to 11-13∘. The physics dictating the
saturation of MRI-driven turbulence and hence the value of the tilt angle is
still not clear, although Pessah (2010) describe such behavior with a model in
which saturation corresponds to disruption of linear-MRI modes by the
parasitic instabilities of Goodman & Xu (1994). Our primary focus with
$\theta_{B}$ is not the final value itself, as this is not a direct measure of
resolution and so may be affected by other factors in the disk, but we instead
look for a stabilization in the value of this parameter.

Figure 2: Convergence metrics: (top and middle) quality factors in the
$\theta$\- and $\phi$-directions, and (bottom) magnetic tilt angle,
$\theta_{B}$. The dotted line indicates the minimum value required to resolve
the linear MRI, while the dashed lines represent the requirements for
resolving the non-linear MRI.  Figure 3: Three-dimensional side-on volume
rendering of the simulation domain for one of our fiducial runs, b200, at a
late time in the disk. Color denotes density.

The time-evolution of these convergence metrics for all of our runs are shown
in Figure 2, and the late-time average is reported in Table 1. Through these
metrics we find that our strong-field simulation b5 is well-resolved with
$Q_{\theta}\sim 16$, $Q_{\phi}\sim 38$, and a tilt angle that stabilizes at
$\theta_{B}\approx 12^{\circ}$. The numbers are less comfortable for b200
(Figure 2). The lower values of $Q_{\theta}\sim 9$ indicates that b200 only
marginally resolves the dominant linear MRI modes, although $\theta_{B}$ does
eventually stabilize for $t\gtrsim 2\times 10^{5}~{}GM/c^{3}$, albeit at a
notably lower value ($\theta_{B}\sim 10)$ than the stronger field case.

A more direct way to assess the effects of resolution is to compare our
standard resolution runs (b200 and b5) with their high-resolution counterparts
(b200_hi and b5_hi). As expected, these high-resolution runs achieve much
higher $Q_{\theta}$ and $Q_{\phi}$ values, well above the fiducial levels
needed. The magnetic tilt angle $\theta_{B}$ take longer to stabilize than in
their lower-resolution counterparts but eventually stabilize at comparable
levels to the low-resolution runs.

We conclude that all simulations have achieved a converged state by the end of
their respective runs, although b200 may be only marginally resolving both the
linear and non-linear evolution of the MRI.

##  3 Simulation Results

###  3.1 Disk structure and evolution

Starting from the initial state described in Section 2.3, the disks undergo
rapid evolution driven principally by the growth and eventual saturation of
MRI-driven turbulence. The outer parts of the simulated disks never achieve a
state of inflow equilibrium; we focus our analysis on the inner (and most
highly resolved) regions of the disks within $r=100r_{g}$. Figure 3 shows a
late-time three-dimensional volume rendering of run b200, and Figure 4 shows
the evolution of mass accretion rate across the ISCO as well as the disk-
averaged plasma-$\beta$, geometric scale height $H$, and Shakura-Sunyaev
$\alpha_{\rm SS}$. Both b5 and b200 show a characteristic early peak and slow
drop-off in mass accretion rate $\dot{M}$, as expected for disks with finite
mass (Figure 4). b200 stabilizes to a higher mass accretion rate of
$\dot{M}\sim 0.25$, possibly related to the presence of a density lump that
inhibits early-time accretion. In all models, we find that
$\langle\beta\rangle$ reaches an approximately steady value relatively
quickly. The memory of the strength of the initial magnetic field is retained,
with run b200 stabilizing at $\langle\beta\rangle\sim 20$ and run b5
stabilizing at $\langle\beta\rangle\sim 10$. The initial magnetization also
makes a difference to the disk thickness, with stabilization at $\langle
H\rangle\sim 0.27$ and $\langle H\rangle\sim 0.32$ for b200 and b5
respectively. The dominant pressure dictating the disk thickness in all cases
is gas pressure, with a modest amount of additional magnetic support
thickening the disk in the stronger initial field case b5. The spatially-
averaged $\langle\alpha_{\rm SS}\rangle$ (lower panel) remains relatively
stable around $10^{-2}$-$10^{-1}$, as is expected for accretion disks, with a
small jump in b5_hi at late times which coincides with a thinning of the disk
(i.e. lower $H$).

Figure 4: (Top) Mass flux $\dot{M}$ through the surface at $r=6r_{g}$, (top
middle) spatially-averaged magnetic plasma $\langle\beta\rangle$, (bottom
middle) spatially-averaged geometric scale height $\langle H\rangle$, and
(bottom) spatially-averaged Shakura-Sunyaev alpha $\langle\alpha_{\rm
SS}\rangle$.

For the stronger field case, doubling the $\theta$\- and $\phi$-resolution
(going from b5 to b5_hi) makes very little difference to the overall evolution
of the disk, with $\dot{M}$, $\langle\beta\rangle$ and $\langle H\rangle$
following very similar time-traces. In the case of the weaker field model, the
early-to-intermediate timescale evolution is very similar for the high-
resolution models whereas the late-time magnetization is a little higher
($\langle\beta\rangle$ marginally smaller) leading to a slightly thicker disk.
This further suggests that the fiducial resolution weaker field run b200 only
marginally resolves the non-linear evolution of the MRI.

Figure 5: Late-time-averaged profiles for the density (top), radial velocity
(top middle), temperature (bottom middle), and Shakura-Sunyaev alpha (bottom)
as a function of disk radius. Dashed lines for b200 and b5 are averaged over
the same time windows as b200_hi and b5_hi respectively. Note: Spatial
averages are taken azimuthally and within $3H$.

The radial profiles of density $\langle\langle\rho\rangle\rangle^{*}$, radial
velocity $\langle\langle v_{r}\rangle\rangle^{*}$, temperature $\langle\langle
P/\rho c^{2}\rangle\rangle^{*}$, and effective Shakura-Sunyaev
$\langle\langle\alpha_{\rm SS}\rangle\rangle^{*}$ are shown for all four runs
in Figure 5. These are time-averaged from the point where the run is deemed to
have stabilized (as laid out in Section 2.4) out to the end of the simulation.
Given that the different runs have quite different end-times, some caution
must be exercised in making direct comparisons between profiles from different
runs. To aid in this comparison, we also provide the profiles for the low-
resolution runs (b200 and b5) averaged over the same time range as the high-
resolution counterparts shown as dashed lines. For example, the large
difference in the overall normalization of density between b5 and b5_hi is
almost all due to the fact that the lower-resolution run was evolved much
longer and hence drained much more of its mass; when the average is calculated
over the same time range, this difference largely disappears. However, the
comparison of b200_hi and b5 is particularly interesting. Both runs well-
resolve the MRI and have evolved for a similar time, with the principal
difference being the initial magnetic field strength. We find that both runs
have very similar radial velocity profiles during these late-time averages,
although b5 experienced high-accretion rate transients at early times that
suppresses its density. This early-time behavior is clearly seen in the
profiles for b5_hi which clearly shows significantly higher radial inflow
within the central 30$r_{g}$.

Figure 6: Late-time snapshots of midplane density for b200 (left) and b5
(right), taken at $t=4.05\times 10^{5}~{}GM/c^{3}$ and $t=1.25\times
10^{5}~{}GM/c^{3}$ respectively.

Figure 7: Density along a midplane slice ($\theta=\pi/2$, $\phi=\\{0,\pi\\}$)
over time, for b200 (top) and b200_hi (bottom). The overdensity, or $m=1$
mode, is clearly visible in b200 as striations between radii $100-200r_{g}$.

The density profile of b200 shows a broad hump around $r=100r_{g}$. This
signposts a particularly interesting feature of this disk that, we argue, is a
subtle and pernicious imprint of the SMR de-refinement boundaries. The
overdensity, or $m=1$ mode, develops at a relatively early time of $t\sim
1\times 10^{5}~{}GM/c^{3}$ and is then sustained throughout the subsequent
evolution. This structure is most apparent in the mid-plane density map shown
in Figure 6 (left panel). The time-dependence of this lump can be seen by
forming space-time plots of the density along a midplane line that cuts
through the center, $(\theta=\pi/2,\,\phi=0,\pi)$; see Figure 7 (top panel).
The sloshing motion of this lump as it orbits the central body is seen as the
long-lived “see-saw” feature, causing an excavation of material from the
central region and creating the density maximum at $r\sim 100r_{g}$ (Figure
5). An examination of the mid-plane velocity field reveals this lump to be a
large scale vortex. It has been argued that such features might arise from a
Rossby Wave Instability (RWI; Lovelace et al., 1999) or generalized form of
the Papaloizou-Pringle Instability (PPI; Papaloizou & Pringle, 1984). The
presence of the RWI within disks has already been argued in the case of Tidal
Disruption Events (TDEs) and the spirals it produces may explain the observed
flaring behaviors in Sgr A* (Falanga et al., 2007; Mignon-Risse et al., 2021).
A similar feature was also reported in an MHD simulation by Mishra et al.
(2020), who observed spiral-like inhomogeneities in the density of disks with
initially weak poloidal fields, however these simulations were only run for 23
orbits (versus our $\sim 6500$ for b200) and so the late-time evolution of the
effect was not noted.

Figure 8: Velocity streamlines overlaid onto density for b200 (left) and b5
(right). Relative velocity magnitude is indicated by the linewidth of
streamlines. Stills are taken from late times: $t=4.05\times
10^{5}~{}GM/c^{3}$ for b200 and $t=1.25\times 10^{5}~{}GM/c^{3}$ for b5.

As seen in Mishra et al. (2020), the lump creates a kind of anomalous
accretion where material preferentially accretes from the atmosphere or
ceiling of the disk rather than in the midplane. This is clearly seen in
Figure 8 which shows the poloidal velocity field on a vertical cross-sectional
slice for late-time snapshots from b200 and b5. Run b200’s velocity
streamlines (Figure 8, left panel) show strong inflows in the jet regions of
the b200 simulation and eddies along the disk-jet boundary which take material
out of the disk and deposit it in the jet streams, while the radial velocity
within the disk is reduced (Figure 5).

The b5 disk, by comparison, does not develop such an overdensity and so
retains a significant amount of its mass below its $\rho_{\rm max}$ at
$20r_{g}$ (Figures 5 and 6, right panel). This smoother disk is also
accompanied by a greater inflow velocity $v_{r}$ and a higher temperature
$P/\rho c^{2}$. Thus stronger magnetic fields suppress the formation of the
vortex. More importantly, an examination of the corresponding space-time
diagram (Figure 7, bottom panel) shows that the lump is not observed in the
high-resolution simulation b200_hi. This fact together with the location of
the b200 lump, strongly suggests that the lump is an artifact caused by
vorticity generation at the SMR refinement boundary which affects only the
weaker-field case. This is a cautionary note for any long-duration accretion
disk models that employ SMR.

###  3.2 Magnetic Field Strength, Morphology, and Evolution

A primary motivation for our study is to understand the evolution of the
magnetic field, particularly the growth and morphology of the poloidal field,
in these disks that start with purely toroidal configurations. Furthermore, we
wish to understand whether the strength of the initial toroidal field has any
long lasting imprint on the magnetic evolution of the disks.

Figure 9: (Top) Magnetic flux, $\Phi_{B}=\int{B_{r}\cdot dS}$, through the
ISCO ($r=6r_{g}$). The signed flux is calculated for each hemisphere
individually before the absolute value is taken and the two hemispheres
averaged to give $|\langle\Phi_{B}\rangle|$. The values are re-normalized
according to Tchekhovskoy et al. (2011) (Section 3). (Bottom) Ratio of signed
and unsigned magnetic fluxes, indicating how much of the flux is aligned in a
common direction.

Our initially weak-field simulation, b200, shows rapid growth of the magnetic
field at early times, settling to $\langle\beta\rangle\sim 20$ (Figure 4),
compared to $\langle\beta\rangle\sim 10$ for the initially strong-field case
b5. Interestingly, we see the development of large-scale poloidal field in
both cases. Figure 9 shows the horizon-threading (hence necessarily poloidal)
flux as defined in equation 16. While there is clearly non-zero poloidal flux
in all cases, b200 shows a consistently low flux $\phi\sim 5$-$15$ with little
sustained growth, whereas b5 exceeds $\phi=40$ and nears the transitional
value for a MAD state of $\phi=50$ (Tchekhovskoy et al., 2012; Liska et al.,
2020) for a non-spinning BH, although fails to cross. The high-resolution
simulations b5_hi and b200_hi again match their lower-resolution analogues
remarkably well.

It is not merely the net horizon-crossing flux that matters, though, when
considering features like jet launching. The degree of order of this field as
it crosses the inner boundary is also important to questions such as whether
the jet is “striped” (Giannios & Uzdensky, 2019). Most disk models that start
with initially-poloidal field configurations end up with a well-organized
(unidirectional) horizon-threading flux bundle. However, that is not
necessarily the case for our models which begin with initially-toroidal field.
To examine this, we calculate the ratio of signed to unsigned fluxes threading
the inner boundary of our simulation,

$\phi_{\parallel}=\frac{\iint_{S}\boldsymbol{B}\cdot
d\boldsymbol{S}}{\iint_{S}|\boldsymbol{B}|\cdot d\boldsymbol{S}},$ (21)

where $S$ is a hemisphere surface at $r=6r_{g}$ and the final value given is
the average of the two hemispheres. A value of $\phi_{\parallel}=1$ indicates
agreement in the direction of the field, while a value $\phi_{\parallel}=0$
indicates that (for example) for every positive $B_{r}$ there is an equal
negative, canceling the flux everywhere on the surface integrated.

With this new measure, we see that b5 has a high degree of field order with
$\phi_{\parallel}\sim 1$ at late times, while b200 barely reaches
$\phi_{\parallel}\sim 0.8$ at $t\sim 2.5\times 10^{5}$ before decreasing again
and spending much of its time with $\phi_{\parallel}<0.6$. Again, we see a
high degree of agreement between the high- and low-resolution runs.

Figure 10: As in Figure 8, but with magnetic field streamlines.

These nature of these differences are seen more clearly in renderings of the
poloidal magnetic field morphology at late times in Figure 10. In b200, close
to the black hole, we see a well-defined field direction dominating the polar
regions and then a strong region of oppositely directed field along a lower-
latitude boundary funnel. The mid-plane field is weak and disorganized,
existing within small magnetic loops distributed seemingly at random.
Recalling that this region is flowing inwards, some of this magnetic topology
may be due to large-scale circulations between the disk atmosphere and the
polar regions. By contrast, b5 displays noticeably stronger fields in the
midplane of the disk, and where there are magnetic loops these are larger than
in b200. The polar region contains an organized and sustained flux bundle
along which there is outflowing material. It must be recognized, though, that
while $\phi_{\parallel}\sim 1$ for this case, the flux bundle still contains
interesting structure (field loops and current sheets) downstream, possibly
due to Kelvin-Helmholtz modes excited at the jet-disk boundary. Such large
magnetic field loops, particularly those on the right side of the upper jet
region and the left side of the lower jet region of Figure 10, are reminiscent
of the large-scale poloidal field loops identified by Liska et al. (2020), and
so leads us to tentatively conclude that simulation b5 is capable of
generating a large-scale poloidal field from an initially-strong toroidal
field, however evidence is much weaker for the weak-field case b200.

Our discussion of magnetic evolution so far has focused on the poloidal field,
but we note that there is also interesting dynamics in the toroidal field. We
recover the well-known butterfly patterns – a “feathering” pattern in the
space-time diagram for the azimuthally-averaged $B_{\phi}$ at a particular
radius (Figure 11). These are interpreted as dynamo-cycles, an emergent
phenomenon that arises from the non-linear MRI physics (Rincon et al., 2007;
Gressel & Pessah, 2015; Hogg & Reynolds, 2018b). In b200 (Figure 11, top
panel), the butterfly pattern is rapidly established and is then sustained for
the full duration of the simulation. This plot also explicitly shows that
(after azimuthal averaging) a non-zero but weak toroidal flux exists within
the polar region. These findings completely track over into the higher-
resolution version of this model b200_hi. By contrast, b5 establishes a
dominant toroidal field sign in the polar regions at a relatively early time
($t\sim 5\times 10^{4}~{}GM/c^{3}$), coincident with the creation of the
strong and organized poloidal flux bundle (Figure 9). From this time onward
the butterfly pattern is significantly weakened, showing that the magnetic
configuration in the polar region is affecting the operation of the dynamo
cycle.

While we have attributed its presence to the (artificial) SMR refinement
boundary, it is still interesting to consider the effect of the $m=1$ lump
found in b200 and discussed in Section 3.1. As the inflow of material from
large radius is impeded by the lump, a lower-density cavity forms within the
inner regions of the disk (Figure 5). This increases the magnetization of the
region (decreases the plasma-$\beta$) which may be responsible for weakening
the amplitude of the dynamo cycles in b200 from $t\approx 1.0\times
10^{5}~{}GM/c^{3}$ onward.

Figure 11: Butterfly diagrams of the azimuthally-averaged $B_{\phi}$ at
$r=25r_{g}$ showing signs of a sporadic butterfly pattern. Note there is a
change in output cadence for b5_hi (bottom) for $t>1.25\times
10^{4}~{}GM/c^{3}$.

##  4 Discussion and conclusion

In this work we have presented the results of two 3D MHD accretion disk
simulations with initially toroidal fields. These simulations were identical
in all respects bar their initial field strength, defined through the initial
magnetic plasma beta, $\beta=P_{\rm gas}/P_{\rm mag}$. One simulation, b200,
was initialized with $\beta=200$, while the other, b5, was initialized with
$\beta=5$. Additionally, two more high-resolution runs were done with double
the $\theta$ and $\phi$ resolution for comparison (b200_hi and b5_hi; for
details, see Section 2).

We find that our two fiducial simulations diverge early in their evolution and
remain on separate paths throughout each run, with notable differences in both
hydrodynamic and magnetic quantities. We show that these differences (with one
exception) are not sensitive to resolution for the values studied. As the only
difference between these two simulations is the initial field strength, we
conclude that it is this initial magnetization that determines the long-term
growth (or lack thereof) of the magnetic field. We produce only a weak large-
scale poloidal field in b200 despite evolving for a very long time ($4.0\times
10^{5}~{}GM/c^{3}$). We do, qualitatively, produce a large-scale poloidal
field in b5, and broadly reproduce the results of Liska et al. (2020) whose
simulation features a comparable strong toroidal field. This finding shows
that the results of Liska et al. (2020) can be reproduced in a pseudo-
Newtonian framework and do not rely on the inclusion of General Relativistic
physics or the peculiarities of spinning black holes.

All simulations produced clear–although disordered–butterfly patterns in the
sign of $\langle B_{\phi}\rangle$, indicating quasi-periodic cycles of the MRI
dynamo. We additionally find that all simulations, both weak and strong
fields, establish a dominant signed field in the jet region, however, in b200
and b200_hi this field is much weaker than that of the disk, and so may be
easily missed if quantities are not plotted with care.

Additionally, although neither simulation produces sufficient horizon-crossing
magnetic flux to exceed the MAD transitional value of $\phi=50$ for a non-
spinning BH, our initially-strong-field case b5 comes close at $\phi\sim 40$,
much higher than is typically seen in SANE disks. In disks such as those
considered here where the magnetic flux (and magnetic field strength though
the magnetic plasma-$\beta$) appears saturated but the MRI simultaneously
appears active, the question of SANE vs MAD relies heavily on the definition
of a MAD state, for which there is still heated debate (e.g. Begelman et al.,
2022).

Our fiducial high $\beta$ simulation b200 produced an overdensity at early
times that persisted throughout the run. No comparable overdensity was
observed in the high resolution counterpart b200_hi. The impact of the
overdensity on disk evolution is mixed: while there are clear changes in the
disk scale height, there are no clearly observable differences in the magnetic
field strength or structure according to the metrics used (magnetic plasma
$\beta$ and magnetic flux $\Phi$). Nonetheless, this finding highlights the
need for great care when choosing a simulation mesh, even at high resolutions.

PER thanks the Gates Cambridge Trust for research funding. This work was
supported, in whole or in part, by the Bill & Melinda Gates Foundation
[OPP1144]. Under the grant conditions of the Foundation, a Creative Commons
Attribution 4.0 Generic License has already been assigned to the Author
Accepted Manuscript version that might arise from this submission. CSR thanks
the STFC for support under the Consolidated Grant ST/S000623/1, as well as the
European Research Council (ERC) for support under the European Union’s Horizon
2020 research and innovation programme (grant 834203).

## References

  * Bai & Stone (2013) Bai, X.-N., & Stone, J. M. 2013, ApJ, 767, 30, doi: 10.1088/0004-637X/767/1/30
  * Balbus (2003) Balbus, S. A. 2003, ARA&A, 41, 555, doi: 10.1146/annurev.astro.41.081401.155207
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Balbus & Hawley (1998) —. 1998, Reviews of Modern Physics, 70, 1, doi: 10.1103/RevModPhys.70.1
  * Balbus et al. (1996) Balbus, S. A., Hawley, J. F., & Stone, J. M. 1996, ApJ, 467, 76, doi: 10.1086/177585
  * Begelman et al. (2022) Begelman, M. C., Scepi, N., & Dexter, J. 2022, MNRAS, 511, 2040, doi: 10.1093/mnras/stab3790
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883, doi: 10.1093/mnras/199.4.883
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433, doi: 10.1093/mnras/179.3.433
  * Childs et al. (2012) Childs, H., Brugger, E., Whitlock, B., et al. 2012, in High Performance Visualization–Enabling Extreme-Scale Scientific Insight, 357–372, doi: 10.1201/b12985
  * Davis & Tchekhovskoy (2020) Davis, S. W., & Tchekhovskoy, A. 2020, ARA&A, 58, 407, doi: 10.1146/annurev-astro-081817-051905
  * Evans & Hawley (1988) Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659, doi: 10.1086/166684
  * Event Horizon Telescope Collaboration et al. (2019) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019, ApJ, 875, L5, doi: 10.3847/2041-8213/ab0f43
  * Falanga et al. (2007) Falanga, M., Melia, F., Tagger, M., Goldwurm, A., & Bélanger, G. 2007, ApJ, 662, L15, doi: 10.1086/519278
  * Fishbone & Moncrief (1976) Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962, doi: 10.1086/154565
  * Flock et al. (2010) Flock, M., Dzyurkevich, N., Klahr, H., & Mignone, A. 2010, A&A, 516, A26, doi: 10.1051/0004-6361/200912443
  * Fragile & Sądowski (2017) Fragile, P. C., & Sądowski, A. 2017, MNRAS, 467, 1838, doi: 10.1093/mnras/stx274
  * Giannios & Uzdensky (2019) Giannios, D., & Uzdensky, D. A. 2019, MNRAS, 484, 1378, doi: 10.1093/mnras/stz082
  * Goodman & Xu (1994) Goodman, J., & Xu, G. 1994, ApJ, 432, 213, doi: 10.1086/174562
  * Gressel & Pessah (2015) Gressel, O., & Pessah, M. E. 2015, ApJ, 810, 59, doi: 10.1088/0004-637X/810/1/59
  * Guan et al. (2009) Guan, X., Gammie, C. F., Simon, J. B., & Johnson, B. M. 2009, ApJ, 694, 1010, doi: 10.1088/0004-637X/694/2/1010
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742, doi: 10.1086/175311
  * Hawley et al. (1996) —. 1996, ApJ, 464, 690, doi: 10.1086/177356
  * Hawley et al. (2011) Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84, doi: 10.1088/0004-637X/738/1/84
  * Hawley et al. (2013) Hawley, J. F., Richers, S. A., Guan, X., & Krolik, J. H. 2013, ApJ, 772, 102, doi: 10.1088/0004-637X/772/2/102
  * Hogg & Reynolds (2016) Hogg, J. D., & Reynolds, C. S. 2016, ApJ, 826, 40, doi: 10.3847/0004-637X/826/1/40
  * Hogg & Reynolds (2018a) —. 2018a, ApJ, 854, 6, doi: 10.3847/1538-4357/aaa6c6
  * Hogg & Reynolds (2018b) —. 2018b, ApJ, 861, 24, doi: 10.3847/1538-4357/aac439
  * Jiang et al. (2019) Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2019, ApJ, 880, 67, doi: 10.3847/1538-4357/ab29ff
  * King (2003) King, A. 2003, ApJ, 596, L27, doi: 10.1086/379143
  * Liska et al. (2020) Liska, M., Tchekhovskoy, A., & Quataert, E. 2020, MNRAS, 494, 3656, doi: 10.1093/mnras/staa955
  * Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805, doi: 10.1086/306900
  * Machida et al. (2001) Machida, M., Matsumoto, R., & Mineshige, S. 2001, PASJ, 53, L1, doi: 10.1093/pasj/53.1.L1
  * Mignon-Risse et al. (2021) Mignon-Risse, R., Aimar, N., Varniere, P., Casse, F., & Vincent, F. 2021, in SF2A-2021: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. Eds.: A. Siebert, ed. A. Siebert, K. Baillié, E. Lagadec, N. Lagarde, J. Malzac, J. B. Marquette, M. N’Diaye, J. Richard, & O. Venot, 113–116
  * Mishra et al. (2020) Mishra, B., Begelman, M. C., Armitage, P. J., & Simon, J. B. 2020, MNRAS, 492, 1855, doi: 10.1093/mnras/stz3572
  * Narayan et al. (2003) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69, doi: 10.1093/pasj/55.6.L69
  * Narayan et al. (2012) Narayan, R., Sądowski, A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426, 3241, doi: 10.1111/j.1365-2966.2012.22002.x
  * Noble et al. (2010) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2010, ApJ, 711, 959, doi: 10.1088/0004-637X/711/2/959
  * Papaloizou & Pringle (1984) Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721, doi: 10.1093/mnras/208.4.721
  * Penna et al. (2013) Penna, R. F., Sądowski, A., Kulkarni, A. K., & Narayan, R. 2013, MNRAS, 428, 2255, doi: 10.1093/mnras/sts185
  * Pessah (2010) Pessah, M. E. 2010, ApJ, 716, 1012, doi: 10.1088/0004-637X/716/2/1012
  * Rincon et al. (2007) Rincon, F., Ogilvie, G. I., & Proctor, M. R. E. 2007, Phys. Rev. Lett., 98, 254502, doi: 10.1103/PhysRevLett.98.254502
  * Salvesen et al. (2016) Salvesen, G., Simon, J. B., Armitage, P. J., & Begelman, M. C. 2016, MNRAS, 457, 857, doi: 10.1093/mnras/stw029
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Sorathia et al. (2010) Sorathia, K. A., Reynolds, C. S., & Armitage, P. J. 2010, ApJ, 712, 1241, doi: 10.1088/0004-637X/712/2/1241
  * Sorathia et al. (2012) Sorathia, K. A., Reynolds, C. S., Stone, J. M., & Beckwith, K. 2012, ApJ, 749, 189, doi: 10.1088/0004-637X/749/2/189
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137, doi: 10.1086/588755
  * Stone et al. (2020) Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4, doi: 10.3847/1538-4365/ab929b
  * Tchekhovskoy et al. (2012) Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2012, in Journal of Physics Conference Series, Vol. 372, Journal of Physics Conference Series, 012040, doi: 10.1088/1742-6596/372/1/012040
  * Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79, doi: 10.1111/j.1745-3933.2011.01147.x
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * White et al. (2020) White, C. J., Quataert, E., & Gammie, C. F. 2020, ApJ, 891, 63, doi: 10.3847/1538-4357/ab718e