# Radiation Transport Two-Temperature GRMHD Simulations of Warped Accretion
Disks

M.T.P. Liska  Institute for Theory and Computation, Harvard University, 60
Garden Street, Cambridge, MA 02138, USA  N. Kaaz Center for Interdisciplinary
Exploration & Research in Astrophysics (CIERA), Physics & Astronomy,
Northwestern University, Evanston, IL 60202, USA  G. Musoke  Anton Pannekoek
Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, The Netherlands  A. Tchekhovskoy Center for Interdisciplinary
Exploration & Research in Astrophysics (CIERA), Physics & Astronomy,
Northwestern University, Evanston, IL 60202, USA  O. Porth Anton Pannekoek
Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, The Netherlands

###### Abstract

In many black hole systems, the accretion disk is expected to be misaligned
with respect to the black hole spin axis. If the scale height of the disk is
much smaller than the misalignment angle, the spin of the black hole can tear
the disk into multiple, independently precessing ‘sub-disks’. This is most
likely to happen during outbursts in black hole X-Ray binaries (BHXRBs) and in
active galactic nuclei (AGN) accreting above a few percent of the Eddington
limit, because the disk becomes razor-thin. Disk tearing has the potential to
explain variability phenomena including quasi-periodic oscillations (QPOs) in
BHXRBs and changing-look phenomena in AGN. Here, we present the first
radiative two-temperature GRMHD simulation of a strongly tilted ($65^{\circ}$)
accretion disk around a $M_{BH}=10M_{\odot}$ black hole, which tears and
precesses. This leads to luminosity swings between a few percent and $50\%$ of
the Eddington limit on sub-viscous timescales. Surprisingly, even where the
disk is radiation pressure dominated, the accretion disk is thermally stable
over $t\gtrsim 14,000r_{g}/c$. This suggests warps play an important role in
stabilizing the disk against thermal collapse. The disk forms two nozzle
shocks perpendicular to the line of nodes where the scale height of the disk
decreases $10$-fold and the electron temperature reaches $T_{e}\sim
10^{8}-10^{9}K$. In addition, optically thin gas crossing the tear between the
inner and outer disk gets heated to $T_{e}\sim 10^{8}K$. This suggests that
warped disks may emit a Comptonized spectrum that deviates substantially from
idealized models.

##  1 Introduction

Accretion disks in black hole X-ray binaries (BHXRBs) and active galactic
nuclei (AGN) are typically described by axisymmetric, time-independent models.
A key factor that determines the most suitable model is the accretion rate
with respect to the Eddington limit. At the Eddington limit, radiation
pressure rivals gravity and plays a dynamically important role. The most well-
known models are the geometrically thin standard accretion disk model (Shakura
& Sunyaev, 1973), which describes sources accreting at an appreciable fraction
of their Eddington limit ($L\gtrsim 0.01-0.1\,L_{\rm Edd}$), and the
geometrically thick advection dominated accretion flow model (Narayan & Yi,
1994; Narayan et al., 2003), which describes highly sub-Eddington sources
($L\lesssim 0.01L_{\rm Edd}$). Slim disk models, in which the advection of
trapped photons dominates, work well near and above the Eddington limit (e.g.
Abramowicz et al. 1978). While these models can provide a reasonable
description of the multi-wavelength emission, they – by construction – do not
address the variability of their light curves and spectra. Such variability
can encode important information about the structure of the accretion disk and
corona (e.g. Kara et al. 2019; Uttley et al. 2014).

For example, recent observational evidence indicates that accretion disks in
changing-look AGN can undergo drastic changes in luminosity and spectral shape
over the course of less than a year (e.g. MacLeod et al. 2016; Yang et al.
2018a) or during quasi-periodic eruptions (QPEs) over the course of hours
(e.g. Yang et al. 2018b). In these AGN, the luminosity can change by up to two
orders of magnitude while the spectrum significantly hardens or softens. This
suggest a dramatic change in the accretion flow on timescales which are too
short to be explained by the viscous timescales inferred from the Maxwell
stresses created by magnetorotational instability (MRI, Balbus & Hawley 1991,
1998) induced turbulence (e.g. Lawrence 2018; Dexter & Begelman 2019). Various
physical mechanisms were proposed to reconcile theory with observations,
including instabilities in the accretion disk (Sniegowska et al., 2020),
magnetically elevated accretion disks (Begelman & Pringle, 2007), reprocessing
of point-source radiation by the outer accretion disk (e.g. Clavel et al.
1992), disk tearing (Raj & Nixon, 2021) and an orbiting compact object (e.g.
King 2020; Arcodia et al. 2021).

Similar to their changing-look AGN counterparts, the power spectra of BHXRBs
display a wide range of mysterious features ranging from broad-spectrum
variability to narrow spectral peaks known as quasi-periodic oscillations
(QPOs). Such features can encode unique information about the structure and
dynamics of the accretion disk in addition to the spin and mass of the black
hole. QPOs are typically divided into low- and high-frequency QPOs, which can
sometimes be observed together (e.g. Ingram & Motta 2019). Various physical
mechanisms have been proposed to explain QPOs. These include geometric effects
such as precession of a misaligned disk around the black hole spin axis (e.g.
Stella & Vietri 1998; Ingram et al. 2009) and intrinsic effects such as
parametric resonances (e.g. Kluzniak & Abramowicz 2002; Rezzolla et al. 2003;
Abramowicz & Kluźniak 2004) and disco-seismic modes (e.g. Kato 2004; Dewberry
et al. 2020a, b). All of these models make simplifying assumptions about the
underlying physics. Most significantly, they typically do not include
magnetized turbulence (though see Dewberry et al. 2020b; Wagoner & Tandon
2021), which can potentially dampen oscillatory modes. Thus, general
relativistic magnetohydrodynamic (GRMHD) simulations, which can simulate
accretion from first-principles, are very attractive for addressing the origin
of QPOs.

Musoke et al. (2022) recently demonstrated using GRMHD simulations (first
presented in Liska et al. 2019b) that accretion disks, which are misaligned
with the black hole spin axis, can be torn apart by the Lense-Thirring torques
(Lense & Thirring, 1918) due to the frame dragging of a spinning black hole.
Figure 1 shows that the inner part of an accretion disk tears off the parent
disk, leading to two independently precessing sub-disks. Such tearing events
lead to spikes in the accretion rate (see also Raj et al. 2021). Disk tearing
occurs in very thin disks when the viscosity cannot communicate the
differential Lense-Thirring torques. This is expected to occur when a disk
accretes at an appreciable fraction of the Eddington limit (see e.g. Esin et
al. 1997 for a review) and was observed in hydrodynamic simulations (Nixon et
al., 2012a, b; Nealon et al., 2015; Raj et al., 2021). Disk tearing can lead
to a rapid burst of accretion on timescales shorter than the viscous time and
is a promising mechanism to explain large luminosity and spectral swings
observed in accreting black holes (e.g. Nixon & Salvesen 2013; Raj & Nixon
2021).

Based on the simulations, Musoke et al. (2022) argued that low-frequency and
high-frequency QPOs can be explained by disk tearing. In their model, a disk
tears off at a near-constant radius ($r\sim 10^{1}r_{g}$) and precesses for
$\lesssim 10$ periods before falling into the black hole. While precessing at
$\nu=2.5$ Hz (for $M_{BH}=10M_{\odot}$), the disk emits a periodically
modulated light curve. Musoke et al. (2022) also found strong radial epicyclic
oscillations at $\nu\sim 56$ Hz in the inner disk, which they argue can
potentially explain some high-frequency QPOs. The geometric origin of low-
frequency QPOs is supported by recent observations, which suggest that that
iron line centroid frequency moves from blue-to-redshifted with the same phase
as the low-frequency QPO (Ingram & van der Klis, 2014; Ingram et al., 2016).
On the other hand, observations have been able to constrain neither an
intrinsic nor geometric origin of the much rarer high-frequency QPOs (e.g.
Ingram & Motta 2019).

Thus, the results of Musoke et al. (2022) suggest that disk tearing is a very
promising mechanism to explain the multi-wavelength variability in many
accreting black holes. However, all GRMHD simulations that found tearing
(Liska et al., 2019a, b) have relied on a cooling function (Noble et al.,
2009) to artificially set the temperature and scale height of the disk. In
reality, the temperature of the disk is determined by complex physics
involving viscous heating, radiative cooling and the advection of energy. Such
physics can modify the radial and vertical structure of the disk, such as its
scale height and vertical temperature profile. This makes it extremely
challenging to benchmark such GRMHD models against multi-wavelength
observations because they might not capture the important physics that sets
the disk thermodynamic state. For example, in radiation pressure supported
disks imbalance between the local dissipation rate and cooling rate can lead
to thermal and viscous instabilities that are not captured by a cooling
function (e.g. Lightman & Eardley 1974; Shakura & Sunyaev 1976; Sa̧dowski
2016; Jiang et al. 2019; Liska et al. 2022). In addition, when the inner disk
precesses, it will beam its radiation field periodically towards the outer
disk. This radiation field will scatter off the outer disk, and the radiation
back-reaction force can cause additional warping and dissipation (e.g. Pringle
1997; Wijers & Pringle 1999).

Here, we present the first radiation-transport two-temperature GRMHD
simulation of an accretion disk tilted by $65^{\circ}$ relative to a rapidly
spinning $M_{BH}=10M_{\odot}$ black hole of dimensionless spin of $a=0.9375$
accreting at $35\%$ the Eddington luminosity. We show that the disk undergoes
tearing. The accretion disk does not contain any large scale vertical magnetic
flux to begin with and does not launch a jet: this makes the simulation
applicable to the soft(-intermediate) states of BHXRBs (see Fender et al.
2004). In Sections 2 and 3 we describe the numerical setup and initial
conditions, before presenting and discussing our results in Sections 4 and 5,
and concluding in Section 6.

##  2 Numerical Setup

In this work we utilize the radiative version of our GPU accelerated GRMHD
code H-AMR (Liska et al., 2018a, 2019b). H-AMR evolves the ideal GRMHD
equations and radiation transfer equations in addition to the electron and ion
entropies. We use a (modified) spherical grid in Kerr-Shild foliation with
coordinates $r$, $\theta$ and $\phi$. Spatial reconstruction of primitive
variables is performed using a 1-dimensional PPM method, which guarantees
second order convergence in 3 dimensions. Magnetic fields are evolved on a
staggered grid as described in Gardiner & Stone (2005); White et al. (2016).
Inversion of conserved to primitive variables is performed using a
2-dimensional Newton-Raphson routine (Noble et al., 2006) or Aitken
Acceleration scheme (Newman & Hamlin, 2014) for the energy equation and a
1-dimensional Newton-Raphson routine for the entropy equation. The energy
based inversion is attempted first, and, if it fails, the entropy based
inversion is used as a backup. This dual energy formulation is now standard in
many GRMHD codes (Porth et al., 2019).

