# Bardeen-Petterson Alignment, Jets and Magnetic Truncation in GRMHD
Simulations of Tilted Thin Accretion Discs

M. Liska1, A. Tchekhovskoy2, A. Ingram3, M. van der Klis1
1Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science
Park 904, 1098 XH Amsterdam, The Netherlands
2Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA),
Physics & Astronomy, Northwestern University, Evanston, IL 60208, USA
3 Department of Physics, Astrophysics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford, OX1 3RH, UK E-mail<EMAIL_ADDRESS>
(Accepted. Received; in original form)

###### Abstract

Prevalent around luminous accreting black holes, thin discs are challenging to
resolve in numerical simulations. When the disc and black hole angular
momentum vectors are misaligned, the challenge becomes extreme, requiring
adaptive meshes to follow the disc proper as it moves through the
computational grid. With our new high-performance general relativistic
magnetohydrodynamic (GRMHD) code H-AMR we have simulated the thinnest
accretion disc to date, of aspect ratio $H/R\approx 0.03\approx 1.7^{\circ}$ ,
around a rapidly spinning ($a=0.9375$) black hole, using a cooling function.
Initially tilted at $10^{\circ}$, the disc warps inside $\sim 5r_{g}$ into
alignment with the black hole, where $r_{g}$ is the gravitational radius. This
is the first demonstration of Bardeen-Petterson alignment in MHD with
viscosity self-consistently generated by magnetized turbulence. The disc
develops a low-density high-viscosity ($\alpha_{\rm eff}\sim 1.0$) magnetic-
pressure–dominated inner region at $r\lesssim 25r_{g}$ that rapidly empties
itself into the black hole. This inner region may in reality, due to thermal
decoupling of ions and electrons, evaporate into a radiatively inefficient
accretion flow if, as we propose, the cooling time exceeds the accretion time
set by the order unity effective viscosity. We furthermore find the unexpected
result that even our very thin disc can sustain large-scale vertical magnetic
flux on the black hole, which launches powerful relativistic jets that carry
$20{-}50\%$ of the accretion power along the angular momentum vector of the
outer tilted disc, providing a potential explanation for the origin of jets in
radio-loud quasars.

###### keywords:

accretion, accretion discs – black hole physics – MHD – galaxies: jets –
methods: numerical

††pagerange: Bardeen-Petterson Alignment, Jets and Magnetic Truncation in
GRMHD Simulations of Tilted Thin Accretion Discs–References††pubyear: 2015

##  1 Introduction

Black holes (BHs) in X-Ray binaries (XRB) and possibly active galactic nuclei
(AGN) cycle during their lifetimes through different accretion states
characterized by the total luminosity and spectral hardness. It is widely
believed that the accretion luminosity expressed as a fraction of the
Eddington limit $L_{\rm Edd}$ is an important factor in determining the BH’s
accretion state (Esin et al., 1997; Remillard & McClintock, 2006; McClintock &
Remillard, 2006). At low luminosities ($L\lesssim 0.01L_{\rm Edd}$), in the
low-hard state, a promising model is the advection dominated accretion flow
(ADAF, Ichimaru 1977; Narayan & Yi 1994, 1995a, 1995b). In an ADAF, the disc
surface density is so low that the plasma can decouple into a two-temperature
electron-ion plasma (Shapiro et al., 1976). Since the ions are unable to cool
on the accretion time, most of the dissipated energy is advected inwards in
the disc or expelled in outflows leading to a low radiative efficiency. It is
also known through general relativistic magnetohydrodynamic (GRMHD)
simulations that these thick accretion discs can sustain and advect inwards
large scale poloidal magnetic flux (e.g. De Villiers et al. 2003; McKinney
2006; Beckwith et al. 2008; McKinney & Blandford 2009; Tchekhovskoy et al.
2011; Tchekhovskoy & McKinney 2012; McKinney et al. 2012) which launches
powerful jets when it reaches the central BH (Blandford & Znajek, 1977). The
large scale height makes ADAFs numerically easy to study since they do not
require high resolutions and have short viscous times. They have been studied
extensively in GRMHD and are relatively well understood.

However, the observed emission from X-ray binaries in the high-soft state, and
from high-Eddington fraction AGN, is incompatible with the ADAF solution.
Their thermal emission spectrum requires a geometrically thin, optically thick
accretion disc (Shakura & Sunyaev, 1973; Novikov & Thorne, 1973). Though X-ray
binaries and AGN only spend a rather short amount of time in the thin disc
regime, most of the BH growth and feedback may still occur there since the
accretion rate is several orders of magnitude higher. Indeed, a typical bright
quasar (or XRB) radiates at $\sim 10\%$ of the Eddington rate, $\dot{M}_{\rm
Edd}=10L_{\rm Edd}/c^{2}$, and its disc thickness is extremely small,
$H/R\simeq 0.01\times(\dot{M}/0.1\dot{M}_{\rm Edd})^{0.9}$ (see e.g. Fig. 3 in
Piran et al., 2015). Thus, it is crucial to study these systems to understand
the growth and feedback of supermassive BHs.

Modelling of relativistic iron line profiles suggests that most bright local
AGN harbour rapidly spinning supermassive BHs, with spin parameter $a>0.8$
(Reynolds, 2014). Because the infalling material is unaware of the orientation
of the BH spin, its angular momentum vector is expected to be misaligned with
respect to the black hole spin vector, resulting in a tilted accretion disc.
In fact, there exist several such candidates for both X-ray binaries
(Hjellming & Rupen, 1995; Greene et al., 2001; Maccarone, 2002) and AGN
(Caproni et al., 2006, 2007). Analytic theory predicts that the inner parts of
a thin disc would align with the BH midplane due to the Bardeen and Petterson
effect (hereafter BP, Bardeen & Petterson 1975; Papaloizou & Pringle 1983;
Kumar & Pringle 1985; Pringle 1992; Ogilvie 1999) out to a Bardeen-Petterson
radius, $r_{\rm BP}$. The alignment can have profound consequences for the
growth and feedback of supermassive BHs (Rees, 1978; Scheuer & Feiler, 1996;
Natarajan & Pringle, 1998) since a larger $r_{\rm BP}$ implies a larger
alignment torque on the BH, forcing it into alignment with the outer disc on
very short timescales. Aligned accretion subsequently leads to rapid BH spinup
by the disc, even if most of the infalling material is initially misaligned
(and even counter-aligned in some cases, see King et al. 2005).

The BP effect is caused by the interplay between general relativistic Lense-
Thirring precession (Lense & Thirring, 1918) and viscosity. BP alignment has
been observed in pioneering smoothed particle hydrodynamics (SPH) simulations
(Nelson & Papaloizou, 2000; Lodato & Pringle, 2007; Lodato & Price, 2010).
However, because these simulations are non-relativistic and hydrodynamic, they
are unable to include the full effects of GR and treat the anisotropy of
magnetized turbulence accurately. In fact, GRMHD is a powerful way to model
the non-linear nature of the magnetorotational instability (MRI) driven
turbulence (Balbus & Hawley, 1991) responsible for the viscosity whilst
including the full effects of GR. However, GRMHD work at a moderate thickness
of $H/R=0.08$ did not find Bardeen-Petterson alignment (Morales Teixeira et
al., 2014; Zhuravlev et al., 2014). Moreover, there is growing evidence that
the interaction between the disc, magnetized corona and the jets should be
taken into account. Jets can torque the inner accretion disc into alignment
(McKinney et al., 2013) before they align with the outer accretion flow (Liska
et al., 2018). The corona, which we define as the hot bloated magnetic
pressure supported flow surrounding the thin disc, on the other hand, is not
expected to align since thick flows cannot exhibit Bardeen-Petterson alignment
(Ivanov & Illarionov, 1997; Papaloizou & Pringle, 1983). Thus, full GRMHD
simulations, which describe the entire _thin_ -disc–corona–jet system, are
uniquely positioned to address the more than $40$ year old fundamental problem
whether BP alignment occurs in the inner parts of thin tilted discs.

Indications are that jets in systems that contain a thin disc are rare: only
$10$% of quasars, or luminous AGN, are observed to produce relativistic jets
and the associated radio emission (e.g Sikora et al. 2007), while there are no
convincing observations in soft state X-ray binaries (though see Rushton et
al. 2012). It is crucial to understand what factors are responsible for the
formation and destruction of jets since they can be the dominant feedback mode
in galaxy clusters (see e.g. Fabian 2012). Early theoretical work suggested
that thin accretion discs are not expected to have jets since the poloidal
magnetic flux may diffuse out before it can advect inwards (Lubow et al.,
1994). However, the non-uniform vertical structure of accretion discs and
their turbulence may aid in the inward advection of poloidal magnetic flux
(Rothstein & Lovelace, 2008; Guilet & Ogilvie, 2012, 2013).

The small vertical extent of thin discs makes them very difficult to study
numerically, with the computational cost scaling as $(H/R)^{-5}$ per accretion
time. Because of this high cost, numerical work studying the physics of such
discs has been mostly limited to shearing box simulations and semi-analytical
studies. In fact, there are no 3D GRMHD simulations available for thin discs
of aspect ratios $H/R<0.05$, and the thinnest discs so far simulated in 3D
GRMHD ($0.05<H/R<0.1$) were aligned, which enabled the vertical wavelength of
the MRI to be resolved by a grid focused on the equatorial plane, leading to
cells compressed in the $\theta$-direction and elongated in the $r$\- and
$\phi$-directions (e.g. Shafee et al. 2008; Noble et al. 2009, 2010; Penna et
al. 2010; Morales Teixeira et al. 2014; Avara et al. 2016; Morales Teixeira et
al. 2017). Studying tilted discs, whose orbital motion does not conform to the
main directions of the grid, is more difficult than aligned ones because one
can assume neither axisymmetry nor use elongated cells to speed up the
simulations.