The radiative transfer equations are closed with a two-moment M1-closure
(Levermore, 1984) whose specific implementation is described in Sadowski et
al. (2013); McKinney et al. (2013). We additionally evolve (see e.g. Noble et
al. 2009; Ressler et al. 2015) the electron and ion entropy tracers
($\kappa_{e,i}=p_{e,i}/\rho^{\gamma_{e,i}}$) and include Coulomb collisions
(Stepney, 1983). We use adiabatic indices $\gamma_{e}=\gamma_{i}=5/3$ for the
electrons and ions (and $\gamma=5/3$ for the gas). This is appropriate for the
ion and electron temperatures in the accretion disk. We use a reconnection
heating model (Rowan et al., 2017) to distribute the dissipative heating
between ions and electrons. This leads to roughly $20-40\%$ of the dissipation
going into the electrons with the remainder heating the ions.

We include Planck-averaged bound-free, free-free, and cyclo-Synchrotron
absorption ($\kappa_{abs}$), emission ($\kappa_{em}$) and electron scattering
($\kappa_{es}$) opacities as given in McKinney et al. (2017) for solar
abundances $X=0.7$, $Y=0.28$, and $Z=0.02$. In addition, we account for
Comptonization through a local blackbody approximation (Sa̧dowski & Narayan,
2015). This works reasonably well when emission, absorption and scattering are
localized, but becomes inaccurate when the radiation field is anisotropic
and/or the energy spectrum deviates from a blackbody. Namely, the radiation
temperature used to calculate the absorption opacity is approximated as a
blackbody with $T_{r}=(\hat{E}_{rad}/a)^{0.25}$ with $\hat{E}_{rad}$ the fluid
frame radiation energy density and $a$ the radiation density constant. If the
radiation deviates from a blackbody this approximation can severely under or
overestimate the radiation temperature. We assume that only electrons can
absorb and emit radiation and set the ion opacity to $\kappa_{i}=0$. This is a
reasonable approximation in the accretion disk where bound-free processes,
which involve energy exchange between ions/electrons/radiation, only become
important when the density is high enough for Coulomb collisions to
equilibrate the ion and electron temperature. The limitations of the M1
closure relevant to our work are further addressed in Section 5.

To resolve the accretion disk in our radiative GRMHD models we use a base grid
resolution of $N_{r}\times N_{\theta}\times N_{\phi}=840\times 432\times 288$.
We then use beyond $r\gtrsim 4r_{g}$ up to 3 layers of adaptive mesh
refinement (AMR) to progressively increase this resolution in the disk from
the base resolution at $r\sim 4r_{g}$ to $N_{r}\times N_{\theta}\times
N_{\phi}=6720\times 2304\times 4096$ at $r\sim 10r_{g}$. (In addition, we use
static mesh refinement Liska et al. 2019b to reduce the $\phi$-resolution to
$N_{\phi}=16$ at the polar axis.) This is exactly half of the resolution used
in Liska et al. (2019b) and Musoke et al. (2022). Nevertheless, the turbulence
is still very well resolved with the number of cells per MRI-wavelength,
$Q^{i}=\lambda_{MRI}^{i}/N^{i}$, exceeding
$(Q^{\theta},Q^{\phi})\gtrsim(20,125)$ within the entire disk. We place the
radial boundaries at $R_{in}=1.13r_{g}\approx 0.84r_{\rm H}$ and
$R_{out}=10^{5}r_{g}$, which ensures that both boundaries are causally
disconnected from the flow; here $r_{\rm H}=1+(1-a^{2})^{1/2}r_{g}$ is the
event horizon radius. The local adaptive time-stepping routine in H-AMR
increases the timestep of each mesh block independent of the refinement level
in factors of 2 (up to a factor of 16) based on the local Courant condition
(Courant et al., 1928). This increases the effective speed of the simulations
several-fold while reducing the numerical noise (Chatterjee et al., 2019).

Since ideal GRMHD is unable to describe the physical processes responsible for
injecting gas in the jet funnel (such as pair creation) we use density floors
in the drift frame (Ressler et al., 2017) to ensure that ${\rho
c^{2}}\geq{p_{mag}}/12.5$. However, due to the absence of any poloidal
magnetic flux, no jets are launched in this work and thus the corresponding
density floors are never activated. Instead, the ambient medium is floored at
the density $\rho=10^{-7}\times r^{-2}$ and internal energy
$u_{tot}=10^{-9}\times r^{-2.5}$, where $u_{tot}=u_{gas}+u_{rad}$ is the sum
of the gas, $u_{gas}$, and radiation, $u_{rad}$, internal energies.

Figure 1: The first demonstration of disk tearing in a radiation pressure
supported accretion disk, as seen through iso-countour 3D renderings of
density at three different times. (Left panel) The disk at $t=28,955r_{g}/c$
forms a rigid body that is warped, but shows no differential precession.
(Middle panel) The disk tears off a precessing inner disk at $t=31,000r_{g}$
that precesses for two periods between $t=35,000r_{g}/c$ and
$t=37,000r_{g}/c$. While precessing, the inner disk aligns with the BH,
shrinking in size until it fully disappears. (Right panel) The disk forms a
new tear at $r\sim 25r_{g}/c$ around $t\sim 45,000r_{g}/c$, and the tearing
process repeats.

##  3 Physical Setup

To study the effect of disk tearing in a highly warped accretion disk we
evolve our model in two stages. In both stages we assume a rapidly spinning
black hole of mass $M_{BH}=10M_{\odot}$ and dimensionless spin $a=0.9375$. The
initial radial ($r$) and vertical ($z$) density profile of the disk is
$\rho(r,z)\propto r^{-1}\exp{(-{z^{2}}/{2h^{2}})}$. The disk extends from an
inner radius, $r_{in}=6r_{g}$, to the outer radius, $r_{out}=76r_{g}$ at a
constant scale height $h/r=0.02$. The covariant magnetic vector potential is
$A_{\theta}\propto(\rho-0.0005)r^{2}$ and normalized to give an approximately
uniform $\beta=p_{b}/p_{tot}\sim 7$ in the initial conditions. Here $p_{b}$ is
the magnetic pressure and $p_{tot}$ is the sum of the electron/ion pressures,
$p_{e,i}$, and the radiation pressure, $p_{rad}$. The whole setup is
subsequently rotated by 65 degrees with respect to the black hole spin axis,
which itself is aligned with polar axis of the grid. Physical equivalence
between tilting the disk (e.g. Liska et al. 2018a; White et al. 2019) and
tilting the black hole (e.g. Fragile & Anninos 2005; Fragile et al. 2007) is
demonstrated in the Appendix of Liska et al. (2018a).

In the first stage (model T65), which is described in Liska et al. (2019b), we
evolve the disk in ideal GRMHD for $\Delta t=150,000r_{g}/c$ with a preset
scale height of $h/r=0.02$ set by a cooling function (Noble et al., 2009). The
disk undergoes a large tearing event around $t\sim 45,000r_{g}/c$, which is
analysed in Musoke et al. (2022). During this tearing event a sub-disk of size
$\Delta r\sim 10-15r_{g}$ tears off and precesses for $\sim 6$ periods between
$t\sim 45,000r_{g}/c$ and $t\sim 80,000r_{g}/c$. While precessing, the inner
disk exhibits a very prominent radial epicyclic oscillation at $\nu=56$ Hz for
a 10 solar mass BH. Another tearing event between $t\sim 120,000r_{g}/c$ and
$t\sim 145,000r_{g}/c$ produces a similar oscillation albeit at a slightly
higher frequency of $\nu=69$ Hz. These frequencies are consistent with
observed high frequency QPOs (e.g. Ingram & Motta 2019).

In the second stage (model RADT65), which is the focus of this work, we
restart this simulation in full radiative GRMHD at $t=29,889r_{g}/c$ as
described above (after downscaling the resolution by a factor $2$). This is
well before the occurrence of any large tearing event in stage 1. We assume
thermodynamic equilibrium between gas and radiation by setting
$T_{r}=T_{e}=T_{i}$ in the initial conditions and run this model (RADT65) for
$\Delta t\sim 14,000r_{g}/c$ until $t\sim 45,000rg/c$. The density scaling is
set such that the average accretion rate corresponds to approximately
$\dot{M}/\dot{M}_{\rm Edd}\sim 0.35$ with $M_{\rm
Edd}=\frac{1}{\eta_{NT}}L_{\rm Edd}/c^{2}$ the Eddington accretion rate,
$L_{\rm Edd}=\frac{4\pi GM_{BH}}{ck_{es}}$ the Eddington luminosity, and
$\eta_{NT}=0.179$ the Novikov & Thorne (1973) efficiency. According to Piran
et al. (2015) this choice of parameters is roughly consistent with a thin disk
solution (Sadowski, 2011) of scale height of $h/r=0.02-0.06$ between
$r=r_{isco}$ and $r=100r_{g}$.

Figure 2:  Panels (a, b) Space-time diagrams of density, $\rho$, and warp
amplitude, $\psi$. The warp amplitude peaks where the disk is about to tear or
has already torn. Since dissipation is enhanced in warps, maxima in warp
amplitude correspond to minima in density. The main tearing event is denoted
by black dotted lines. Panels (c, d) Space-time diagrams of the electron,
$T_{e}$, and ion, $T_{i}$, temperatures. The ion and electron temperature
correlate with the warp amplitude and anti-correlate with the density. Close
to the black hole and in the tear the ions are significantly hotter than the
electrons since Coulomb collisions are unable to equilibrate the ion and
eletron temperatures. Panels (e, f) Space-time diagrams of the tilt,
$\mathcal{T}$, and precession, $\mathcal{P}$, angle of the disk. Around $t\sim
36,500r_{g}/c$ a disk tears off and precesses for $\lesssim 2$ periods, during
which it gradually aligns with the black hole spin axis.

##  4 Results

###  4.1 Tearing Process

The disk in RADT65 undergoes a tearing event between $t=31,000r_{g}/c$ and
$t=40,000r_{g}/c$, as visualised using a 3D density isocontour rendering at
three different times in Fig. 1). Figure 2 shows the space-time diagram, in
which we calculate the average density of the disk
$\bar{\rho}=\frac{\int^{2\pi}_{0}\int^{\pi}_{0}\rho^{2}\sqrt{-g}d\theta
d\phi}{\int^{2\pi}_{0}\int^{\pi}_{0}\rho\sqrt{-g}d\theta d\phi}$ (Fig. 2a),
the warp amplitude, $\Psi=r\frac{\partial\vec{l}}{\partial r}$ with $\vec{l}$
being the angular momentum (Fig. 2b), the average electron and ion
temperatures,
$T_{e,i}=\frac{\int^{2\pi}_{0}\int_{0}^{\pi}p_{e,i}\sqrt{-g}d\theta
d\phi}{\int^{2\pi}_{0}\int^{\pi}_{0}\rho\sqrt{-g}d\theta d\phi}$ (Fig. 2c,d),
the tilt and precession angles, $T_{disk}$ and $P_{disk}$ (Fig. 2e,f as
calculated in Fragile & Anninos 2005).

During this tearing event the disk precesses for $\lesssim 2$ periods between
$36,000<t<38,000r_{g}/c$ (Fig. 2e) before it aligns with the spin axis of the
black hole (Fig. 2f). Alignment happens through warp driven dissipation of
misaligned angular momentum on the accretion time of the inner disk, which
differs from the Bardeen-Petterson alignment (Bardeen & Petterson, 1975)
mechanism that manifests itself on much shorter timescales during which the
structure and density of the disk remains in a steady state (see discussion in
Liska et al. 2018b). In addition to this ‘main’ tearing event there are
indications, such as the formation of an inwards moving ring-like structure
and an increase in the warp amplitude around $r\sim 30r_{g}$, that a much
larger tearing event will occur after the end of this simulation. Due to the
large computational expense associated with extending the duration of RADT65,
we leave the analysis of this (potentially) much longer tearing event to
future work.

Whether a disk is able to tear depends on how the internal torques in the disk
react to differential Lense-Thirring precession (e.g. Nixon & King, 2012). The
Lense-Thirring precession rate of a point particle around a rapidly spinning
black hole follows a steep radial dependence of $\nu_{LT}\propto a/r^{3}$.
However, an accretion disk will precess with a single frequency set by the
integrated LT torque in the disk if the viscous torques are able to counteract
the differential precession rate (e.g. Fragile & Anninos 2005; Fragile et al.
2007). Naively, one expects that the warp amplitude decreases monotonically
with the increasing distance from the black hole. This is (roughly) the case
in the inner and outer disk of RADT65 for $t\lesssim 32,000r_{g}/c$.

As Fig. 2(a),(b) shows, the warp amplitude and surface density in RADT65 are
anti-correlated for the entire runtime. At $r\lesssim 20r_{g}$ the warp
amplitude is rather high and a cavity of low density gas forms at both early
and late times. In this cavity the plasma thermally decouples into hot
electrons with $T_{e}\sim 10^{8}-10^{9}K$ and extremely hot ions with
$T_{i}\sim 10^{9}-10^{10}K$. The strong anti-correlation between warp
amplitude and density suggests that accretion is driven by dissipation in
warps (e.g. Papaloizou & Pringle 1983; Ogilvie 1999; Nelson & Papaloizou
2000). In a companion paper (Kaaz et al 2022, submitted) we demonstrate that
warp-driven dissipation is responsible for the bulk of the dissipation and
that magneto-rotational instability driven turbulence (e.g. Balbus & Hawley
1991, 1998) plays only a very subdominant role.

However, after $t\gtrsim 32,000r_{g}/c$ the warp amplitude is not a
monotonically decreasing function of radius. Instead, the warp amplitude forms
local minima and maxima, and concentric rings of higher density gas centered
at the (local) minima of the warp amplitude form. It is known that when the
warp amplitude exceeds a critical value, it can become unstable (Doǧan et al.,
2018; Doğan & Nixon, 2020). When this happens, the viscous torques in the warp
drop to zero and the disk tears. In RADT65 the tearing process starts at
$r\sim 19r_{g}$ around $33,000r_{g}/c$ and the disk fully detaches at $r\sim
10r_{g}$ around $t\sim 36,500r_{g}/c$. We see later that during this process
warp-driven dissipation leads to a temporary burst in the BH accretion rate
and luminosity around $t\sim 36,500r_{g}/c$ (Sec. 4.3). The disk even
(briefly) enters the slim disk regime where $\dot{M}\sim\dot{M}_{edd}$ and
$L\sim 0.5L_{edd}$, suggesting the luminosity is slightly suppressed due to
photon trapping (e.g. Abramowicz et al. 1978).

Figure 3:  A volume rendering of the density during the main tearing event at
$t=37,026r_{g}/c$. The red plane cuts the disk at the $\phi=0$ surface.
Streamers from the outer disk (green) deposit low angular momentum gas onto
the inner disk (blue). This causes the inner accretion disk to shrink in size
(Kaaz et al. 2022, submitted).

As the gas from the inner disk falls into the black hole it transports angular
momentum outwards. Angular momentum conservation dictates that the radius of
the inner disk should increase. This process, called viscous spreading,
affects all finite size accretion disks that are not resupplied externally
with gas (e.g. Liska et al. 2018a; Porth et al. 2019). However, the inner disk
in RADT65 does not spread viscously, but, in fact, decreases in size. This
discrepancy could potentially be explained by cancellation of misaligned
angular momentum (e.g. Nixon & King 2012; Hawley & Krolik 2018). Namely, when
gas from the outer disk falls onto the inner disk the misaligned components of
angular momentum cancel out leading to a net decrease in the inner disk
radius. In addition, it is possible that the inner disk transfers some of its
excess angular momentum to the outer disk. In a companion paper (Kaaz et al.
2022, submitted) we quantify the relative contributions of the angular
momentum cancellation and outward angular momentum transport to the evolution
of the tearing radius.

Figure 4: A top down view of the accretion disk during the main tearing event
at $t=37,026r_{g}/c$ corresponding to the volume rendering in Figure 3. (Panel
a) The surface density, $\Sigma$, drops sharply in the tear due to additional
dissipation as the gas undergoes a radical plane change. (Panel b) The
vertically integrated scattering optical depth ($\tau_{es}$) drops below unity
in the tear suggesting optically thin emission. (Panels c, d) The electron,
$T_{e}$, and ion, $T_{i}$, temperatures peak perpendicular to the line of
nodes (due to a ‘nozzle’ shock, see Sec. 4.2) and in the tear at $r\sim
10r_{g}$. (Panel e) The density scale height, $\theta_{d}$, of the accretion
disk oscillates in azimuth. It peaks along the x-axis and reaches a minimum
along the y-axis. (Panel f) The bolometric luminosity, $\lambda$, forms a non-
axisymmetric pattern ranging from $\lambda\sim 0.15L_{\rm Edd}$ to
$\lambda\sim 0.5L_{\rm Edd}$. (Panels g, h) The gas entropy entropy,
$s_{gas}$, does not increase substantially when gas crosses the scaleheight
compression along the vertical y-axis, but the radiation entropy, $s_{rad}$,
does. This suggests the shock heating is radiatively efficient.

###  4.2 Dissipation and Radiative Signatures

To better understand the structure of the disk during the tearing event at
$t=37,026r_{g}/c$ we show a volume rendering of the density in Figure 3 and a
vertical projection of several quantities in Figure 4. In this projection the
line of nodes, where the disk crosses the equatorial plane, is aligned with
the horizontal x-axis. The projected quantities include the surface density
$\Sigma=\int_{0}^{\pi}\rho\sqrt{g_{\theta\theta}}d\theta$ (panel a), the
scattering optical depth $\tau_{es}\sim 0.34\Sigma$ (panel b), electron and
ion temperature
$T_{e,i}=\frac{\int_{0}^{\pi}p_{e,i}\sqrt{g_{\theta\theta}}d\theta}{\Sigma}$
(panels c and d), density scale height
$\theta_{d}=\frac{\int_{0}^{\pi}\rho(\theta-\pi/2)\sqrt{g_{\theta\theta}}d\theta}{\Sigma}$
(panel e), effective bolometric luminosity
$\lambda=\frac{\int_{0}^{\pi}R^{r}_{t}\sqrt{g_{\theta\theta}}2\pi
d\theta}{L_{edd}}$ (panel f), gas entropy per unit mass
$s_{gas}=\frac{\int_{0}^{\pi}p_{e,i}/\rho^{\gamma}\rho
u^{t}\sqrt{g_{\theta\theta}}d\theta}{\Sigma}$ (panel e), and radiation entropy
per unit mass
$s_{rad}=\frac{\int_{0}^{\pi}E_{rad}^{3/4}\sqrt{g_{\theta\theta}}d\theta}{\Sigma}$
(panel h). Integration is performed in a spherical coordinate system aligned
with the local angular momentum vector of the disk. An animation of Figure 3
and panels a, c and f of Figure 4 can be found on our YouTube channel (movie).

Interestingly, the electron and ion temperatures do not decline smoothly with
distance from the black hole (Fig. 4b,c) as predicted by idealized models of
thin accretion disks Novikov & Thorne (1973). At the tearing radius of
$10r_{g}$, where the warp amplitude reaches a maximum, the electron
temperature reaches $T_{e}\gtrsim 10^{8}K$ and the plasma becomes optically
thin. Here streamers of gas get thrown onto highly eccentric orbits and
subsequently rain down on the inner disk (Fig. 3). The rapid rise in
temperature in the tear suggests that gas crossing the tear is subject to
additional dissipation as it undergoes a rapid orbital plane change. If we
assume the gas has a temperature of $T_{i}=0K$ before crossing the tear at
$r\sim 7r_{g}$, and the gas dissipates $\epsilon=10^{-2}$ of its orbital
kinetic energy $U_{kin}\approx\frac{1}{2r}$, we find that dissipation can
easily heat the gas to $T_{i}\sim 10^{10}K$. In reality, $\epsilon$ can be
much higher if radiative cooling is efficient.

The scale height of the disk exhibits a prominent $m=2$ azimuthal oscillation
(Fig. 4e). This can also be observed in Figure 5 where we show a $\theta-\phi$
slice at $r=25r_{g}$ of the density, electron temperature, and gas entropy. At
the line of nodes the scale height reaches a maximum of $h/r\sim 0.1$ while
perpendicular to the line of nodes the scale height drops to $h/r\sim 0.005$.
Around this nozzle-like ‘compression’, the gas reaches a temperature ranging
from $T_{e}\sim 2\times 10^{7}K$ at $r\sim 20r_{g}$ to $T_{e}\gtrsim 10^{8}K$
within $r\lesssim 5r_{g}$. At select times $T_{e}$ can even reach $T_{e}\sim
10^{9}K$ closer to the black hole (see YouTube movie). The radial and
azimuthal fluctuations in the temperature lead to a highly non-axisymmetric
emission pattern (Fig. 4f). At this specific snapshot ($t=37,026r_{g}/c$) the
effective bolometric luminosity varies azimuthally between $\lambda\sim
0.15L_{edd}$ and $\lambda\sim 0.5L_{edd}$ at $r=200r_{g}$.

Figure 5: A $\theta-\phi$ slice through the disk at $r=15r_{g}$ in tilted
coordinates with the ascending node at $\phi=1/2\pi$ and the descending node
at $\phi=3/2\pi$. The pink and black lines correspond to the $\tau_{es}=1$ and
$\tau_{bf}=1$ surfaces respectively. (Upper panel) The gas density, $\rho$,
peaks in the ‘nozzle’, which is located perpendicular to the local line of
nodes at $\phi=0$ and $\phi=\pi$. (Middle panel) The electron temperature,
$T_{e}$, increases by a factor $\sim 5$ in the nozzle, suggesting the
radiative emission will be hardened. (Lower panel) The gas entropy per unit
mass, $\kappa_{gas}$, does not increase in the nozzle, suggesting the rise in
gas temperature is adiabatic (Sec. 4.2)

The specific gas entropy does not increase in the nozzle (Fig. 4g), suggesting
heating of gas in the nozzle is adiabatic. This might seem surprising since
nozzles in tidal disruption events (TDEs) are typically associated with shock
heating and steep gradients in the gas entropy (e.g. Rees 1988; Kochanek 1994;
Andalman et al. 2022). However, if the radiative cooling timescale is very
short, shocks do not automatically lead to an increase of the gas entropy.
Instead, they will lead to an increase in the radiation entropy. This increase
in radiation entropy is visible along the vertical y-axis in Fig. 4(g). In a
companion paper (Kaaz et al 2022, submitted), we demonstrate that a shock
forms in the nozzle that dissipates roughly $\sim 1.8\%$ of the orbital energy
each time the gas passes through it.