In this work we present the thinnest global GRMHD accretion disc simulations
to date and study Bardeen-Petterson alignment and jet launching. We describe
our setup in Sec. 2, present the results and discussion in Secs. 3 and 4, and
conclude in Sec. 5.

Figure 1: Vertical slice at $t=18,\\!091r_{g}/c$ through density of our tilted
disc with an aspect ratio $H/R=0.03$, the thinnest to date in a GRMHD
simulation, shows that the inner parts of the disc ($r\lesssim 5r_{g}$) align
with the rapidly spinning BH’s equator (horizontal in the figure; see the
inset for a zoom-in). The warp is accompanied by a $\sim 180^{\circ}$ sweep in
the precession angle (see Fig. 2b), creating a ‘negative’ tilt around $r\sim
8r_{g}$. Thick black curve shows that the disc tilt $\mathcal{T}$ smoothly
decreases toward small radii before flattening at $0^{\circ}$ for $r\lesssim
5r_{g}$ (see also Fig. 2a). This is the first demonstration of Bardeen &
Petterson (1975) alignment in a GRMHD simulation, i.e., in GR and in the
presence of large- and small-scale turbulent magnetic stresses. Magenta lines
show the corona-jet boundary, $p_{\rm b}=1.5\rho c^{2}$, and cyan lines the
disc-corona boundary, $\rho=1.2$.

##  2 Numerical models

We use for this work our state-of-the-art GPU-accelerated GRMHD code H-AMR
(Liska et al., 2018), which finds its heritage in the HARM2D code (Gammie et
al., 2003; Noble et al., 2006) but has been significantly expanded to take
advantage of vectorization and SIMD instructions, achieving $1.2\times 10^{6}$
zone-cycles per second on a single Intel Skylake 3.3GHz CPU core, and include
advanced features discussed below. We have developed a CUDA version of the
code that reaches $10^{8}$ zone-cycles/s on an NVIDIA Tesla V100 GPU. It also
features a staggered grid for constrained transport of magnetic fields
(Gardiner & Stone, 2005), uses an Harten-Lax-van Leer (HLL) Riemann solver
(Harten, 1983) and advanced features such as adaptive mesh refinement (AMR)
and local adaptive time-stepping (LAT) that bring down the cost of the
simulation described in this work by $2$ extra orders of magnitude (in
comparison to using a uniform grid and a global, fixed timestep).

Our thin disc model considers a spinning BH of spin parameter $a=0.9375$. We
insert an initial Fishbone & Moncrief torus (Fishbone & Moncrief, 1976), which
is maintained in hydrostatic equilibrium by the BH’s vertical component of
gravity counterbalancing the pressure forces in the disc. This torus has an
inner radius $r_{\rm in}=12.5r_{g}$, where $r_{g}=GM_{\rm BH}/c^{2}$ is the
gravitational radius, with the pressure maximum at $r_{\rm max}=25r_{g}$. We
use an ideal gas law equation of state, $p_{\rm g}=(\Gamma-1)u_{\rm g}$, with
$u_{\rm g}$ the gas internal energy density, $p_{g}$ the gas pressure and
$\Gamma=5/3$ the adiabatic index. We insert a poloidal magnetic field in the
torus described by a covariant vector potential
$A_{\phi}=(\rho-0.05)^{2}r^{3}$, with $\rho$ the gas density of the torus. The
magnetic field is subsequently normalized by setting $\bar{\beta}={\max p_{\rm
g}}/{\max p_{\rm b}}=30$, where $p_{\rm b}$ is the magnetic pressure and both
maxima are taken over the torus separately. The disc is tilted by
$\mathcal{T}_{0}=10^{\circ}$ with respect to the BH equator (see Liska et al.
2018 for details). Since the code is scale-free, we set the initial gas
density maximum to $\rho_{\rm max}=1$. To maintain the desired value of disc
thickness, $(H/R)_{\rm target}=0.03$, we cool the disc towards its target
temperature on the Keplerian timescale using a prescribed source term (Noble
et al. 2009). We disable this cooling function when $p_{b}/(\rho c^{2})\gtrsim
5$ at $r\gtrsim 10r_{g}$ to avoid cooling the jets.

In this work we for the first time use the full AMR capability of H-AMR, in
order to ensure sufficient resolution within the thin disc. For the refinement
criterion we use a density cutoff equal to $\sim 4$ percent of the maximum
disc density. The lower refinement levels are set such that jumps in spatial
resolution are limited to a factor $2$. To avoid noise from sporadic
refinement and derefinement, we only derefine at a two times lower density
than set as the refinement criterion. Typically, we use $3$ levels of AMR and
attain a speedup by a factor 32–60 in comparison to an equivalent uniform
grid. By evolving lower AMR levels or parts of the grid further from the black
hole, which have larger cell sizes, at a larger timestep, LAT gives an
additional speed-up of factor $\sim 5$ while reducing inversion errors in
relativistic regions (Chatterjee et al 2018, in prep). The effective
resolution of $2880\times 864\times 1200$ cells in $r$-, $\theta$-, and
$\phi$-directions, respectively, resolves the target disc thickness,
$(H/R)_{\rm target}=0.03$, by approximately $8$ cells in all 3 dimensions, and
the base grid of $720\times 216\times 300$ cells guarantees that the jets and
corona are also sufficiently resolved. As we will see in Sec. 3.5, the fastest
growing MRI wavelength is resolved by $\gtrsim 10$ cells over most of the
disc.

Operating in Kerr-Schild spherical polar coordinates, we place the inner
$r$-boundary just inside the event horizon and the outer $r$-boundary at
$10^{5}r_{g}$, so that the flow is unaffected by the boundaries. We use
transmissive polar boundary conditions in the $\theta$-direction (Liska et
al., 2018), and periodic boundary conditions in the $\phi$-direction.

In BH powered jets, the gas either drains off the field lines into the BH or
gets flung away along the field lines into the jet. This leads to a runaway
drop in density around the jet’s stagnation surface, at which the outflow
velocity vanishes. To avoid the development of vacuum regions and the
breakdown of ideal MHD, we replenish the density in the regions where the
density drops too low. For this, we follow the approach of Ressler et al.
(2017) and approximate physical processes that mass-load relativistic jets at
their base by applying a density floor of $\rho_{\rm floor}c^{2}=p_{\rm b}/10$
throughout the jet. This adds a small amount of density on the field lines and
does not noticeably affect the energetics of the jets.

##  3 Results

We start our analysis at $t=10^{4}r_{g}/c$ when the disc has cooled to its
equilibrium thickness, reaching a highly turbulent state, as seen in Figure 1.
We measure throughout this work vector quantities for the disc, corona and jet
in tilted spherical polar coordinates $r,\tilde{\theta},\tilde{\phi}$ that are
aligned with the disc’s angular momentum at each radius. We define the jet-
corona boundary at $p_{\rm b}=1.5\rho c^{2}$ and the corona-disc boundary at
$\rho=1.2$; this is approximately $10^{-2}$ times the maximum density in the
disc after the initial cooling completes.

###  3.1 Bardeen-Petterson Alignment

Figure 1 shows that the inner parts of the accretion disc ($r\lesssim 5r_{g}$)
align with the BH equator. This is the first demonstration of BP alignment in
GRMHD. This alignment happens relatively soon, already $\sim(200-500)r_{g}/c$
after accretion starts, and persists throughout the simulation. The outer
parts remain tilted, and the disc develops a smooth warp in between.

We can see this more quantitatively by introducing the tilt, $\mathcal{T}$,
and precession, $\mathcal{P}$, angles. We define them as polar and azimuthal
angles between the disc’s angular momentum and the BH spin vector (e.g.
Fragile & Anninos 2005; Fragile et al. 2007). While this definition works well
for the disc and the corona, it produces large fluctuations when applied to
jets since they can undergo strong kinks and pinches, which violently change
the jet angular momentum direction (by more than $10^{\circ}$) on very short
time and length scales. For this reason, we adopt a more robust way of
calculating the tilt and precession angles for the jets, as described in
Appendix B of Liska et al. (2018). Namely, instead of considering the angular
momentum vector orientation, we calculate the jet tilt and precession angles
at each radius from the jet’s average position, ($X,Y,Z$), weighted by the
magnetic pressure $p_{b}$:

$\epsilon=p_{\rm b}\times\operatorname{H}\left({p_{\rm b}}-1.5{\rho
c^{2}}\right),$ (1)

where $\operatorname{H}$ is the Heaviside step function that zeroes out the
weight outside of the magnetized jet body.

Figure 2:  Radial profiles of quantities in our thin disc system, averaged
over the time interval, $2\times 10^{4}<t<2.25\times 10^{4}r_{g}/c$. The
shaded regions represent the time variability over this interval. (a) Tilt
angles of the disc, corona and jet. Whereas the outer disc is tilted by
$\mathcal{T}\approx 5^{\circ}$, the inner disc at $r\lesssim 5r_{g}$ is
aligned with the BH equator, indicating the presence of Bardeen & Petterson
(1975) alignment. While the disc aligns rapidly near the BH, the corona and
jet show a single radial tilt oscillation peaking at $r\sim 10r_{g}$, which is
characteristic behaviour of thick accretion discs that do not show Bardeen-
Petterson alignment (e.g., Fragile et al. 2007; Liska et al. 2018). That the
corona and the jets at smaller radii align less readily with the BH than the
disc, suggests that the alignment is driven by the disc dynamics rather than
that of the corona or jet. (b) The precession angle for the disc, corona, and
jet are roughly consistent with each other at large radii, suggesting they
become co-aligned. As the disc aligns at small radii, $r\lesssim 5r_{g}$, the
precession angles become ill-defined, and we do not show them there. (c) The
half-opening angle of our jets (red curve) exceeds the opening angle of the
jets in a thick $H/R\sim 0.3$ disc model (blue dotted curve, Liska et al.,
2018), suggesting that a smaller disc thickness in this work provides less
pressure support for the jet and causes its opening angle to widen.

Figure 2(a) shows that the inner parts of the accretion disc, at $r\lesssim
5r_{g}$, are well-aligned with the BH equator: The tilt angle
$\mathcal{T}^{\rm inner}_{\rm disc}\ll 1^{\circ}$ of the inner disc
($r\lesssim 5r_{g}$) is much smaller than the tilt angle $\mathcal{T}^{\rm
outer}_{\rm disc}\gtrsim 5^{\circ}$ of the outer disc ($r\gtrsim 10r_{g}$).
The alignment radius remains steady at $r_{\rm BP}\approx 5r_{g}$ during the
course of our simulation. Note that $\mathcal{T}^{\rm outer}_{\rm disc}\simeq
5^{\circ}$ is smaller than the initial tilt of the disc,
$\mathcal{T}_{0}=10^{\circ}$. This is because the disc as a whole undergoes
global alignment, which is distinctly different than the BP effect (Sec. 3.3).

Since the corona is relatively thick, spanning an opening angle of around
$30-70^{\circ}$ (e.g. Figure 1), it is not expected to show BP alignment (Sec.
1). For such thick structures as the corona, the analytic theory predicts that
BP alignment is suppressed and accretion instead occurs in a misaligned
fashion through radial tilt oscillations (e.g. Ivanov & Illarionov 1997; Lubow
et al. 2002). Indeed, Fig. 2(a) shows a single radial tilt oscillation in the
corona and jet peaking around $r\sim 10r_{g}$. However, the corona still
manages to align within $r\lesssim 2.5r_{g}$, possibly due to torque from the
disc. Since the angular momentum of the corona is negligible, the disc can
easily affects its alignment.

While the inner jets are relatively closely aligned with the BH, the outer
jets are torqued into misalignment by the corona. This is consistent with the
previous work (Liska et al., 2018), which in the context of thick $H/R\sim
0.3$ discs found that the outer disc-corona system is responsible for
reorienting and collimating the jets. Similarly the precession angle of the
disc, corona and jet are closely related at large radii, as seen in Fig. 2(b).
At small radii, as the system becomes aligned, the precession angle becomes
ill-defined, and we do not show it.

Since the disc shows much ($\sim 2\times$) better alignment with the BH than
the jet and corona, it is unlikely that the jet can be responsible for
torquing the inner disc into (partial) alignment as has been demonstrated for
thicker discs (e.g., McKinney et al. 2013). If this were the case, one would
expect the jet tilt to be smaller than the disc tilt and, since the jet
actually has to transmit its torque through the corona to the disc, one would
expect better alignment for the corona as well.

###  3.2 Global precession and alignment

Figure 3: The tilt angle (a) and precession angle (b), as functions of time
for the entire disc (blue), coronal wind within $r<250r_{g}$ (black) and jet
(between $10r_{g}<r<100r_{g}$, red). The shaded regions show the standard
deviations within $10r_{g}<r<40r_{g}$. The whole system precesses and over
time aligns with the BH. The analytic estimate for the precession rate (green)
is only valid at later times. (c) The jet half-opening angle at $r=200r_{g}$
decreases as function of time, possibly due to a fall in jet power (Figure
4a,b).

In addition to exhibiting the BP effect, which depends on space and is
independent of time and in which the inner part of the disc aligns with the
BH, the entire accretion system additionally undergoes global alignment that
depends on time and is independent of radius. Figure 3(a),(b) shows the time
evolution of the average tilt, $\mathcal{T}$, and precession, $\mathcal{P}$,
angles for the disc, jet and corona. The disc and corona are angular momentum
averaged, the jet is magnetic pressure averaged (Sec. 3.1). Similar to thick
$H/R\sim 0.3$ discs (Liska et al. 2018; Liska et al 2018b, in prep), our thin
disc aligns _as a whole_ with the BH spin. Sorathia et al. (2013) proposed
that this happens by the turbulent mixing between disc annuli with different
precession angles (Fig. 2b shows that the precession angle decreases as
function of radius), which leads to cancellation of misaligned angular
momentum and thus produces net alignment.

As we discussed above, this is _not_ BP alignment for two reasons. First, BP
alignment is characterized by a steady state solution that is established on a
timescale shorter than an accretion time of the disc (our disc lost $\sim
20\%$ of its initial mass by the end of this simulation): otherwise, most of
the disc’s mass would accrete misaligned. Second, if the disc is fed
externally, such as by a larger thin disc or by fallback material in a tidal
disruption event, such global alignment may disappear as the inner disc is
replenished by the gas carrying the misaligned angular momentum on the
accretion time of the inner disc. Indeed, Liska et al. (2018) showed that
global alignment becomes slower as the disc size becomes larger due to viscous
spreading.

The precession period of around $10^{5}r_{g}/c$ is consistent with a Type-C
quasi-periodic oscillation (QPO) frequency of $0.2\rm{Hz}$ for a $10M_{\odot}$
BH (e.g. Ingram et al. 2016). Although Type-C QPOs are not observed for the
clean thin discs thought to be present in the soft state, we note that the
precessing corona in our simulation could give rise to the QPOs typically
observed in the Comptonised radiation from X-ray binaries during the
transition from hard to soft state. Note that the precession angle of the
corona is smaller than that of the disc, as seen in Fig. 3(b). Therefore, the
precession of the corona lags that of the disc. Thus, the hard variations
might lag the soft variations (see Liska et al. 2018b, in prep). Note that,
due to global alignment, which causes the outer disc to align with the black
hole to within $2^{\circ}$ by $6\times 10^{4}r_{g}/c$ (Fig. 3a), the
precession cannot be sustained for more than a single period.

Figure 4: Time-dependence of various quantities in our simulation. (a) Mass
$\dot{M}$ and energy $\dot{E}$ accretion rates, which initially oscillate but
stabilize later on. (b) The total outflow efficiency (black) and its
constituents: jets (blue), coronal wind (red), and radiative (green)
contributions. The total efficiency of $60{-}80\%$ is unprecedentedly high for
such a thin disc as ours and exceeds the standard thin disc efficiency,
$17.9\%$ (Novikov & Thorne, 1973), by a factor of $3$. Particularly surprising
is the high jet efficiency, $\eta_{\rm jet}\sim 20\%$, indicating that thin
discs are capable of producing powerful relativistic jets. (c) The magnetic
flux on the BH (blue) and in the disc (blue). While the disc magnetic flux,
$\Phi_{\rm disc}$ remains roughly constant at late times, $t\gtrsim 4\times
10^{4}r_{g}/c$, the flux on the BH decreases, suggesting that the BH might be
leaking its flux into the disc and that the jet formation by our thin disc
might be a transient phenomenon reminiscent of transient jets in XRBs (see
Sec. 4).

We can analytically estimate the precession period $P_{\rm LT}$ by calculating
the total perpendicular angular momentum $L_{\bot}$ and Lense-Thirring torque
$\tau_{\bot}$ directly from the stress energy tensor $T^{\mu}_{\nu}$ and test-
particle LT precession rate $\Omega_{\rm LT}$:

$\displaystyle L_{\bot}$ $\displaystyle=\iiint_{r_{\rm in}}^{r_{\rm
max}}T^{t}_{\tilde{\phi}}\sin\mathcal{T}\,dV,$ (2) $\displaystyle\Omega_{\rm
LT}$
$\displaystyle=\frac{1}{r^{3/2}+|a|}\left(1-\sqrt{1-4\frac{|a|}{r^{3/2}}+3\frac{a^{2}}{r^{2}}}\right),$
(3) $\displaystyle\tau_{\bot}$ $\displaystyle=\iiint_{r_{\rm in}}^{r_{\rm
max}}\Omega_{LT}\times T^{t}_{\tilde{\phi}}\sin\mathcal{T}\,dV,$ (4)
$\displaystyle P_{\rm LT}$
$\displaystyle=2\pi\times\frac{L_{\bot}}{\tau_{\bot}}.$ (5)

Figure 3(b) shows that the above analytic expression for the precession period
approximately agrees with the simulation at later times if we assume $r_{\rm
in}=1.5\times r_{\rm isco}$. However, at earlier times this overestimates the
precession rate. To remedy this, $r_{\rm in}$ would need to increase in time,
however it is unclear why this would happen.

Figure 5: Transverse slices of logarithm of density $\rho$ and the logarithm
of the magnitude of the proper velocity $\gamma v$, at $t\approx 6\times
10^{4}r_{g}/c$. The field lines (black) that connect to the BH launch
relativistic jets (demarcated by magenta lines at $p_{\rm b}=1.5\rho c^{2}$),
while field lines threading the disc give rise to a non-relativistic coronal
wind (disc-corona boundary is shown with the cyan line at $\rho=1.2$, see Sec.
3). The twin low-density polar jets accelerate to $\gamma v\sim 5c$, whereas
the winds are limited to at most $\gamma v\lesssim c$, or $v\lesssim 0.7c$,
which occurs near the edges of the jets.