###  4.3 Outflows and Variability

To understand the variability associated with disk tearing and the energetics
of the outflows we plot in Figure 6 the mass accretion rate (panel a),
radiative emission rate (panel b) and (radiative) outflow efficiencies (panel
c). The mass accretion rate measured at the event horizon
($\dot{M}=\int_{0}^{2\pi}\int_{0}^{\pi}\rho u^{r}\sqrt{-g}d\theta d\phi$ with
$g$ the metric determinant and $u^{\mu}$ the fluid 4-velocity) exhibits a
factor $\sim 40$ variation in time, ranging from $\dot{M}=2\times
10^{-2}\dot{M}_{\rm Edd}$ to $\dot{M}=\dot{M}_{\rm Edd}$ (defined in Sec. 3).
However, the luminosity
($L=\int_{0}^{2\pi}\int_{0}^{\pi}R^{r}_{t}\sqrt{-g}d\theta d\phi$ with
$R^{r}_{t}$ the $r$-$t$ component of the radiation stress energy tensor) only
exhibits a factor $\sim 5$ variation and peaks at $L\sim 0.5L_{edd}$.

The outflow efficiencies with respect to the BH mass accretion rate are
defined at the event horizon as
$\eta_{wind}=\frac{\dot{M}-\dot{E}}{|\dot{M}|_{t}}$,
$\eta_{rad}=\frac{L_{r=r_{200r_{g}}}}{|\dot{M}|_{t}}$, and
$\eta_{adv}=\frac{-L_{r=r_{BH}}}{|\dot{M}|_{t}}$. Here
$\dot{E}=\int_{0}^{\pi}\int_{0}^{2\pi}T^{r}_{t}\sqrt{-g}d\theta d\phi$ is the
energy accretion rate with $T^{r}_{t}$ the stress energy tensor, and,
$|\dot{M}|_{t}$ is a running average of $\dot{M}$ over a time interval of
$500r_{g}/c$. Unless stated otherwise, we measure $\dot{M}$ and $\dot{E}$ at
the event horizon, $L$ at the event horizon for $\eta_{adv}$ and at
$r=200r_{g}$ for $\eta_{rad}$. As Fig. 6(c) shows, when the accretion rate
reaches the Eddington limit the advective efficiency, $\eta_{adv}$, is of
comparable magnitude to the radiative efficiency, $\eta_{rad}$ (Fig. 6d). This
implies that a significant fraction of the emitted radiation will fall into
the black hole before it reaches the observer suggesting the disk might have
entered the slim disk regime at the peak of luminosity. Nevertheless, the
time-averaged radiative efficiency reaches $|\eta_{rad}|=14.7\%$, which is
only slightly below the canonical Novikov & Thorne (1973) efficiency of
$\eta_{NT}=17.8\%$.

We do not observe any evidence for a radiation or magnetic pressure driven
wind. While the outflow efficiency measured at the event horizon exceeds
$\eta_{wind}\gtrsim 15\%$ this drops to $\eta_{wind}\lesssim 0.1\%$ when
measured at $r=200r_{g}$, which demonstrates that both the magnetic pressure
(e.g. Liska et al. 2018b) and radiation pressure (e.g. Kitaki et al. 2021) are
insufficient to accelerate the wind to escape velocities. This suggests that
poloidal magnetic fields are a key ingredient to accelerate sub-relativistic
outflows, even when the accretion rate approaches the Eddington limit.

Figure 6: Time evolution of RADT65 with the main tearing event denoted by
black dotted lines. (Panel a) The mass accretion rate increases from
$\dot{M}\sim 10^{-2}\dot{M}_{\rm Edd}$ to $\dot{M}\sim 1\dot{M}_{\rm Edd}$.
(Panel b) The luminosity increases from $L\sim 0.1L_{\rm Edd}$ to $L\sim
0.5L_{\rm Edd}$. (Panel c) The tilt $\mathcal{T}$ of the accretion disk. The
inner accretion disk ($r\lesssim 10r_{g}$) aligns with the black hole just
after the main tearing event while the outer accretion disk remains
misaligned. (Panel d) The wind efficiency, $\eta_{wind}$, radiative
efficiency, $\eta_{rad}$, advective radiative efficiency, $\eta_{adv}$, and
Novikov & Thorne (1973) efficiency, $\eta_{NT}$, and total efficiency,
$\eta_{tot}$, as defined in Sec. 4.3. When the accretion rate peaks, advection
of radiation becomes important and suppresses the total luminosity of the
system. (Panel e) The thermal scale height of the disk, $\theta_{t}$, stays
stable or increases throughout the disk, which suggests the disk is thermally
stable (Sec. 4.5).

###  4.4 Radial Structure

To better understand the structure of the accretion disk we show in Figure 7
time-averaged (between $t=[43,000,44,000]r_{g}/c$) radial profiles of the
density, $\bar{\rho}$ (panel a), luminosity, $L$ (panel b), and radiation,
ion, electron, and magnetic pressure, $p_{rad}$, $p_{i}$, $p_{e}$, and
$p_{b}$, respectively (panel c). During this time interval a new tear starts
to form around $r\sim 30r_{g}$, which manifests itself as a drop in density
and increase in the radiative flux. While the magnetic and radiation pressure
are in approximate equipartition within $r\lesssim 20r_{g}$, the accretion
disk becomes fully radiation pressure dominated between $r=[20,40]r_{g}$. The
sum of the ion and electron pressure is a factor $\sim 10-100$ lower than the
radiation pressure, suggesting that gas pressure does not play a dynamically
important role.

We have already mentioned in Section 4.1 that warp driven dissipation plays an
important role in RADT65. To illustrate how much warp-driven dissipation
accelerates the inflow of gas in RADT65, we translate the radial infall speed
of the gas into an effective viscosity, $\alpha_{\rm eff}=\frac{\langle
v^{r}v^{\phi}\rangle_{\rho}}{\langle c_{s}^{2}\rangle_{\rho}}$ with
$v^{r}=u^{r}/u^{t}$ and $v^{\phi}=u^{\phi}/u^{t}$. Here, $\langle
x\rangle_{y}=\frac{\int_{0}^{2\pi}\int_{0}^{\pi}xy\sqrt{-g}d\theta d\phi
d\phi}{\int_{0}^{2\pi}\int_{0}^{\pi}y\sqrt{-g}d\theta d\phi}$ gives the
$y-$weighted average of quantity $x$. In Fig. 7(c), we find $\alpha_{\rm
eff}\sim 1-5\times 10^{1}$ throughout the disk. The Maxwell stresses
$\alpha_{\rm m}=\frac{\langle b^{r}b^{\phi}\rangle_{\rho}}{\langle
p_{tot}\rangle_{\rho}}\sim 10^{-2}-10^{-1}$ induced by magnetized turbulence
cannot account for this effective viscosity, which suggests that other
physical processes such as shocks (e.g. Fragile & Blaes 2008) drive accretion
in warped accretion disks.

###  4.5 Thermal Stability

Interestingly, Liska et al. (2022) found that an aligned version of the
accretion disk in RADT65 collapsed within $t\lesssim 5,000r_{g}/c$ to an
exceedingly thin slab that could not be resolved numerically. This runaway
cooling (see also Sa̧dowski 2016; Fragile et al. 2018; Jiang et al. 2019;
Mishra et al. 2022) is a manifestation of the thermal instability (Shakura &
Sunyaev, 1976). However, we do not observe any evidence for the thermal
instability in RADT65. In fact, as Fig. 6(e) demonstrates, the thermal scale
height, $\theta_{t}=\frac{\langle c_{s}\rangle_{\rho}}{\langle
v^{\phi}\rangle_{\rho}}$, remains stable or even increases depending on the
radius.

This unexpected result could potentially be explained within $r\lesssim
20r_{g}$ by magnetic fields, which can stabilize a thermally unstable
accretion disk if they reach equipartition with $p_{b}\sim p_{rad}$ (e.g.
Begelman & Pringle 2007; Sa̧dowski 2016; Jiang et al. 2019; Liska et al. 2022;
Mishra et al. 2022). However, magnetic fields will not be able to stabilize
the disk between $r=[20,40]r_{g}$ where the disk is strongly radiation
pressure dominated with $p_{rad}/p_{b}\sim 10^{1}-10^{2}$.

We found in Sec. 4.3 that advection of radiation becomes significant at the
event horizon. If the radiation diffusion timescale becomes longer than the
accretion timescale, photon trapping might stabilize the disk against runaway
cooling (e.g. Abramowicz et al. 1978). To test if photon trapping can explain
the apparent thermal stability of the disk we plot the radiation diffusion
time, $t_{\rm diff}\sim\frac{h}{v_{rad}}\sim 3h\tau_{es}$ and the accretion
time, $t_{\rm acc}=\frac{r}{\langle v_{r}\rangle_{\rho}}$, in Fig. 7(e). Here
we use for $v_{rad}$ the optically thick radiation diffusion speed $v_{\rm
rad}\sim\frac{1}{3\tau_{es}}c$ (e.g. McKinney et al. 2013). As is evident from
Fig. 7(e), while the radiation diffusion and accretion timescales are
comparable in the inner disk, the radiation diffusion timescale becomes much
shorter than the accretion timescale in the outer disk. This suggests that
photon trapping will not be able to stabilize the outer disk against runaway
cooling or heating.

Instead, we propose in Section 5 that nozzle shock driven accretion disks are
inherently thermally stable.

Figure 7: Radial structure of the disk time-averaged between
$t=[43,000;44,000]r_{g}/c$. This is the end state of our simulation at which
point a new tear starts to form. (Panel a) The gas density, $\bar{\rho}$,
peaks around $r\sim 25r_{g}$. We expect the new tear to form just behind the
peak in density. (Panel b) The cumulative luminosity of the disk increases
until $35r_{g}$, partially driven by enhanced dissipation in the forming tear
between $r=20r_{g}$ and $r=30r_{g}$. (Panel c) The disk is radiation pressure,
$p_{rad}$, dominated between $r=[15,40]r_{g}$ and magnetic pressure, $p_{b}$,
dominated everywhere else. The ion, $p_{i}$), and electron, $p_{e}$, pressures
are very weak compared to either the radiation or magnetic pressure and thus
do not play an dynamically important role. (Panel d) The effective viscosity,
$\alpha_{\rm eff}$, is much larger than the Maxwell, $\alpha_{\rm m}$, stress.
This suggests that dissipation is not driven by magnetized turbulence induced
by magnetic stresses. Instead, we demonstrate in a companion paper that shocks
are driving accretion (Kaaz et al 2022, submitted). (Panel e) The photon
diffusion timescale, $t_{\rm diff}$, is shorter than the accretion timescale,
$t_{\rm acc}$, except very close to the black hole. This suggest advection of
radiation internal energy is only dynamically important very close to the
black hole.

##  5 Discussion

###  5.1 Applicability to Spectral States

The accretion disk in RADT65 does not launch any relativistic jet. This
suggests RADT65 is applicable to the high-soft and potentially soft-
intermediate spectral states. In these spectral states the spectrum is
dominated by thermal blackbody radiation and no radio jets are detected
(though see Russell et al. 2020 for detection of compact jets in the
infrared).

Nevertheless, the emission profile of RADT65 presents a radical departure from
the Novikov & Thorne (1973) model of a thin accretion disk. Instead of an
axisymmetric power-law emission profile, emission in warped accretion disks
can become hardened around the nozzle shock or at the tearing radius. Here the
temperature of plasma can reach $T_{e}\sim 10^{8}-10^{9}K$, which might lead
to Comptonization of cold accretion disk photons. This could potentially
contribute to the high-energy power law emissions tail observed in the spectra
of soft state BHXRBS (e.g. Remillard & McClintock 2006) in addition to the
soft X-Ray excess in the spectra of high-luminosity AGN (e.g. Gierliński &
Done 2004, 2006). Note that part of this high-energy emission can come from
within the ISCO as suggested by recent GRMHD simulations (Zhu et al., 2012)
and semi-analytical models (Hankla et al., 2022). Future general relativistic
ray-tracing (GRRT) calculations will need to quantify how much the spectrum
deviates from a blackbody spectrum and in which wavelength bands this
optically thin shock-heated gas can be detected. The energy spectrum of the
accelerated electrons will depend on the microphysics of the shocks, which
will need to be addressed by particle-in-cell simulations.

In the hard-intermediate state (HIMS) jetted ejections imply that the
accretion disks might be saturated by vertical magnetic flux (e.g. McKinney et
al. 2012; Tchekhovskoy et al. 2011; Chatterjee et al. 2020). When an accretion
disk is saturated by magnetic flux it enters the magnetically arrested disk
(MAD) regime (Narayan et al., 2003) in which magnetic pressure overcomes
gravity and accretion proceeds through non-axisymmetric instabilities. Recent
two-temperature radiative GRMHD simulations (Liska et al., 2022) have
demonstrated that magnetic flux saturation can lead to ‘magnetic’ truncation
during which the inner disk decouples into a magnetic pressure supported
‘corona’ of hot electrons ($T_{e}\gtrsim 5\times 10^{8}K$) and extremely hot
ions ($T_{i}\gtrsim 10^{10}K$). This happens within the magnetic truncation
radius, which is determined by the amount and location of the excess poloidal
magnetic flux in the system. It is conceivable that the combined effects of
magnetic reconnection in current sheets and shocks in warps can heat the
plasma to temperatures far above what was observed in this work or in Liska et
al. (2022). Future work will need to address if, where and how such misaligned
truncated disks tear, and, if their spectral and timing signatures are
consistent with sources in the HIMS.

###  5.2 Quasi-periodic Variability

This work and Musoke et al. (2022) have demonstrated that disk tearing can
lead to short luminosity outbursts of length $\Delta t\sim 1000-3000r_{g}/c$
every $\Delta t_{tear}\sim 10,000-50,000r_{g}/c$. Coincident with a tearing
event the disk precesses on a timescale of $\Delta t_{prec}\sim
1000-6000r_{g}/c$ suggesting tearing events might be accompanied by a
sinusoidial modulation of the light curve. In parallel, the mass-weighted
radius of the precessing inner disk oscillates on a timescale of $\Delta
t_{osc}\sim 150-400r_{g}/c$ suggesting another sinusoidal modulation of the
lightcurve albeit at a higher frequency (Fig. 8 and Musoke et al. 2022). We
summarize the timescales of these variability phenomena in Table 1 for a BHXRB
($M_{BH}=10M_{\odot}$), a low-mass AGN ($M_{BH}=10^{6}M_{\odot}$), and a high-
mass AGN ($M_{BH}=10^{9}M_{\odot}$). Please note that these timescales might
vary by at least an order of magnitude depending the location of the tearing
radius, which we expect will depend on e.g. the accretion rate, black hole
spin, misalignment angle, and amount of large scale magnetic flux.

Figure 8: The mass weighted radius,
$r_{M}=\frac{\int_{0}^{2\pi}\int_{0}^{\pi}\Sigma r^{3}\sqrt{-g}d\theta
d\phi}{\int_{0}^{2\pi}\int_{0}^{\pi}\Sigma r^{2}\sqrt{-g}d\theta d\phi}$, of
the inner disk (after tearing) in RADT65 (blue) features a radial epicyclic
oscillation consistent with the radial epicyclic frequency of an oscillating
ring of gas at $r=6.9r_{g}$ (orange). Such oscillations may be associated with
HFQPOs detected in BHXRB lightcurves.

Table 1 suggests that the timescales associated with disk tearing are
consistent with variety of astrophysical phenomena. The tearing events
themselves might be associated with the heartbeat modes (on $1-100$ s
timescales) detected in the BHXRBs GRS 1915+105 and IGR J17091–3624 (e.g.
Belloni et al. 2000; Neilsen et al. 2011; Weng et al. 2018; Altamirano et al.
2011) in addition to QPEs detected in various AGN (e.g. Miniutti et al. 2019;
Giustini et al. 2020; Arcodia et al. 2021; Miniutti et al. 2022). If the
tearing events in some spectral states are spaced further apart, they could
potentially also explain some changing-look phenomena observed in AGN with
masses of $10^{7}-10^{9}M_{\odot}$ that occur on timescales of years to
decades (e.g. Sniegowska et al. 2020). In all of these cases the luminosity
fluctuates by a factor $10^{0}-10^{2}$ suggesting some significant flaring
event in the accretion disk reminiscent of a tearing event. The precession of
the disk that follows each tearing event can potentially explain low frequency
QPOs detected in numerous BHXRB with $\nu\sim 10^{-1}-10^{1}$ Hz (e.g. Ingram
& Motta 2019). Precession is also consistent with the $44$ day LFQPO detected
in the $2\times 10^{8}M_{\odot}$ AGN KIC 9650712 (Smith et al., 2018) and the
$3.8H$ LFQPO detected in the $10^{5}-10^{6}M_{\odot}$ AGN GJ1231+1106 (Lin et
al., 2014). In addition, the oscillation of the mass-weighted radius might be
associated with HFQPOs of $\nu\sim 10^{1}-10^{2}$ Hz in BHXRBs and some of the
short duration QPOs detected in various AGN. These include the $\sim 1H$ QPO
in the lightcurve of the $10^{6}-10^{7}M_{\odot}$ AGN RE J1034+396 (Gierlinski
et al., 2008; Czerny et al., 2016), the $1.8H$ QPO in the $5\times
10^{6}M_{\odot}$ MRK 766 (Zhang et al., 2017), and the $2H$ QPO in the
$4\times 10^{6}M_{\odot}$ AGN MS 2254.9-3712 (Alston et al., 2015).

The heated gas at the tearing radius suggests the QPO emission will be
substantially hardened through Comptonization. This is consistent with
frequency-resolved spectroscopy of BHXRBs which has demonstrated that the
Type-C QPO waveform is dominated by hardened emission (Sobolewska & Życki,
2006) and recent observations of AGN that suggest an increase in the soft-X-
ray excess during QPEs (Giustini et al., 2020). Interestingly, our simulations
suggest that tearing events are accompanied by a low-frequency sinusoidal
signal caused by precession of the inner disk and a high-frequency sinusoidial
signal caused by a radial epicyclic oscillation at the tearing radius. This
can potentially be tested using observational data for both BHXRBs and AGN.

Table 1: Time period between large tearing events, $\Delta t_{tear}$, disk precession period, $\Delta t_{prec}$, and period of radial epicyclic oscillations, $\Delta t_{osc}$, for various black hole masses. These timescales can be associated with various variability phenomena in BHXRBs and AGN such as QPEs, QPOs, and heartbeat oscillations (Sec. 5).  $M_{BH}$ | $\Delta t_{tear}$ | $\Delta t_{prec}$ | $\Delta t_{osc}$
---|---|---|---
N/A | $50,000r_{g}/c$ | $6000r_{g}/c$ | $400r_{g}/c$
$10M_{\odot}$ | $2.5s$ | $0.3s$ | $20ms$
$10^{6}M_{\odot}$ | $70h$ | $9h$ | $0.5h$
$10^{9}M_{\odot}$ | $8yr$ | $1yr$ | $20d$

###  5.3 Comparison to Semi-analytical Models

Semi-analytical models typically divide warped accretion disks into two
categories (e.g. Papaloizou & Pringle 1983) depending on the scale height of
the disk relative to the $\alpha-$viscosity parameter. In the first category
bending waves excited by oscillating radial and vertical pressure gradients in
the warp are damped and the warp is propagated through viscous diffusion only.
This happens when the scaleheight is smaller than the viscosity parameter,
e.g. $h/r<\alpha$. In this regime instability criteria were derived that
predict the tearing of accretion disks (e.g. Ogilvie & Dubus 2001; Doǧan et
al. 2018; Doğan & Nixon 2020). In the second category the warp is propagated
through bending waves that travel at a fraction of the sound speed (Papaloizou
& Lin, 1995). These bending waves are excited by pressure gradients in a warp.
Recent work has made substantial progress in describing warp propagation in
the bending wave regime (e.g. Ogilvie & Latter 2013; Paardekooper & Ogilvie
2018; Deng et al. 2020; Ogilvie 2022), but has not yet come up with criteria
for disk tearing.

Interestingly, Deng & Ogilvie (2022) predicted using the affine model of
Ogilvie (2018) that oscillations in the vertical pressure gradient lead to an
azimuthal $m=2$ oscillation in the scale height of the disk similar to what is
observed in our work. These authors also demonstrated that a parametric
instability can induce turbulence and create ring-like structures in warped
accretion disks. Since $h/r\sim\alpha_{m}$ (Fig. 7d) and the warp is highly
non-linear it is unclear in which regime of warp propagation our disk falls
(see also Sorathia et al. 2013; Hawley & Krolik 2018). We do note though that
there is some evidence for wave-like structures in the disk (Fig. 10), which
suggests the disk might be in the bending wave regime of warp propagation.
Nevertheless, future work will need to isolate wave-like from diffusive
(misaligned) angular momentum transport to make more direct comparison to
semi-analytical models possible.

Figure 9: A scatterplot of the effective viscosity, $\alpha_{\rm eff}$, with
respect to the thermal scaleheight, $\theta_{t}$ at $r=5r_{g}$ in green,
$r=10r_{g}$ in black, and $r=20r_{g}$ in blue. We show power-law fits to the
data as solid lines with fitting function $\alpha_{\rm
eff}\propto(h/r)^{-\gamma}$. The effective viscosity increases as the disk
becomes thinner, which suggests warp-driven dissipation becomes more important
for thinner disks. We argue in Section 5.4 that warp-driven dissipation can
stabilize the disk against runaway thermal collapse (e.g. Shakura & Sunyaev
1976).

###  5.4 Thermal and Viscous Stability of Warped Disks

In Section 4.5 we demonstrated that warped, radiation pressure supported,
accretion disks remain thermally stable for much longer than similar non-
warped, radiation pressure supported, accretion disks considered in previous
work (Liska et al., 2022). We concluded that neither magnetic pressure nor
advection of radiation internal energy can explain this thermal stability.
Here we propose an alternative explanation.