###  3.3 Outflow and radiative efficiency

Figure 4(a) shows that even after the disc has approached the equilibrium
thickness at $t\sim 10^{4}r_{g}/c$, both the mass $\dot{M}$ and energy
$\dot{E}$ accretion rates at the event horizon, $r_{\rm
H}=r_{g}[1+(1-a^{2})^{1/2}]=1.35r_{g}$, show several peaks. They stabilise to
within a factor of a few at $t\gtrsim 4\times 10^{4}r_{g}/c$, suggesting that
by then the disc has reached a quasi-steady state. That $\dot{M}c^{2}>\dot{E}$
shows that the BH accretes more energy than it ejects out. To quantify this,
we introduce dimensionless energy outflow efficiencies for the relativistic
jets,

$\eta_{\rm jet}=-\frac{\dot{E}_{jet}}{\langle\dot{M}c^{2}\rangle_{t}},$ (6)

where $\dot{E}_{jet}$ is the energy accretion rate in the jets at $r=20r_{g}$
(chosen to be larger than $r_{\rm H}$ for robustness of identification of the
jets), and the disc winds,

$\eta_{\rm wind}=\frac{\langle\dot{M}c^{2}\rangle_{t}-\dot{\langle
E\rangle_{t}}+\dot{E}_{jet}}{\langle\dot{M}c^{2}\rangle_{t}}.$ (7)

In the above, e.g. $\langle\dot{M}\rangle_{t}$ is the running time-average of
the mass accretion rate over an interval of $\pm 500r_{g}/c$. Note that, from
these equations, any of the rest mass energy of disc material not accreted
onto the BH horizon, and not ejected in the jets, is by definition lost in a
wind. Figure 4(b) shows that the jet efficiency resides in the range
$\eta_{\rm jet}=20{-}50\%$, implying that as much as $50\%$ of the mass-energy
accreted by the BH can be extracted by large-scale magnetic flux from the
black hole spin energy and carried out in the form of Poynting-flux dominated
jets. Figure 5 shows that the field lines responsible for the jet launching
are anchored in the BH: they extract the energy via the Blandford & Znajek
(1977) effect. That the jets are so efficient is a surprising result: the
standard expectation is that thin discs are incapable of holding on to large-
scale poloidal magnetic flux and hence are not expected to have powerful jets
(Lubow et al., 1994).

In addition to the jets, the disc also launches a sub-relativistic disc wind,
which is magnetic pressure dominated and thus similar to a corona, with energy
efficiency $\eta_{wind}\sim 20\%$. The field lines threading this outflow are
anchored in the disc, suggesting it may be magneto-centrifugally driven by the
Blandford & Payne (1982) mechanism. However, the precise nature of the
launching mechanism is not entirely clear, since time-dependence, (magnetic)
pressure gradient forces, and buoyancy forces, as well as energy extraction
from the BH, may also play significant roles in launching this outflow. The
outflow has a radial velocity of around $0.02-0.2c$ between
$10r_{g}<r<500r_{g}$. This compares favorably to AGN ultra-fast outflows,
whose velocities lie in the range of $0.1{-}0.4c$ (e.g., Tombesi et al., 2010,
2011), but is on the (very) high end of winds detected in the soft state of
XRBs (e.g., Ponti et al., 2012, 2016; Miller et al., 2016).

To characterize the radial distribution of the magnetic flux, we define a 1D
poloidal magnetic flux function, $\Phi(B,r)$, as the maximum of the poloidal
magnetic flux at each radius of magnetic field $B$,

$\Phi(B,r)=\max_{\tilde{\theta}_{\rm m}}\iint_{0}^{\tilde{\theta}_{\rm
m}}B(r,\tilde{\theta},\tilde{\phi})dA_{\tilde{\theta}\tilde{\phi}},$ (8)

where $dA_{\tilde{\theta}\tilde{\phi}}=\sqrt{-g}d\tilde{\theta}d\tilde{\phi}$
is the area element and $g$ is the metric determinant. This gives us the
magnetic flux on the BH, $\Phi_{\rm BH}=0.5\Phi(|B_{r}|,r_{\rm H})$, and the
the magnetic flux content of the disc, $\Phi_{\rm
disc}=\max_{r}\Phi(B_{r},r)$. Figure 4(c) shows that at early times $\Phi_{\rm
BH}\simeq\Phi_{\rm disc}$, implying that most of the positively-oriented
magnetic flux resides on the BH and very little in the disc. Over time, both
initially decline roughly in the same proportion, but at $t\gtrsim 4\times
10^{4}r_{g}/c$ the value of $\Phi_{\rm disc}$ flattens out whereas $\Phi_{\rm
BH}$ continues to decline. This suggests that the magnetic flux diffuses out
of the BH and into the disc. The stability of $\Phi_{\rm BH}$ on short
timescales suggests that the disc is not in the magnetically arrested disc
regime (MAD, Narayan et al., 2003; Igumenshchev et al., 2003). However, it is
not yet known whether the magnetic flux expulsions and associated magnetic
flux variations, as characteristic of thick MADs (Tchekhovskoy et al., 2011),
are also present at the disc thickness considered here. Assuming that the
thinner the disc the easier it is to saturate the BH with magnetic flux
(Tchekhovskoy et al., 2014), it is likely that our disc is close to the
saturation of the MAD state, since the dimensionless magnetic flux ($\Phi_{\rm
BH}/(\dot{M}r_{g}^{2}c)^{1/2}\approx 25$) is only a factor of $2$ smaller than
the saturation value for thick discs (Tchekhovskoy et al., 2011; McKinney et
al., 2012).

Figure 6: Late-time vertical slice through density in our GRMHD simulation of
a thin disc (compare to Fig. 1 at earlier times). The disc reaches a steady
state where poloidal magnetic flux (black lines) gets trapped in the inner
disc. The density maximum in the disc at $\left|x\right|\sim 25r_{g}$ marks
the transition from the inner low-density, high-magnetization disc to the
high-density, low-magnetization outer disc (see also Fig. 7a). In Sec 4.2 we
argue that the high effective viscosity $\alpha_{\rm eff}\sim 1$ in the inner
disc may cause radiative cooling to become inefficient and the flow to
transition into an ADAF.

So far we have considered mechanical outflows. To compute the full efficiency,
we also need to account for radiative losses in the disc. A rough estimate for
the radiative efficiency is obtained by integrating the cooling rate,
$-\dot{U}$, used to keep our disc thin (at the equilibrium thickness), from
the photon orbit at $r_{\rm photon}=1.43r_{g}$ to a large enough radius beyond
which the disc luminosity is negligible, e.g., $r_{\rm max}=100r_{g}$, over
volume:

$\eta_{\rm rad}=\frac{-\iiint_{r_{\rm photon}}^{{\rm
r_{max}}}\dot{U}dV}{\langle\dot{M}c^{2}\rangle_{t}},$ (9)

where $dV=\sqrt{-g}drd\tilde{\theta}d\tilde{\phi}$ is the volume element. This
gives a radiative efficiency consistent with the Novikov & Thorne (1973) value
of $\sim 18\%$ (Fig. 4(b)), which is encouraging for the continuum fitting
method used to determine BH spin (see e.g. McClintock et al. 2014). However,
as also seen in Fig. 4(b), the total efficiency $\eta_{\rm tot}=\eta_{\rm
jet}+\eta_{\rm wind}+\eta_{\rm rad}$ of our disc reaches $60{-}80$%.

###  3.4 Jet geometry

A major surprise of this work is the finding of powerful jets, even in our
thin disc accretion system. How do the jets from thin discs compare to the
more familiar jets from thick discs? To carry out the comparison, we estimate
the cross-sectional area $A$ of the jet, defined by $p_{b}>1.5\rho c^{2}$, at
each radius r:

$A=\iint\operatorname{H}(p_{b}-1.5\rho c^{2})dA_{\tilde{\theta}\tilde{\phi}},$
(10)

Assuming the jet cross-section is circular, we obtain its effective half-
opening angle as

$\Delta\theta=\sqrt{\frac{A}{\pi r^{2}}}.$ (11)

Fig. 2(c) compares the radial dependence of our jet opening angle to that of a
jet from a thick $H/R\sim 0.3$ disc, as found in non-radiative GRMHD
simulations described in Liska et al. (2018). The jet opening angle in the
present work is much wider, presumably due to lack of the support in
collimation pressure from the much thinner and cooler ambient medium
consisting out of the disc and corona. This wider opening angle associated
with thin discs may lead to differences in the radiative flux (Fragile et al.
2012), possibly explaining the absence of any observations of soft state jets
in XRBs. Figure 3(c) shows that $\Delta\theta$ decreases over time, possibly
due to a decrease in jet power, as implied by the decreasing mass accretion
rate and approximately constant jet efficiency, as seen in Fig. 4(a),(b). If
the power of the jet is lower, its pressure is also lower, so the ambient
medium would compress the jet into a smaller opening angle.

###  3.5 Radial structure