For a disk to remain thermally stable the derivative of the dissipation rate,
$Q_{vis}$, with respect to the temperature needs to be smaller than the
derivative of the cooling rate, $Q_{rad}$, with respect to the temperature
such that $\frac{Q_{vis}}{dT}<\frac{Q_{rad}}{dT}$. For a radiation pressure
supported $\alpha-$disk $Q_{vis}\propto T^{r}_{\phi}H\propto\alpha
p_{tot}H\propto\alpha p_{tot}^{2}\propto\alpha T^{8}$ and $Q_{rad}\propto
T^{4}$ (Shakura & Sunyaev, 1973; Meier, 2012). This suggests that with a
constant $\alpha$-viscosity radiation pressure supported accretion disks are
thermally unstable. However, if the dissipation rate follows $Q_{vis}\propto
T^{\zeta}$ with $\zeta\lesssim 4$ due to warping, the disk will remain
thermally stable.

To determine $\zeta$ in our simulation we present in Figure 9 a scatter plot
of the effective viscosity, $\alpha_{\rm eff}$, with respect to the thermal
scale height, $\theta_{t}$. We only use data after $32,000r_{g}/c$ to allow
the disk to adjust after seeding it with radiation and exclude data where the
tilt angle is smaller than $\mathcal{T}<60^{\circ}$. We visually fit the data
to find that $\alpha_{\rm eff}\propto H^{-\gamma}$ with $\gamma\sim
1.25-1.75$. From the relation $H\propto p_{tot}\propto T^{4}$ we find that
$\zeta=8-4\gamma\sim 1-3$. This suggests that even though our accretion disk
is radiation pressure dominated at $r=20r_{g}$, warp-driven dissipation can
potentially prevent the disk from becoming thermally unstable.

Following a similar argument, one can show that the warp-driven dissipation
also avoids the viscous Lightman & Eardley (1974) instability of radiation
pressure dominated disks. If a disk is subject to the viscous inability, it
will disintegrate into rings of high density gas (e.g. Mishra et al. 2016).
Even though high-density rings do form in our simulation (Figs. 2 and 4), we
conclude that this cannot be a manifestation of the viscous instability since
similar rings form in the non-radiative analogue to RADT65 presented by Musoke
et al. (2022). This suggests that the ring formation in RADT65 is warp-
induced.

Figure 10: A transverse slice through the density along the line of nodes at
$50r_{g}$. This snapshot is taken beyond the ‘official’ duration of our
simulation (Sec. 3), just after the onset of the next tearing event. Outwards
propagating wave-like structures (visible as vertical ‘stripes’) form in the
outer disk suggesting that wave-like angular momentum transport might be
present. The vertical extent of the outer disk is artificially inflated
because the frame is sliced at an angle with respect to the disk. Figure 11:
Radiation field warps the inner disk, as seen in a transverse slice of density
through the $x-z$ plane of our initial conditions (left) and after $\Delta
t=2,000r_{g}/c$ (right) with radiation streamlines shown in red. Even though
our M1 radiation scheme is unable to properly model multi-beam radiation (Sec.
5.5), these results suggest that radiation warping might be an important
effect.

###  5.5 Radiation Warping

When the inner disk precesses it will periodically beam radiation towards the
outer disk. This can potentially warp the outer disk and induce additional
dissipation (similar to Pringle 1997; Wijers & Pringle 1999). In a similar
context Liska et al. (2019a) has demonstrated that when a precessing jet
collides with the outer accretion disk it can inject energy and angular
momentum into the disk, leading to an increase of particle orbits in the outer
disk. Since the energy efficiency of the radiative outflows in RADT65 is non-
negligible ($|\eta_{rad}|\sim 14.7\%$) we expect that radiation feedback can
have a profound effect on the structure and dynamics of an accretion disk.

Radiative feedback is an interesting avenue of future research that is hard to
address in this article due to limitations of our M1 radiation closure
(Levermore, 1984). M1 works well in modeling radiative cooling and energy
transport in optically thick regions, and in optically thin regions for a
single radiation bundle. However, M1 is unable to properly model crossing
radiation bundles. The reason for this is that M1 treats radiation as a highly
collisional fluid and is thus unable to model crossing radiation bundles in
optically thin media. Instead the energy and momentum of the two radiation
bundles will be summed together. While the M1 closure conserves the total
amount of energy and momentum in the photon-fluid, it can lead to a
nonphysical redistribution of energy and angular momentum.

Despite these limitations we demonstrate in Figure 11 that the structure of
the inner accretion disk in RADT65 differs significantly after seeding it with
radiation. In Figure 11 the radiation streamlines get deflected by the inner
disk presumably leading to an increase in the warp amplitude compared to the
start of the simulation. In this snapshot the tearing process has not started
yet, which suggests that radiation warping can be important even in the
absence of disk tearing. While these preliminary results are encouraging, more
advanced numerical simulations will be necessary to reliably address radiation
warping.

##  6 Conclusion

In this article we have presented the first radiative two-temperature GRMHD
simulation of an $65^{\circ}$ misaligned disk accreting at $\dot{M}\sim
0.35\dot{M}_{edd}$. Similar to the idealized GRMHD models presented in Liska
et al. (2019a, b); Musoke et al. (2022) the radiation pressure supported
accretion disk tears apart and forms a precessing disk. During the tearing
process the mass accretion rate increases by a factor $\sim 10-40$ and the
luminosity increases by a factor $\sim 5$, almost exceeding the Eddington
limit at the peak of the outburst. The timescales and amplitudes of these
luminosity swings are roughly consistent with QPEs in AGN (e.g. Miniutti et
al. 2019; Giustini et al. 2020; Arcodia et al. 2021; Miniutti et al. 2022) and
several heartbeat modes in BHXRBs (e.g. Belloni et al. 2000; Neilsen et al.
2011; Altamirano et al. 2011; Weng et al. 2018). Following a tearing event,
the inner disk precesses for several periods with a frequency that is
consistent with low-frequency QPOs (e.g. Ingram & Motta 2019). While
precessing, the inner disk also exhibits a radial epicyclic oscillation of its
center of mass whose frequency is consistent with high frequency QPOs (e.g.
Ingram & Motta 2019). Future long-duration GRMHD simulations combined with
radiative transfer calculations will need to test if disk tearing indeed
produces QPO-like variability. This will require the tearing radius to remain
stable over many cycles of tearing.

Warping of the disk forms two nozzle shocks directed perpendicular to the line
of nodes. At the nozzle shock, the scale height of the disk decreases by a
factor of $\sim 5$ and the temperature increases up to $T_{e}\sim
10^{8}-10^{9}K$ and $T_{i}\sim 10^{9}-10^{10}K$. In a companion paper (Kaaz et
al 2022, submitted), we demonstrate that dissipation in the nozzle shock leads
to a much shorter accretion time scale than predicted by $\alpha$-viscosity
based models (e.g. Shakura & Sunyaev 1973) seeded by magnetized turbulence
(e.g. Balbus & Hawley 1991). This can potentially explain the timescales
associated with some changing look phenomena in AGN including QPEs. In
addition, we find during tearing events that gas crossing the tear heats up to
$T_{e}\sim 10^{8}K$ and $T_{i}\sim 10^{10}K$. This is caused by an increase in
the dissipation rate as the gas undergoes a rapid orbital plane change. These
results imply that the emission profile from warped accretion disks will
deviate substantially from the Novikov & Thorne (1973) model of a
geometrically thin accretion disk. This can, pending future radiative transfer
calculations, have far-reaching consequences for spectral fitting and black
hole spin measurements (e.g. Zhang et al. 1997; Kulkarni et al. 2011).

While radiation pressure supported accretion disks in GRMHD simulations
typically collapse (e.g. Sa̧dowski 2016; Fragile et al. 2018; Jiang et al.
2019; Mishra et al. 2020; Liska et al. 2022) into an infinitely thin slab due
to the thermal instability (Shakura & Sunyaev, 1976), this did not happen in
our simulation. Here we argue that the warp and associated nozzle shock
stabilize the disk against thermal collapse (Sec. 5). These shocks differ from
spiral density wave induced shocks (e.g. Arzamasskiy & Rafikov 2018). In a
companion paper (Kaaz et al 2022, submitted), we demonstrate that nozzle shock
dissipation drives accretion in geometrically thin warped accretion disks
within at least $r\lesssim 50r_{g}$. This upends the standard paradigm that
invokes magneto-rotational instability induced turbulence for the dissipation
of orbital kinetic energy (e.g. Balbus & Hawley 1991, 1998). Future work will
need to address the range of spectral states, accretion rates, and
misalignment angles where nozzle shock dissipation becomes important.

Past work almost exclusively invoked magnetic reconnection in coronae to
explain high energy (non-)thermal emission from accretion disks (e.g. Sironi &
Spitkovsky 2014; Beloborodov 2017; Ripperda et al. 2021). However, this work
suggests that heated plasma formed in the nozzle and at the tearing radius can
potentially be another source of high energy emission. This is especially
appealing to explain non-thermal emission in the high-soft and soft-
intermediate states of BHXRBs. Those states do not launch any large-scale jets
and are thus unlikely to be saturated by magnetic flux (e.g. McKinney et al.
2012), which might be a key ingredient to form a corona where the prime source
of dissipation comes from magnetic reconnection (e.g. Liska et al. 2022;
Ripperda et al. 2021). Test-particle methods (e.g. Bacchini et al. 2019)
informed by first principles particle-in-cell simulations are an interesting
avenue for future research to quantify the contribution of magnetic
reconnection and warp-driven shocks to the acceleration of non-thermal
particles and corona-like emission.

##  7 Acknowledgements

We thank Michiel van der Klis, Adam Ingram, Sera Markoff and Ramesh Narayan
for insightful discussions. An award of computer time was provided by the
Innovative and Novel Computational Impact on Theory and Experiment (INCITE),
OLCF Director’s Discretionary Allocation, and ASCR Leadership Computing
Challenge (ALCC) programs under award PHY129. This research used resources of
the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science
User Facility supported under Contract DE-AC05-00OR22725.This research used
resources of the National Energy Research Scientific Computing Center, a DOE
Office of Science User Facility supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC02-05CH11231 using NERSC award
ALCC-ERCAP0022634. We acknowledge PRACE for awarding us access to JUWELS
Booster at GCS@JSC, Germany. ML was supported by the John Harvard
Distinguished Science Fellowship, GM is supported by a Netherlands Research
School for Astronomy (NOVA), Virtual Institute of Accretion (VIA) postdoctoral
fellowship, NK by an NSF Graduate Research Fellowship, and AT by the National
Science Foundation grants AST-2206471, AST-2009884, AST-2107839, AST-1815304,
OAC-2031997, and AST-1911080.