Figure 7: Radial profiles averaged over $6\times 10^{4}<t<6.25\times
10^{4}r_{g}/c$. Radii within the ISCO ($r_{\rm ISCO}\approx 2.04r_{g}$) are
shaded grey. (a) The accretion disc transitions from an inner low-$\beta$ low-
density flow to an outer high-$\beta$ high-density flow at $r\approx 25r_{g}$.
(b) The radial dependence of the viscous stresses. Since the effective
viscosity exceeds the Reynolds and Maxwell stresses by one order of magnitude,
$\alpha_{\rm eff}\gg\alpha_{M,R}$, it is likely that large scale torques
transport angular momentum outwards. The dotted lines denote negative values
of $\alpha_{\rm eff}$. The disc is outflowing beyond the stagnation surface at
$r\approx 23r_{g}$. (c) Within $25r_{g}$ the density scale height
$(H/R)_{\rho}$ increases above the thermal scale height $(H/R)_{\rm thermal}$
due to magnetic pressure support in the inner disc. (d) We maintain more than
$10$ cells per MRI wavelength, $Q$, in all $3$ dimensions,
$r,\tilde{\theta},\tilde{\phi}$ over most of the disc.

Figure 6 shows a vertical slice through a late-time state of the system. The
accretion disc exhibits a sharp transition in density around $|x|\sim
25r_{g}$: at smaller radii, the accretion disc is lower-density and contains a
larger magnetic flux (larger number of magnetic field lines) than at larger
radii. We can see the same more quantitatively in Fig. 7, which shows disc
radial profiles at late time. The disc reaches a quasi-steady state with a
low-density, low plasma-$\beta$ inner disc truncated at $r\sim 25r_{g}$
coupled to a high density, high plasma-$\beta$ outer disc at larger radii.
Here $\beta$ is calculated by taking the ratio of the density weighted average
gas pressure $p_{g}$ and the density weighted average magnetic pressure
$p_{b}$. Note the presence of a density bump before the ISCO due to mass
accumulation, as is characteristic for thin discs (Penna et al. 2010). To
understand how this density and magnetic structure affects the disc dynamics,
we compute the Maxwell viscosity, $\alpha_{M}$, Reynolds viscosity,
$\alpha_{R}$, and effective viscosity, $\alpha_{\rm eff}$, parameters as
follows,

$\displaystyle\alpha_{M}$ $\displaystyle=\frac{\langle
b_{r}b_{\tilde{\phi}}\rangle_{\rho}}{\langle p_{\rm b}+p_{\rm
g}\rangle_{\rho}},$ (12) $\displaystyle\alpha_{R}$
$\displaystyle=\frac{\langle(\rho c^{2}+u_{\rm g}+p_{\rm g})\delta u_{r}\delta
u_{\tilde{\phi}}\rangle_{\rho}}{\langle p_{\rm b}+p_{\rm g}\rangle_{\rho}}$
(13) $\displaystyle\alpha_{\rm eff}$ $\displaystyle=\frac{\langle
u_{\tilde{\phi}}u_{r}\rangle_{\rho}}{\langle c_{i}^{2}\rangle_{\rho}},$ (14)
$\displaystyle c_{i}$ $\displaystyle=\sqrt{\frac{p_{\rm g}}{\rho c^{2}}}.$
(15)

Here $c_{i}$ is the isothermal sound speed and $\delta u=u-\langle
u\rangle_{\rho}$ is the velocity deviation from the mean. The density-weighted
average of a given quantity $\rm x$ is denoted by $\langle\rm
x\rangle_{\rho}$. The averaging is performed within $\sim 10$ scaleheights of
the disc’s midplane, more specifically for $\rho>0.01$.

Figure 7(b) shows that the Reynolds stress exceeds the Maxwell stress by a
factor of $\sim 2$ in the inner disc while in the outer disc the Maxwell
stress dominates. The effective viscosity exceeds the sum of Reynolds and
Maxwell viscosity contribution by an order of magnitude. This might be
indicative of large-scale magnetic torques removing the angular momentum from
the inner disc and/or large-scale torques between the disc and corona. In
fact, deviations of (GR)MHD results from the $\alpha$-viscosity description
are rather typical (Avara et al., 2016; Morales Teixeira et al., 2014;
McKinney et al., 2012; Penna et al., 2010; Sorathia et al., 2010) and
indicative of the intrinsic limitations of the $\alpha$-disc (Shakura &
Sunyaev, 1973), especially at greater field strengths. It is an interesting
question whether the reported discrepancy (King et al., 2007) between the
values for the $\alpha$-viscosity found in local MHD simulations (derived from
Maxwell and Reynolds stresses) and those constrained from observations
(derived from the effective viscosity set by accretion timescale) might be due
to the limitations of the $\alpha$-disc description rather than a mismatch
between the simulations and observations.

Interestingly, while the Maxwell and Reynolds stresses remain positive
throughout the disc, the effective viscosity becomes negative for $r>23r_{g}$,
causing the disc to spread out (Fig. 7(b)). Viscous spreading is
characteristic for all finite-size accretion discs and is caused by an
outwards flux of conserved angular momentum as mass moves inwards. However, in
contrast to the thick $H/R\sim 0.3$ discs considered in our previous work
(Liska et al., 2018), where viscous spreading caused the precession to stall,
it is unlikely to significantly affect the internal disc dynamics in this work
since the timescale for viscous spreading is more than an order of magnitude
longer than the simulation runtime.

Due to the global alignment of the disc over time, the tilt angle evolves from
the initial $10^{\circ}$ to $\sim 2^{\circ}$ at late times, and therefore our
simulation sweeps through a wide range of tilt angles. Throughout the
evolution, we do not see any evidence of significant tilt-related effects on
the internal disc dynamics, such as standing shocks/sharp entropy gradients
aligned with the lines of nodes (Fragile & Blaes, 2008).

We define the thermal and density disc scale heights, respectively, as

$\displaystyle(H/R)_{\rm thermal}$ $\displaystyle=\frac{\langle
c_{i}\rangle_{\rho^{2}}}{\langle v_{\tilde{\phi}}\rangle_{\rho^{2}}},$ (16)
$\displaystyle(H/R)_{\rm\rho}$
$\displaystyle=\langle|\tilde{\theta}-\tilde{\theta}_{avg}|\rangle_{\rho}$
(17)

Here $c_{i}$ is the isothermal sound speed and $v$ is 3-velocity, and
$\tilde{\theta}_{avg}$ is the average $\theta$-position of the disc’s midplane
in tilted coordinates. Figure 7(c) shows that while the thermal scale height
remains approximately constant for $r>3r_{g}$, $(H/R)_{\rm thermal}\approx
0.03$, the density scale height develops a bump, $(H/R)_{\rho}\sim 0.08$ at
$r\sim 10r_{g}$, due to the excess magnetic pressure support in a strongly
magnetized, $\beta\ll 1$, disc. Around the ISCO the scaleheight increases,
because the disc becomes super Keplerian, while the cooling function assumes a
Keplerian disc.

To quantify the degree to which our simulation resolves the magnetized
turbulence in the disc, we compute the quality factors
$Q_{r,\tilde{\theta},\tilde{\phi}}$, which give the number of cells per MRI
wavelength in each of the three directions:

$\displaystyle Q_{r,\tilde{\theta},\tilde{\phi}}$
$\displaystyle=\frac{2\pi}{\Delta_{r,\tilde{\theta},\tilde{\phi}}}\frac{\langle
v^{\rm
A}_{r,\tilde{\theta},\tilde{\phi}}\rangle_{w}}{\langle\Omega\rangle_{w}},$
(18) $\displaystyle v^{A}_{r,\tilde{\theta},\tilde{\phi}}$
$\displaystyle=\sqrt{\frac{b_{r,\tilde{\theta},\tilde{\phi}}^{2}}{(\rho
c^{2}+u_{\rm g}+p_{\rm g}+b^{2})}},$ (19) $\displaystyle\Omega$
$\displaystyle=\frac{v_{\tilde{\phi}}}{r},$ (20)

where $\Delta_{r,\tilde{\theta},\tilde{\phi}}$ is the size of the cell. The
Alfven speed, $v_{\rm A}$, and angular frequency, $\Omega$, are volume-
averaged with $w=\sqrt{b^{2}\rho}$ as the weight since a large fraction of the
mass in our thin discs resides in equatorial current sheets where the magnetic
field vanishes. Figure 7(d) shows that the whole disc is well resolved with
$Q_{r,\tilde{\theta}}>30$, $Q_{\tilde{\phi}}>200$ in the inner disc and
$Q_{r,\tilde{\theta}}>10$, $Q_{\tilde{\phi}}>70$ over most of the outer disc.
This satisfies the numerical convergence criteria for MRI turbulence,
$Q_{\theta}>10$ and $Q_{\phi}>20$ (see e.g. Hawley et al. 2011).

##  4 Discussion

###  4.1 Bardeen-Petterson Alignment

We found that the inner $\sim 5r_{g}$ of a thin disc, $H/R\simeq 0.03$,
initially tilted by $10^{\circ}$ relative to the central spinning BH,
undergoes alignment with the BH equator (Sec. 3.1). This is the first
demonstration of the Bardeen & Petterson (1975) effect in a GRMHD simulation,
in the presence of non-local and anisotropic turbulent MHD stresses.