## References

  * Abramowicz et al. (1978) Abramowicz, M., Jaroszynski, M., & Sikora, M. 1978, A&A, 63, 221
  * Abramowicz & Kluźniak (2004) Abramowicz, M. A., & Kluźniak, W. 2004, in American Institute of Physics Conference Series, Vol. 714, X-ray Timing 2003: Rossi and Beyond, ed. P. Kaaret, F. K. Lamb, & J. H. Swank, 21–28, doi: 10.1063/1.1780993
  * Alston et al. (2015) Alston, W. N., Parker, M. L., Markevičiūtė, J., et al. 2015, MNRAS, 449, 467, doi: 10.1093/mnras/stv351
  * Altamirano et al. (2011) Altamirano, D., Belloni, T., Linares, M., et al. 2011, ApJ, 742, L17, doi: 10.1088/2041-8205/742/2/L17
  * Andalman et al. (2022) Andalman, Z. L., Liska, M. T. P., Tchekhovskoy, A., Coughlin, E. R., & Stone, N. 2022, MNRAS, 510, 1627, doi: 10.1093/mnras/stab3444
  * Arcodia et al. (2021) Arcodia, R., Merloni, A., Nandra, K., et al. 2021, Nature, 592, 704, doi: 10.1038/s41586-021-03394-6
  * Arzamasskiy & Rafikov (2018) Arzamasskiy, L., & Rafikov, R. R. 2018, ApJ, 854, 84, doi: 10.3847/1538-4357/aaa8e8
  * Bacchini et al. (2019) Bacchini, F., Ripperda, B., Porth, O., & Sironi, L. 2019, ApJS, 240, 40, doi: 10.3847/1538-4365/aafcb3
  * 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
  * Bardeen & Petterson (1975) Bardeen, J. M., & Petterson, J. A. 1975, ApJ, 195, L65, doi: 10.1086/181711
  * Begelman & Pringle (2007) Begelman, M. C., & Pringle, J. E. 2007, MNRAS, 375, 1070, doi: 10.1111/j.1365-2966.2006.11372.x
  * Belloni et al. (2000) Belloni, T., Klein-Wolt, M., Méndez, M., van der Klis, M., & van Paradijs, J. 2000, A&A, 355, 271.  https://arxiv.org/abs/astro-ph/0001103
  * Beloborodov (2017) Beloborodov, A. M. 2017, ApJ, 850, 141, doi: 10.3847/1538-4357/aa8f4f
  * Chatterjee et al. (2019) Chatterjee, K., Liska, M., Tchekhovskoy, A., & Markoff, S. B. 2019, MNRAS, 490, 2200, doi: 10.1093/mnras/stz2626
  * Chatterjee et al. (2020) Chatterjee, K., Younsi, Z., Liska, M., et al. 2020, MNRAS, 499, 362, doi: 10.1093/mnras/staa2718
  * Clavel et al. (1992) Clavel, J., Nandra, K., Makino, F., et al. 1992, ApJ, 393, 113, doi: 10.1086/171490
  * Courant et al. (1928) Courant, R., Friedrichs, K., & Lewy, H. 1928, Mathematische Annalen, 100, 32, doi: 10.1007/BF01448839
  * Czerny et al. (2016) Czerny, B., You, B., Kurcz, A., et al. 2016, A&A, 594, A102, doi: 10.1051/0004-6361/201628103
  * Deng & Ogilvie (2022) Deng, H., & Ogilvie, G. I. 2022, Monthly Notices of the Royal Astronomical Society, 512, 6078, doi: 10.1093/mnras/stac858
  * Deng et al. (2020) Deng, H., Ogilvie, G. I., & Mayer, L. 2020, Monthly Notices of the Royal Astronomical Society, 500, 4248, doi: 10.1093/mnras/staa3504
  * Dewberry et al. (2020a) Dewberry, J. W., Latter, H. N., Ogilvie, G. I., & Fromang, S. 2020a, MNRAS, 497, 435, doi: 10.1093/mnras/staa1897
  * Dewberry et al. (2020b) —. 2020b, MNRAS, 497, 451, doi: 10.1093/mnras/staa1898
  * Dexter & Begelman (2019) Dexter, J., & Begelman, M. C. 2019, MNRAS, 483, L17, doi: 10.1093/mnrasl/sly213
  * Doğan & Nixon (2020) Doğan, S., & Nixon, C. J. 2020, MNRAS, 495, 1148, doi: 10.1093/mnras/staa1239
  * Doǧan et al. (2018) Doǧan, S., Nixon, C. J., King, A. R., & Pringle, J. E. 2018, MNRAS, 476, 1519, doi: 10.1093/mnras/sty155
  * Esin et al. (1997) Esin, A. A., McClintock, J. E., & Narayan, R. 1997, ApJ, 489, 865, doi: 10.1086/304829
  * Fender et al. (2004) Fender, R. P., Belloni, T. M., & Gallo, E. 2004, MNRAS, 355, 1105, doi: 10.1111/j.1365-2966.2004.08384.x
  * Fragile & Anninos (2005) Fragile, P. C., & Anninos, P. 2005, ApJ, 623, 347, doi: 10.1086/428433
  * Fragile & Blaes (2008) Fragile, P. C., & Blaes, O. M. 2008, ApJ, 687, 757, doi: 10.1086/591936
  * Fragile et al. (2007) Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007, ApJ, 668, 417, doi: 10.1086/521092
  * Fragile et al. (2018) Fragile, P. C., Etheridge, S. M., Anninos, P., Mishra, B., & Kluźniak, W. 2018, ApJ, 857, 1, doi: 10.3847/1538-4357/aab788
  * Gardiner & Stone (2005) Gardiner, T. A., & Stone, J. M. 2005, Journal of Computational Physics, 205, 509, doi: 10.1016/j.jcp.2004.11.016
  * Gierliński & Done (2004) Gierliński, M., & Done, C. 2004, MNRAS, 349, L7, doi: 10.1111/j.1365-2966.2004.07687.x
  * Gierliński & Done (2006) —. 2006, MNRAS, 371, L16, doi: 10.1111/j.1745-3933.2006.00199.x
  * Gierlinski et al. (2008) Gierlinski, M., Middleton, M., Ward, M., & Done, C. 2008, arXiv e-prints, arXiv:0807.1899.  https://arxiv.org/abs/0807.1899
  * Giustini et al. (2020) Giustini, M., Miniutti, G., & Saxton, R. D. 2020, A&A, 636, L2, doi: 10.1051/0004-6361/202037610
  * Hankla et al. (2022) Hankla, A. M., Scepi, N., & Dexter, J. 2022, MNRAS, 515, 775, doi: 10.1093/mnras/stac1785
  * Hawley & Krolik (2018) Hawley, J. F., & Krolik, J. H. 2018, ApJ, 866, 5, doi: 10.3847/1538-4357/aadf90
  * Ingram et al. (2009) Ingram, A., Done, C., & Fragile, P. C. 2009, MNRAS, 397, L101, doi: 10.1111/j.1745-3933.2009.00693.x
  * Ingram & van der Klis (2014) Ingram, A., & van der Klis, M. 2014, Monthly Notices of the Royal Astronomical Society, 446, 3516, doi: 10.1093/mnras/stu2373
  * Ingram et al. (2016) Ingram, A., van der Klis, M., Middleton, M., et al. 2016, MNRAS, 461, 1967, doi: 10.1093/mnras/stw1245
  * Ingram & Motta (2019) Ingram, A. R., & Motta, S. E. 2019, New A Rev., 85, 101524, doi: 10.1016/j.newar.2020.101524
  * Jiang et al. (2019) Jiang, Y.-F., Blaes, O., Stone, J. M., & Davis, S. W. 2019, ApJ, 885, 144, doi: 10.3847/1538-4357/ab4a00
  * Kara et al. (2019) Kara, E., Steiner, J. F., Fabian, A. C., et al. 2019, Nature, 565, 198, doi: 10.1038/s41586-018-0803-x
  * Kato (2004) Kato, S. 2004, PASJ, 56, 905, doi: 10.1093/pasj/56.5.905
  * King (2020) King, A. 2020, MNRAS, 493, L120, doi: 10.1093/mnrasl/slaa020
  * Kitaki et al. (2021) Kitaki, T., Mineshige, S., Ohsuga, K., & Kawashima, T. 2021, arXiv e-prints, arXiv:2101.11028.  https://arxiv.org/abs/2101.11028
  * Kluzniak & Abramowicz (2002) Kluzniak, W., & Abramowicz, M. A. 2002, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0203314
  * Kochanek (1994) Kochanek, C. S. 1994, ApJ, 422, 508, doi: 10.1086/173745
  * Kulkarni et al. (2011) Kulkarni, A. K., Penna, R. F., Shcherbakov, R. V., et al. 2011, MNRAS, 414, 1183, doi: 10.1111/j.1365-2966.2011.18446.x
  * Lawrence (2018) Lawrence, A. 2018, Nature Astronomy, 2, 102, doi: 10.1038/s41550-017-0372-1
  * Lense & Thirring (1918) Lense, J., & Thirring, H. 1918, Physikalische Zeitschrift, 19
  * Levermore (1984) Levermore, C. 1984, J. Quant. Spec. Radiat. Transf., 31, 149, doi: 10.1016/0022-4073(84)90112-2
  * Lightman & Eardley (1974) Lightman, A. P., & Eardley, D. M. 1974, ApJ, 187, L1, doi: 10.1086/181377
  * Lightman & Eardley (1974) Lightman, A. P., & Eardley, D. M. 1974, The Astrophysical Journal, 187, L1, doi: 10.1086/181377
  * Lin et al. (2014) Lin, D., Webb, N. A., & Barret, D. 2014, ApJ, 780, 39, doi: 10.1088/0004-637X/780/1/39
  * Liska et al. (2018a) Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018a, MNRAS, 474, L81, doi: 10.1093/mnrasl/slx174
  * Liska et al. (2019a) —. 2019a, arXiv e-prints, arXiv:1904.08428.  https://arxiv.org/abs/1904.08428
  * Liska et al. (2018b) Liska, M., Tchekhovskoy, A., Ingram, A., & van der Klis, M. 2018b, ArXiv e-prints.  https://arxiv.org/abs/1810.00883
  * Liska et al. (2019b) Liska, M., Chatterjee, K., Tchekhovskoy, A., et al. 2019b, arXiv e-prints, arXiv:1912.10192.  https://arxiv.org/abs/1912.10192
  * Liska et al. (2022) Liska, M. T. P., Musoke, G., Tchekhovskoy, A., Porth, O., & Beloborodov, A. M. 2022, arXiv e-prints, arXiv:2201.03526.  https://arxiv.org/abs/2201.03526
  * MacLeod et al. (2016) MacLeod, C. L., Ross, N. P., Lawrence, A., et al. 2016, MNRAS, 457, 389, doi: 10.1093/mnras/stv2997
  * McKinney et al. (2017) McKinney, J. C., Chluba, J., Wielgus, M., Narayan, R., & Sadowski, A. 2017, MNRAS, 467, 2241, doi: 10.1093/mnras/stx227
  * McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083, doi: 10.1111/j.1365-2966.2012.21074.x
  * McKinney et al. (2013) —. 2013, Science, 339, 49, doi: 10.1126/science.1230811
  * Meier (2012) Meier, D. L. 2012, Black Hole Astrophysics: The Engine Paradigm
  * Miniutti et al. (2022) Miniutti, G., Giustini, M., Arcodia, R., et al. 2022, arXiv e-prints, arXiv:2207.07511.  https://arxiv.org/abs/2207.07511
  * Miniutti et al. (2019) Miniutti, G., Saxton, R. D., Giustini, M., et al. 2019, Nature, 573, 381, doi: 10.1038/s41586-019-1556-x
  * Mishra et al. (2022) Mishra, B., Fragile, P. C., Anderson, J., et al. 2022, arXiv e-prints, arXiv:2209.03317.  https://arxiv.org/abs/2209.03317
  * Mishra et al. (2016) Mishra, B., Fragile, P. C., Johnson, L. C., & Kluźniak, W. 2016, MNRAS, 463, 3437, doi: 10.1093/mnras/stw2245
  * Mishra et al. (2020) Mishra, B., Kluźniak, W., & Fragile, P. C. 2020, MNRAS, 497, 1066, doi: 10.1093/mnras/staa1848
  * Musoke et al. (2022) Musoke, G., Liska, M., Porth, O., van der Klis, M., & Ingram, A. 2022, arXiv e-prints, arXiv:2201.03085.  https://arxiv.org/abs/2201.03085
  * Narayan et al. (2003) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69, doi: 10.1093/pasj/55.6.L69
  * Narayan & Yi (1994) Narayan, R., & Yi, I. 1994, ApJ, 428, L13, doi: 10.1086/187381
  * Nealon et al. (2015) Nealon, R., Price, D. J., & Nixon, C. J. 2015, MNRAS, 448, 1526, doi: 10.1093/mnras/stv014
  * Neilsen et al. (2011) Neilsen, J., Remillard, R. A., & Lee, J. C. 2011, ApJ, 737, 69, doi: 10.1088/0004-637X/737/2/69
  * Nelson & Papaloizou (2000) Nelson, R. P., & Papaloizou, J. C. B. 2000, MNRAS, 315, 570, doi: 10.1046/j.1365-8711.2000.03478.x
  * Newman & Hamlin (2014) Newman, W. I., & Hamlin, N. D. 2014, SIAM Journal on Scientific Computing, 36, B661, doi: 10.1137/140956749
  * Nixon et al. (2012a) Nixon, C., King, A., Price, D., & Frank, J. 2012a, ApJ, 757, L24, doi: 10.1088/2041-8205/757/2/L24
  * Nixon & Salvesen (2013) Nixon, C., & Salvesen, G. 2013, Monthly Notices of the Royal Astronomical Society, 437, 3994, doi: 10.1093/mnras/stt2215
  * Nixon & King (2012) Nixon, C. J., & King, A. R. 2012, MNRAS, 421, 1201, doi: 10.1111/j.1365-2966.2011.20377.x
  * Nixon et al. (2012b) Nixon, C. J., King, A. R., & Price, D. J. 2012b, MNRAS, 422, 2547, doi: 10.1111/j.1365-2966.2012.20814.x
  * Noble et al. (2006) Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626, doi: 10.1086/500349
  * Noble et al. (2009) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, ApJ, 692, 411, doi: 10.1088/0004-637X/692/1/411
  * Novikov & Thorne (1973) Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus), ed. C. Dewitt & B. S. Dewitt, 343–450
  * Ogilvie (1999) Ogilvie, G. I. 1999, MNRAS, 304, 557, doi: 10.1046/j.1365-8711.1999.02340.x
  * Ogilvie (2018) —. 2018, MNRAS, 477, 1744, doi: 10.1093/mnras/sty588
  * Ogilvie (2022) —. 2022, MNRAS, 513, 1701, doi: 10.1093/mnras/stac939
  * Ogilvie & Dubus (2001) Ogilvie, G. I., & Dubus, G. 2001, Monthly Notices of the Royal Astronomical Society, 320, 485, doi: 10.1046/j.1365-8711.2001.04011.x
  * Ogilvie & Latter (2013) Ogilvie, G. I., & Latter, H. N. 2013, MNRAS, 433, 2420, doi: 10.1093/mnras/stt917
  * Paardekooper & Ogilvie (2018) Paardekooper, S.-J., & Ogilvie, G. I. 2018, Monthly Notices of the Royal Astronomical Society, 483, 3738, doi: 10.1093/mnras/sty3349
  * Papaloizou & Lin (1995) Papaloizou, J. C. B., & Lin, D. N. C. 1995, ARA&A, 33, 505, doi: 10.1146/annurev.aa.33.090195.002445
  * Papaloizou & Pringle (1983) Papaloizou, J. C. B., & Pringle, J. E. 1983, MNRAS, 202, 1181, doi: 10.1093/mnras/202.4.1181
  * Piran et al. (2015) Piran, T., Sa̧dowski, A., & Tchekhovskoy, A. 2015, MNRAS, 453, 157, doi: 10.1093/mnras/stv1547
  * Porth et al. (2019) Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26, doi: 10.3847/1538-4365/ab29fd
  * Pringle (1997) Pringle, J. E. 1997, Monthly Notices of the Royal Astronomical Society, 292, 136, doi: 10.1093/mnras/292.1.136
  * Raj & Nixon (2021) Raj, A., & Nixon, C. J. 2021, ApJ, 909, 82, doi: 10.3847/1538-4357/abdc25
  * Raj et al. (2021) Raj, A., Nixon, C. J., & Doğan, S. 2021, ApJ, 909, 81, doi: 10.3847/1538-4357/abdc24
  * Rees (1988) Rees, M. J. 1988, Nature, 333, 523, doi: 10.1038/333523a0
  * Remillard & McClintock (2006) Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49, doi: 10.1146/annurev.astro.44.051905.092532PDF: http://arjournals.annualreviews.org/doi/pdf/10.1146/annurev.astro.44.051905.092532
  * Ressler et al. (2015) Ressler, S. M., Tchekhovskoy, A., Quataert, E., Chandra, M., & Gammie, C. F. 2015, MNRAS, 454, 1848, doi: 10.1093/mnras/stv2084
  * Ressler et al. (2017) Ressler, S. M., Tchekhovskoy, A., Quataert, E., & Gammie, C. F. 2017, MNRAS, 467, 3604, doi: 10.1093/mnras/stx364
  * Rezzolla et al. (2003) Rezzolla, L., Yoshida, S., & Zanotti, O. 2003, MNRAS, 344, 978, doi: 10.1046/j.1365-8711.2003.07023.x
  * Ripperda et al. (2021) Ripperda, B., Liska, M., Chatterjee, K., et al. 2021, arXiv:2109.15115 [astro-ph, physics:gr-qc, physics:physics].  https://arxiv.org/abs/2109.15115
  * Rowan et al. (2017) Rowan, M. E., Sironi, L., & Narayan, R. 2017, ApJ, 850, 29, doi: 10.3847/1538-4357/aa9380
  * Russell et al. (2020) Russell, D. M., Casella, P., Kalemci, E., et al. 2020, MNRAS, 495, 182, doi: 10.1093/mnras/staa1182
  * Sadowski (2011) Sadowski, A. 2011, arXiv e-prints, arXiv:1108.0396.  https://arxiv.org/abs/1108.0396
  * Sa̧dowski (2016) Sa̧dowski, A. 2016, MNRAS, 459, 4397, doi: 10.1093/mnras/stw913
  * Sa̧dowski & Narayan (2015) Sa̧dowski, A., & Narayan, R. 2015, MNRAS, 454, 2372, doi: 10.1093/mnras/stv2022
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Shakura & Sunyaev (1976) —. 1976, MNRAS, 175, 613, doi: 10.1093/mnras/175.3.613
  * Sironi & Spitkovsky (2014) Sironi, L., & Spitkovsky, A. 2014, ApJ, 783, L21, doi: 10.1088/2041-8205/783/1/l21
  * Sadowski et al. (2013) Sadowski, A., Narayan, R., Tchekhovskoy, A., & Zhu, Y. 2013, MNRAS, 429, 3533, doi: 10.1093/mnras/sts632
  * Smith et al. (2018) Smith, K. L., Mushotzky, R. F., Boyd, P. T., & Wagoner, R. V. 2018, ApJ, 860, L10, doi: 10.3847/2041-8213/aac88c
  * Sniegowska et al. (2020) Sniegowska, M., Czerny, B., Bon, E., & Bon, N. 2020, A&A, 641, A167, doi: 10.1051/0004-6361/202038575
  * Sobolewska & Życki (2006) Sobolewska, M. A., & Życki, P. T. 2006, MNRAS, 370, 405, doi: 10.1111/j.1365-2966.2006.10489.x
  * Sorathia et al. (2013) Sorathia, K. A., Krolik, J. H., & Hawley, J. F. 2013, ApJ, 768, 133, doi: 10.1088/0004-637X/768/2/133
  * Stella & Vietri (1998) Stella, L., & Vietri, M. 1998, ApJ, 492, L59, doi: 10.1086/311075
  * Stepney (1983) Stepney, S. 1983, MNRAS, 202, 467, doi: 10.1093/mnras/202.2.467
  * Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79, doi: 10.1111/j.1745-3933.2011.01147.x
  * Uttley et al. (2014) Uttley, P., Cackett, E. M., Fabian, A. C., Kara, E., & Wilkins, D. R. 2014, 22, 72, doi: 10.1007/s00159-014-0072-0
  * Wagoner & Tandon (2021) Wagoner, R. V., & Tandon, C. R. 2021, ApJ, 920, 114, doi: 10.3847/1538-4357/ac1868
  * Weng et al. (2018) Weng, S.-S., Wang, T.-T., Cai, J.-P., Yuan, Q.-R., & Gu, W.-M. 2018, ApJ, 865, 19, doi: 10.3847/1538-4357/aad9a6
  * White et al. (2019) White, C. J., Quataert, E., & Blaes, O. 2019, arXiv e-prints, arXiv:1902.09662.  https://arxiv.org/abs/1902.09662
  * White et al. (2016) White, C. J., Stone, J. M., & Gammie, C. F. 2016, ApJS, 225, 22, doi: 10.3847/0067-0049/225/2/22
  * Wijers & Pringle (1999) Wijers, R. A. M. J., & Pringle, J. E. 1999, MNRAS, 308, 207, doi: 10.1046/j.1365-8711.1999.02720.x
  * Yang et al. (2018a) Yang, Q., Wu, X.-B., Fan, X., et al. 2018a, ApJ, 862, 109, doi: 10.3847/1538-4357/aaca3a
  * Yang et al. (2018b) —. 2018b, ApJ, 862, 109, doi: 10.3847/1538-4357/aaca3a
  * Zhang et al. (2017) Zhang, P., Zhang, P.-f., Yan, J.-z., Fan, Y.-z., & Liu, Q.-z. 2017, ApJ, 849, 9, doi: 10.3847/1538-4357/aa8d6e
  * Zhang et al. (1997) Zhang, S. N., Cui, W., & Chen, W. 1997, ApJ, 482, L155, doi: 10.1086/310705
  * Zhu et al. (2012) Zhu, Y., Davis, S. W., Narayan, R., et al. 2012, Monthly Notices of the Royal Astronomical Society, 424, 2504, doi: 10.1111/j.1365-2966.2012.21181.x