This confirmation of the BP effect has profound consequences for the growth
and spin evolution of supermassive BHs (SMBHs), since BP alignment is a
crucial ingredient that has been assumed to take place for misaligned
accretion episodes (e.g. Volonteri et al. 2005; King & Pringle 2006; Fanidakis
et al. 2011). Because the alignment radius acts as lever arm helping to torque
the BH, the BP effect can torque the BH and align its spin vector with the
outer misaligned accretion flow on much shorter timescales than otherwise
(Rees, 1978; Scheuer & Feiler, 1996; Natarajan & Pringle, 1998). This rapid
reorientation of BH spin has the potential to create the right conditions for
rapid BH spin-up. If BH spin reorientation occurs on a shorter timescale than
the timescale of a single accretion episode in a chaotic accretion scenario
(in which the direction of the supplied gas angular momentum randomly changes
between different accretion episodes, see e.g. Volonteri et al. 2005), then
supermassive BHs can be efficiently spun up (e.g. Natarajan & Pringle 1998).
In the opposite case, the accretion-supplied angular momenta would tend to
cancel out, and the central BHs would be on average spun down.

Due to the high cost of GRMHD simulations, it is appealing to use them to
calibrate both computationally cheaper SPH simulations and analytic theory.
Our simulations are in the viscosity-dominated, $H/R<\alpha$ regime
(Papaloizou & Pringle, 1983) since both our effective viscosity parameter,
$\alpha_{\rm eff}\simeq 1$, and local Maxwell plus Reynolds stress related
viscosity parameter, $\alpha_{M}+\alpha_{R}\lesssim 0.1$, exceed $H/R=0.03$
(see Figure 7). However, we did not find Bardeen-Petterson alignment in all of
our simulations featuring thicker discs ($0.05<H/R<0.1$; Liska et al 2018b in
prep; Liska et al 2018c in prep). Because all of these simulations have
$\alpha_{M,R}\lesssim H/R<\alpha_{\rm eff}$, it appears likely that values of
local Maxwell and Reynolds stresses determine the transition from the wave-
($H/R>\alpha$) to the viscosity- ($H/R<\alpha$) dominated regime, where BP
alignment is expected.

However, the alignment radius substantially differs from both the analytical
theory and SPH simulations. Specifically, Kumar & Pringle (1985) found
analytically using the corrected Bardeen & Petterson (1975) equations in
Papaloizou & Pringle (1983) a radius $r_{\rm BP}\approx 300$ $r_{g}$, while
Nelson & Papaloizou (2000) find $r_{\rm BP}\approx 100$ $r_{g}$ for
$\alpha=0.05$ and $h/r=0.03$. SPH simulations of Nelson & Papaloizou (2000)
find $r_{\rm BP}\approx 30r_{g}$ for $\alpha=0.1$ and $h/r=0.03$. Other more
recent SPH work considers a smaller disc thickness $h/r\approx 0.013$ and
smaller viscosity $\alpha\approx 0.03$ (Lodato & Price, 2010), which makes
comparison to our work difficult. However, assuming a scaling relation,
$r_{\rm BP}\sim[(h/r)^{-2}/\alpha]^{4/7}$ (Kumar & Pringle, 1985), these
results are in similar disagreement.

One contributing factor to this disagreement with SPH simulations may be the
$\sim 2$ times larger density scale height, caused by the buildup of magnetic
pressure, between $5r_{g}<r<20r_{g}$ in our simulations (see Sec. 3.5). We
indeed see that $r_{\rm BP}$ (Fig. 2a) coincides with the radius where the
disc becomes thicker (Fig. 7(c)). Naively, analytically one would only predict
(Kumar & Pringle, 1985) an $\approx 0.4$ times smaller alignment radius for a
similar increase in disc thickness, insufficient to account for the full
extent of the discrepancy. However since $H/R\approx\alpha_{R}+\alpha_{M}$ the
disc may not be in the $H/R<\alpha$ viscosity dominated regime required for BP
alignment. Another contributing factor might be the presence of large-scale
magnetic torques in the system, which can affect the alignment radius in at
least two ways. First, as discussed in Sec. 3.5, large-scale magnetic torques
can remove angular momentum from the disc and increase the inflow velocity.
Because BP alignment is expected to occur more rapidly when the radial inflow
velocity is smaller (Bardeen & Petterson, 1975), this might reduce the
efficiency of BP effect (Nealon et al. 2015). Second, these torques might
induce coupling between the inner aligned disc and outer misaligned disc-
corona-jet system moving $r_{\rm BP}$ inwards.

A smaller $r_{\rm BP}$ could have a significant effect on the predictions of
SMBH growth models because smaller values of $r_{\rm BP}$ lead to less rapid
alignment between BH and outer disc and perhaps consequently less rapid spin
up. For the same reason, our result implies that initially misaligned X-ray
binary systems will take even longer to align than previously predicted (King
& Nixon, 2016), indicating that there could be many misaligned X-ray binaries
today, as implied by the Lense-Thirring precession QPO model of Ingram et al.
(2009). Future work should study the effect of BP alignment and disc warp on
the measured values of BH spin (e.g., McClintock et al., 2014).

The inclusion of a gas pressure dominated equation of state with adiabatic
index $\Gamma=5/3$ as in this work is only applicable in the outer accretion
disc of X-ray binaries (see e.g. Zhu & Narayan 2013), while the inner part
could be radiation pressure dominated with $\Gamma=4/3$. We find that changing
the adiabatic index to $\Gamma=4/3$ for a thin disc tilted by $45^{\circ}$
behaves qualitatively similar to $\Gamma=5/3$ (Liska et al 2018c, in prep).
This is not unexpected since the most dominant effect of a softer equation of
state is that the disc becomes thinner for a given specific internal energy
$u_{g}/(\rho c^{2})$. However, the cooling function (Noble et al., 2009)
automatically adapts to the equation of state in order to maintain the desired
thermal scale height and thus the absence of any strong dependence on the
adiabatic indices is not unexpected.

###  4.2 Disc Evaporation?

The transition from a low-viscosity, high-density outer disc into a high-
viscosity, low-density inner disc (See Sec. 3.5 and Fig. 6/7) might provide
clues into a long-standing puzzle in accretion physics: How do cool thin discs
(Shakura & Sunyaev, 1973) evaporate into hot thick radiatively-inefficient
accretion flows (Narayan & Yi, 1994)? At low mass accretion rates, the ions
and electron in the accretion flow become weakly coupled, and the cooling
timescale becomes limited by the timescale for Coulomb collisions to
equilibrate the temperature of the hot non-radiative ions with the
radiatively-cooled electrons (Shapiro et al., 1976). While our cooling
function does not include this weak ion-electron coupling, the shorter
accretion time due to the high, order unity, $\alpha$-viscosity in the inner
disc may cause the inflow timescale to become shorter than the electron-ion
coupling timescale. This would prevent the inner disc from cooling and could
conceivably lead to a radiatively inefficient thick accretion flow at a radius
$r\lesssim 25r_{g}$. Indeed, the finite timescale for electron-ion coupling
implies that one would generally expect such a disc to form for
$\dot{M}<\alpha^{2}\dot{M}_{\rm Edd}$ (e.g. Esin et al. 1997). This
evaporation mechanism, through the elevated $\alpha$-viscosity, is attractive
in that it does not require conduction of heat from the corona to the disc
(e.g. Meyer & Meyer-Hofmeister 1994; Liu et al. 1999; Czerny et al. 2000; Qian
et al. 2007).

How does the high-viscosity, low-$\beta$ inner disc, seen in Figure 7(a),
form? In our simulation it may have formed due to the rapid cooling of an
initial torus threaded with poloidal magnetic flux (Sikora & Begelman, 2013;
Begelman & Armitage, 2014): the cooling causes the thermal pressure to
decrease, but – due to vertical magnetic flux conservation – the magnetic flux
stays about the same. This causes the disc to become more strongly magnetized
and plasma $\beta$ to drop. Shearing box simulations seeded with strong
vertical magnetic flux appear to develop a similarly highly magnetized
accretion state with strong outflows (Salvesen et al., 2016; Bai & Stone,
2013). If this scenario is indeed the case, it would require the presence of
large scale magnetic flux in the accretion disc prior to the disc becoming
thin, which would limit the applicability of this simulation to the
intermediate states for X-ray binaries.

Future work will investigate the effect of the different initial magnetic
field geometries (e.g., Liska et al., 2018), exploring if large scale poloidal
magnetic flux is indeed a necessary ingredient for the high $\alpha$-viscosity
inner disc. It will also include electron-ion coupling, and on-the-fly
radiation transfer, to accurately model the cooling of the disc.

###  4.3 Jet Launching

This work shows that thin discs down to at least $H/R=0.03$ can efficiently
launch relativistic Blandford & Znajek (1977) jets of substantial power,
carrying out $\gtrsim 20\%$ of the accretion power, over timescales comparable
to the accretion time (Sec. 3.4). This suggests that even such thin discs as
considered in this work are capable of retaining for their accretion time
large-scale poloidal (vertical) magnetic flux on the BH, a necessary
ingredient for launching relativistic jets (Blandford & Znajek, 1977). This is
particularly interesting given that simple analytical arguments suggest that
thin discs should lose their large-scale magnetic flux to outward diffusion
(Lubow et al., 1994). It is possible that the large effective viscosity caused
by large scale external torques may overcome this problem by dragging flux
inwards before it has time to diffuse out (see also Guilet & Ogilvie 2012,
2013).

How can we reconcile the formation of powerful jets from thin discs with
observations? There are no observations that have convincingly detected jets
from thin discs in the soft state of X-ray binaries (though see Rushton et al.
2012), however, about $10\%$ of quasars are radio loud and form jets (Sikora
et al., 2007). Because our simulated jets have wide opening angles,
$\Delta\theta\sim 20^{\circ}$ (Fig. 2c), they might become less optically
thick and more difficult to detect (see also Russell et al. 2011; Fragile et
al. 2012). Another possible explanation is that our simulations do not apply
to the soft state of XRBs (Sec. 4.4).

###  4.4 A Transitional Disc?

An interesting possibility is that our simulations apply to _transitional_
discs, in the middle of the hard-to-soft state transition (e.g. Fender et al.
2004). In fact, we set up our simulations in a very similar way: the initial
thick torus rapidly cools down to the target thickness, $H/R=0.03$, which is
much smaller than the initial thickness, $H/R\sim 0.3$. Since thick accretion
discs may be able to generate and advect large scale poloidal magnetic flux
through large scale dynamo action (e.g., Liska et al., 2018), they are
expected to retain a substantial amount of it after their collapse into a thin
disc. This can lead both to a highly viscous inner disc that evaporates into
an ADAF (Sec. 4.3) and sustains a strong jet (Sec. 4.3, see also Sikora &
Begelman 2013; Begelman & Armitage 2014). Indeed XRBs in the hard-to-soft
state transition are known to produce jets (e.g. Fender et al. 2004).
Similarly, it is possible that radio-loud quasars contain such transitional
discs (Tchekhovskoy, 2015).

However as proposed in Lubow et al. (1994) this flux may slowly diffuse out
and cause the jet to shut down. In addition, if this large scale poloidal
magnetic field indeed leads to evaporation of the inner disc into an ADAF
(Sec. 4.2), the truncation radius between the inner thick and outer thin disc
will move inwards. This is consistent with observational evidence of the
truncation radius moving in during the evolution towards the soft state (Esin
et al., 1997; Done et al., 2007; Ingram & Done, 2011). Figure 4(c) indeed
shows signs of magnetic flux diffusing out of the BH: the flux in the disc
($\Phi_{\rm disc}$) stays roughly constant while the flux on the BH
($\Phi_{\rm BH}$) drops. However, the drop is small and appears to be leveling
off. Several mechanisms have been suggested that can prevent the poloidal
magnetic flux from diffusing out in thin discs (Rothstein & Lovelace, 2008;
Guilet & Ogilvie, 2012, 2013). Future simulations spanning much longer
runtimes can probe if thin discs are able to retain poloidal magnetic flux for
a more extended time period, or are always transitional.

##  5 Conclusions

In this work we have performed the thinnest disc GRMHD simulations to date. We
started with an $H/R=0.03$ accretion disc tilted by $10^{\circ}$ relative to a
rapidly spinning $a=0.9375$ BH. Using 3 AMR levels, we carried out GRMHD
simulations at sufficiently high effective resolution, $2880\times 864\times
1200$, which for the first time resolved the MRI turbulence in a thin disc in
all 3 dimensions with near-cubical cells (of order unity aspect ratio). Our
results can be summarized in 3 key points.

First, we have confirmed for the first time that the inner parts of tilted
thin discs can align with the BH equatorial plane as theorized 40 years ago by
Bardeen & Petterson (1975), even when the full effects of GR, anisotropic MRI
turbulence and torquing of the disc by magnetized corona and jets are
included. The disc aligns with the BH within the BP radius, $r_{\rm BP}\simeq
5r_{g}$, whose value is expected to increase for thinner discs (e.g. Kumar &
Pringle 1985). The development of a BP configuration can have profound
consequences for the evolution of BH spins in AGN, as the large lever arm of
$r_{\rm BP}$ out to which the disc is aligned can torque the BH into alignment
with the outer, tilted disc on a much shorter timescale than without the BP
effect (e.g. Scheuer & Feiler 1996; Natarajan & Pringle 1998).

Second, we have shown that an accretion disc can develop an inner low-density,
high-viscosity disc coupled to an outer high-density, low-viscosity disc at
$r\lesssim 25r_{g}$. We suggested that the order unity viscosity of the inner
disc we find might lead to it evaporating into a radiatively inefficient
accretion flow when the electron-ion coupling time exceeds the accretion time
(e.g. Esin et al. 1997). This high viscosity may be caused by the presence in
the initial conditions of large-scale poloidal magnetic flux, which removes
the angular momentum through large-scale outflows. Large scale poloidal
magnetic flux may be present in thin discs during hard-to-soft state
transitions (e.g. Sikora & Begelman 2013; Begelman & Armitage 2014), as
discussed in Sec. 4.4.

Third, we have shown that BH accretion systems with thin discs, if initially
threaded with large scale poloidal magnetic flux, can launch powerful
Blandford & Znajek (1977) jets on the viscous timescale, with their power
reaching $20{-}50\%$ of the accretion power. This challenges the standard
paradigm that thin discs in the soft state cannot advect inwards poloidal
magnetic flux needed to launch jets (Lubow et al., 1994) and is seemingly in
tension with the lack of any clear detection of jets in X-ray binaries.
However the morphology of our jets, specifically their twice as large opening
angle as of those produced by thick discs (e.g. Liska et al. 2018; McKinney
2006; Chatterjee et al 2018a, in prep), may make them more optically thin and
thus more difficult to detect (see also Fragile et al. 2012). Another
possibility is that our simulations describe transitional discs in the hard-
to-soft state transition which are known to produce powerful jets (e.g. Fender
et al. 2004) and, like our simulations, may naturally harbor large scale
poloidal magnetic flux (Sikora & Begelman, 2013), which is required to produce
powerful jets (Blandford & Znajek, 1977). This flux may eventually diffuse out
(e.g. Begelman & Armitage 2014) causing the jets to shut down. Outwards flux
diffusion might indeed be present in our simulation (Sec. 4.4).

##  6 Acknowledgments

We thank Chris Fragile and Cole Miller for useful suggestions. AI thanks James
Matthews for useful discussions. This research was made possible by NSF PRAC
award no. 1615281 and OAC-1811605 at the Blue Waters sustained-petascale
computing project and supported in part under grant no. NSF PHY-1125915. ML
and MK were supported by the Netherlands Organisation for Scientific Research
(NWO) Spinoza Prize, AI by the Royal Society URF, AT by Northwestern
University and the TAC and NASA Einstein (grant no. PF3-140131) postdoctoral
fellowships.

##  7 Supporting Information

Additional Supporting Information may be found in the online version of this
article: movie file (link).

## References

  * Avara et al. (2016) Avara M. J., McKinney J. C., Reynolds C. S., 2016, MNRAS, 462, 636
  * Bai & Stone (2013) Bai X.-N., Stone J. M., 2013, ApJ, 767, 30
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Bardeen & Petterson (1975) Bardeen J. M., Petterson J. A., 1975, ApJ, 195, L65
  * Beckwith et al. (2008) Beckwith K., Hawley J. F., Krolik J. H., 2008, ApJ, 678, 1180
  * Begelman & Armitage (2014) Begelman M. C., Armitage P. J., 2014, ApJ, 782, L18
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Caproni et al. (2007) Caproni A., Abraham Z., Livio M., Mosquera Cuesta H. J., 2007, MNRAS, 379, 135
  * Caproni et al. (2006) Caproni A., Abraham Z., Mosquera Cuesta H. J., 2006, ApJ, 638, 120
  * Czerny et al. (2000) Czerny B., Różańska A., Janiuk A., Życki P. T., 2000, New Astronomy Reviews, 44, 439
  * De Villiers et al. (2003) De Villiers J.-P., Hawley J. F., Krolik J. H., 2003, ApJ, 599, 1238
  * Done et al. (2007) Done C., Gierliński M., Kubota A., 2007, A&A Rev., 15, 1
  * Esin et al. (1997) Esin A. A., McClintock J. E., Narayan R., 1997, ApJ, 489, 865
  * Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  * Fanidakis et al. (2011) Fanidakis N., Baugh C. M., Benson A. J., Bower R. G., Cole S., Done C., Frenk C. S., 2011, MNRAS, 410, 53
  * Fender et al. (2004) Fender R. P., Belloni T. M., Gallo E., 2004, Monthly Notices of the Royal Astronomical Society, 355, 1105
  * Fishbone & Moncrief (1976) Fishbone L. G., Moncrief V., 1976, ApJ, 207, 962
  * Fragile & Anninos (2005) Fragile P. C., Anninos P., 2005, ApJ, 623, 347
  * Fragile & Blaes (2008) Fragile P. C., Blaes O. M., 2008, The Astrophysical Journal, 687, 757
  * Fragile et al. (2007) Fragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007, ApJ, 668, 417
  * Fragile et al. (2012) Fragile P. C., Wilson J., Rodriguez M., 2012, MNRAS, 424, 524
  * Gammie et al. (2003) Gammie C. F., McKinney J. C., Tóth G., 2003, ApJ, 589, 444
  * Gardiner & Stone (2005) Gardiner T. A., Stone J. M., 2005, Journal of Computational Physics, 205, 509
  * Greene et al. (2001) Greene J., Bailyn C. D., Orosz J. A., 2001, ApJ, 554, 1290
  * Guilet & Ogilvie (2012) Guilet J., Ogilvie G. I., 2012, MNRAS, 424, 2097
  * Guilet & Ogilvie (2013) Guilet J., Ogilvie G. I., 2013, MNRAS, 430, 822
  * Harten (1983) Harten A., 1983, Journal of Computational Physics, 49, 357
  * Hawley et al. (2011) Hawley J. F., Guan X., Krolik J. H., 2011, ApJ, 738, 84
  * Hjellming & Rupen (1995) Hjellming R. M., Rupen M. P., 1995, Nature, 375, 464
  * Ichimaru (1977) Ichimaru S., 1977, ApJ, 214, 840
  * Igumenshchev et al. (2003) Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ, 592, 1042
  * Ingram & Done (2011) Ingram A., Done C., 2011, MNRAS, 415, 2323
  * Ingram et al. (2009) Ingram A., Done C., Fragile P. C., 2009, MNRAS, 397, L101
  * Ingram et al. (2016) Ingram A., van der Klis M., Middleton M., Done C., Altamirano D., Heil L., Uttley P., Axelsson M., 2016, MNRAS, 461, 1967
  * Ivanov & Illarionov (1997) Ivanov P. B., Illarionov A. F., 1997, MNRAS, 285, 394
  * King & Nixon (2016) King A., Nixon C., 2016, MNRAS, 462, 464
  * King et al. (2005) King A. R., Lubow S. H., Ogilvie G. I., Pringle J. E., 2005, MNRAS, 363, 49
  * King & Pringle (2006) King A. R., Pringle J. E., 2006, MNRAS, 373, L90
  * King et al. (2007) King A. R., Pringle J. E., Livio M., 2007, MNRAS, 376, 1740
  * Kumar & Pringle (1985) Kumar S., Pringle J. E., 1985, MNRAS, 213, 435
  * Lense & Thirring (1918) Lense J., Thirring H., 1918, Physikalische Zeitschrift, 19
  * Liska et al. (2018) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2018, MNRAS, 474, L81
  * Liska et al. (2018) Liska M. T. P., Tchekhovskoy A., Quataert E., 2018, ArXiv:1809.04608
  * Liu et al. (1999) Liu B. F., Yuan W., Meyer F., Meyer-Hofmeister E., Xie G. Z., 1999, ApJ, 527, L17
  * Lodato & Price (2010) Lodato G., Price D. J., 2010, MNRAS, 405, 1212
  * Lodato & Pringle (2007) Lodato G., Pringle J. E., 2007, MNRAS, 381, 1287
  * Lubow et al. (2002) Lubow S. H., Ogilvie G. I., Pringle J. E., 2002, MNRAS, 337, 706
  * Lubow et al. (1994) Lubow S. H., Papaloizou J. C. B., Pringle J. E., 1994, Monthly Notices of the Royal Astronomical Society, 267, 235
  * Maccarone (2002) Maccarone T. J., 2002, MNRAS, 336, 1371
  * McClintock et al. (2014) McClintock J. E., Narayan R., Steiner J. F., 2014, Space Sci. Rev., 183, 295
  * McClintock & Remillard (2006) McClintock J. E., Remillard R. A., 2006, Black hole binaries. pp 157–213
  * McKinney (2006) McKinney J. C., 2006, MNRAS, 368, 1561
  * McKinney & Blandford (2009) McKinney J. C., Blandford R. D., 2009, MNRAS, 394, L126
  * McKinney et al. (2012) McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, MNRAS, 423, 3083
  * McKinney et al. (2013) McKinney J. C., Tchekhovskoy A., Blandford R. D., 2013, Science, 339, 49
  * Meyer & Meyer-Hofmeister (1994) Meyer F., Meyer-Hofmeister E., 1994, A&A, 288, 175
  * Miller et al. (2016) Miller J. M., Raymond J., Cackett E., Grinberg V., Nowak M., 2016, ApJ, 822, L18
  * Morales Teixeira et al. (2017) Morales Teixeira D., Avara M. J., McKinney J. C., 2017, ArXiv e-prints
  * Morales Teixeira et al. (2014) Morales Teixeira D., Fragile P. C., Zhuravlev V. V., Ivanov P. B., 2014, ApJ, 796, 103
  * Narayan et al. (2003) Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ, 55, L69
  * Narayan & Yi (1994) Narayan R., Yi I., 1994, ApJ, 428, L13
  * Narayan & Yi (1995a) Narayan R., Yi I., 1995a, ApJ, 444, 231
  * Narayan & Yi (1995b) Narayan R., Yi I., 1995b, ApJ, 452, 710
  * Natarajan & Pringle (1998) Natarajan P., Pringle J. E., 1998, ApJ, 506, L97
  * Nealon et al. (2015) Nealon R., Price D. J., Nixon C. J., 2015, MNRAS, 448, 1526
  * Nelson & Papaloizou (2000) Nelson R. P., Papaloizou J. C. B., 2000, MNRAS, 315, 570
  * Noble et al. (2006) Noble S. C., Gammie C. F., McKinney J. C., Del Zanna L., 2006, ApJ, 641, 626
  * Noble et al. (2009) Noble S. C., Krolik J. H., Hawley J. F., 2009, ApJ, 692, 411
  * Noble et al. (2010) Noble S. C., Krolik J. H., Hawley J. F., 2010, ApJ, 711, 959
  * Novikov & Thorne (1973) Novikov I. D., Thorne K. S., 1973, in Dewitt C., Dewitt B. S., eds, Black Holes (Les Astres Occlus) Astrophysics of black holes.. pp 343–450
  * Ogilvie (1999) Ogilvie G. I., 1999, MNRAS, 304, 557
  * Papaloizou & Pringle (1983) Papaloizou J. C. B., Pringle J. E., 1983, MNRAS, 202, 1181
  * Penna et al. (2010) Penna R. F., McKinney J. C., Narayan R., Tchekhovskoy A., Shafee R., McClintock J. E., 2010, MNRAS, 408, 752
  * Piran et al. (2015) Piran T., Sa̧dowski A., Tchekhovskoy A., 2015, MNRAS, 453, 157
  * Ponti et al. (2016) Ponti G., Bianchi S., Muñoz-Darias T., De K., Fender R., Merloni A., 2016, Astronomische Nachrichten, 337, 512
  * Ponti et al. (2012) Ponti G., Fender R. P., Begelman M. C., Dunn R. J. H., Neilsen J., Coriat M., 2012, MNRAS, 422, L11
  * Pringle (1992) Pringle J. E., 1992, MNRAS, 258, 811
  * Qian et al. (2007) Qian L., Liu B. F., Wu X.-B., 2007, The Astrophysical Journal, 668, 1145
  * Rees (1978) Rees M. J., 1978, Nature, 275, 516
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
  * Ressler et al. (2017) Ressler S. M., Tchekhovskoy A., Quataert E., Gammie C. F., 2017, MNRAS, 467, 3604
  * Reynolds (2014) Reynolds C. S., 2014, Space Sci. Rev., 183, 277
  * Rothstein & Lovelace (2008) Rothstein D. M., Lovelace R. V. E., 2008, ApJ, 677, 1221
  * Rushton et al. (2012) Rushton A., Miller-Jones J. C. A., Campana R., Evangelista Y., Paragi Z., Maccarone T. J., Pooley G. G., Tudose V., Fender R. P., Spencer R. E., Dhawan V., 2012, MNRAS, 419, 3194
  * Russell et al. (2011) Russell D. M., Miller-Jones J. C. A., Maccarone T. J., Yang Y. J., Fender R. P., Lewis F., 2011, ApJ, 739, L19
  * Salvesen et al. (2016) Salvesen G., Simon J. B., Armitage P. J., Begelman M. C., 2016, MNRAS, 457, 857
  * Scheuer & Feiler (1996) Scheuer P. A. G., Feiler R., 1996, MNRAS, 282, 291
  * Shafee et al. (2008) Shafee R., McKinney J. C., Narayan R., Tchekhovskoy A., Gammie C. F., McClintock J. E., 2008, ApJ, 687, L25
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shapiro et al. (1976) Shapiro S. L., Lightman A. P., Eardley D. M., 1976, ApJ, 204, 187
  * Sikora & Begelman (2013) Sikora M., Begelman M. C., 2013, ApJ, 764, L24
  * Sikora et al. (2007) Sikora M., Łukasz Stawarz Lasota J.-P., 2007, The Astrophysical Journal, 658, 815
  * Sorathia et al. (2013) Sorathia K. A., Krolik J. H., Hawley J. F., 2013, ApJ, 777, 21
  * Sorathia et al. (2010) Sorathia K. A., Reynolds C. S., Armitage P. J., 2010, ApJ, 712, 1241
  * Tchekhovskoy (2015) Tchekhovskoy A., 2015, in Contopoulos I., Gabuzda D., Kylafis N., eds, The Formation and Disruption of Black Hole Jets Vol. 414 of Astrophysics and Space Science Library, Launching of Active Galactic Nuclei Jets. p. 45
  * Tchekhovskoy & McKinney (2012) Tchekhovskoy A., McKinney J. C., 2012, MNRAS, 423, L55
  * Tchekhovskoy et al. (2014) Tchekhovskoy A., Metzger B. D., Giannios D., Kelley L. Z., 2014, MNRAS, 437, 2744
  * Tchekhovskoy et al. (2011) Tchekhovskoy A., Narayan R., McKinney J. C., 2011, MNRAS, 418, L79
  * Tombesi et al. (2011) Tombesi F., Cappi M., Reeves J. N., Palumbo G. G. C., Braito V., Dadina M., 2011, ApJ, 742, 44
  * Tombesi et al. (2010) Tombesi F., Cappi M., Reeves J. N., Palumbo G. G. C., Yaqoob T., Braito V., Dadina M., 2010, A&A, 521, A57
  * Volonteri et al. (2005) Volonteri M., Madau P., Quataert E., Rees M. J., 2005, The Astrophysical Journal, 620, 69
  * Zhu & Narayan (2013) Zhu Y., Narayan R., 2013, MNRAS, 434, 2262
  * Zhuravlev et al. (2014) Zhuravlev V. V., Ivanov P. B., Fragile P. C., Morales Teixeira D., 2014, ApJ, 796, 104