ABSTRACT
Title of dissertation: MAGNETOHYDRODYNAMIC SIMULATIONS OF
BLACK HOLE ACCRETION
Mark J. Avara, Doctor of Philosophy, 2017
Dissertation directed by: Professors Jonathan McKinney and
Christopher S. Reynolds
University of Maryland, College Park
Black holes embody one of the few, simple, solutions to the Einstein field
equations that describe our modern understanding of gravitation. In isolation they
are small, dark, and elusive. However, when a gas cloud or star wanders too close,
they light up our universe in a way no other cosmic object can. The processes of
magnetohydrodynamics which describe the accretion inflow and outflows of plasma
around black holes are highly coupled and nonlinear and so require numerical ex-
periments for elucidation. These processes are at the heart of astrophysics since
black holes, once they somehow reach super-massive status, influence the evolution
of the largest structures in the universe. It has been my goal, with the body of work
comprising this thesis, to explore the ways in which the influence of black holes on
their surroundings differs from the predictions of standard accretion models. I have
especially focused on how magnetization of the greater black hole environment can
impact accretion systems.
MAGNETOHYDRODYNAMIC SIMULATIONS OF BLACK HOLE
ACCRETION
by
Mark J. Avara
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2017
Advisory Committee:
Professor Jonathan McKinney, Advisor
Professor Christopher Reynolds, Chair
Professor David Levermore
Professor M. Coleman Miller
Professor Richard Mushotzky
Dr. Jeremy Schnittman
c© Copyright by
Mark James Avara
2017
Preface
The body of this thesis contained in Chapters 2-5 are comprised of four papers
I have written or substantially contributed to. Chapter 2 is published as Avara,
McKinney, and Reynolds (2016) and I carried out the writing of the paper and the
bulk of all work entailed in the research presented.
Chapter 3 is Morales, McKinney, and Avara (2017), currently in preparation
for submission to Monthly Notices of the Royal Astronomical Society. While I will
not be first author on that paper, the version presented in this thesis is presented in
predominately my own words, and I have contributed substantially to the analysis.
Many of the interpretations presented here are my own contribution, and I made all
of the figures aside from Fig. 3.4.
Chapter 4 is published as McKinney, Dai, and Avara (2015b). The paper
is included here in its published form, with only minor revisions. I aided in the
interpretation of the results of this work and made some contributions to the written
text.
Chapter 5 is published as Avara, Reynolds, and Bogdanovic´ (2013). I carried
out the writing of the paper and the bulk of all work entailed in this research as
well, under the supervision of Prof. Christopher Reynolds with assistance by Dr.
Tamara Bogdanovic´.
ii
Dedication
For my parents, to whom I owe opportunity, support, and my thirst for knowledge.
‘I may not always have the answer, but I know where to look to find it.’
- James E. Avara
iii
Acknowledgments
Without the generous support of many, I would not have reached my goal of
obtaining a doctorate degree in astrophysics. The tireless dedication of my parents
and many talented teachers, and the support and strength of friends and loved ones
has been essential, and I am eternally grateful.
I would like to specifically thank some of the teachers that helped me achieve
important milestones without which I would not have gotten here: Mrs. Allen,
who brought joy to mathematics and learning in general, Mr. Wightman, who’s
portrayal of history and literature made me a better writer and deeper thinker,
and who’s recommendation help get me into TAMS. At TAMS I owe a growing
appreciation and understanding of mathematics and the sciences to Neal Brand and
Paul Braterman, especially, though many others had an impact on my scientific
understanding as well. At MIT I thank those tireless in the art of teaching and
physics, especially Max Tegmark, Peter Dourmashkin, and David Kaiser. It is there
that I first learned what scientific research entails under the instruction of Herman
Marshall.
Without the opportunity of attending graduate school at UMD this thesis
would never have come to pass, and for that opportunity I thank those on the
admissions committee, especially Lee Mundy, and those scientists who recommended
me, including Belinda Wilkes, who generously offered me a place in her group when
I was trying to get into a good graduate program.
It is at UMD that I have grown mature as a scientist, though there is still a
long road ahead. I have been lucky to work under the wise tutelage of Christopher
Reynolds and Jonathan McKinney. Working with two scientists of such depth of
understanding and appreciation for their work has been a worthwhile challenge, and
a true pleasure.
These are the shoulders I stand on, and I hope to make a contribution to
science worthy of the time and devotion of those who have trained me.
I would also like to thank Manuela Campanelli for offering me an opportunity
to continue research at RIT as a post-doc, and for being flexible during a difficult
iv
time.
Finally, I especially owe so much to my parents. JoAnn Avara and James E.
Avara helped me find my path, time after time, with support and love beyond what
I could have asked for. I wish he’d grown up with the opportunities I’ve been given,
but in place of that, both my parents, have helped me achieve my dreams.
v
Table of Contents
List of Tables viii
List of Figures ix
List of Abbreviations xvii
1 Background and Motivation 1
1.1 Accretion in astrophysics . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Observations of AGN and XRBs . . . . . . . . . . . . . . . . . 3
1.1.2 Geometrically thin, radiatively efficient accretion . . . . . . . 9
1.2 Magnetic fields in disks . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.1 Magnetorotational instability . . . . . . . . . . . . . . . . . . 12
1.2.2 Magnetic topology and field strength . . . . . . . . . . . . . . 14
1.3 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3.1 Finite volume approach . . . . . . . . . . . . . . . . . . . . . 16
1.3.2 HARM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.3.3 Radiative GRMHD . . . . . . . . . . . . . . . . . . . . . . . . 20
1.4 Outline - Description of chapters . . . . . . . . . . . . . . . . . . . . 21
2 Efficiency of Thin MADs 23
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2 Fully 3D GRMHD Thin MAD Model . . . . . . . . . . . . . . . . . . 25
2.2.1 Physical Model . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.2.2 Numerical Grid and Density Floors . . . . . . . . . . . . . . . 31
2.2.3 Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.4 Diagnostics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.1 Demonstration of flux accumulation to MAD state . . . . . . 41
2.3.2 Long-term evolution of MAD state . . . . . . . . . . . . . . . 42
2.3.3 Radiative Efficiency . . . . . . . . . . . . . . . . . . . . . . . . 45
2.3.4 Jet Efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.3.5 Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
vi
2.3.6 Power-law Fits for Radial Dependence . . . . . . . . . . . . . 54
2.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . 57
3 3D GRRMHD Simulations of a Thin Accretion Disk in the MAD state 59
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.2 3D GRRMHD thin disk model . . . . . . . . . . . . . . . . . . . . . . 62
3.2.1 Numerical grid . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.2.2 Diagnostics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.3.1 Accretion rate . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.3.2 Scale Height . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.3.3 Efficiencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.3.4 Magnetic flux . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3.3.5 Luminosities . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.3.6 Resolution effects . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.4 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 81
4 Efficiency of Super-Eddington MADs 83
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.2 Fully 3D GR Radiative MHD MAD Model . . . . . . . . . . . . . . . 85
4.2.1 Opacities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.2.2 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.2.3 Numerical Grid and Density Floors . . . . . . . . . . . . . . . 88
4.2.4 Diagnostics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5 Role of Magnetic Field Strength and Numerical Resolution in Simulations of
the Heat-flux Driven Buoyancy Instability 99
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.2 MHD of Weakly Magnetized Plasmas . . . . . . . . . . . . . . . . . . 102
5.2.1 Background Equilibrium Conditions . . . . . . . . . . . . . . . 104
5.2.2 Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.3 Numerical Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.4 Plasma Behavior as a Function of Resolution and β . . . . . . . . . . 108
5.4.1 Plasma β . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.4.2 Grid Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.5 Dynamics and Persistence of Filaments . . . . . . . . . . . . . . . . . 115
5.6 3-D Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
5.7 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 122
6 Conclusions and a look forward 126
vii
List of Tables
3.1 Grid resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.2 Polar grid Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . 67
viii
List of Figures
1.1 This figure reproduced from Remillard and McClintock (2006) shows
in red dots the instantaneous measure of hardness and luminosity
of a hypothetical LMXRB. The disk spends most of its time in the
quiescent state on the lower right corner. During an outburst it goes
through disk states i-iv. The diagrams around the plot show the
dominate components of the accretion/jet inflow/outflow during each
accretion state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1 Snapshots of the initially sub-MAD (MADfHR - left panel) and nearly-
MAD (MADiHR - right panel) models at t=0, showing ρ in color
(with legend) with magnetic field lines (black lines). In MADfHR
there is no magnetic field initially present within the ISCO or on the
horizon and the field is weak througout the disk. On the other hand,
in MADiHR, a sufficient magnetic flux is present throughout the disk
and on the horizon so that the disk will quickly become MAD out to
large radius. The higher density surrounding the horizon in MADiHR
at t=0 is the result of the density floor activated due to the presence
of an initially strong magnetic field. . . . . . . . . . . . . . . . . . . . 26
2.2 (Caption on next page.) . . . . . . . . . . . . . . . . . . . . . . . . . 27
ix
2.2 (Previous page.) Evolved snapshot of the MADfHR model at t ≈
70, 000rg/c showing log of rest-mass density in color (with legend)
in both the z − x plane at y = 0 (top-left panel) and y − x plane
at z = 0 (top-right panel). Black lines with arrows trace field lines,
where thicker black lines show where field is lightly mass-loaded. The
bottom panel has 3 subpanels. The top subpanel shows the mass
accretion rate through the BH (M˙H). The middle subpanel shows
the magnetic flux threading the BH (ΥH), normalized so that order
unity is a dynamically-substantial amount of magnetic flux. The
bottom subpanel shows the gas efficiency measured on the horizon
(ηH,g) and jet efficiency measured at r = 50rg (ηj). In this model, the
disk initially has a weak large-scale poloidal magnetic field, which is
transported by the MRI and other magnetic stresses toward the BH.
Eventually, by t ∼ 55000rg/c, the magnetic flux saturates on the BH
and cannot grow despite plentiful magnetic flux outside. The MAD
state that develops consists of a highly non-axisymmetric clumpy thin
disk with a weak (ηj ∼ 1%) jet. . . . . . . . . . . . . . . . . . . . . . 28
2.3 (Previous page.) Plotted from left to right are snapshots from MADfHR
at t=[0,10000,16320,40000]rg/c. Black contours (with a single colored
white) trace the same values on all snapshots and show the φ-averaged
magnetic field. The contours demonstrate the transport of magnetic
flux radially toward the horizon. All magnetic flux is initially thread-
ing the disk only, then accretes most rapidly through coronal pro-
cesses along the surface of the disks where it piles up at these latitudes
before it does in the mid-plane of the accretion flow. Reconnection
near the horizon where the disk is very thin leads to magnetic loops
(seen in middle panel) in the inner disk out to ∼ 20rg, where fur-
ther reconnection and decoupling from the magnetically driven wind
cause a drop in the effective disk viscosity. Ultimately the magnetic
loops are also accreted and some flux originally threading the disk
then threads only the horizon. . . . . . . . . . . . . . . . . . . . . . . 29
2.4 (Caption on next page.) . . . . . . . . . . . . . . . . . . . . . . . . . 37
x
2.4 (Figure on previous page.) Time and φ-averaged quantities are plot-
ted for MADfHR. The time average is taken over 55000-68000 rg/c,
during which time the horizon has reached a MAD saturation of
flux. Rest mass gas flow is traced by black streamlines. The colored
(green and blue) thick lines correspond to time-averages of quantities
Y = {ur, b2/ρ}, respectively, at values [Y ]t = {0, 1}. While near the
BH the flow has [b2/ρ]t & 1 as averaged directly, the dense inflow
has [b2/ρ]t . 1 at all other radii. The time and spatial variability
of the jet causes the time-φ average to blend jet and disk material,
so the average [b2/ρ]t & 1 bleeds into the disk. Similarly, [b2]t/[ρ]t
identifies a jet boundary in the opposite way, causing too much of
the polar/jet region to appear to have [b2/ρ]t . 1. No particular
average is correct, but our choice includes the part of the wind and
disk that experiences extra acceleration through interaction with the
jet. Streamlines nearest the polar axes are calculated using veloc-
ity structure five cells away from the axis for better visualization of
the qualitative geometry of the jet flow. The flow exhibits equatorial
asymmetry over many inflow times as found in thicker MAD simu-
lations, likely a result of bulk flows associated with the very initial
evolution of the disk. The red magnetic field line contours show how
the magnetic flux threads the black hole and inner disk. . . . . . . . . 38
2.5 Time (from 30000-70000rg/c) and φ-averaged quantities are plotted
for MADiHR. Identical layout as Fig. 2.4. The disk is magnetically
arrested out to a transition radius of ∼ 35rg, which leads to a short
inflow time-scale compared to the MADfHR disk at those radii, so
the time-averaged velocity fields and other disk structure is smoother
and more symmetric. . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.6 Evolved snapshot of the MADiHR model at t ≈ 70, 000rg/c, with
identical layout to Fig. 2.2. The MAD state is highly non-axisymmetric
and clumpy with a weak (ηj ∼ 1%) jet. The rest-mass accretion rate
and magnetic flux on the horizon reach the MAD state quickly, which
was as desired so that we can study a quasi-steady MAD out to large
radii over long time periods. The dimensionless magnetic flux ΥH on
the horizon is in agreement with our MADfHR model, showing a disk
initial condition with a relatively stronger poloidal field strength still
reaches the same MAD state near the BH. . . . . . . . . . . . . . . . 40
2.7 (Caption on next page.) . . . . . . . . . . . . . . . . . . . . . . . . . 43
xi
2.7 (Figure on previous page.) Time-averaged quantities, averaged over
10000-70000rg/c, as a function of radius. From top to bottom, panels
are: Total mass accretion rate (M˙0, blue line), energy efficiency (η),
specific angular momentum (j), magnetic flux (Ψ), and disk thick-
ness (H/R, n=1 in Eq. 2.4, blue line). Plots for η and j show total
gas+radiative (black solid), electromagnetic (green), matter (red),
and radiation (black dashed) components. Remaining plots show
their initial value (magenta dot-dash). The plot of Ψ shows the cu-
mulative net vertical flux including BH and disk contributions (blue).
The constancy of M˙ , η, j show where the flow has achieved inflow
equilibrium. Ψ shows the MAD to non-MAD transition at r ∼ 35rg,
while throughout the disk H/R ≈ 0.1 (bound gas), as maintained
by cooling. The primary result of this plot is the efficiency panel
showing the radiative contribution of ηr ∼ 15% with most of emission
beyond the standard NT disk model coming from quite close to the
photon orbit where the MAD state differs significantly in structure
and behavior from the NT disk. . . . . . . . . . . . . . . . . . . . . . 44
2.8 Shows temporal evolution of a few quantities for our MADiHR model.
The disk scale height H/R is shown in the top panel. We show the
temporally-smoothed H/R as plotted for radii RISCO (solid-black),
15rg (short-dashed), 30rg (dots), and 60rg (long-dashed). The red
line is the un-smoothed evolution at 15rg to demonstrate the level of
variability. The H/R value is regulated by cooling to be H/R ∼ 0.1,
which is what is achieved for radii just outside the ISCO and out to
the inflow equilibrium radius of r ∼ 35rg. The middle panel shows
the αb viscosity parameter spatially-averaged over the disk at a radius
10rg. The relatively high αb ∼ 0.5 is indicative of MADs, implying
efficient transport of mass and angular momentum. The third panel
shows the ratio of the total magnetic flux threading the BH to the
total magnetic flux available at the start of the simulation inside the
radius out to which the simulation ultimately reaches at least one
inflow time. This shows that the horizon retains only a fraction of
the flux available during the course of the simulation, so the quasi-
steady state has a horizon-threading flux that has reached saturation. 47
2.9 Shows temporal evolution of same quantities as Figure 2.8 for our
MADfHR model. The disk scale height H/R is shown in the top
panel. We show the temporally-smoothed H/R as plotted for radii
RISCO (solid-black), 15rg (short-dashed), 30rg (dots), and 60rg (long-
dashed). The red line is the un-smoothed evolution at 10rg to demon-
strate the level of variability and show the evolution of the under-
resolved dense torus present at early times. αb ∼ 0.5 averaged around
8rg drops to low values during that time as expected. . . . . . . . . . 48
xii
2.10 Plotted as in the top frames of Fig. 2.6 is a snapshot of MADiHR
at t = 31568rg/c when the disk experiences an extreme eruption
of magnetic flux from the BH horizon. A large magnetic flux tube
creates a low-density cavity in the disk before it is partially blended
with disk material by the magnetic Rayleigh-Taylor instability. Such
prominences occur frequently in the MAD state. . . . . . . . . . . . . 52
2.11 MRI quality factors indicating time-φ averaged resolution, viscosity
of the disk, and vertical MRI half-wavelengths per disk scale height
(Sd, blue line) are plotted as a function of radius. Q factors indicated
the number of grid cells per critical MRI wavelength, averaged with
weight
√
ρb2 (solid lines) and ρ (short dashed) across the disk where
b2/ρ < 0.5. The third panel shows the viscosity parameter αb (solid
line) and the effective viscosity calculated via the mass accretion rate
(short dashed). The magnetic evolution of the disk is well resolved
and leads to a high viscosity as measured by local stresses, αb, and
the rate of mass inflow αb,eff . The magenta dot-dash line represents
the value of Sd at t=0 and the time-average is blue-solid. . . . . . . 55
2.12 Several quantities characterizing the accretion flow in MADiHR are
plotted as a function of radius. The top panel shows the density
profile ρ (solid), the internal gas energy ug (dashed), and magnetic
energy density ub (dotted). The inner disk, that has reached the
MAD state, has ug ∼ ub indicating equipartition is achieved even
within the dense parts of the disk. The middle panel includes the
disk radial velocity, azimuthal velocity, and Keplerian velocity (vr,
solid; vφ, black-dashed; vK,blue-dashed, respectively). The final panel
shows the radial profiles of each component of the comoving magnetic
field, br,θ,φ, in solid, dashed, and dotted lines respectively. The density
is shallower than thicker disk solutions, the disk is fairly Keplerian,
and the magnetic field falls off as expected for accretion disks. . . . . 56
3.1 Snapshot of our simulation RADHR at the time t=40.000Rg/c where
the top two panels represents the X −Z at y=0 and X − Y at z = 0
slices of the density with the red line showing where b2/ρ=1 and the
black lines with arrows trace field lines where thicker black lines show
where field is lightly mass loaded. The top sub-panel shows the mass
accretion rate at the horizon (red line) and the luminosity (green
line). The middle sub-panel shows the magnetic flux in the horizon
normalized so that order unity is a dynamically substantial amount of
magnetic flux. The bottom sub-panels shows the accretion efficiencies
in the horizon (red line), jet efficiency (green line) measured at 15Rg
and radiation (blue line) measure at 50Rg. . . . . . . . . . . . . . . . 70
3.2 (Caption on next page.) . . . . . . . . . . . . . . . . . . . . . . . . . 71
xiii
3.2 (Figure on previous page.) Characteristic quantities as a function
of radius, time-averaged over the period t =[30,000-40,000]Rg/c for
the simulation RADHR. Panels consecutively show, from top to bot-
tom: the mass accretion rate (M˙), energy efficiencies (η) for different
parts of the stress-energy tensor, magnetic flux (Ψ), disk scale-height
(H/R), and Sd. In efficiencies plot the black line corresponds to the
total gas+radiative (black line), electromagnetic (green line), matter
(red) and radiation (black dashed) components. The magenta lines
in the plots of Ψ, H/R, and Sd correspond to their initial values. The
thick lines show values for RADHR and the thin line MADiHR from
(Avara et al., 2016) time-averaged as in that paper. . . . . . . . . . . 72
3.3 Red is for RADLR, black for RADHR. Top panel: Disk scale-height
as function of time. The solid line corresponds to the values for 10Rg
(lines for all radii are smoothed, with only the 10Rg line duplicated
without smoothing in gray), the dotted and dashed lines correspond
to the values of H/R at the horizon and 20Rg respectively. Middle:
Magnetic α measured at 10Rg (full, un-smoothed RADHR data in
gray). Bottom: Amount of magnetic flux on the horizon normalized
by the amount of magnetic flux contained in the portion of the disk
which reaches one inflow time by the end of the simulation, i.e. the
total magnetic flux available for the inner disk and horizon. . . . . . . 77
3.4 Top: Electromagnetic luminosity per unit angle (∂θLEM)/(M˙Hc
2)
with time-φ averaged magnetic flux (translucent gray lines); blue lines
show the b2/ρ=1 boundary; the green line shows the ur=0 boundary;
yellow is the surface where optical depth reaches unity to electron
scattering, integrating inward from a large, but at least moderately
evolved, radius. Bottom: lab-frame radiation flux and radiation lu-
minosity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.5 Plot showing the efficiencies for the simulation RADHR (dot-dashed
line), RADvHR (dashed line), RADLR (dotted line) and MADiHR
(continuous line). Black: total, Red: PAKE, Green: EM, Cyan:
Radiation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.1 The initial (t = 0) state consists of a weakly magnetized radially-
extended vertically-Gaussian disk with H/R ∼ 0.2 around a spinning
(a/M = 0.8) BH. Rest-mass density is shown as color with legend,
while green lines show magnetic field lines. Density values are made
dimensionless using M˙Edd, length rg, and time rg/c. The initial disk
is threaded by weak (but ordered) magnetic flux, capable of accumu-
lating onto the BH and leading to the MAD state. . . . . . . . . . . . 91
xiv
4.2 The evolved (t = 31, 200rg/c) super-Eddington MAD state, with the
figure split at x = 0 as two panels (x < 0 on left and x > 0 on
right). Left panel shows radiation-frame radiation energy density
(color, with legend), radiation velocity lines (black, fixed line width),
and optical depth of unity away from each polar axis, τθ = 1 (cyan
lines). Right panel shows fluid-frame rest-mass density (color, same
legend), magnetic field lines (black, thicker lines for more magnetized
gas), where magnetic energy is equal to rest-mass energy density (red
lines), and same optical depth (cyan lines). The MAD state reaches a
quasi-steady state out to r ∼ 20rg. Radiation advects inward within
the equatorial disk, outward through the low-density jet channel, and
outward along with the optically thick wind. . . . . . . . . . . . . . . 92
4.3 The time-averaged angle-integrated mass flux, efficiency, disk thick-
ness, and magnetic flux. From top to bottom, panels are: Total mass
accretion rate (M˙0) per unit M˙Edd out to where there is inflow and
thermal equilibrium (r ∼ 20rg, beyond which the line is truncated),
energy efficiency η (total: solid black line, electromagnetic: blue
dashed, matter without rest-mass (i.e. kinetic+gravitational+thermal):
dark red dotted, radiation (negative within 3rg): red long-dashed,
radiative luminosity L: magenta dot-dashed), disk thickness H/R
(evolved state: black solid, initial state: blue dashed), and magnetic
flux (total on BH and threading disk: black solid, only threading disk:
dark red dotted). The disk is magnetically compressed near the BH,
and the thin disk generates a significant radiative flux that effectively
releases starting at r ∼ 300rg. The initial disk with H/R ≈ 0.2 has
puffed up to H/R ≈ 0.3 from r ∼ 20–70rg after evolving for several
thermal times. By r ∼ 400rg, the luminosity reaches L ∼ 50LEdd. . . 93
4.4 The evolved (t = 31, 200rg/c) super-Eddington MAD state at large
radii showing the jet, wind, and photosphere. Shows fluid-frame rest-
mass density (color, with legend), where magnetic energy equal to
rest-mass energy density (red lines), radiative photosphere at τr = 1
(black line), and the wide-angle wind component’s magnetic field line,
within which (toward the polar axis) the outflow has achieved more
than a single flow time along the line extending ∼ 1000rg (blue line).
Quantities are φ-averaged. The jet extends to r ∼ 30, 000rg with
Lorentz factor γ ∼ 5 by r ∼ 103rg. The wind extends to r ∼ 3000rg,
with most of the wind’s trapped radiation having an opening angle
of 30◦ at r ∼ 400rg. The disk at R  100rg has not evolved much.
Radiation within r . 400rg can more freely stream within the jet’s
low-density channel. By t ∼ 30, 000rg/c, at r ∼ 400rg within 15◦
around the polar axis, the outflow become conical and optically thin. 95
5.1 H2d128 7 and H2d128 9 at four times (in units of dynamical times,
ω−1dyn) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.2 Horizontally averaged vertical heat flux as a function of time. . . . . 110
xv
5.3 The contribution to kinetic energy from vertical motion spatially av-
eraged across the entire anisotropically conducting region is shown as
a function of time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.4 Evolution of energy components. . . . . . . . . . . . . . . . . . . . . 113
5.5 Horizontally averaged vertical heat flux plotted in the β ∼ 107 column
of Fig. 5.2 versus time in dynamical units. . . . . . . . . . . . . . . . 115
5.6 H2d512 5 (upper) and H2d512 7 (lower) filaments. Stability of the
filaments is demonstrated by plotting the three stability diagnostics
of Eqn. (5.19) measured for each filament. . . . . . . . . . . . . . . . 117
5.7 Black contours (dash-dot) represent filaments from H2d128 7PS; red
(solid) represent those from H2d128 7. Contour levels are set the
same as in Fig. 5.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
5.8 A frame from the saturated state of H3d128 7PS. . . . . . . . . . . . 122
5.9 Vertical Heat Flux integrated over the horizontal plane slicing through
the center of the active domain in H3d128 7PS. . . . . . . . . . . . . 123
xvi
List of Abbreviations
BH Black hole
MHD Magnetohydrodynamics
GRRMHD General relativistic - Radiation - MHD
AGN Active galactic nucleus
XRB X-ray binary
MRI Magneto-rotational instability
HBI Heat-flux driven buoyancy instability
xvii
Chapter 1: Background and Motivation
1.1 Accretion in astrophysics
Astrophysics, as a mature field of science, truly encompasses the full breadth
of physical law, from quantum theory to the extremes of gravitation. It generously
provides environments inaccessible in the laboratory but necessary in order to test
many modern physical theories, and discover new ones. Although much about the
formation and behavior of galaxies, stars, planets, and even black holes can be
understood through the complex coupling of predominantly known laws of physics,
it is by seeking an understanding of that complexity that we hope to use astronomical
observations to gain knowledge of new physical theories.
Accretion is the process a gravitationally self-bound object undergoes, whether
it be a gas cloud, a star, a galaxy, a large-scale cosmological structure, or other
object, as it accumulates more mass by gravitational attraction of the gas, or plasma
if it has become ionized. The accumulating mass generally carries with it angular
momentum leading to rotation as it contracts toward the central accretor.
There are two fundamental modes of accretion, separated by whether the
specific angular momentum of the accreting matter (the matter can be composed
of dust grains, molecules, a mixture of elements or ionization fractions, etc., but
will hereon be referred to simply as ‘gas’) defines a plunging or a circulating orbit.
Neglecting gas pressure forces, etc., for sufficiently high angular momentum, matter
in an accretion inflow will reach a minimum radius where the centrifugal force equals
the gravitational attractive force before reaching the surface of the central mass. A
gas flow with this specific angular momentum or more is thus centrifugally choked,
and as gas streams on different crossing orbits interact they form a disk.
1
On the other hand, the lower angular momentum case describes Bondi accre-
tion (Bondi, 1952), the zero angular momentum limit of which provides the simplest
model for accretion of gas onto a point-like gravitating object. The Bondi accretion
rate onto the central object, M˙ only depends on the mass of the central object and
the sound speed of the surrounding gas reservoir, and is determined by a critical
point condition. Because of the simplicity of the Bondi model for accretion from
homogeneous and rotationally invariant non-self gravitating gas, it is widely used
to estimate the growth rate of gravitationally bound structures, such as supermas-
sive BHs (Sakurai et al., 2016; Inayoshi et al., 2016), below the resolving power of
simulations of galaxies and stars.
The disadvantage of the Bondi model, especially in its application to compact
objects like black holes and neutron stars, where the mass to surface radius ratio
is very high, is that most of these systems in the universe are accreting gas with
significant angular momentum. Therefore for accretion to occur, the system must
transport angular momentum away from the gas. Any process that transports the
angular momentum outward dissipates energy and the mode of accretion thereby
branches again depending on whether the accreting gas is efficient or inefficient at
radiating the released energy. The outward transport of angular momentum and
the energy released by that process give rise to a variety of theoretical models of
accretion which we believe apply to the zoo of seemingly different accreting systems
astronomers observe. A second limitation of the Bondi model is that it fails to take
into account radiative transport (McKinney et al., 2014).
The total energy liberated as matter falls from an effectively infinite radius to
the central object is often written as a conversion of mass flux to energy flux
L = ηM˙c2. (1.1)
The power of stars originates from the fusion of hydrogen into helium, a process
with η ≈ 0.007. In contrast, as I will describe below, a maximally rotating accreting
black hole (BH) can reach η ≈ 40% even using classical accretion models! This
extreme efficiency in liberating gravitational potential energy makes some accreting
2
BHs the most luminous objects in the universe in electromagnetic radiation.
With the advent of X-ray telescopes in the 1960’s, the discovery of many
highly energetic and luminous astronomical objects summoned a renewed interested
in astrophysical manifestations of extreme solutions to Einstein’s field equations.
Because of their compact nature, non-rotating (Schwarzschild, 1916) and rotating
(Kerr and Schild, 1965; see also Kerr and Schild, 2009) BHs, time-independent
solutions to the field equations, can host a disk reaching far enough into the gravita-
tional potential that the opportunity for conversion of gravitational potential energy
to other forms is maximized. Accreting black holes, due to their high luminosity and
compact extent (most energy is liberated close to the BH), are an obvious choice to
explain the high efficiency of the unresolved massive objects in active galactic nuclei
(AGN) and X-ray binaries (XRBs).
1.1.1 Observations of AGN and XRBs
The most luminous AGN are historically referred to as quasars, or quasi-
stellar-objects. Seyfert (1943) was the first to characterize galaxies that contained
stellar-like cores using their optical spectra, and soon after the detection of Seyfert
galaxies NGC 1068 and NGC 1275 as radio sources in 1955, Woltjer (1959) tried to
understand the physics of these unusual galactic centers. He used virial arguments to
infer large central masses and constrained the unresolved nuclei to less than 100pc.
Thus something was required to release a huge amount of energy in a tiny volume. It
did not take long before BHs were suggested as the central accretor (e.g., Zeldovich,
1964).
Accreting supermassive BHs (SMBH), i.e., > 106 solar masses (M), are now
believed to lie in the center of nearly all galaxies at low redshift. A major argument
for this is that of Soltan (1982). Even assuming a moderate radiative efficiency of
accretion of %10, the observed characteristics of typical quasars leads to a simple
argument for the average accreted central mass per galaxy (for those on the bright
end of the galaxy luminosity function, i.e., those capable of supporting a massive
central object) of greater than 107M.
3
These massive AGN exhibit a continuum of emission across the electromag-
netic spectrum, from radio to γ-rays (Peterson 1997). Observations of AGN across
these wavelengths has led to unified models centered around the principle of ac-
cretion onto a central SMBH, where the variety of AGN spectra and morphologies
observed are determined by only a few simple parameters. The first unified model
centered around inclination of the disk to the line of sight, and a more modern pic-
ture typically includes spin of the central BH and the mass accretion rate. These
newer additions to the unified model, including the most recent, large-scale mag-
netic field presence in the galactic environment of the AGN, couple the BH more
to the history of accretion and the effective boundary condition that is the hosting
galactic environment.
The mass accretion rate onto BHs is typically scaled by the so-called Eddington
rate,
M˙E =
LE
ηc2
≈ 1.4× 1018(M/M)gs−1 (1.2)
where a typical (radiation) efficiency factor of η ≈ 0.1 is taken and LE is the Edding-
ton luminosity, the maximum luminosity allowed by spherically symmetric accretion
of hydrogen gas. At this accretion rate, and assuming solar gas composition (i.e.,
metallicity), the outward radiation pressure equals the gravitational force acting on
the accreting matter. Of course, a break of the symmetry of the system and other
realities of BH accretion have called into question just how limiting LE really is.
Nevertheless, it is a useful parameter for characterizing different accretion modes.
Most broadly, there are Low-Luminosity AGN (LLAGN), where M˙/M˙E is
smallest, Seyfert galaxies, typically undergoing thin-disk accretion with intermediate
M˙ , and quasars, the brightest sources which may exhibit accretion rates above M˙E.
Quasars are active galaxies where the central source outshines the galaxy and they
are typically high-redshift objects. Seyfert galaxies fall into two categories, Seyfert
Is with a line of sight to the central accretion disk, and Seyfert IIs, where the central
disk is obscured by the outermost part of the disk, known as the ‘dusty torus’. These
AGN are further characterized by broad and narrow emission line regions, varying
levels of absorption, presence or absence of highly collimated relativistic outflows
4
called jets, varying radio emission rates and morphologies potentially linked to jet
power, and lower velocity outflows referred to as winds.
Though the process is far from fully understood, there is ample observational
evidence (see Fabian, 1994, for a review) that the high radiative output and espe-
cially the high-velocity outflows can influence the host galaxies or galaxy clusters
of AGN. Powerful radio emission is detected to emanate from the narrow body and
base of many jets, but they can be especially luminous at the ‘radio lobe’ struc-
tures where highly collimated relativistic jets shock and couple to the intra-cluster
medium in some galaxy clusters.
In so-called ‘cool-core’ galaxy clusters, X-ray brightness measured from the
central regions implies a cooling time-scale much shorter than what would be con-
sistent with observations. Early explanations invoked conduction of heat into the
cores from the large heat reservoir of the surrounding intra-cluster medium (ICM)
(Zakamska and Narayan, 2003, etc). However, without some sort of self-regulation
of the conduction by turbulence or another process, a conduction-based explanation
constitutes a fine-tuning problem.
Arguably a more natural explanation for this cooling flow problem is that the
AGN that are clearly interacting with the ICM via their jets are heating it either
by driving weak shocks, driving turbulence (Yang and Reynolds, 2016), or even
by introducing cosmic rays (Ruszkowski et al., 2016). Clearly there are significant
open questions in how AGN can use jets or winds to transport energy to their
large-scale environment, but as we will show, there is still much work to be done
in understanding how much energy an AGN can convert into outflows and what
determines its typical accretion mode.
As already mentioned, the large variety of AGN observed is likely due in part to
the different properties of the gas funneled to the central SMBHs by the surrounding
galactic environments. Complementary to the study of AGN is their low-mass (a
few to 10s of solar masses) counterpart, the XRB (see Remillard and McClintock,
2006, for a review). In an XRB, a BH or neutron star accretes through Roche Lobe
overflow (low mass XRB) or stellar wind (high mass XRB where the star is high
5
enough in mass that it undergoes significant mass-loss via. wind) from an orbiting
companion that is thought to provide boundary conditions for the accreting system
with fewer free parameters.
Near-Keplerian angular momentum in the material feeding onto the disk from
the companion star likely means that these systems are always disk-like (see Begel-
man, 2014), and in most cases the spin of the BH is expected to be in near alignment
with the disk and larger angular momentum of the binary (Fragos et al., 2010). On
the other hand, SMBHs could be misaligned with the net angular momentum of the
accreting gas clouds raining from the surrounding ISM, effectively introducing a free
parameter to the system which has relatively little numerical exploration thus far.
For nearly all XRBs M˙ << M˙E most of the time, and only approaches a
significant fraction of M˙E during ‘outburst’ events. At present, only two systems
are observed to likely exhibit super-Eddington states: SS 433, and GRS 1915. They
spend most of their time in the ‘low’ or quiescent states where the disk is expected
to be a radiatively inefficient accretion flow (RIAF). In such flows the density, and
hence collisionality, is too low for most of the energy converted during accretion to
radiate away. One type of analytic RIAF commonly applied to model the low/hard
and quiescent states is the advection dominated accretion flow (ADAF; Narayan and
Yi, 1994; Abramowicz et al., 1995), where heat is advected across the event horizon
with the gas rather than escaping the system. However, on timescales humans
can observe, sections of the disks in low-mass XRBs (LMXRB) are believed to fall
subject to an instability caused by varying hydrogen ionization levels, which causes
periodic bursts in the accretion rate (see Cannizzo, 1993, for a review). Many low-
mass BH XRBs seem to at least roughly follow a specific type of cyclical behavior
through one of the outburst events illustrated in Figure 1.1.
Starting in a low state of accretion where M˙ << M˙E, XRBs are referred to
as ‘quiescent’ if they are extremely faint, which may be the limiting case of the
low/hard state. The accretion flow of a RIAF, is hot because the heat liberated by
accretion processes is trapped in the low-density gas. This high temperature causes
gas pressure to support a geometrically thick, optically thin disk. Thick disks are
6
jet
 lin
e
HS LS
VHS/IS
Soft           Hard
Γ > 2
i
ii
iii
jet
Je
t L
or
en
tz
 fa
ct
or
iiiiv
iv
iii
in
te
ns
ity
hardnessX−
ray
D
is
c 
in
ne
r r
ad
iu
s
no
Figure 1.1: This figure reproduced from Remillard and McClintock (2006) shows
in red dots the instantaneous measure of hardness and luminosity of a hypothetical
LMXRB. The disk spends most of its time in the quiescent state on the lower right
corner. During an outburst it goes through disk states i-iv. The diagrams around
the plot show the dominate components of the accretion/jet inflow/outflow during
each accretion state.
7
thought to be good at collimating outflows either driven by accretion processes in
the inner disk or accelerated by the BZ process near the horizon. Power from the
latter has long been thought to be amplified by the thick disk because it is effective
at advecting poloidal magnetic flux with the flow, increasing flux ‘threading’ the
horizon.
Most emission in the low state is attributed to non-thermal processes and is
likely dominated by inverse Compton scattering of lower energy thermal photons
by high-velocity electrons in a corona. This emission has a strong component that
extends into the hard part of the X-ray spectrum, 10-100keV. That, along with low
X-ray flux from any thermally dominated disk component in the soft X-ray band of
0.1-10keV leads to a naming convention of ‘low/hard’ for this disk state. Hardness
is defined as the ratio of the hard spectral component to the soft.
During an outburst event as M˙ rises, the thermal and non-thermal compo-
nents of the spectrum increase along with the overall luminosity of the system. The
hardness of the spectrum remains roughly constant as the flux increases by as much
as several orders of magnitude. The epitome of this is called the very high state,
where there is, temporarily, a power-law component dominating across most of the
X-ray and even up to γ-ray energies in the MeV. Physically, this is usually inter-
preted as a high-M˙ intermediate state leading towards geometrical collapse of the
disk as it becomes radiative efficient, but there is still the presence of a hot corona.
As the corona, or disk collapses, the jet usually is thought to disappear, i.e., radio
emission associated with the jet drops significantly. The intermediate states between
the very high state and the next state are broken into two, hard and soft, as the
coronal spectral components give way to a thermal multi-temperature black body.
The multi-temperature accretion disk model is thought to be the closest state
to the standard SS disk. It is thin, so is traditionally thought to not support a jet
(though we shall demonstrate this could be incorrect) and most flux is in thermal
components of the spectrum, which is ‘soft’ rather than hard. This high/soft state
remains as M˙ drops. Once the accretion rate falls low enough there are two branches
in the optical depth-M˙ parameter space, a very thin optically thick branch and the
8
ADAF branch. Through an unknown mechanism it is believed that the thin disk
evaporates into the ADAF state. This returns the system to the low/hard state and
the jet appears to return, along with the associated radio emission.
Since the transition from hard to soft states happens at very different total
luminosities depending on the direction taken in the luminosity-hardness diagram,
there seems to be some sort of hysteresis that implies the presence of another funda-
mental parameter of the system besides M˙ and spin. Just as the magnetic field can
help to explain the complexities beyond the simplest AGN unified models, there are
several ways in which magnetic fields could play a role in explaining the hysteresis
(Begelman and Armitage, 2014) and other complicated behavior of XRBs.
Even during the high soft state there is often a strong non-thermal component
and the interplay between large-scale magnetic field and the corona of disks has
become the subject of increasing numerical investigation. How much material is
strongly magnetized? How does magnetization effect the presence of a corona in
the thin-disk state? How does it affect the present of a jet? These are questions
intimately tied into our understanding of the states of XRBs that we investigate in
this thesis.
There is a rich complexity of variability phenomena associated with the study
of XRBs that can also be used to gain insight into these accreting systems, but I
do not summarize it here because it is not the focus on any of the Chapters that
follow.
1.1.2 Geometrically thin, radiatively efficient accretion
The accretion mode associated with the high-soft state in XRB outbursts is
most simply modeled using the so-called α-disk model attributed to Shakura and
Sunyaev (1973). In their model they assume that energy dissipated by the disk is
radiated locally as a blackbody, and angular momentum transport is mediated by
an anomalously parameterized torque of the form
νT = αcsH, (1.3)
9
where cs is the sound speed and H is the scale height of the disk. Reasoning that any
hydrodynamic stress can be scaled using the total pressure P , as αP , Shakura and
Sunyaev arrived at the above formula starting with a simplified form of the stress
energy tensor in a Keplerian disk. However, one can arrive at the same equation
considering a simple physical description: consider momentum transport through
turbulent eddies communicating angular momentum between radii within the disk.
A natural scale for the eddies is the height of the disk with a natural bulk motion
moving with the speed of sound.
The anomalous viscosity model of SS was generalized to the Kerr space-time
of a spinning BH by Novikov and Thorne (1973a) (NT73). As a simplified picture,
consider the potential energy U = GMm/r of a mass m at radius r in a disk sur-
rounding a central mass M. The rate of energy conversion, from the Virial theorem,
is
L ≈ 1
2
dU
dt
=
GM
2r
dm
dt
=
GMM˙
2r
. (1.4)
This combined with Eqn. 1.1 shows that η ∝M/r. This is a measure of compactness
of the system at its terminal radius, which for a BH is naturally given by the
innermost stable circular orbit (ISCO), inside of which gas plunges into the horizon
and there is no further dissipation (in these idealized models) of energy. The total
energy available from a mass falling from infinite radius to the ISCO is
U
2
=
GMm
2 · 3Rs =
GMm
12GM/c2
≈ 0.08mc2, (1.5)
where Rs is the Schwarzschild radius. This gives an approximation to the actual
value for a non-spinning BH in NT73 that is surprisingly close to the expected value.
While the SS or NT thin disk, parameterized by a constant average viscosity
parameter at all radii, α, is an extraordinarily simple picture of accretion when the
disk material is able to efficiently radiate dissipated energy, that simplicity has made
it valuable when applied to constraining observations.
The fundamental physical assumptions in the standard disk model are the
form of the viscosity, zero stress inside the ISCO, a near-Keplerian rotation profile,
and lack of any wind or deviation from black body radiation of energy from the disk,
10
radiation of dissipated energy locally to where its dissipated, and a shared axis be-
tween the BH and the disk angular momenta. Also, there is an inherent assumption
of locality in the sense that different disk radii are only causally connected through
the viscous torque in the disk. It should be noted also that the NT and SS solutions
reveal the disk is radiation-pressure dominated close to the BH.
The focus of the next two chapters in this thesis concern the relaxation of
several of these assumptions while remaining in the optically thick, geometrically
thin accretion flow mode. Our primary concern will be with the effects of large-scale
and dynamically relevant magnetic fields.
1.2 Magnetic fields in disks
The simple models for accretion discussed so far have assumed only the unmag-
netized fluid equations for conservation of mass, angular momentum, and energy.
Before considering the effects of magnetic fields of varying strengths and topology
on accretion disks around BHs, first consider the equations of ideal MHD in conser-
vative form,
∂ρ
∂t
+∇·[ρv] = 0, (1.6)
∂ρv
∂t
+∇· [ρvv −BB+ P∗] = 0, (1.7)
∂E
∂t
+∇ · [(E + P ∗)v −B(B · v)] = 0, (1.8)
∂B
∂t
−∇× (v ×B) = 0, (1.9)
where P∗ is a diagonal tensor with components P ∗ = P +B2/2 (P is the gas thermal
pressure component), E is the total energy density
E =
P
γ − 1 +
1
2
ρv2 +
B2
2
, (1.10)
and B2 = B ·B. In the units chosen, the magnetic permeability µ = 1.
Closure is provided by the equation of state for an ideal gas, P = (γ − 1)e
(where γ is the ratio of specific heats, and e is the internal energy density). Equa-
11
tions 1.6, 1.7, and 1.8 are the non-relativistic conservation equations for rest mass,
momentum, and energy, respectively with no external sources in Minkowski (i.e.,
flat) space-time, for simplicity. We could add to these equations some source terms
on the right hand side, and will use numerical solutions of those more complex forms
of the equations of magnetohydrodynamics (MHD) in the chapters representing the
work of my thesis below. For relativistic accretion flows around BHs it is also neces-
sary to formulate the equations of MHD in a covariant way. For now, let us continue
with the simplified form.
1.2.1 Magnetorotational instability
In a landmark invocation of a property of rotational shearing flows first ex-
plored in the context of Taylor-Coulette flows by Velikhov (1959) and later again by
Chandrasekhar (1961), Balbus and Hawley (1991) showed that a weak-field MHD
instability, the magneto-rotational instability (MRI), can drive vigorous turbulence
in accretion disks. Thus, the MRI provides a natural process for self-consistent
angular momentum transport. Turbulence was thought to be necessary to trans-
port angular momentum because shear viscosity of the gas alone in disks would
provide torques many orders of magnitude too small to explain observations, if the
empirically determined α ∼ 0.1 is typical.
The MRI manifests in its simplest form for a vertical (z) (poloidal) magnetic
field with a radial (r) perturbation, where there is shear in the φ direction, such as
that found in a differentially rotating accretion disk.
Using a standard perturbation analysis on the linearized MHD equations in
a locally co-rotating frame in a disk with a Keplerian radial profile for vφ(r), one
reduces the linearized conservation equations to get the linear dispersion relation
for perturbation ~k = kzˆ with frequency ω,
ω2
Ω2
= q2 +
1
2
(
1±
√
1 + 16q2
)
, (1.11)
where the reduced wavenumber q = kvA/Ω, and vA = B/
√
4piρ is the Alfve´n speed.
k must not be smaller than ∼ 1/H, otherwise the perturbation would not fit within
12
the disk, but otherwise, for small allowed q, the growth rate is set by the imaginary
component of ω for the minus version of Eqn. 1.11, which is Im(ω) ' √3qΩ. The
maximum growth rate (fastest growing and therefore dominating MRI mode) is with
q =
√
15/16, where one e-folding happens at the frequency (3/4)Ω. This means the
MRI can amplify the magnetic field on a sub-orbital timescale! If the magnetic field
becomes strong, the Alfven velocity increases and when q ≥ √3 this mode is stable.
If all modes that fit within the scale-height H of the disk are stable, then the MRI
is quenched and some other process must drive accretion.
Physically one can understand the MRI using the following Gedanken exper-
iment: consider two parcels of gas at adjacent radially locations, connected by the
radially perturbed magnetic field oriented in the z direction as above. You can think
of the field as a spring connecting the two masses (Gammie, 2004). As the gas mass
m2 at the outer radius falls behind gas mass m1 with a faster rotational velocity,
the magnetic field connecting them deforms further in the ‘S’-like shape sinusoidal
shape of the perturbation, storing a bit of energy itself in magnetic tension, while
also transferring a bit of angular momentum and kinetic energy from m1 to m2.
With less angular momentum m1 falls to a lower orbit and m2 similarly moves out-
ward to a larger orbit, thereby further stretching the connecting field. Of course, if
the field is too strong, the magnetic tension will resist the initial radial displacement
altogether.
The non-linear saturation of the MRI and the transition to turbulence has a
rich history of exploration (Balbus and Hawley, 1991, 1998a; Hawley et al., 1995;
Stone et al., 1996; Pessah et al., 2007; Sorathia et al., 2010, 2012). With the discovery
of the MRI, simulations of accretion disks could be achieved with fully self-consistent
angular momentum transport!
As you can see, even the inclusion of weak (plasma β parameter, ratio of
magnetic to thermal pressure, is  1) magnetic fields represents a significant de-
parture from the idealized thin disk framework. Though the MRI does supply
self-consistency to the idea of high α angular momentum transport, the magneto-
dynamo action of the saturated-state MRI continually produces magnetic patches
13
which can be transported to the surface of the disk by convection and drive a disk
wind. The nature of the magnetic field can also further enhance the asymmetry of
stresses within the disk.
As simulations of BH accretion, typically using a pseudo-Newtonian approxi-
mation to the gravitational potential, developed from 2D to 3D, and from shearing
box to vertically stratified shearing box, to semi-global and ultimately to fully global
simulations, a wealth of non-linear behavior has been uncovered which may explain
the variability and complex spectral energy distributions seen in real systems. Even
with the sub-sample of arguments presented so far in this thesis, it is clear that in-
cluding the event horizon of the BH as a realistic inner boundary, and including the
full influence of the Kerr space-time is absolutely key in understanding the highest
energy emission and most rapid variability from BH systems.
I have discussed the MRI in the context of weak magnetic fields, but there is
no a priori reason the field should be weak! If the poloidal magnetic field is strong
enough for the smallest unstable mode to not fit vertically within the disk body,
then the axisymmetric MRI is quenched, and accretion proceeds through magnetic
Raleigh-Taylor interchange modes instead. This is the so-called ‘magnetically ar-
rested’ disk state that is explored in various experiments in the next three chapters.
1.2.2 Magnetic topology and field strength
The first GRMHD simulations (starting with Wilson (1977)) with full time-
independent GR of global disks have typically included weak initial magnetic fields
of either a single set of large poloidal loops tracing pressure contours of a hydro-
statically stable torus, or a series of poloidal loops with alternating field direction
to avoid introducing surfaces of artificial reconnection in the initial conditions.
These magnetic field initializations are useful in studying the magneto-dynamo
and carrying out simulations with high resolution to study disk turbulence. A lack
of large-scale magnetic field threading the disk minimizes activity very far off the
equatorial plane, so one can even save computing resources by cutting out large
polar conical zones without significantly effecting disk behavior.
14
However, many real BH systems at any given time have access to large amounts
of large-scale poloidal magnetic flux (Zamaninasab et al., 2014; Mocz and Guo, 2015)
either from the surrounding environment or from remnant structures from a prior
epoch of accretion, and so the inclusion of these large-scale fields in simulations
with focused self-consistency is key. For instance, large scale magnetic field can
be advected inward and accumulate on the horizon affecting many aspects of the
accretion inflow/outflow behavior, as we will soon discuss. Also, large-scale magnetic
fields launch powerful mass-loaded winds to high latitudes and so no conical excision
of the grid around the poles will suffice if we want to capture this part of the system.
In addition to magnetic topology, the initial field strength can significantly
affect the accretion disk as well (Penna et al., 2010a; Bai and Stone, 2013; Salvesen
et al., 2016). Strong magnetic fields can enhance magnetic torques within the disk
and affect the efficiency of magnetically driven winds and the strength of the jet
driven by the spin of the BH via the Blandford-Znajek mechanism (BZ; Blandford
and Znajek, 1977), Though simulations have yet to show if toroidally dominated
magnetic field configurations, i.e., disks with only small patches of local net vertical
flux, can lead to powerful jets, inclusion of large scale poloidal flux threading through
a disk can drive powerful jets in simulations when either the disk is very thick and
confining, or the spin of the BH is close to maximal.
Of fundamental importance is whether large scale magnetic fields can be trans-
ported to the inner disk and the horizon. If they can, the impact of magneto-
centrifugally driven winds (Blandford and Payne, 1982; Bai and Stone, 2013, etc)
on disk stability and structure, contribution of large-scale poloidal flux to jet power,
and overall accretion efficiency of BH systems, could mean that the large-scale mag-
netic field is as fundamental a free parameter governing the observed behavior of
BH systems as M˙ and inclination of the observer.
Most of the body of work presented in this thesis concerns the transport of
large-scale magnetic flux and the associated influence on accretion flows.
15
1.3 Numerical methods
There are a plethora of techniques invoked by numerical experimentalists to
solve the equations of MHD in various astrophysical contexts. I only attempt here
to provide an overview of the particular choice of techniques implemented in the
codes Athena (Stone et al., 2008) and HARM (Gammie et al., 2003; McKinney and
Blandford, 2009; McKinney et al., 2012), used to carry out the work of this thesis,
with a focus on the implementation in the latter. Teyssier (2015) and the Living
Review of Font (2008) provide nice reviews of techniques used to solve the fluid and
MHD equations in astrophysics.
1.3.1 Finite volume approach
Both Athena and HARM (and the extension HARMRAD (McKinney et al.,
2014) which includes self-consistent radiative transport) are built upon a funda-
mentally conservative approach of the finite volume form. Rather than evolving
quantities like ρ, ρ~v, ~B, etc. at discrete points and approximating fluxes, by track-
ing a set of “conserved” quantities that are volume averaged over grid cells and
accounting for transport of those quantities via fluxes across cell faces, one can pre-
cisely conserve their volume integrated totals to machine precision. It is the simple
principle that whatever conserved quantity is lost by one cell is exactly captured by
another.
While the simpler implementation of schemes that were not total-energy-
conserving, such as ZEUS (Stone and Norman, 1992), makes them somewhat easier
to write and run faster on computers, they can be less accurate, for instance, when
used to study systems where the internal energy is a significant fraction of the total.
Some non-energy-conservative schemes with artificial viscosity have trouble with
relativistic shocks (Norman and Winkler, 1986). A major advantage of conservative
schemes is that in one dimension, those methods that are total variation stable are
guaranteed to converge to a weak solution of the equations (Lax and Wendroff, 1960;
LeVeque, 1998). Though this does not strictly apply for higher spatial dimensions,
16
it is a promising starting point.
The basic method for solving a set of hyperbolic equations (like those of MHD)
upon which the higher-order schemes of Athena and HARMRAD grew is called
Godunov’s scheme, which in its simplest form is first order accurate in space and
time. It makes use of updates to conserved volume-integrate quantities using fluxes
calculated by solving the Riemann problem at cell interfaces. A Riemann problem
is an initial value problem associated with a conservation equation and piecewise
constant values, such as those on the grid considered here.
Without sources, a system of hyperbolic conservation laws can be written in
equations of the form
∂q
∂t
+∇ · Fq = 0, (1.12)
where q is one conserved quantity of the system transported via flux Fq.
Consider discretization of these quantities onto a grid (of 1 dimension for
purpose of illustration) in the x direction, where at each discrete time step the
spatial average of q(tn, x) is associated with cell center i, and is averaged over the
volume between evenly spaced cell interfaces at i− 1 and i+ 1, so that the volume
averages are
Qni =
1
∆x
∫ xi+1/2
xi−1/2
q(tn, x)dx. (1.13)
Let’s re-write Eqn. 1.12 as
qt + (F (q))x = 0, (1.14)
where the sub-scripts represent partial derivatives.
If we integrated this over the volume of the cell we obtain
∂
∂t
Qi(t) = − 1
∆x
(
F (q(t, xi+1/2))− F (q(t, xi−1/2))
)
. (1.15)
Now we time integrate through one step tn to tn+1, which yields the exact update
formula,
Qn+1i = Q
n
i −
1
∆x
∫ tn+1
tn
(
F (q(t, xi+1/2))− F (q(t, xi−1/2))
)
dt. (1.16)
17
In this low-order formulation, Godunov’s method uses the forward Euler method
to replace each of the time integrated terms so that, for example,∫ tn+1
tn
F (q(t, xi+1/2))dt ≈ ∆tf(Qni−1, Qni ), (1.17)
where f is found by solving the Riemann problem for flux at the interface with left
and right values given by the left and right states Qni−1 and Q
n
i .
Therefore, the final update for the unknown Qi from t
n to tn+1 is given by
Qn+1i = Q
n
i −
∆t
∆x
(
fni+1/2 − fni−1/2
)
. (1.18)
You can identify immediately from inspection which of the above steps in the
method might be improved to higher order accuracy. For instance, the left and right
states used to find fluxes by solving the Riemann problem at cell face i− 1/2, could
be found using an interpolation rather than the strict values Qni−1 and Q
n
i . Similarly,
a higher order time-step could be taken than the first-order Euler approximation
in Eqn. 1.17. As a final example, approximate Riemann solvers can be used in
circumstances where it is computationally expensive to iterate to an exact solution.
1.3.2 HARM
Though I will not attempt to fully cover the methodology of the numerical
scheme called HARM, for high-accuracy relativistic magnetohydrodynamics (Gam-
mie et al., 2003), I will attempt to outline the ways it extends beyond the basic
introduction above.
Certainly the most significant extension is that of a general relativistic for-
mulation of the MHD equations accounting for an arbitrary space-time metric, and
which allows HARM to be used to simulate BH accretion disks following the flow
self-consistently all the way past the event horizon. This, for instance, allows the ex-
traction of rotational energy from the spinning black hole by surrounding magnetic
field which has been shown to launch powerful jets via the BH mechanism.
The conservative scheme in HARM updates a set of conserved variables
U ≡ √−g(ρut, T tt , T ti , Bi), (1.19)
18
where T µν is the stress energy tensor subject to the four equations expressing con-
servation of energy and momentum,
T µν;µ = 0, (1.20)
and Bi is the magnetic field 3-vector.
The MHD conservation equations are solved in a coordinate basis and if we
write this out as such,
∂t(
√−gT tν) = −∂i(
√−gT iν) +
√−gT κλΓλνκ, (1.21)
one can see how the curvature of space effective adds a source term to the conser-
vations equations.
Along with the vector of conserved variables, a vector of “primitive” variables
is needed to model the flow within zones,
P = (ρ, u, vi, Bi). (1.22)
vi = ui/ut is the 3-velocity. Both the primitive and conservative variables are
used within the numerical scheme of HARM, but only the conserved quantities
are updated directly. Therefore, HARM uses a multi-dimensional Newton-Raphson
routine to solve for P(U) at the end of each time-step.
The original version of HARM in the 2003 paper evaluates fluxes using an
MUSCL-type scheme, the HLL approach (Harten et al., 1983), with a monotized
central “Woodward” limiter (other limiters available). More sophisticated formu-
lations such as the HLLC scheme are used for the work presented in the following
chapters.
The version of HARM used in the work of this thesis implements choices of
constrained transport scheme based on weighing trade-offs of stability and speed
in each problem studied. In general, the method of constrained transport, which
manifestly conserves ∇· ~B = 0, follows a similar procedure to the Gudunov method
outline above. Electromotive forces, integrated along the edges of grid cell faces are
used to update the face-centered (i.e., staggered) magnetic field vector quantities.
19
HARM has gone through multiple generations of code verification tests with
each new code addition, and scales very well for > 104 cores on supercomputers,
making it an ideal code to test global accretion problems onto BHs and to study
the effects of the extreme BH space-time.
1.3.3 Radiative GRMHD
The work in Chapters 3 and 4 involve a new code HARMRAD (McKinney
et al., 2014), which solves the equations of GRMHD + radiative transport using a
closure scheme for the radiation called M1 (Levermore, 1984). This closure formu-
lation has the form
Pij = ETij(E ,Fi), (1.23)
where P is the radiation pressure tensor, E is the energy density of radiation in the
lab frame, and T is the Eddington tensor. In the M1 closure, the Eddington tensor
is taken to be only a function of the first two moments of the specific intensity, the
radiative energy density and the radiative flux vector, respectively, and is further
assumed to be a linear combination of the isotropic unit tensor and the directional
tensor taken in the direction of the radiative flux. In the rest frame of the radiation
fluid, the specific intensity is isotropic. In any other frame (rest frame of fluid or
laboratory, for instance) the radiation field is Dirac-distributed in the direction of
the radiation ‘flow’.
The M1 formulation has several advantages in its treatment of radiation. One
key feature is the ability of the M1 approach to handle radiation fields that have
both optically thick and thin regions, with a natural scaling for all intermediate
states! It deals with shadows much better than flux limited diffusion, where the
radiation pressure tensor is assumed to be proportional to the local gradient of
E alone, and M1 does not have the spurious angular affects when free-streaming
over large distances that the variable Eddington tensor (VET) formulation has (see
Jiang et al., 2012, and references therein). However, M1 does not deal well with
converging radiative zones, where the two incoming radiative fluids can merge in
20
unphysical ways rather than free-streaming through each other.
In any case, M1 offers a good choice when trying to understand the radiation
produced by optically thick, geometrically thin flows, where radiation is processed
by the optically thick disk before being radiated away (and, as I will show, into a
wind that can substantially re-processes the radiation).
The numerical techniques used to solve the equations in HARMRAD and for
GRRMHD (GRMHD + radiation) in general are similar to those outlined above
and those implemented in HARM. The radiation adds sometimes stiff source terms
to the conservation equation, and so they are usually dealt with using a hybrid
explicit/implicit method where fluxes, conserved quantities, and source terms are
updated iteratively, solving the radiative MHD equations to 2nd order or more in
time and space.
1.4 Outline - Description of chapters
As introduced above, there are several fundamental problems in BH accretion
that can be addressed using a code that solves the equations of GR(R)MHD:
(1) Are geometrically thin, radiatively efficient disks able to advect large-scale
poloidal magnetic flux inward? If so, can it reach the MAD state, where the
MRI is quenched the the fundamental behavior of the accretion mode could
be very different than our standard model of thin-disk accretion? Is MAD an
inevitable outcome?
(2) What is the behavior of a MAD thin disk?
(3) What is the radiative efficiency of a MAD thin disk? Is it enhanced compared
to that predicted by the standard NT model of accretion?
(4) What is the radiative efficiency of super-Eddington accretion in the MAD
mode?
Chapter 2 addresses the first three problems listed, but uses an ad hoc cooling
function to keep the disk thin and to estimate radiative efficiency of the MAD
21
thin disk. Chapter 3 is very similar in setup but uses HARMRAD to address the
same problems with self-consistent radiative transport. Chapter 4 describes the
use of HARMRAD in simulating thicker, super-Eddington MAD accretion flows,
where characteristics of the flow associated with the MAD state are responsible
for significant enhancement of radiative efficiency compared to traditional super-
Eddington, or slim-disk models of accretion.
In addition to my work studying BH accretion flows using ideal MHD, as part
of my early graduate work I studied an astrophysical problem in which non-ideal
MHD effects, namely anisotropic conduction, must be accounted for. Chapter 5 is
comprised of this study of properties of heat transport in the intracluster medium
of galaxy clusters. The efficiency of accreting SMBHs in these systems, and the
coupling of that energy output to the ICM, is still under dispute under the title
of the ‘cooling flow problem’. While Chapters 2-4 address the energy output of
AGN in various states, Chapter 5 addresses a related issue: is efficient coupling of
AGN power to the ICM necessary to provide the required heating, or can regulated
heat conduction solve the problem without the invocation of AGN feedback. Future
studies of BH accretion flows will necessitate similar departures from ideal MHD.
In fact, simulations of disks with weakly collisional and collision-less plasmas are
already under way.
22
Chapter 2: Efficiency of Thin Magnetically-Arrested Disks Around
Black Holes
Abstract
The radiative and jet efficiencies of thin magnetized accretion disks around
black holes (BHs) are affected by BH spin and the presence of a magnetic
field that, when strong, could lead to large deviations from Novikov-Thorne
(NT) thin disk theory. To seek the maximum deviations, we perform gen-
eral relativistic magnetohydrodynamic (GRMHD) simulations of radiatively
efficient thin (half-height H to radius R of H/R ≈ 0.10) disks around mod-
erately rotating BHs with a/M = 0.5. First, our simulations, each evolved
for more than 70, 000rg/c (gravitational radius rg and speed of light c), show
that large-scale magnetic field readily accretes inward even through our thin
disk and builds-up to the magnetically-arrested disk (MAD) state. Second,
our simulations of thin MADs show the disk achieves a radiative efficiency
of ηr ≈ 15% (after estimating photon capture), which is about twice the NT
value of ηr ∼ 8% for a/M = 0.5 and gives the same luminosity as a NT disk
with a/M ≈ 0.9. Compared to prior simulations with . 10% deviations, our
result of an ≈ 80% deviation sets a new benchmark. Building on prior work,
we are now able to complete an important scaling law which suggest that ob-
served jet quenching in the high-soft state in BH X-ray binaries is consistent
with an ever-present MAD state with a weak yet sustained jet.
2.1 Introduction
Black hole (BH) accretion systems can operate as efficient engines converting
gravitational potential energy and BH spin energy into radiation. Quasars, active
23
galactic nuclei (AGN), X-ray binaries, gamma-ray bursts, and other BH accretion
systems enter a thin radiatively efficient mode when the luminosity is between ∼
10−2 and ∼ 0.5 times the Eddington luminosity (Abramowicz et al., 1995; Narayan
and Yi, 1995). The Novikov and Thorne (1973b), NT, calculation for the thin disk
radiative efficiency assumes emission terminates at the inner-most stable circular
orbit (ISCO). However, accretion is driven by magnetic stresses via the magneto-
rotational instability (MRI) (Balbus and Hawley, 1991, 1998a) or magnetic braking
by large scale field threading the disk (Blandford and Payne, 1982) or BH (Blandford
and Znajek, 1977), hereafter BZ. A fundamental question in accretion theory has
been whether the NT model is applicable to actual magnetized disks where magnetic
flux couples to a rotating BH and increases the radiative efficiency near the ISCO
(Gammie, 1999; Krolik, 1999). Any extra radiative efficiency could, in principle,
affect BH spin estimates obtained from Soltan’s argument (e.g., see Shapiro 2005)
or from continuum emission spectrum observations (McClintock et al., 2011).
An upper limit to the amount of magnetic flux a disk and BH in steady-state
can support is achieved when an accretion flow reaches the so-called magnetically-
arrested disk (MAD) state, when magnetic flux accretes until magnetic forces push-
ing out balance gas forces pushing in (Igumenshchev et al., 2003; Narayan et al.,
2003; Igumenshchev, 2008; Tchekhovskoy et al., 2011; McKinney et al., 2012). If jet
power scales as Pj ∝ B2a2 for magnetic field B and BH spin a/M (Blandford and
Znajek, 1977), then the MAD state leads to the maximum jet power. The existence
of such a maximum well-defined magnetic flux for a given spin is assumed when
using transient jets in BH X-ray binaries to test the BZ prediction of power vs. spin
(McClintock et al., 2014). Destruction of such a strong magnetic flux could lead
to such transient jets without the need for BH spin (Igumenshchev, 2009; Dexter
et al., 2014), while BH spin measurements may be more reliable for thin MADs that
could have angular momenta aligned with the BH spin axis near the BH (Polko and
McKinney, 2015).
General relativistic magnetohydrodynamic (GRMHD) simulations of thin disks
have found different degrees of deviations from NT theory, but in these cases there
24
are no more than 10% fractional deviations for the radiative efficiency (Shafee et al.,
2008; Noble et al., 2009; Penna et al., 2010a; Noble et al., 2010, 2011). Differences
in results are likely due to the amount of magnetic flux that threads the BH and
disk (Penna et al., 2010a). Unlike other choices for initial conditions, the amount of
magnetic flux threading the BH and disk in a MAD state is independent of the initial
magnetic field strength or geometry as long as there is a plentiful supply of mag-
netic flux (McKinney et al., 2012; Tchekhovskoy and McKinney, 2012). Numerous
non-MAD GRMHD simulations of relatively thick disks have been performed (De
Villiers and Hawley, 2003; McKinney and Gammie, 2004; De Villiers et al., 2005;
McKinney, 2005, 2006; Komissarov and McKinney, 2007; McKinney and Narayan,
2007; Fragile et al., 2007; McKinney and Blandford, 2009; Narayan et al., 2012).
Simulations of MADs have so-far involved only relatively thick disks with no cool-
ing (Tchekhovskoy et al., 2011; McKinney et al., 2012) and thick disks with radiative
transfer (McKinney et al., 2015a).
We perform fully 3D GRMHD simulations of radiatively efficient thin (half-
height H to radius R of H/R ≈ 0.1) disks around rotating BHs (dimensionless
spin, a/M = 0.5) that reach the MAD state. We expect the MAD to maximize the
radiative efficiency of thin disks and lead to the maximum deviations from the NT
model. First, we test how easily magnetic flux advects inward by threading the disk
with weak magnetic field. Second, we setup a disk that is initially nearly MAD, let
it become MAD to large radius, and then measure the radiative efficiency.
We describe the physical and numerical setup in §2.2, present results and
provide discussions in §4.3, and summarize in §4.4.
2.2 Fully 3D GRMHD Thin MAD Model
For this study we use the HARM GRMHD code (Gammie et al., 2003) with
various improvements (McKinney and Blandford, 2009; McKinney et al., 2012) with
two physical setups, both with spin a/M = 0.5. Each was run with a low resolution,
medium resolution, and high resolution for convergence testing and the high resolu-
25
tion was used for detailed analysis. We name the simulations in this work MADxy
where x=i,f signifies whether the disk starts with enough flux to be nearly MAD in
the initial, i, or only MAD in the final, f, state. y=HR,MR,LR distinguishes the
grid as high, medium, or low resolution.
2.2.1 Physical Model
Figure 2.1: Snapshots of the initially sub-MAD (MADfHR - left panel) and nearly-
MAD (MADiHR - right panel) models at t=0, showing ρ in color (with legend) with
magnetic field lines (black lines). In MADfHR there is no magnetic field initially
present within the ISCO or on the horizon and the field is weak througout the disk.
On the other hand, in MADiHR, a sufficient magnetic flux is present throughout
the disk and on the horizon so that the disk will quickly become MAD out to large
radius. The higher density surrounding the horizon in MADiHR at t=0 is the result
of the density floor activated due to the presence of an initially strong magnetic
field.
The initial disk is Keplerian with rest-mass density that is Gaussian in angle
with a half-height-to-radius ratio H/R ≈ 0.1 and a power-law radial distribution
ρ ∝ r−0.6. Since this is not an equilibrium solution near and inside the ISCO at rISCO,
the density is tapered as ρ → ρ(1 − .95(rISCO/R)3) and truncated at rISCO. The
26
Figure 2.2: (Caption on next page.)
27
Figure 2.2: (Previous page.) Evolved snapshot of the MADfHR model at t ≈
70, 000rg/c showing log of rest-mass density in color (with legend) in both the z−x
plane at y = 0 (top-left panel) and y − x plane at z = 0 (top-right panel). Black
lines with arrows trace field lines, where thicker black lines show where field is lightly
mass-loaded. The bottom panel has 3 subpanels. The top subpanel shows the mass
accretion rate through the BH (M˙H). The middle subpanel shows the magnetic flux
threading the BH (ΥH), normalized so that order unity is a dynamically-substantial
amount of magnetic flux. The bottom subpanel shows the gas efficiency measured
on the horizon (ηH,g) and jet efficiency measured at r = 50rg (ηj). In this model, the
disk initially has a weak large-scale poloidal magnetic field, which is transported by
the MRI and other magnetic stresses toward the BH. Eventually, by t ∼ 55000rg/c,
the magnetic flux saturates on the BH and cannot grow despite plentiful magnetic
flux outside. The MAD state that develops consists of a highly non-axisymmetric
clumpy thin disk with a weak (ηj ∼ 1%) jet.
gas internal energy density egas is obtained from the vertical hydrostatic equilibrium
condition of H/R ≈ h¯ ≡ cs/vrot for sound speed cs ≈
√
ΓPgas/ρ and rotational
speed vrot chosen to be Keplerian with velocity vK = (r/rg)/((r/rg)
3/2 + a/M). In
order to seed the MRI we use a 10% random perturbation to the pressure that obeys
the ideal gas law Pgas = (Γ−1)egas. We choose Γ = 4/3 because the disk is expected
to be radiation dominated near the BH. We set an atmosphere for the disk with
ρ = 10−4(r/rg)−3/2 and egas = 10−6(r/rg)−5/2.
We choose an initially purely poloidal large-scale magnetic field. With other
low resolution (MADxLR) simulations already performed, we know ahead of time
the magnetic flux saturation level reached on the BH and throughout the disk in
the quasi-equilibrium MAD state, so we then tune our initial conditions to be either
nearly MAD (for MADi) or sufficiently sub-MAD to be MRI dominated initially
(for MADf). Net magnetic flux on the BH and through the disk out to some radius
is conserved and only changes by accretion or expulsion through that radius. We
refer to the simulation as in ‘quasi-equilibrium’ when characteristic diagnostics of
28
Figure 2.3: (Previous page.) Plotted from left to right are snapshots from MADfHR
at t=[0,10000,16320,40000]rg/c. Black contours (with a single colored white) trace
the same values on all snapshots and show the φ-averaged magnetic field. The
contours demonstrate the transport of magnetic flux radially toward the horizon.
All magnetic flux is initially threading the disk only, then accretes most rapidly
through coronal processes along the surface of the disks where it piles up at these
latitudes before it does in the mid-plane of the accretion flow. Reconnection near the
horizon where the disk is very thin leads to magnetic loops (seen in middle panel)
in the inner disk out to ∼ 20rg, where further reconnection and decoupling from the
magnetically driven wind cause a drop in the effective disk viscosity. Ultimately the
magnetic loops are also accreted and some flux originally threading the disk then
threads only the horizon.
29
the accretion state remain steady, rather than saying it is in ‘inflow equilibrium’,
since the magnetic flux beyond that MAD/non-MAD transition radius is still slowly
evolving. We consider the simulation in quasi-equilibrium out to the radius at which
quantities like M˙ are constant, and at least one inflow timescale has been reached,
tinflow = r/ |vr|.
In both simulation types, MADi and MADf, our initial conditions, shown in
Figure 2.1, provide an effective reservoir of magnetic flux at larger radii that can
be accreted by the inner disk and onto the black hole to reach the maximal amount
of magnetic flux allowed. The supply reservoir is the same for our two setups to
ensure that our two disk states are not the result of initial ‘background solution’
conditions. The supply is vast compared to what ends up on the BH or in the disk
near the BH, but not so large that outward diffusion of flux would be prevented by
our choice.
The coordinate basis φ-component of the vector potential is matched to the
radial gas pressure profile so that the plasma parameter β = Pgas/Pb (Pb as magnetic
pressure), and thus the number of critical MRI wavelengths spanning a full scale-
height of the disk, is a tunable constant for all radii. For MADiy runs, we found
it useful to transition from MAD-flux levels to a weaker sub-MAD field beyond a
transition radius, r0,MAD = 30, with a transition function F± = 0.5± 1pi arctan((r −
r0,MAD)/4). This keeps the evolution of the disks at larger radii the same and allows
us to better capture the ‘pile-up’ effect of flux capture producing the MAD state on
the BH and inner disk. The vector potential is then calculated by radially integrating
Aφ,r ∝ (sF− + F+) sin θν ·
B
θ
0 if r > rtr
Bθtr,0r
(−0.6/2+3/2)/µ if r <= rtr
(2.1)
where Bθ0 is the field threading the disk mid-plane at time t=0 necessary to start
with s = 2 and β = 200/s2 = 50 inside r0,MAD, and with β = 200 elsewhere.
Thus, the disk within r0,MAD initially has Sd ≈ 0.7, where Sd is the number of MRI
wavelengths within two disk scale heights, Sd . 0.5 where the MRI is suppressed
(Balbus and Hawley, 1998a), and MADs typically have had Sd . 0.25 (McKinney
30
et al., 2012). This means the MRI is marginally suppressed in MADs. The field
threading the BH horizon and within the ISCO matches to the disk values at a
transition radius of rtr = 12. µ = 4 ensures the initial flux residing within rtr is
small enough (≈ 2.5 times smaller than what ends up on the BH plus out to rtr)
that the BH and the region between the horizon and rtr can slowly build-up to a
MAD level.
To ensure a MAD can readily develop for a thin disk, our MADf model starts
with sub-MAD flux at t=0, with β = 200 outside rtr. The solution for Aφ is the
same as the r > rtr branch of Eq. 2.1 but with r0,MAD → −∞ and a linear taper in
the field from rtr to the ISCO so there is no initial magnetic flux any closer.
While there is only a factor of 2 in field strength in the inner disk separating
our two initial conditions, these set-ups straddle the boundary of MRI stability in the
disk. While MADfy simulations would still be considered to have a strong magnetic
field, with the associated large effective viscosity, they are clearly distinguished from
the limiting MAD state. Over the course of our MADfHR simulation, the disk and
BH inside r = 35rg accumulate an amount of flux equivalent to that contained in
the disk initially out to r ∼ 110rg. The disk thus evolves significantly away from its
initial conditions toward the MAD limit.
2.2.2 Numerical Grid and Density Floors
We use similar numerical grid mapping and identical boundary conditions as
McKinney et al. (2012), except a smooth arctan transition from exponential to
hyper-exponential radial grid spacing is used. For our high (mid) resolution grid
we choose resolution Nr × Nθ × Nφ = 192x96x208 (168x64x152) with radial grid
spanning Rin ≈ 0.75rH (horizon radius rH) to Rout = 104rg.
The θ-grid spans from 0 to pi with the mapping of McKinney et al. (2012),
but with njet = 0 used for our entire disk since it is thin at all radii. We adjust the
constants of the mapping so that at r/rg = 6, 30, 100 we resolve the disk vertically
with 56,62,60 cells within ±2H= 2 · 0.1R. We flare the inner radial grid to avoid
unnecessarily small limiting time-steps. The polar boundary conditions are trans-
31
missive (see appendix in McKinney et al. 2012), and the φ-grid is equally spaced
from 0 to 2pi with periodic boundary conditions.
Apart from our direct resolution studies, we also measure convergence quality
factors for the MRI and turbulence correlations lengths. Our values demonstrate
good resolution and are reported in Section 2.3.5.
2.2.3 Cooling
To keep the disk thin and radiatively efficient, we implement an ad hoc cooling
function similar to Noble et al. (2010) with
Sν =
(
dU
dτ
)
uν , (2.2)
as a source term to the conservation equations, ∇µT µν = Sν , for gas stress-energy
tensor T and 4-velocity uν . This assumes an isotropic comoving rate of thermal
energy loss dU/dτ . dU is computed so that cooling decreases internal energy toward
a target temperature set by the desired H/R ≈ 0.1 and local density and radius.
Radiative cooling is efficient for a NT thin disk with a cooling timescale on the
order of the orbital timescale (Novikov and Thorne, 1973b). The NT model assumes
radiation gets out only by diffusion, while our disk is clumpy and can be optically
thin near the ISCO (Zhu et al., 2012). Radiative GRMHD MAD simulations are
required for an accurate cooling timescale, which could be faster than for literal NT
disks (McKinney et al., 2015a). We only try to keep the disk thin, and we settled
upon a graded cooling timescale dτ so that
dU
dτ
= −ρTtarget
dτ
= −ρTtarget ·

ΩK/τcool
(
H
νAlf
)
> τcool/ΩK
νAlf
H
τcool/ΩK >
(
H
νAlf
)
> 2∆t
1
2∆t
(
H
νAlf
)
< 2∆t
(2.3)
where νAlf is the Alfvenic frequency and thus H/νAlf is the Alfven time across one
disk scale height, 2∆t is the code time-step with a prefactor 2 for stability, and
ΩK = vK/r. The target temperature is from thin disk theory and is isotropic across
32
cylindrical radii R,
Ttarget ≡
(
0.1 ·R
a+ r3/2
)2
where the denominator uses spherical radius r to control the temperature along the
polar axes. If Tgas = Pgas/ρ < Ttarget then dU/dτ = 0. To ensure there is always
rapid enough cooling, we set τcool = 0.1 · 2pi that enforces cooling on 1/10th the
orbital timescale.
We also use a temperature ceiling that acts as another cooling term, operating
after each timestep. We limit the temperature where β < 0.1 and Tgas > Ttarget.
Unfortunately these losses were not tracked in the version of HARM used. When
we sum the radiative losses from the cooling function during our analysis, in order
to compute the contribution lost to the temperature ceiling, we restarted MADiHR
at a late time with the ceiling turned off. The additional heat is then captured
by the cooling function which we can measure for the radii near the BH where the
temperature ceiling mattered. Accounting for all cooling effects gives us a nearly
constant radial profile for total efficiency, and so the procedure correctly recovers
how a steady-state flow should behave.
Rest-mass density is injected whenever 2Pb/ρ > 200 in order to ensure the code
remains stable, and internal energy density is added when egas/ρ > 100 (though this
limit is never reached).
The cooling keeps the disk very close to the target scale height. After a very
short (< 5000rg/c) initial period of disk evolution, the disk height is very stable.
2.2.4 Diagnostics
The disk’s geometric half-angular thickness (H) per radius (R) is
H
R
(r, φ) ≡
(∫
θ
ρ(θ − θ0)ndAθφ
)1/n(∫
θ
ρdAθφ
)1/n , (2.4)
where
θ0(r, φ) ≡ pi
2
+
(∫
θ
ρ(θ − pi/2)ndAθφ
)1/n(∫
θ
ρdAθφ
)1/n , (2.5)
33
and with surface differential dAθφ. Typically we use n = 1 but also compare to
using other n.
The radiative energy and angular momentum lost are
Lν(r) =
∫
Rrν(r)dAθφ = −
1
Dt
∫
Dt
∫ r
Rcap
SνdtdV, (2.6)
respectively, for ν = t, φ for some period Dt over all angles with photon capture ra-
dius Rcap (see §4.3). The radiation contribution from beyond the inflow equilibrium
radius (∼ 30rg for our MADiHR model) is conservatively estimated by using NT’s
radial dependence.
The mass accretion rate, gas+radiative energy efficiency, and specific angular
momentum accreted are, respectively,
M˙0 =
∣∣∣∣∫ ρurdAθφ∣∣∣∣ , (2.7)
η = −
∫
(T rt + ρu
r +Rrt )dAθφ
[M˙0]H
, (2.8)
 =
∫
(T rφ +R
r
φ)dAθφ
[M˙0]H
(2.9)
where R is the radiation stress-energy tensor and [M˙0]H is the time-averaged M˙ on
the horizon. These diagnostic quantities are calculated as totals and for the jet and
wind. The jet is defined as where 2Pb/ρ > 1, and the wind is where 2Pb/ρ < 1
with an outgoing flow. The radiative efficiency and specific angular momentum vs.
radius are, respectively, ηr = −Lt(r)/[M˙0]H and r = Lφ(r)/[M˙0]H .
Time averages for the MADfHR simulation are made over 55000-68000rg/c,
the time period during which the horizon flux has reached the MAD saturation.
For MADiHR, M˙ is average from 10000-70000rg/c, and all other quantities from
30000-70000rg/c except where explicitly stated.
With immediate proper normalization (rather than for only final quantities in
the flux definitions in McKinney et al. 2012), the half-hemisphere horizon flux is
ΨH ≡ 0.7
∫ θ=pi
θ=pi/2
dAθφB
r√
[M˙0]H
, (2.10)
34
which is negative compared to the integral from θ = 0 to pi/2, for radial magnetic
field strength Br in Heaviside-Lorentz units.
The θ magnetic flux vs. radius at angle θ is
Ψθ(r, θ) ≡ 0.7
∫ r
rH
√−gdx(1)dx(3)Bx(2)√
[M˙0]H
, (2.11)
integrated in internal coordinates, xα ≡ (t, x(1), x(2), x(3)). Thus, the vertical mag-
netic flux threading the equator is
Ψeq(r) ≡ Ψθ(r, θ = pi/2), (2.12)
and we define the total magnetic flux along the equator as
Ψ(r) ≡ ΨH + Ψeq(r). (2.13)
For all forms of Ψ, all φ-angles are integrated over.
The absolute magnetic flux (Φ) is computed similarly to Ψ, except one 1)
inserts absolute values around the field (e.g. Br and Bθ in the integrals); 2) puts
absolute values around the integral; 3) integrates over all θ then divides by 2; and
4) does not immediately normalize. For example, Φr(r, θ) = (1/2)
∣∣∫ dAθφ|Br|∣∣.
We introduce the Gammie (1999) model version of Υ for consistency in com-
paring to previous literature, normalized like our version of Ψ,
Υ ≈ 0.7 Φr√
[M˙ ]t
, (2.14)
which accounts for Φr being in Heaviside-Lorentz units (Penna et al., 2010a). Com-
pared to Gaussian units version of φH ≡ ΦH/
√
M˙r2gc defined in Tchekhovskoy et al.
(2011), Υ ≈ 0.2φH.
Viscous theory gives a GR α-viscosity estimate for vr of vr,visc ∼ −Gα(H/R)2 |vφ|,
as defined in McKinney et al. (2012), leading to a definition for the measure of ef-
fective α-viscosity
αb,eff ≡ vr
vvisc/α
, (2.15)
35
whereas viscosity from local shear-stress in the disk is
αb ≈ αmag = −brbφ
pb
, (2.16)
where we have dropped other (e.g. Reynolds) terms from the full α because they are
negligible compared to the magnetic term in the MAD regions of our simulations.
These α’s are obtained by separately volume-averaging the numerator and denomi-
nator in θ, φ for each r with a weight of ρ and only including material with b2/ρ < 1
in order to focus the measurement on the disk material and not the lower density
corona. The quantity αb measures the magnetic stress within the dense regions,
while the αb,eff measures the actual effective α based upon the radial velocity of the
dense regions. For example, very little magnetic flux may thread the dense regions,
but external magnetic torques may push and drive angular momentum transport
and inflow of the dense regions. αb,eff is only accurate far outside rISCO because of
the G term.
2.3 Results and Discussion
With our fully global 3D GRMHD simulations we can demonstrate whether the
non-linear MRI and other MHD physics in thin accretion flows is able to transport
large-scale magnetic flux towards the BH and build-up the magnetic field to the
MAD state. We first consider a disk with a weak magnetic field using our MADfHR
setup to see if magnetic flux transports inward in a thin, MRI-dominated disk. Then,
having confirmed thin disks transport magnetic flux efficiently to the BH and inner-
radial disk, we use our second setup MADiHR to study a MAD where magnetic
saturation is reached to large radii. This two-step procedure is necessary due to
the extreme computational resources that would otherwise be necessary to build-up
large-scale flux to the MAD state out to large radii in an initially MRI-dominated
disk. In the thin-MAD state, we find the horizon contains an equivalent amount of
flux as a significant radial portion of the disk, so our simulation with no flux initially
inside the inner edge of the disk takes a long time to accrete enough flux to reach
the MAD state on the horizon.
36
Figure 2.4: (Caption on next page.)
37
Figure 2.4: (Figure on previous page.) Time and φ-averaged quantities are plotted
for MADfHR. The time average is taken over 55000-68000 rg/c, during which time
the horizon has reached a MAD saturation of flux. Rest mass gas flow is traced
by black streamlines. The colored (green and blue) thick lines correspond to time-
averages of quantities Y = {ur, b2/ρ}, respectively, at values [Y ]t = {0, 1}. While
near the BH the flow has [b2/ρ]t & 1 as averaged directly, the dense inflow has
[b2/ρ]t . 1 at all other radii. The time and spatial variability of the jet causes the
time-φ average to blend jet and disk material, so the average [b2/ρ]t & 1 bleeds into
the disk. Similarly, [b2]t/[ρ]t identifies a jet boundary in the opposite way, causing
too much of the polar/jet region to appear to have [b2/ρ]t . 1. No particular average
is correct, but our choice includes the part of the wind and disk that experiences
extra acceleration through interaction with the jet. Streamlines nearest the polar
axes are calculated using velocity structure five cells away from the axis for better
visualization of the qualitative geometry of the jet flow. The flow exhibits equatorial
asymmetry over many inflow times as found in thicker MAD simulations, likely a
result of bulk flows associated with the very initial evolution of the disk. The red
magnetic field line contours show how the magnetic flux threads the black hole and
inner disk.
38
Figure 2.5: Time (from 30000-70000rg/c) and φ-averaged quantities are plotted for
MADiHR. Identical layout as Fig. 2.4. The disk is magnetically arrested out to a
transition radius of ∼ 35rg, which leads to a short inflow time-scale compared to
the MADfHR disk at those radii, so the time-averaged velocity fields and other disk
structure is smoother and more symmetric.
39
Figure 2.6: Evolved snapshot of the MADiHR model at t ≈ 70, 000rg/c, with iden-
tical layout to Fig. 2.2. The MAD state is highly non-axisymmetric and clumpy
with a weak (ηj ∼ 1%) jet. The rest-mass accretion rate and magnetic flux on the
horizon reach the MAD state quickly, which was as desired so that we can study
a quasi-steady MAD out to large radii over long time periods. The dimensionless
magnetic flux ΥH on the horizon is in agreement with our MADfHR model, showing
a disk initial condition with a relatively stronger poloidal field strength still reaches
the same MAD state near the BH.
40
2.3.1 Demonstration of flux accumulation to MAD state
First, we consider a thin radiatively efficient disk threaded by weak large-scale
magnetic field. This model, MADfHR/MR, started in the sub-MAD state with
β ≈ 200 (a field strength that is ∼ 7 times weaker than when the MRI becomes
suppressed and ∼ 40 times weaker than the final MAD state) so that the MRI
operates over long times.
In MADfHR, disk and coronal processes lead to accretion of magnetic flux
inward, such that over t ≈ 70, 000rg/c the BH and disk become MAD out to r ∼ 15rg
in the main disk body. Figure 2.2 shows a steady accumulation of flux onto the
horizon during the first ∼ 30000rg/c. Movies of the phi-averaged vector potential
and the time-r-θ dependence of flux accumulation reveal a faster pile-up of flux at
higher disk latitudes via coronal processes, and a slower steady inward transport
of flux near the equator by some combination of MRI processes and wind driving.
Figure 2.2.1 shows the different stages of the inward flux transport. For detailed
physical descriptions in prior explorations of this behavior see, e.g. Rothstein and
Lovelace 2008; Beckwith et al. 2009.
Accumulation onto the BH is dominated in the first 5000rg/c by large-scale
MRI ‘channel’-type modes typical of these initial conditions that last until φ-symmetry
of the poloidal field is lost to parasitic modes, as well as coronal processes in the
inner-most disk. Even during this very early evolution there is some separation of
net flux passing through the disk from the most dense disk material, development
of a two-phase flow, and during the initial transient evolution inside r ∼ 15rg the
field transported from the disk to the horizon leaves behind some extra mass. This
separation, as well as reconnection of field near the equator and on the horizon, leads
to the remaining matter of the inner disk having a weaker net field and the MRI
is under-resolved in this region for about ∼ 20000rg/c. As flux is accumulated into
this under-resolved region the effective viscosity rises until the over-dense region is
fully accreted and the field in the inner disk is strong enough for the MRI to be
resolved once again. Since the accretion rate drops and the strong field in the disk
41
near the horizon now provides low-order asymmetries previously suppressed by the
under-resolved material, the field on the horizon, now also above the MAD limit,
can reconnect and push sporadically into the disk. That sporadic reconnection is a
behavioral signature of the MAD state near/on the horizon.
In the much longer run, MADfMR becomes MAD out to a radius of r ∼ 18rg
by the time t ≈ 183, 000rg/c. Plenty more magnetic flux exists at larger radii in both
MADfMR and MADiMR, and it would presumably contribute to the MAD region
of the disk at later times if we continued running the simulations. It is unclear if
this slower transport of net flux inward through the main disk body is primarily due
to the MRI processes or wind driving.
Though the physics behind the saturation of magnetic flux needs more ex-
ploration in future work, it has reached a maximum far below the supply, yet far
beyond the initial value. This demonstrates that large-scale magnetic flux can read-
ily accumulate even with H/R ≈ 0.1 and radiative cooling, and that it saturates at
a maximum value like seen in prior MAD simulations.
Fig. 2.4 shows the time and φ-averaged rest-mass velocity structure and dis-
tinguishes several key features of the accretion flow. This includes showing the jet
boundary and the flow stagnation surface (where material has zero radial velocity,
moves outward at higher latitudes in the disk, and moves inward at lower latitudes).
Most importantly, the plot indicates the time-φ-averaged structure of the magnetic
field once the inner disk and BH have reached magnetic flux saturation. As with the
thicker MAD model plotted similarly in Fig. 6 of McKinney et al. (2012), the evo-
lution of the φ-averaged field structure, which we studied in movies of these frames,
demonstrate ordered poloidal magnetic flux is easily transported through the disk
from large radii through a combination of coronal and disk processes.
2.3.2 Long-term evolution of MAD state
In order to study the thin MAD out to large radii to ensure an accurate quasi-
steady radiative efficiency, we next considered the disk having an initially nearly
MAD-level of magnetic field. Figure 2.1, and Figures 2.5 and 2.6 show the initial and
42
evolved quasi-equilibrium state of our high-resolution model that is initially nearly
MAD inside r = 30rg. As the accretion flow evolves from its initial state, rotation
drives outward angular momentum transport by turbulent magnetic stresses via the
MRI as well as by a magneto-centrifugal wind. Mass, energy, angular momentum,
and magnetic flux are accreted and ejected at large latitudes, and by ∼ 5000rg/c
magnetic stresses balance incoming gas forces so the accretion is MAD out to r ∼
30rg. By the end of the simulation, t ∼ 70000rg/c for MADiHR, the MAD extends
to r ∼ 35rg.
Figure 2.7: (Caption on next page.)
The evolved state is evaluated for specific reported values and time-φ-averaged
quantities using the period t = [30000, 70000]rg/c. However, a slightly longer period
43
Figure 2.7: (Figure on previous page.) Time-averaged quantities, averaged over
10000-70000rg/c, as a function of radius. From top to bottom, panels are: Total
mass accretion rate (M˙0, blue line), energy efficiency (η), specific angular momentum
(j), magnetic flux (Ψ), and disk thickness (H/R, n=1 in Eq. 2.4, blue line). Plots
for η and j show total gas+radiative (black solid), electromagnetic (green), matter
(red), and radiation (black dashed) components. Remaining plots show their initial
value (magenta dot-dash). The plot of Ψ shows the cumulative net vertical flux
including BH and disk contributions (blue). The constancy of M˙ , η, j show where
the flow has achieved inflow equilibrium. Ψ shows the MAD to non-MAD transition
at r ∼ 35rg, while throughout the disk H/R ≈ 0.1 (bound gas), as maintained by
cooling. The primary result of this plot is the efficiency panel showing the radiative
contribution of ηr ∼ 15% with most of emission beyond the standard NT disk model
coming from quite close to the photon orbit where the MAD state differs significantly
in structure and behavior from the NT disk.
t = [10000, 70000]rg/c is used for the time-averaged functions of radius plotted in
Figure 2.7 to minimize the effects of long-term trends predominantly at larger radii.
The time-φ-averaged structure is show in Figure 2.5 and shows qualitatively similar
behavior to MADfHR regarding the structure of inwardly accumulating magnetic
flux and rest mass flow.
The gas+radiative efficiency η ∼ 20%, and specific angular momentum  ≈
−0.4 are constant within the MAD/non-MAD transition r ∼ 35rg, indicating a
quasi-steady MAD state has been reached out to a relatively large radius.
The disk thickness at r ∼ 20rg is regulated by cooling to be H/R ≈ 0.11 for
n=1 in Eq. 2.4, H/R ≈ 0.14 for n = 2, density-weighted normalized H/R ≈ 0.13
(Penna et al., 2012), and thermal thickness of h¯ ≡ arctan (cs/vrot) ≈ 0.13. An
accurate measurement of H/R is complicated by the fact that magnetic compression
leads to a geometric H/R ≈ 0.05 near the BH, mass is launched into a broad wind,
and vertical disk fluctuations are the same order as H/R. For the φ-averaged density
profile and including only bound material (i.e. −ut(ρ+egas +Pg +2Pb)/ρ < 1), valid
44
when comparing with NT, we get H/R ≈ 0.10, which is our target value we desired
for the disk proper. Figure 2.8 shows the time evolution of H/R at four radii for
n = 1. The same quantities are plotted for MADfHR for comparison in Figure 2.9.
While it is less physically informative than Ψ in evaluating the quantity of
magnetic flux present, we refer to values of Υ for comparison to prior work. When
only considering the highly ordered dipole magnetic field on the horizon in the
MAD state, these quantities are usually very similar in magnitude anyway. Our
thin MAD achieves a magnetic flux on the BH of ΥH ≈ 5. The stability of ΥH with
time and the presence of plentiful magnetic flux remaining in the disk (see Fig. 2.8,
lower panel) indicate that the disk has achieved a balance between magnetic forces
and gas forces with a strong magnetic field (Υ  1). Thus one expects maximum
magnetically-induced deviations from the standard NT thin disk solution.
Large temporal deviations are driven by magnetic flux repeatedly building up
on the BH horizon to a level the disk cannot sustain. As soon as a weak point in the
disk exists, a non-axisymmetric mode in in-falling matter, gas pressure distribution,
or rest-mass, then there is a sudden re-emergence of magnetic flux back into the
disk, a process that creates large low-density bubbles in the form of a magnetic
prominence. The magnetic flux that is ejected is then processed by turbulence
by the outer parts of the disk until being accreted again at which point the process
repeats. Figure 2.10 shows one of the most extreme magnetic prominences occurring
in the MADiHR simulation.
2.3.3 Radiative Efficiency
We find that nearly as much radiation is released near and inside the ISCO
as in the entire rest of the disk, so we calculate the radiative efficiency assuming
100% photon capture within a choice among three radii: the horizon, the photon
orbit radius r ≈ 2.4, and the critical impact parameter radius, r ≈ 4.1, for which
geodesics with any smaller impact parameter intersect the horizon. An estimate
of photon capture using ray-tracing from a thin disk gives as much as 30% escape
fraction from the photon orbit for equivalent BH spin (see Beckwith et al. 2006),
45
and more than 70% for r ≈ 4.1. This suggests using the photon orbit may give the
best choice for our simple approximation.
The radiative efficiency is ηr ∼ 15% for MADiHR (< 0.5% lower for MADiMR)
when including only radiation from outside the photon orbit. This is about twice
larger than the ηr ≈ 8% efficiency in standard NT thin disk theory for BH spin
a/M ≈ 0.5. The impact parameter estimate is an extreme, albeit unphysical, limit
giving efficiency ηr ∼ 9.4%. Also unphysical, including all radiation released outside
the horizon gives ηr ∼ 19%. Unbound gas only contributes ∼ 1% to ηr, a smaller
relative contribution compared to prior simulations (Penna et al., 2010a; Zhu et al.,
2012).
Thus, our thin MAD with H/R ≈ 0.05–0.1 and a radiative efficiency of ηr ≈
15%, demonstrates an ≈ 80% deviation from NT. Prior GRMHD simulations of thin
disks have had H/R ≈ 0.1 with 6% deviations (Noble et al., 2009), H/R ≈ 0.06
with 6% to 10% deviations (Noble et al., 2011), and H/R ≈ 0.07–0.13 with 4% to
5% deviations (Penna et al., 2010a). So despite having similar H/R, our thin MAD
has vastly higher deviations from NT.
2.3.4 Jet Efficiency
Some of the horizon magnetic flux reaches into the polar regions leading to a
jet efficiency of ηj ≈ 1% at r = 50rg for both MADiHR and MADfHR. A caveat is
that the jet within b2/ρ > 1 at r = 50rg is only resolved by about 7 grid cells. So we
check the jet efficiency and resolution on the horizon, where the jet is resolved with
more than 12 grid points per hemisphere and still has a jet power of η ≈ 1%. These
horizon jet grid cells are sufficient resolution to capture the force-free behavior of
the jet near the polar axes, according to prior resolution studies (McKinney and
Gammie, 2004). So the jet efficiency measurement is robust to resolution changes,
slight model changes, and what radius it is measured at.
The BH spin and disk rotation drive a magnetized wind with efficiency ηw ≈
4%, though at least an additional 1% wind efficiency may be lost to radiation via
the very efficient cooling we use as well as lost to the comparatively low resolution
46
Figure 2.8: Shows temporal evolution of a few quantities for our MADiHR model.
The disk scale height H/R is shown in the top panel. We show the temporally-
smoothed H/R as plotted for radii RISCO (solid-black), 15rg (short-dashed), 30rg
(dots), and 60rg (long-dashed). The red line is the un-smoothed evolution at 15rg
to demonstrate the level of variability. The H/R value is regulated by cooling to
be H/R ∼ 0.1, which is what is achieved for radii just outside the ISCO and out to
the inflow equilibrium radius of r ∼ 35rg. The middle panel shows the αb viscosity
parameter spatially-averaged over the disk at a radius 10rg. The relatively high
αb ∼ 0.5 is indicative of MADs, implying efficient transport of mass and angular
momentum. The third panel shows the ratio of the total magnetic flux threading
the BH to the total magnetic flux available at the start of the simulation inside the
radius out to which the simulation ultimately reaches at least one inflow time. This
shows that the horizon retains only a fraction of the flux available during the course
of the simulation, so the quasi-steady state has a horizon-threading flux that has
reached saturation.
47
Figure 2.9: Shows temporal evolution of same quantities as Figure 2.8 for our
MADfHR model. The disk scale height H/R is shown in the top panel. We show
the temporally-smoothed H/R as plotted for radii RISCO (solid-black), 15rg (short-
dashed), 30rg (dots), and 60rg (long-dashed). The red line is the un-smoothed
evolution at 10rg to demonstrate the level of variability and show the evolution of
the under-resolved dense torus present at early times. αb ∼ 0.5 averaged around 8rg
drops to low values during that time as expected.
48
in the jet and highly magnetized wind regions.
The jet and wind contribute only a bit to the total η ≈ 20%. However, the jet
and wind carry nearly as much angular momentum outward as the dense Keplerian
disk carries inward. This results in the BH spin-up parameter s = −−2(a/M)(1−
η) ≈ −0.4, i.e. the BH is somewhat spinning down (McKinney et al., 2012).
Further, the jet efficiency as a function of H/R is constrained by our thin
MAD result. Magnetic compression leads to varying geometric H/R vs. radius, so
the more constant thermal thickness h¯ ≡ arctan (cs/vrot) is used (McKinney et al.,
2012). Our thin MAD models have h¯ ≈ 0.13 and older MADs had h¯ ≈ 0.5 (e.g.
with geometric H/R ≈ 0.3 at r = 20rg) or h¯ ≈ 1.5.
Prior thick MAD simulations found that the saturated magnetic flux on the BH
given by ΥH is roughly constant vs. BH spin. That data showed a slight tendency
to peak at a/M ∼ 0.5, but this is probably intrinsic to the initial disk conditions or
due to intrinsic errors that suggest a/M ∼ 0 could be the peak as well (McKinney
et al., 2012; Tchekhovskoy et al., 2012). The thick MAD simulations of moderate
thickness show a clear asymmetry in the ΥH vs. spin for positive vs. negative spin
that is consistent with the asymmetry in jet efficiency vs. spin. So we do not try
to fit the peak in ΥH being at a/M ∼ 0.5, but we capture the asymmetry. We also
capture the fact that this spin asymmetry is seen to go away as the disc becomes
even thicker.
Combining more than 25 different thicker MAD models with our new thin
MAD model, a rough fit of the jet efficiency is given by
ηj ≈ 400%ω2H
(
1 +
0.3ωH
1 + 2h¯4
)2
h¯2, (2.17)
where ωH ≡ a/rH, which assumes steady activation of a Blandford and Znajek
(1977) jet. The ηj ∝ h¯2 scaling is required (compared to ∝ h¯) to reach our thin
MAD with ηj ≈ 1%. That ηj ∝ (H/R)2 is consistent with S7.1 of Penna et al.
(2010a) who analytically found ΥH ∝ (H/R), and with Meier (2001) who argued
that ηj ∝ (H/R)2.
In our model we fit to jet efficiency rather than jet power since this is easily
49
compared to observations of stellar and super-massive BH systems, and because
rest-mass density can be scaled to any value in these simulations. However, our
formula is roughly consistent with what is expected based upon the BZ formula
for jet power ∝ Υ2HΩ2H. Models with h¯ ≈ 0.5 and a/M = 0.5 gave ηj ∼ 20% and
ΥH ≈ 13. Our thin MAD has ΥH ≈ 5, so the BZ formula suggests we should get
ηj ≈ 3% while we got ηj ≈ 1%. Why the factor of three difference?
First, Υ measured on the BH is not relevant to the BZ jet. One cannot use ΥH
with ηj ∝ Υ2H, because in general (and as found in thicker MAD models) magnetic
flux in the jet is a fraction of ΥH and the jet only partially covers the rotating BH’s
surface. Instead, the magnetic flux passing through the jet is what is relevant to
the jet power. The h¯ ≈ 0.5 model had Υj ≈ 10, while our model has Υj ≈ 4. These
are significantly different from the magnetic flux over the entire horizon. This still
implies we should have ηj ≈ 3% for consistency with the BZ expectation.
Second, the BZ formula fails to account for the disk covering a significant part
of the potential power output (McKinney, 2005). The disk at the equator gets in
the way, even for MADs (especially at lower spin when less magnetic compression
occurs). However, most of the BZ power would be at the equator where the Poynting
flux peaks. The disk chops-off this peaked part of the Poynting flux for the jet, giving
some of it to the wind. Some of that BZ power is never produced in the first place
due to the heavy mass-loading by the disc. A better model of that disk part that
connects to the horizon is the Gammie (1999) model of a thin magnetized disc. In
that case, Υ locally in the disk is much less than over the entire horizon, so the
Gammie (1999) model predicts much less than BZ for extraction of energy. One
cannot just fit ηj vs. ΥH and expect it to match BZ.
Third, we have defined the jet arbitrarily as where b2/ρ > 1. The wind is
defined as the rest of the non-radiative flux that adds up to the total efficiency. For
our thin MAD, this is 4% (i.e. 20%-1%-15%). The wind efficiency is partially directly
due to the disk rotation, but some of the wind efficiency derives from magnetic flux
threading the BH that passes into the disk leading to a BH-spin-driven wind. Some
of this BH-spin-driven wind absorbs what would have otherwise (without a disk)
50
become the jet and can account for another missing 2%, giving a total BZ-like
efficiency matching expectations. If this were the sole effect in comparison to the
thicker MAD models, this would suggest that thin MADs may have more powerful
winds per unit jet power relative to thick MADs.
Fourth, the jet efficiency for the h¯ ≈ 0.5 a/M = 0.5 model is a peak among all
the otherwise identical models with different BH spins. Given how much the results
vary for slight changes in the model setup (like initial magnetic field strength), the
ηj ≈ 20% could easily have been ηj ≈ 10% within systematic errors. Then there
would be no discrepancy between our formula and the BZ approximation.
So our Eq. (2.17) is not supposed to be derived from BZ because the BZ-
driven process creates both a “jet” and a “wind”, while we have defined the jet
by that part which would be capable of becoming relativistic. So in that sense,
our formula is useful in application to real jets instead of a test of whether the
theoretical BZ extraction efficiency applies to the simulation data. However, only
direct observational comparisons without ad hoc definitions will test the validity
of these simulations being applicable to real systems. The fitting formula is just a
useful guide.
2.3.5 Resolution
We test the resolution of our simulations using convergence quality factors
for the MRI (Qθφ,MRI; number of grid cells per MRI wavelength in each direction)
and turbulence (Qrθφ,corr; number of cells per turbulence correlation length in each
direction) (Shiokawa et al., 2012; Hawley et al., 2011; McKinney et al., 2012; Hawley
et al., 2013; McKinney et al., 2013, 2014). For our high resolution runs, MADiHR
(MADfHR), at r = 12, 30rg centered on the equatorial flow, we have Qθ,MRI =
160, 55 (Qθ,MRI = 88, 18) with a weighting of
√
ρb2 and we have Qφ,MRI = 120, 100
(Qφ,MRI = 79, 77) with a weight of ρ. Weighted averages are computed for the
numerator and denominator separately before taking the ratio to get Q. The purely
density-weighted values represent a stricter constraint on whether we have sufficient
resolution (because MADs have strong fields mostly outside the dense region, and
51
Figure 2.10: Plotted as in the top frames of Fig. 2.6 is a snapshot of MADiHR at
t = 31568rg/c when the disk experiences an extreme eruption of magnetic flux from
the BH horizon. A large magnetic flux tube creates a low-density cavity in the disk
before it is partially blended with disk material by the magnetic Rayleigh-Taylor
instability. Such prominences occur frequently in the MAD state.
52
the external stresses are important to whether transport occurs). These are plotted
as a function of radius in Figure 2.11. Resolution of turbulent structures in rest-
mass and magnetic energy density are evaluated at four radii r = {rH/rg, 4, 8, 30}rg.
For MADiHR, rest-mass correlation values are, at each radius respectively, Qr,corr ∼
6, 10, 14, 14, Qθ,corr ∼ 40, 35, 31, 34, and Qφ,corr ∼ 6, 7, 8, 11. Correlation values for
magnetic energy density are, similarly, Qr,corr ∼ 6, 12, 17, 15, Qθ,corr ∼ 34, 33, 29, 36,
and Qφ,corr ∼ 8, 9, 12, 16.
For MADfHR, calculated for the disk once the BH has reached the MAD state,
rest-mass correlation values are, at each radius respectively, Qr,corr ∼ 6, 10, 14, 14,
Qθ,corr ∼ 3, 29, 28, 32, and Qφ,corr ∼ 7, 9, 11, 13. Correlation values for magnetic
energy density are, similarly, Qr,corr ∼ 6, 12, 16, 15, Qθ,corr ∼ 5, 28, 30, 31, and
Qφ,corr ∼ 7, 9, 12, 14.
These MRI and turbulent quality factors show convergence according to strict
application of conventional tests. Of course there could be behavior characteristic of
the MAD state that eludes tests for convergence in accretion physics so far devised.
We therefore directly convergence tested our results by performing simulations of
each setup with three resolutions. The qualitative structure of the disk in the MR
and HR simulations are the same, but the lower resolution in MR allowed us to run
these simulations for significantly longer. The nature of magnetic flux accumulation
is the same, but the longer MR simulations ultimately become magnetically arrested
to larger radius, r ∼ 18rg in MADfMR compared to only r ∼ 15rg in MADfHR. A
difference of only ∼ 0.5% in radiative efficiency between MADiMR and MADiHR
indicates very good convergence in our quantity of interest. We find only small
(typically ∼ 10% or less) differences between our MR and HR simulations for the
other quantities reported which characterize the accretion flow. For instance, the jet
power differs by only ∼ 50% between MADiMR and MADiHR, unsurprising since
ΥH ≈ 5 for both. Low resolution ‘testing’ simulations were largely in qualitative
agreement, but with some Q’s indicating lack of convergence.
Figures 2.8 and 2.11 show values of the αb viscosity parameter, dominated by
the magnetic stress. We find αb ∼ 0.5, as consistent with resolved disks in prior
53
MAD simulations. The effective viscosity αb,eff , corresponding to a measurement of
the actual inflow rate indicated by an α parameter, is much larger than αb due to
enhanced inflow caused by external torques driven by low-density highly-magnetized
material that introduces large-scale stresses on the more dense inflow. αb does not
reflect the effective viscosity in large part because the ρ weighting in calculating αb
predominantly excludes the flux tubes that drive the Blandford-Payne-type trans-
port of angular momentum. The dense parts of the disk lose angular momentum by
running into these flux tubes on their orbit. One can think of the MAD, then, as a
two-phase flow where MRI drives turbulence, still, in small dense clumps, but the
large-scale evolution of the disk, at scales where the MRI is quenched by magnetic
field strength, is dominated by magnetic-Raleigh-Taylor and wind/jet processes.
2.3.6 Power-law Fits for Radial Dependence
We now consider radial power-law fits (of the form f = f0(r/r0)
n) for quantities
associated with the disk structure and flow, plotted in Figure 2.12, for MADiHR.
The fits to disk quantities are performed as in McKinney et al. (2012) between
radii r = 12rg and r = 30rg for both our high resolution models. The profiles for
MADiHR (MADfHR) are ρ ∝ r−0.4±0.1 (∝ r−0.8±0.2), pg ∝ r−1.5±0.1 (∝ r−1.8±0.1),
|vr| ∝ r−1.6±0.06 (∝ r−0.8±0.2), |vφ| ∝ r−0.38±0.01 (∝ r−0.49±0.01) , |br| ∝ r−1.4±0.07 (∝
r−1±0.1), and |b| ∝ r−1.2±0.05 (∝ r−0.89±0.02).
The MADiHR fits are more representative of a steady-state MAD solution,
because the disk in MADiHR is actually magnetically arrested over all radii used
for fitting and the model is in quasi-equilibrium out to large radii. Our thin disk
has a somewhat shallower radial density profile than prior MAD models, but are
otherwise consistent. To obtain a density profile accurate to large radii, much longer-
run simulations would be required. The overall fits are consistent with pg ∝ r−1.5
and |b| ∝ r−1 expected for accretion disks, while the shallow density profile could
be connected with the stronger wind (relative to the jet) compared to thicker MAD
simulations.
54
Figure 2.11: MRI quality factors indicating time-φ averaged resolution, viscosity
of the disk, and vertical MRI half-wavelengths per disk scale height (Sd, blue line)
are plotted as a function of radius. Q factors indicated the number of grid cells
per critical MRI wavelength, averaged with weight
√
ρb2 (solid lines) and ρ (short
dashed) across the disk where b2/ρ < 0.5. The third panel shows the viscosity
parameter αb (solid line) and the effective viscosity calculated via the mass accretion
rate (short dashed). The magnetic evolution of the disk is well resolved and leads to
a high viscosity as measured by local stresses, αb, and the rate of mass inflow αb,eff .
The magenta dot-dash line represents the value of Sd at t=0 and the time-average
is blue-solid.
55
Figure 2.12: Several quantities characterizing the accretion flow in MADiHR are
plotted as a function of radius. The top panel shows the density profile ρ (solid),
the internal gas energy ug (dashed), and magnetic energy density ub (dotted). The
inner disk, that has reached the MAD state, has ug ∼ ub indicating equipartition
is achieved even within the dense parts of the disk. The middle panel includes the
disk radial velocity, azimuthal velocity, and Keplerian velocity (vr, solid; vφ, black-
dashed; vK,blue-dashed, respectively). The final panel shows the radial profiles of
each component of the comoving magnetic field, br,θ,φ, in solid, dashed, and dotted
lines respectively. The density is shallower than thicker disk solutions, the disk is
fairly Keplerian, and the magnetic field falls off as expected for accretion disks.
56
2.4 Summary and Conclusions
We have performed fully 3D GRMHD simulations of radiatively efficient thin
accretion disks seeking to maximize the magnetically induced deviations from NT
thin disk theory. We first demonstrated that even in thin disks, after a time
183, 000rg/c and out to r ∼ 18rg, magnetic flux readily advects inward and builds-up
to a MAD level (for which accumulation results in magnetic forces pushing out bal-
ancing against gas forces pushing in). This occurs despite possible outward magnetic
diffusion through the disk.
Our simulation that is MAD out to r ∼ 35rg after 70, 000rg/c has a radiative
efficiency of ηr ≈ 15%. This is ≈ 80% higher than the standard NT thin disk
value of ηNT ≈ 8% and is as if the disk were a standard NT thin disk but with
a/M ≈ 0.9. BH X-ray binaries have L/LEdd ∼ 0.1 in the high-soft state (Kulkarni
et al., 2011), for which the NT model givesH/R . 0.02 near the BH. Prior non-MAD
simulations suggest deviations scale with H/R (Penna et al., 2010a), so typical thin
MAD deviations might be conservatively ∼ 20%. Such a magnitude of deviations
could significantly impact BH spin measurements, such as in Noble et al. (2011) who
found 6% deviations lead to a BH with a/M = 0 giving a spectrum like NT but with
a/M ≈ 0.3. However, it is unclear into what form the energy would be dissipated.
Zhu et al. (2012) found that emission from near the ISCO is mostly high-energy,
but more physics (e.g., dissipation-radiation energy balance, optical depth effects,
BH photon capture) is required for BH spin fitting (Kulkarni et al., 2011). We plan
to obtain spectra via post-processing (Zhu et al., 2015; Narayan et al., 2016). Also,
clumpy rest-mass distribution in thin MAD flows may help to explain spectra from
AGN (Dexter and Agol, 2011) and BH X-ray binaries (Dexter and Quataert, 2012).
It is not a new idea that X-ray binaries in the low-hard (LH) state supply
enough magnetic flux to the BH to enable BZ-type powering of an observable jet
by being geometrically thick, and that the jet quenches when the disk collapses to
become thin in the high-soft state. For instance, plunging material inside the ISCO
may enhance the flux on the BH in the low-hard state only, even in sub-MAD disks
57
(Reynolds et al., 2006). However, these ideas generally propose transition mecha-
nisms in which the jet switches between being either powered or not, depending on
disk state. Our model does not require a special mechanism for the LH state.
Radiatively inefficient flows have H/R ≈ 0.4 (Narayan et al., 2012). Our
fitting formula Eq. (2.17) suggests that jet power in the MAD state can drop
by & (0.4/0.02)2 ∼ 400 (or drop by & (0.4/0.01)2 ∼ 1600 for 4U 1957+11 if at
L/LEdd ≈ 0.06) when undergoing hard-to-soft state transitions. This implies that
the observed jet quenching in BH X-ray binaries (Russell et al., 2011) can occur
despite the presence of a large-scale MAD-level magnetic field in both disk states.
Such hard-to-soft transitions may also be applicable to tidal disruption events at late
times (Tchekhovskoy et al., 2014). The MAD state could be necessary to explain
powerful jet systems (Zamaninasab et al., 2014; Mocz and Guo, 2015), while the
thermodynamically determined disk thickness may also play an important role in
setting just how much magnetic flux a MAD has. Even thinner disks than studied
here would be implied to have ΥH ∝ H/R and so could have much weaker magnetic
flux and much weaker jets rather than being independent of the disk state.
Acknowledgments
We thank Ramesh Narayan, Jack Steiner, Roman Gold, Peter Polko, Alexan-
der Tchekhovskoy, and Eliot Quataert for discussions. We acknowledge NASA/NSF/TCAN
(NNX14AB46G), NSF/XSEDE/TACC/Stampede (TG-PHY120005), NASA/Pleiades
(SMD-14-5451), and UMD Deepthought2.
58
Chapter 3: General Relativistic Radiation Magnetohydrodynamic Sim-
ulations of Thin Magnetically Arrested Disks
Abstract
The classical, relativistic thin-disk theory of Novikov and Thorne (NT)
predicts a maximum accretion efficiency of 42% for an optically thick, radia-
tively efficient accretion disk around a maximally spinning black hole (BH).
However, when a strong magnetic field is introduced to numerical simulations
of thin disks, large deviations in efficiency are observed, in part due to mass
and energy carried by jets and winds launched by the disk or BH spin. The to-
tal efficiency of accretion can be significantly enhanced beyond that predicted
by NT but it has remained unclear how the radiative component is affected.
In order to study the affect of a dynamically relevant large-scale magnetic field
on radiatively efficient accretion, we have performed numerical 3D GRRMHD
simulations of a disk with scale height to radius ratio of H/R ∼ 0.1 around a
moderately spinning BH (a = 0.5) using the code HARMRAD. We use fully
global simulations with different numerical grids and resolutions to measure
the jet, wind, and radiative properties of a magnetically arrested disk (MAD)
that is kept thin via self-consistent radiative transport using the M1 closure
scheme. Our fiducial disk is MAD out to a radius of ∼ 17Rg and the majority
of the total ∼ 13% efficiency of the accretion flow is carried by a magneti-
cally driven wind. We find that the radiative efficiency is slightly suppressed
compared to NT, contrary to prior GRMHD simulations with an ad hoc cool-
ing function, but it is unclear how much of the radiation and thermal energy
trapped in the outflows could ultimately escape.
59
3.1 Introduction
Black hole accretion disks are present in Active Galactic Nuclei (AGNs), X-
ray binaries, gamma-ray bursts (GRBs) and tidal disruption events (TDEs) and are
able to convert gravitational potential energy and BH spin energy into radiation, jets
and winds through the stresses induced via the magnetorotational instability (MRI)
(Balbus and Hawley, 1991, 1998a), magnetic field threading the disks (Blandford
and Payne, 1982) and magnetic field threading the black hole (BH) (Blandford and
Znajek, 1977). Novikov and Thorne (1973b) developed a general relativistic model
for thin disks which has been successfully applied to understand many astrophysical
systems. In this model it is assumed that all energy carried away from the disk
is radiation. In addition, it is assumed that radiative emission ceases inside the
inner-most stable circular orbit (ISCO) where stresses in the NT model disappear.
However, Gammie (1999) and Krolik (1999) demonstrated that additional stress
inside the ISCO can increase the gravitational potential energy converted by the
disk.
Another way to increase the conversion efficiency of a BH accretion disk is to
increase the magnetic stress within the disk by making the field dynamically impor-
tant. If the disk piles up magnetic flux on the BH horizon and in the disk until the
MRI is quenched and the magnetic forces pushing out balance gas forces pushing in,
no more magnetic flux can be transported inward. This is the so called magnetically
arrested disk (MAD) state (Narayan et al., 2003), which has been shown to maxi-
mize jet power. GRMHD simulations of thick disks in the MAD state have found
accretion efficiencies as high as 30% for a black hole with a/M=0.5 and 140% for
a black hole with a/M=0.99 (Tchekhovskoy et al., 2011). Tchekhovskoy and McK-
inney (2012) studied the differences in the jet power of prograde and retrograde
black holes and found accretion efficiencies of up to ∼300% for the prograde case
and ∼80% for the retrograde case, for a black hole with |a/M |=0.9375. These sim-
ulations were the first to demonstrate that the amount of magnetic flux threading
the horizon and inner disk in the MAD state is determined by the accretion state
60
alone, not the initial magnetic conditions of the simulations, so long as sufficient
large-scale flux is available to the disk.
GRMHD simulations of non-MAD thin disks (untilted and titled) have found
roughly 10% of deviation of the accretion efficiency compared with the expected
value from the NT73 model (Shafee et al., 2008; Noble et al., 2009, 2010; Penna
et al., 2010a; Morales Teixeira et al., 2014) where this difference has been associated
with the amount of vertical magnetic flux that threads the black hole and disk.
Chapter 2 (Avara et al., 2016) built upon prior numerical work to study MAD
disks in the thin, radiatively efficient Sub-Eddington regime. We were able to
demonstrate the self-consistent inward transport of magnetic flux by a thin disk
thereby saturating into the MAD state, and characterized the radiation, wind, and
jet efficiencies, disk structure, and jet power in such a system. We found that a
sub-Eddington, moderately thin disk can be as much as twice as radiatively efficient
as the NT model would predict, with an even greater total efficiency. However,
in Chapter 2 we made use of an ad hoc cooling function to maintain a particular
disk state, and large unknowns in converting the cooling in those simulations to the
dissipation of a real system called for more realistic numerical experiments which
include self-consistent radiative transport and proper ray-tracing. The present work
provides that next step.
In this chapter we study the accretion efficiency of a thin disk in the sub-
Eddington limit with half-height (H) to radius (R) of H/R ≈ 0.1 around a black
hole with dimensionless spin (a/M)=0.5 that reaches the MAD state, potentially
maximizing accretion efficiency. We compare our values with the results from the
simulation MADiHR from Chapter 2 (Avara et al., 2016), that used an ad-hoc
cooling function to control the scale-height, in order to understand how radiation
affects the jet and wind efficiencies. Here we alternatively make use of the GRRMHD
code HARMRAD (McKinney et al., 2014) which includes radiative transport self-
consistently and solves the GRRMHD equations using the M1 closure for radiation.
The structure of the chapter is as follows: the methodology and numerical
setup are presented in §3.2, results are presented in §??, and in §?? we present our
61
conclusions and look forward.
3.2 3D GRRMHD thin disk model
We start our simulations with an accretion disk around a black hole with
MBH=1M and spin a/M=0.5 where the initial conditions of our disk follows the
NT73 solution and the density was tuned in oder to have a quasi steady-state ac-
cretion flow with M˙=0.4M˙Edd where M˙Edd is the Eddington accretion rate defined
as M˙Edd = (1/ηNT )LEdd/c
2. ηNT is the NT73 accretion efficiency and LEdd is the
Eddington luminosity defined as
LEdd =
4piGMc
κes
≈ 1.3× 1046 MBH
108M
ergs−1, (3.1)
where for our case we have ηNT ≈8.6% and M˙Edd ≈1.68×1018g s−1.
We initialized our simulations with a Keplerian disk that has a density profile
given by:
ρ(r) = ρ0e
−z2/(2H2), (3.2)
with
ρ0 = Σ/(2H), (3.3)
where Σ is the surface density, H the height, and their expressions are given by
NT73 solution for the inner (radiation-pressure dominated) region (equation 5.9.10
in NT73).
The initial radial velocity profile of the disk is given by:
Vr = αvisc(H/R)
2r|Ω|, (3.4)
where Ω is the angular velocity and αvisc was set to 0.5.
Due to the poor resolution in the radial direction at large radii we have trun-
cated the disk at Rtr1=120M , tapered off by an exponential profile for the density,
62
ρ′(r) = ρ(r)e(−r/Rtr1)+1. (3.5)
The total ideal pressure is given by Ptot = (Γtot − 1)utot, where Γtot=4/3,
utot is the internal energy density, and the gas pressure was randomly perturbed
at the 10% level to seed the MRI. The disc is initialized with an atmosphere of
density ρ = 10−5(r/Rg)−1.1, gas internal energy density egas = 10−6(r/Rg)−5/2,
and the radiation energy density and flux were set by the local thermodynamic
equilibrium (LTE) and flux-limited diffusion approximation (McKinney et al. 2014)
with a negligible radiation atmosphere. We set the rest mass and internal energy
densities to zero near the black hole within the jet and near the axis, corrected by
the numerical ceilings of b2/ρ=300, b2/egas=10
9 and egas/ρ=10
10.
The magnetic field in the disk is large scale and poloidal only, with open field
lines and a transition to monopolar at Rtr2. For r < Rtr2 the φ component of the
vector potential is given by:
Aφ = MAX((r + 5)
ν − 1040 − 0.02, 0)(sin θ)1+h, (3.6)
with ν=1 and h=4. For r ≥ Rtr2 the field transitions to monopolar using a vector
potential given by
Aφ = MAX((Rtr2 + 5)10
40 − 0.02)(sin θ)1+h(Rtr2/r). (3.7)
Our field is normalized inside one MRI wavelength per half-height H giving a ratio
of gas+radiation pressure to average magnetic pressure of β ≈10 (β = (pgas +
prad)/pmag). In order to investigate the effects of the transition to monopolar field, we
have choosen Rtr2=120Rg for our simulations RADHR and RADLR and Rtr2=200Rg
for our simulation RADvHR. The motivation for the changing the field transition
radius was to try to reach a MAD state out to large radii.
For the radiation transport processes, we have assumed that our initial disk
has solar abundance with mass fractions of hydrogen, helium and metals equal to
X=0.7, Y=0.28 and Y=0.02, which gives an electron fraction of Ye = (1 + X)/2
63
Table 3.1: Grid resolution
Simulation Polar grid Nr Nθ Nφ Tstop(Rg/c)
MADiHR Avara et al. (2016) 192 96 208 70,000
RADLR McKinney et al. (2012) 128 64 32 40,000
RADHR McKinney et al. (2012) 128 128 64 40,000
RADvHR This work 128 128 64 18,000
and mean molecular weight of µ¯ ≈0.62. The expressions for the electron scattering,
absorption-mean energy, bound-free, and free-free opacities are the same used by
McKinney et al. (2015b).
3.2.1 Numerical grid
In this section we will describe the two grids used in our simulations RADHR,
RADLR, and RADvHR. For simulations RADHR and RADLR we used the same
grid as ub McKinney et al. (2012). Our grid for RADvHR is similar, but with
important modifications that will be described below. In all of our simulations
we have adopted the same boundary conditions as McKinney et al. (2012), where
the radial grid has Nr cells spanning from Rin ∼1.5Rg ∼0.084Rh (where Rh is the
horizon radius) to Rout=10
5Rg with cell size increasing exponentially until rbreak
and then increasing hyper-exponentially. Outflow radial boundary conditions were
adopted. The θ grid has Nθ cells spanning from 0 to pi and the grid of our simulation
RADvHR has been properly adjusted for a disk with H/R=0.1 concentrating mostly
of the cells in the disk. The φ grid has Nφ cells spanning uniformly from 0 to 2pi with
periodic boundary conditions. The resolutions for the radial, θ and φ directions and
the grid type for each of our simulations are given in Table 3.1.
In the results section we will show that in our simulation RADHR the MAD
state has only built up until R=17Rg and for this reason we have modified the θ
grid in order to better resolve the MRI and raise the efficiency of the transport of
64
magnetic flux to try to reach MAD state beyond 17Rg. Our simulation RADvHR
has only been run up to a time 18,000Rg/c due to the limitations of computer
resources and algorithmic improvements in HARMRAD will aid in feasibly running
this simulation substantially longer.
Before we introduce the new equations for the new θ grid, we must define some
new variables:
Tr(x) =
e−1/x
e−1/x + e−1/(1−x)
, (3.8)
and
Trans(X,L,R) =

0.0 if x ≥ L,
1.0 if x ≥ R,
Tr
(
x−L
R−L
)
if L < x < R.
(3.9)
To concentrate most of the cells in the disk, the expressions for S0 and S2 to
control the cells at small, middle and large radii were changed. In our case S0 was
changed to:
S0(r, ra, rb) = Trans(log(r), log(ra), log(rb)), (3.10)
where the objective of S0 is to control the cell size at middle and large radii. The
expression for S2 has been changed and in our case is given by:
S2(r, ra, rb) = 1− S0(r, ra, rb), (3.11)
where this expression for S2 controls the cell size at small and middle radii. Then
the new expression for h2 is given by:
h2 = h3S2(r, 40, 200) + S0(r, 40, 200)(h
′
2(r)S2(r, 200, 500) + h
′
2(500)S0(r, 200, 500)),
(3.12)
where
65
h′2(r) = h3 +
(
r − rsj3
r0j3
)nj
, (3.13)
is the same equation for h2 in McKinney et al. (2012). The original equation for θ2
has also been modified and in or case is given by
θ′2 = θ
a
2S0(r, 20, 40) + θ
b
2S0(r, 20, 40), (3.14)
where θa2 and θ
b
2 are the original expression for θ2 and T2 presented in McKinney
et al. (2012).
Similar to the other expressions, the original expression for h0 was modified to
h0(r) = h3S2(r, 40, 200) +S0(r, 40, 200)(h
′(r)0S2(r, 200, 500) +h′(500)S0(r, 40, 200)),
(3.15)
where the expressions for h′ is now given by
h′(r) = 2−Qj(r/r1j)−nj2[(1/2)+(1/pi) arctan(r/r0j−rsj/r0j)]. (3.16)
The expression for θ1 in our case is given by
θ1 = T0S2(r, 20, 200) + θ
′
2S0(r, 20, 200), (3.17)
where T0 is the original expression from McKinney et al. (2012).
Finally the expression for θ is
θ = θa2S2(r, 1, 10) + S0(r, 1, 10)θ1. (3.18)
The parameters for the polar grid of each of ours simulations are listed in 3.2.
As in McKinney, Tchekhovskoy & Blandford (2012, 2013), McKinney et al.
(2014) and McKinney, Dai & Avara (2015), we tested the convergence quality factors
for the MRI in the θ and φ directions (Qθ,MRI , Qφ,MRI) and turbulence (Qnlm,corr)
measuring, respectively, the number of grid cells per MRI wavelength in the θ and
φ directions and the number of correlation wavelengths in the radial, θ and φ di-
rections. The values were averaged over the period 30,000-40,000Rg/c and over θ
66
and φ at r=10Rg, which represents the non-choked part of the accretion flow. Our
simulation RADHR has Qθ,MRI ∼110, Qφ,MRI ∼13 and Qnlm,corr ∼11, 13 and 4
respectively. The value of Qnlm,corr in the φ direction shows that the turbulence
is not well resolved since the previous works mentioned above have shown that 6
or more cells are required to resolve the turbulence, but even being slightly under
resolved we can still identify some physical effects as we will discuss in the next
sections. We have also measured the disk thickness per unit MRI wavelength, Sd
(where for Sd <0.5 the MRI is suppressed). Initially for our simulation RADHR,
Sd ∼0.9 and the time averaged flow has Sd ∼0.26 out to r ∼17Rg. The simula-
tion RADLR has Qθ,MRI ∼21, Qφ,MRI ∼8 and Qnlm,corr ∼12, 8, 3 where the values
were time averaged over 30,000-40,000Rg/c. Our simulation RADvHR at its end
has Qθ,MRI ∼56, Qφ,MRI ∼13 and Qnlm,corr ∼15, 22, 7 and the MAD state built up
to 13Rg. As we can see our simulation RADLR is underesolved and will be used
for reference to study resolution effects of the flow at late times and check how the
MAD state builds up while the turbulence in our simulation RADvHR was much
better resolved.
Table 3.2: Polar grid Parameters
Simulation n0 c2 n2 rs r0 h3 r0j3 rsj3 nj1 r1j roj rsj Qj nθ hθ rsj2 r0j2 rbreak
RAD(HR,vHR) 1 1 6 40 40 0.1 40 0 0.5 30 30 40 1.8 5 0.1 8 3 200
RADLR 1 1 6 100 60 0.1 40 0 0.7 30 30 40 1.7 5 0.1 5 2 300
3.2.2 Diagnostics
The diagnostics evaluated in this study are identical to those used in Chapter
2 save for a few cases, so we do not reproduce the majority of the definitions here.
However, unless otherwise stated, the half-angular thickness H normalized by the
radius R is calculated with equation 2.4 in this study as before, but with n = 2
unless otherwise specified.
With the addition of the radiation fluid, we introduce several new characteristic
67
quantities. In order to obtain the optical depth in the radial direction for each instant
in time we compute:
τ =
∫
ρκtotdl, (3.19)
where dl = −fγdr, fγ = ut(1− (v/c) cos θ), where for large radii (v/c) ∼ 1−1/(ut)2,
θ=0 and the integral is carried out from r0=3,000 (to consider only transient material
that would contribute to the optical depth, a radius where the disk wind has reached)
to r in order to obtain τt(r). In the angular direction dl = fγrdθ, θ = pi/2 and the
integral is carried out from each polar axis towards the equator to obtain τθ(θ). The
radiation photosphere is defines as τr=1, which is an upper limit to the radius of
the photosphere for an observer since radiation can escape.
Finally the radiative luminosity is given by:
L = −
∫
dAθφR
r
t , (3.20)
where the luminosity is measured at r =400Rg where only the gas has τt(r) <1 is
included.
3.3 Results
In this section we will use our three new simulations which are magnetically
arrested on the BH and in the innermost accretion flow, evolving from highly mag-
netized NT-type initial conditions. Our discussion focuses on results from RADHR
and RADvHR to study how self-consistent radiative transport affects the accretion
efficiencies when compared with the simulation MADiHR from Avara et al. (2016).
We will also show how resolution affects disk behavior.
RADHR and RADLR have been run until t=40,000Rg/c, by which time the
inner disk accretion flow reaches a MAD quasi-equilibrium, where magnetic flux
continues to accumulate at the edge of the MAD region, and the disk is slightly
sub-MAD beyond that transition radius. Figure 3.1 shows three main regimes of
behavior in RADHR, culminating in the long-term evolution of the MAD inner disk
68
illustrated in the disk cross-sections. The first stage of evolution involves the rapid
evolution of magnetic flux threading the horizon to balance the inner-most accretion
flow and dynamical effects of the plunging region inside the innermost stable circular
orbit (ISCO), as the system evolves away from equilibrium NT conditions, subject to
the dynamically-relevant large-scale polloidal magnetic flux. The second behavioral
stage involves mass accretion rate and magnetic flux growing until t ∼17,000Rg/c.
At that point there is a weak point in the inner disk and a large magnetic disruption
occurs which redistributes a lot of polloidal flux into the disk, increasing the size of
the MAD region. The final stage, covering roughly the second half of the simulation
involves a quasi-steady state MAD accretion flow, order unity oscillations of the
horizon threading flux (consistent with prior numerical findings), and a general
accretion structure very similar to the MADfHR simulation in Avara et al. (2016),
where only the horizon and inner disk reach the MAD state over the time they were
able to run their simulation.
In Figure 3.2 we show characteristic quantities that describe the quasi-steady
state of the disk. We time averaged these quantities over the last quarter of our high
resolution simulation that we were able to run for a longer time, RADHR (RADvHR
was too expensive to run longer for this work). We average over the period t=30,000-
40,000Rg/c, when the disk evolution is less subject to initial conditions different than
the unknown true thin-disk steady-state structure, for which there is no detailed
analytic radiative model. In the next sub-sections we will make a discussion of
each of these quantities that we believe best characterize the accretion flow in our
simulations, and provide a direct comparison to prior work with ad hoc cooling. We
will also describe the time dependence of some quantities, will present the results of
our highest resolution simulation, RADvHR, and describe resolution effects.
3.3.1 Accretion rate
Analyzing the mass accretion rate profile of our simulation we can observe that
the simulation has reached a quasi-steady state out to a radius of r ∼ 15Rg, but
it is clear from the radial dependence of M˙ that going further out there are larger
69
Figure 3.1: Snapshot of our simulation RADHR at the time t=40.000Rg/c where
the top two panels represents the X − Z at y=0 and X − Y at z = 0 slices of the
density with the red line showing where b2/ρ=1 and the black lines with arrows
trace field lines where thicker black lines show where field is lightly mass loaded.
The top sub-panel shows the mass accretion rate at the horizon (red line) and the
luminosity (green line). The middle sub-panel shows the magnetic flux in the horizon
normalized so that order unity is a dynamically substantial amount of magnetic flux.
The bottom sub-panels shows the accretion efficiencies in the horizon (red line), jet
efficiency (green line) measured at 15Rg and radiation (blue line) measure at 50Rg.
70
Figure 3.2: (Caption on next page.)
71
Figure 3.2: (Figure on previous page.) Characteristic quantities as a function of
radius, time-averaged over the period t =[30,000-40,000]Rg/c for the simulation
RADHR. Panels consecutively show, from top to bottom: the mass accretion rate
(M˙), energy efficiencies (η) for different parts of the stress-energy tensor, magnetic
flux (Ψ), disk scale-height (H/R), and Sd. In efficiencies plot the black line corre-
sponds to the total gas+radiative (black line), electromagnetic (green line), matter
(red) and radiation (black dashed) components. The magenta lines in the plots of
Ψ, H/R, and Sd correspond to their initial values. The thick lines show values for
RADHR and the thin line MADiHR from (Avara et al., 2016) time-averaged as in
that paper.
deviations where the disk is still evolving from initial conditions. The large-scale
magnetic flux in these types of simulations enhances these effects due to large-scale
coronal flows, large mass transfer from MRI channel modes, and initial radiative
evolution since the disk in this simulation does not start in an equilibrium state due
to the radiation. Inside the ISCO there is also a deviation in the radial dependence
of M˙ due to mass injection associated with the floors which demands a careful
analysis. It is worth noting that in these and the simulations of Avara et al. (2016)
the radial dependence of M˙ near the BH flattens when the number of snapshots
over the averaging time period increases, though it is unclear if it would improve
beyond the average shown here if more snapshots were available from the simulation
output.
The additional mass added by the floor in the simulation is removed from
the calculation of M˙ where cells have rest-mass values of b2/ρ > 1. The profile of
Avara et al. (2016) is much flatter because they time average over a later and longer
period, 30, 000 − 70, 000Rg/c in MADiHR, an ad hoc cooling function evolves the
disk toward equilibrium faster than the self-consistent processes we capture with full
radiative transport, which demands more time to reach thermal equilibrium, and
because the higher resolution allows large-scale initial transient modes to break up
faster thus having a smaller impact on later evolution (e.g. parasitic modes of the
72
MRI at intermediate and large radii).
3.3.2 Scale Height
In these simulations with self-consistent radiative transport, we start with a
NT disk with a scale height of 0.1, but there is no a priori reason to expect that the
MAD disk resulting from these initial conditions should persist at this disk width.
There was some tuning of initial conditions using test simulations to try to attain a
scale height of 0.1 but this result is entirely self-consistent with the mass accretion
rate and radiative transport alone. On the other hand, the disk scale height was fixed
directly in Avara et al. (2016). This is one of the central advantages of including self-
consistent radiative transport, which would allow an inherently thermally unstable
disk to collapse if that were the true physical behavior.
Unlike the prediction from thin-disk theory, and similar to the findings of
Mishra et al. (2016); Sa¸dowski et al. (2016), our disk is thermally stable in the
inner parts. Over the 30, 000 − 40, 000Rg/c time period we average RADHR we
find that the disk has reached thermal equilibrium only in the inner parts, not
having enough time to reach an equilibrium state further out. We determine that a
given radius has reached thermal equilibrium when that part of the disk has reached
inflow/outflow quasi-equilibrium and the scale height, in part governed by vertical
radiative transport and pressures support evolving under temperature evolution,
stabilizes to roughly a fixed value. The time evolution of H/R at a sampling of
radii is shown in Figure 3.3. By 25,000M the values of H/R in RADHR at H/R =
RH , 10Rg, 20Rg have converged to the expectation for a magnetically choked flow
(standard inner disk structure in the MAD state) where the disk has H/R ≈0.1.
The radial profile of the disk shown in 3.2 shows the inner region ‘choked’
by magnetic compression, where the disk scale height is reduced to ∼ 0.04. In the
inner edge of the disk the scale height grows to 0.1 at r=5Rg and to 0.1 by r=10Rg.
Beyond that radius we don’t show the values because the disk hasn’t reached thermal
equilibrium in this locations. For comparison, the MADiHR simulation in Avara
et al. (2016) had H/R ∼ 0.05 near the black hole and 0.1 throughout the rest
73
of the disk. Thus there is significant similarity in disk structure as a function of
radius in simulations with radiative transfer and those using only an ad hoc cooling
function. This gives additional credibility to the findings taken from the much longer
simulations with ad hoc cooling since artificial cooling and the associated choices
for cooling timescale are arguably the least trustworthy in the plunging region and
very inner disk.
In a future paper we will provide an extensive study of stability of the disk and
describe the processes by which thin disks, formally thermally unstable in standard
thin-disk theory, appear to be stabilized in the MAD state.
3.3.3 Efficiencies
Figures 3.2 and 3.5 show the radial profile of the total disk efficiency and the
different contributions to the total including all radiation (in both optically thick
and thin regions), energetics of the matter, and electromagnetic (EM) flux. Total
disk efficiency, efficiency of the accretion flow in converting gravitational potential
energy to combined forms escaping beyond the domain studied (transport to infinite
radius opens questions addressed in the conclusions), is roughly 17% out to r=20Rg
in RADvHR, ∼ 13% in RADHR, and is more constant for RADHR in the region
inside 10Rg where the disk is more securely in inflow quasi-equilibrium. While the
conversion of energy into radiation is significantly less in this simulation than was
estimated using the ad hoc cooling procedure for MADiHR, there is some additional
contribution by the matter component, and possibly the EM component beyond
the choked region of the inflow. RADvHR has significantly more outward going
EM energy flux, but it is unclear if this enhancement is due to better resolution
in RADvHR or if it is a transient that would disappear if the simulation were ran
the same length as RADHR and RADLR. If one considers only the radiation that
reaches the optically thin region and would travel outward unimpeded, the radiative
efficiency of our MAD thin disk is even less, only ∼ 2%. Thus, the total efficiency
of this type of disk is significantly enhanced beyond the expected NT value but the
energy is transported outward in the form of winds and jets predominantly rather
74
than radiation.
For our system with spin a/M = 0.5, the total efficiency of RADHR is 18.6%
and the corresponding value from NT73 is 8.2%, a deviation of 125%. At the hori-
zon, the efficiency in matter, electromagnetic and radiation forms are ηEMH =1.5%,
ηMAKEH =22.5% and η
Rad
H =-5.5%, respectively. At r=50Rg and we found that the jet
is carrying 4.28%, the wind 7.39%, and the radiation 2.94%, where for the jet effi-
ciency 4.09% is in the electromagnetic form and 0.20% in the matter, while for the
wind 4% is being carried by the electromagnetic part and 3.39% in the matter. The
summed values of the efficiencies carried by the jet, wind and radiation is less than
the total efficiency and this difference is because even at R=50Rg it’s not trivial
to guarantee that the wind material will reach infinity and for this reason we have
assumed that the wind efficiency is ηwind = ηBH − ηjet − ηRad; using this condition
we found that the wind has an efficiency of 11.4%.
3.3.4 Magnetic flux
Measuring the absolute flux threading the horizon as a function of time,
demonstrated in Figure 3.3, we find that the black hole gains a significant quantity
of polloidal magnetic flux across the first half of the simulation, in the initial transi-
tory phase due to large scale MRI channel modes that survive until the φ-symmetry
of the field falls subject to parasitic modes and non-linearity. Then, similar to the
MADiHR in Chapter 2, during the early evolution we observe the two modes of flux
accumulation as it is transported both along the disk surfaces through a coronal-
type mechanism and through disk advection. Similar to the build-up of flux from
sub-MAD to MAD conditions on the horizon seen in MADfHR there is a steady
build-up of flux on the horizon in RADHR until at ∼ 17, 000Rg/c there is more flux
on the horizon than the disk can contain, and a significant portion of the flux on the
horizon finds a low pressure, ‘weak’ spot in the previously rather axisymmetric disk
and reconnects, pushing out into the disk through that point. Future work will be
needed to determine if the velocity profile of the disk over some time-periods renders
it stable to the magnetic Rayleigh-Taylor instability, which could also contribute to
75
determining when the field will break free from the horizon. Whatever the process,
from that moment on the more highly magnetized disk inside r ∼ 10Rg is MAD.
Due to constraints imposed by having NT-type initial conditions necessary
to study the disk with self-consistent radiative transport, and having attainable
resolution for the outer disk in our simulations, a significant portion of the polloidal
magnetic flux in the outer disk is lost as the system evolves through initial transients.
Despite this loss of flux in the disk, analyzing the suppression factor (Sd) of RADHR
reveals that at late times the disk reaches the MAD state out to r = 15Rg. As
in Avara et al. (2016) this radius grows over time as disk and coronal processes
transport more magnetic flux inward.
When the net vertical magnetic flux threading any given radius in the disk
increases as flux is transported inward (or alternatively decreases if it is diffusing
outward), we expect the value of αmag to increase at that radius. Once a section
of the disk has reached the MAD state, however, the value of αmag is expected to
saturate at a maximum and this can be seen in the evolution of αmag at 10Rg in
Figure 3.3, where it’s possible to see that by the time t=20,000Rg/c the value of
effective αmag rises and saturates to be roughly constant for the remainder of the
simulation.
3.3.5 Luminosities
In Figure 3.4 we use the time-φ averaged structure of the radiative and EM
luminosities to demonstrate the angular distribution of this outgoing energy, the
dissipate locations where radiation originates in the accretion inflow/outflow, and
the distribution of reprocessing. Most of the EM flux comes from the jet in the
south which is less isotropically distributed in θ, while the jet pointing to the north
direction shows a more uniform distribution. The electron scattering surface has a
large opening angle due to the substantial quantity of material blocking the radiation
at large distances.
The radiation luminosity along jets in both directions also exhibits a non-
uniform distribution having a more enhanced emission in the north side in a opening
76
Figure 3.3: Red is for RADLR, black for RADHR. Top panel: Disk scale-height as
function of time. The solid line corresponds to the values for 10Rg (lines for all radii
are smoothed, with only the 10Rg line duplicated without smoothing in gray), the
dotted and dashed lines correspond to the values of H/R at the horizon and 20Rg
respectively. Middle: Magnetic α measured at 10Rg (full, un-smoothed RADHR
data in gray). Bottom: Amount of magnetic flux on the horizon normalized by
the amount of magnetic flux contained in the portion of the disk which reaches one
inflow time by the end of the simulation, i.e. the total magnetic flux available for
the inner disk and horizon.
77
Figure 3.4: Top: Electromagnetic luminosity per unit angle (∂θLEM)/(M˙Hc
2) with
time-φ averaged magnetic flux (translucent gray lines); blue lines show the b2/ρ=1
boundary; the green line shows the ur=0 boundary; yellow is the surface where
optical depth reaches unity to electron scattering, integrating inward from a large,
but at least moderately evolved, radius. Bottom: lab-frame radiation flux and
radiation luminosity.
78
angle roughly similar to the EM emission. The disk seems to spontaneously break
the north-south symmetry of the emission and then the larger radiation flux to
the north may maintain that asymmetry for long periods of time. It is not clear
if the asymmetry would persist if the RADHR simulation was run longer, but we
see this present to some degree in each of our GRRMHD simulations and a weaker
asymmetry was seen in some simulations of Avara et al. (2016) as well, suggesting
the magnetic wind may be playing a role in the observed behavior.
For both EM and radiation streaming, it is possible to see that the emission
largely originates from the region inside r=15Rg.
3.3.6 Resolution effects
Figure 3.5: Plot showing the efficiencies for the simulation RADHR (dot-dashed
line), RADvHR (dashed line), RADLR (dotted line) and MADiHR (continuous
line). Black: total, Red: PAKE, Green: EM, Cyan: Radiation.
In this section we will demonstrate how the resolution affects the disk thick-
ness, efficiencies, accretion rate, stress and the MAD state.
The disk thickness at the horizon for simulation RADHR(LR) is 0.04 (0.03),
at r=5Rg is 0.1 (0.08), and at r=20Rg is 0.12 ( 0.15). When comparing the values at
79
the horizon it’s possible to see that the disk thickness may be sensitive to resolution:
low resolution may lead to enhanced dissipation through radiation, with enhanced
escape potential, thereby cooling the disk further and into a thinner distribution.
where the disk hasn’t reached thermal equilibrium at large radii.
When comparing Sd for our simulations RADHR and RADLR, we found that
RADHR has reached the MAD state out to r=16Rg while RADLR is MAD out
to r=13Rg. This difference can be explained considering resolution of the fastest
growing MRI modes, which are important to the processes of magnetic field dynamo
and transport in the part of the disk that is sub-MAD. Our simulation RADHR has
Qθ,φ=110,13 while RADLR has Qθ,φ=21,8 showing that the MRI is not resolved for
our simulation S2, affecting this transport of magnetic flux.
In Figure 3.3 one can see that the higher resolution of RADHR makes the time
evolution of the disk scale height smoother and more magnetic flux is held by the
horizon and inner disk.
Since the numerical resolution seems to affect the transport of magnetic flux
we might expect our low resolution simulation RADLR to have lower efficiencies in
the horizon and jet. For the jet we expect a lower efficiency since less rotational
energy is extracted from the black hole while for the jet and radiation at r=50Rg
we must expect the efficiencies to be larger since the wind and radiation in this case
absorb more energy. As expected the jet efficiency is 3.06%, compared to 4.28% for
RADHR, the wind has 15.8% (7.39% for RADHR) and the radiation 4.58% (2.94%).
The most robust measure of resolution effects in our simulation suite is the
monotonically increasing total efficiency with increasing resolution seen in Figure
3.5. While we expect the EM efficiency component of RADvHR to likely drop to
a value consistent with other simulations if we were to run it longer, even then, it
would have a higher total efficiency compared to RADHR and RADLR respectively.
Thus, the true total efficiency of a MAD thin disk is likely underestimated by our
simulations. While the estimate for radiative efficiency in MADiHR is found to be
an overestimate in light of this work, it is worth noting that the resolution was
substantially higher, which would be consistent with the trend seen in these new
80
simulations.
3.4 Discussion and Conclusions
We have shown the results of the first known simulations of radiatively efficient
thin disks, with H/R ∼0.1, around a spinning black hole with a/M=0.5, that are
in the MAD state and include self-consistent radiative transport. Our longer high-
resolution simulation RADHR reached the MAD state out to 17Rg and most energy
extracted by the disk is carried away by the wind consistent with the findings in
Chapter 2 but with a smaller total efficiency due to a significantly smaller radiation
component. This is because our calculations include self-consistent absorption and
scattering of the radiation, advection by plunging material, and true ray-tracing. We
measured the electromagnetic and radiation luminosities and we found anisotropic
emission in both cases. The electromagnetic emission is dominated by the southern
jet/wind, with the jet in the north direction having more uniform distribution.
Radiation, on the other hand, prefers emission to the north.
If one considers the radius at which radiation dissipated in the central disk
no longer has a time-averaged path leading to the infinite observer, and is instead
captured by the horizon (whether it streams away from the disk material first or
not), one can make a more accurate estimate for a radiation inflow/outflow tran-
sition radius than the simple photon orbit chosen in the previous chapter. For
RADHR this is ∼ 4Rg. If only radiation outside that radius is considered from
MADiHR, the disk has a total radiative efficiency of ∼ 6 − 7%. Allowing for some
additional dissipation of radiative energy by the wind component, this could easily
match the 8% NT prediction. The additional 4% discrepancy between RADHR and
MADiHR is primarily from contributions between the ISCO and ∼ 10Rg. This is
where additional behavior of the energetics of the flow, advection, choking/magnetic
compression, and the atypically fast cooling rate chosen in Chapter 2 could account
for the difference. Because the geometry, magnetization, and resolution of the ac-
cretion flow in this region is not identical in RADHR and MADiHR, we are not yet
81
able to determine if the 4% difference is inherent to radiative transport, or due to
different disk initial conditions and simulation parameters. Future simulations will
shed light on this difference.
One should note that the radiative and magnetic energy converted from grav-
itational potential energy of the accretion flow has above-NT contributions every-
where the disk is MAD in RADHR and MADiHR. This could imply that a disk
MAD out to larger radii than studied here would have increasing efficiency with in-
creasing MAD/non-MAD transition radius. However, this is somewhat speculative
and requires further numerical exploration.
A key result of this work that will also be explored in the future is the thermal
stability of the disk, likely a result of the very strong magnetic fields intrinsic to
the MAD state. It remains to be shown if there is a qualitative difference in disk
stabilization by a strong magnetic field in a standard disk that is MRI dominated,
compared to a disk in the MAD state where the MRI no longer dominates the
accretion process.
82
Chapter 4: Efficiency of Super-Eddington Magnetically-Arrested Ac-
cretion
Abstract
The radiative efficiency of super-Eddington accreting black holes (BHs) is
explored for magnetically-arrested disks (MADs), where magnetic flux builds-
up to saturation near the BH. Our three-dimensional general relativistic radi-
ation magnetohydrodynamic (GRRMHD) simulation of a spinning BH (spin
a/M = 0.8) accreting at ∼ 50 times Eddington shows a total efficiency ∼ 50%
when time-averaged and total efficiency & 100% in moments. Magnetic com-
pression by the magnetic flux near the rotating BH leads to a thin disk,
whose radiation escapes via advection by a magnetized wind and via trans-
port through a low-density channel created by a Blandford-Znajek (BZ) jet.
The BZ efficiency is sub-optimal due to inertial loading of field lines by opti-
cally thick radiation, leading to BZ efficiency ∼ 40% on the horizon and BZ
efficiency ∼ 5% by r ∼ 400rg (gravitational radii) via absorption by the wind.
Importantly, radiation escapes at r ∼ 400rg with efficiency η ≈ 15% (lumi-
nosity L ∼ 50LEdd), similar to η ≈ 12% for a Novikov-Thorne thin disk and
beyond η . 1% seen in prior GRRMHD simulations or slim disk theory. Our
simulations show how BH spin, magnetic field, and jet mass-loading affect the
radiative and jet efficiencies of super-Eddington accretion.
4.1 Introduction
Black hole (BH) accretion drives a broad range of phenomena, including
quasars, active galactic nuclei, BH X-ray binaries, some gamma-ray bursts, and tidal
disruption events (TDEs). The disk’s gravitational potential energy and BH spin
83
energy are converted into radiation, winds, and relativistic jets via magnetic stresses
generated by the magneto-rotational instability (MRI) (Balbus and Hawley, 1998a),
magnetic field threading the disk (Blandford and Payne, 1982), or magnetic field
threading the BH (Blandford and Znajek, 1977). Magnetic torques are maximized
when magnetic flux advects inward and piles-up until magnetic stresses balance in-
coming gas forces – the so-called magnetically-arrested disk (MAD) state (Narayan
et al., 2003). Simulations of non-radiative MADs show up to ∼ 300% total efficiency
(Tchekhovskoy et al., 2011; McKinney et al., 2012), corresponding to 3 times the
energy out as going into the BH. A fundamental question is: What form does this
energy take at large distances and what fraction is radiative?
The radiative efficiency of BH systems is expected to be dependent upon the
accretion rate M˙0, which controls the density, disk thickness, and dynamical im-
portance of radiation (from gas-pressure dominated at low M˙0 to radiation-pressure
dominated at high M˙0). To scale M˙0c
2 or the radiative luminosity L, for a black
hole mass M , speed of light c, gravitational constant G (giving gravitational ra-
dius rg ≡ GM/c2), and Thomson electron scattering opacity κes, one can use the
Eddington luminosity
LEdd =
4piGMc
κes
≈ 1.3× 1046 M
108M
erg s−1. (4.1)
One can also choose to normalize M˙0 by M˙Edd = (1/ηNT)LEdd/c
2, where ηNT is
the nominal accretion efficiency for the Novikov-Thorne thin disk solution (Novikov
and Thorne, 1973a) (commonly, a fixed ηNT = 0.1 is used, but we include the spin
dependence).
For luminosities L & 0.3LEdd, the accretion flow is expected to become ge-
ometrically thick and optically thick, and in this regime the photons can remain
trapped within the flow as in the “slim disk” model (which includes no magnetic
field) (Abramowicz et al., 1988). The super-Eddington accretion regime where
M˙ & M˙Edd or L & LEdd may help explain ultra-luminous X-ray sources as highly
super-Eddington stellar-mass BHs (Watarai et al., 2001; Miller and Colbert, 2004).
Also, some black hole X-ray binaries can have L & LEdd (e.g., SS433, Margon et al.
84
1979; Takeuchi et al. 2010; GRS1915+105, Fender and Belloni 2004), while tidal
disruption events seem to require M˙  M˙Edd (Bloom et al., 2011).
In this paper, we seek to test whether super-Eddington MADs around spinning
BHs are radiatively efficient. We use the general relativistic radiation magnetohy-
drodynamic (GRRMHD) code HARMRAD, which uses the M1 closure for radiation
(McKinney et al., 2014). The slim disk model and recent GRRMHD simulations of
non-MAD (i.e. lower magnetic flux than MAD on the BH and in the disk) super-
Eddington flows are radiatively inefficient with L . LEdd even for M˙ ∼ 100M˙Edd
(Sa¸dowski et al., 2014; McKinney et al., 2014), far below the thin disk Novikov-
Thorne (NT) efficiency (e.g. about 6% efficient for BH spin of a/M = 0 to 12%
efficient for a/M = 0.8). However, MADs (not yet studied with radiation) can
maximize the efficiency and may lead to a radiatively efficient state. This is also
plausible for rapidly spinning BHs because MADs are then compressed into a thin
disk (McKinney et al., 2012), which can lose radiation more rapidly.
We discuss the physical and numerical setup in §4.2. Results and discussions
are presented in §4.3. We summarize in §4.4.
4.2 Fully 3D GR Radiative MHD MAD Model
This model with black hole mass of M = 10M and dimensionless spin a/M =
0.8 has nominal thin disk efficiency ηNT ≈ 12.2%, so that M˙Edd ≈ 1.2× 1019g/s (see
Eq. 4.1). While we choose a specific BH mass, the flow is strongly electron scattering
dominated and so the results might apply roughly equally to all BH masses.
4.2.1 Opacities
We assume solar abundances (mass fractions of Hydrogen, Helium, and “met-
als”, respectively, X = 0.7, Y = 0.28, Z = 0.02) giving electron fraction Ye =
(1 +X)/2 and mean molecular weight µ¯ ≈ 0.62, entering gas entropy, pressure, and
temperature Tg [Kelvin].
We use a frequency (ν) mean of the opacity αν to get an absorption-mean ab-
85
sorption opacity (units of [cm2/g]) of κ = (
∫
ν
dνανJν)/(
∫
ν
dνJν), where the absorbed
radiation Jν is assumed to be a Planck distribution at a radiation temperature of
Tr = (Eˆ/arad)
1/4 [Kelvin], where Eˆ = uµuνRµν , u
µ is fluid 4-velocity, R is radiation
stress-energy tensor, and arad is the radiation constant.
The electron scattering opacity is
κes ≈ 0.2(1 +X)κkn, (4.2)
where the Klein-Nishina correction for thermal electrons is κkn ≈ (1 + (Tg/(4.5 ×
108))0.86)−1 (Buchler and Yueh, 1976).
The absorption-mean energy absorption opacity is
κabs ≈
(
1
κm + κH−
+
1
κChianti + κbf + κff
)−1
, (4.3)
which bridges between different temperature regimes. The molecular opacity is κm ≈
0.1Z, and the H− opacity is κH− ≈ 1.1×10−25Z0.5ρ0.5T 7.7g . Rest-mass density, ρ, is in
cgs units. The Chianti opacity accounts for bound-free at slightly lower temperatures
and is given by κChianti ∼ 4.0 × 1034ρ(Z/Zsolar)YeT−1.7g T−3.0r (this accounts for the
assumed Z = Zsolar = 0.02 for figure 34.1 in Draine (2011), most applicable for
baryon densities of nb ∼ 1cm−3). The bound-free opacity is
κbf ≈ 3× 1025Z(1 +X + 0.75Y )ρT−0.5g T−3.0r ln (1 + 1.6(Tr/Tg)) , (4.4)
(Rybicki and Lightman, 1986), where the 1 + X + 0.75Y term is roughly accurate
near solar abundances, and the ln() term comes from fitting the absorption-mean
integral. The free-free opacity is
κff ≈ 4× 1022(1 +X)(1− Z)ρT−0.5g T−3.0r ln (1 + 1.6(Tr/Tg))
× (1 + 4.4× 10−10Tg) , (4.5)
for thermal electrons and no pairs (see Eq.5.25 in Rybicki and Lightman 1986 and
Shu 1991). Thermal energy Comptonization is included as in Sa¸dowski et al. (2015).
The mean emission opacity κemit is the same as κabs but letting Tr → Tg, such that
Kirchoff’s law gives an energy density emission rate of λ = cρ0κemitaradT
4
g . The total
opacity is κtot = κes +κabs. The low-temperature opacities avoid unphysical opacity
divergences during the simulation.
86
4.2.2 Initial Conditions
The initial disk is Keplerian with a rest-mass density that is Gaussian in angle
with a height-to-radius ratio of H/R ≈ 0.2 and radially follows a power-law of
ρ0 ∝ r−0.6. The solution near and inside the inner-most stable circular orbit (ISCO)
is not an equilibrium, so near the ISCO the solution is tapered to a smaller density
(ρ0 → ρ0(r/15)7, within r = 15rg) and a smaller thickness (H/R → 0.2(r/10)0.5,
within r = 10rg – based upon a low-resolution simulation). The total internal energy
density utot is estimated from vertical equilibrium of H/R ≈ cs/vK for sound speed
cs ≈
√
ΓtotPtot/ρ0 with Γtot ≈ 4/3 and Keplerian speed vK ≈ (r/rg)/((r/rg)3/2 +
a/M). The total ideal pressure Ptot = (Γtot − 1)utot is randomly perturbed by
10% to seed the MRI. The disk gas has Γgas = 5/3. The disk has an atmosphere
with ρ0 = 10
−5(r/rg)−1.1 and gas internal energy density egas = 10−6(r/rg)−5/2.
The disk’s radiation energy density and flux are set by local thermal equilibrium
(LTE) and flux-limited diffusion (McKinney et al., 2014) with a negligible radiation
atmosphere.
We do not use polish donuts (Abramowicz et al., 1978) or equilibrium tori
as initial conditions. The outer parts of tori have an extended column of gas at
high latitudes that contributes significantly to spurious luminosity (see section 6.8
in McKinney et al. 2014), trapping of the disk wind and radiation, and artificial
(instead of the self-consistent wind) collimation of the jet.
The initial magnetic field is large-scale and poloidal. For r < 300rg, the
coordinate basis φ-component of the vector potential is
Aφ = MAX(r
ν1040 − 0.02, 0)(sin θ)1+h, (4.6)
with ν = 1 and h = 4. For r ≥ r0 = 300rg, the field transitions to monpolar using
Aφ = MAX(r
ν
010
40 − 0.02, 0)(sin θ)1+h(r0/r). The field is normalized with ∼ 1 MRI
wavelength per half-height H giving a ratio of average gas+radiation pressure to
average magnetic pressure of β ≈ 24 for r < 100rg.
87
4.2.3 Numerical Grid and Density Floors
The numerical grid mapping equations and boundary conditions used here are
identical to that given in McKinney et al. (2012). The radial grid of Nr = 256
cells spans from Rin ≈ 0.688rH (horizon radius rH) to Rout = 105rg with cell size
increasing exponentially till r ∼ 500rg and then even faster. Radial boundaries
use absorbing conditions. The θ-grid of Nθ = 128 cells spans from 0 to pi with
mapping given in McKinney et al. (2012) but with njet = 0.7 to follow the jet.
Other coefficients are tuned so that the grid aspect ratio at r ∼ 30rg near the
equator is 1 : 2 : 3 in r, rdθ, r sin θdφ. Approaching the horizon the disk thins due to
magnetic compression, so we refine the grid to have 20 points across a half-height
thickness H/R ∼ 0.1 on the horizon. The Poynting-dominated polar jet contains,
respectively, about 80, 60, and 110 θ grid cells at r ∼ rH, r ∼ 20rg, and r & 500rg.
This gives sufficient resolution of the Poynting-dominated jet. The φ-grid of Nφ = 64
cells spans uniformly from 0 to 2pi with periodic boundary conditions.
As in our other papers (McKinney et al., 2012, 2013, 2014), we test the so-
called convergence quality factors for the MRI in the θ and φ directions (Qθ,MRI
and Qφ,MRI) and turbulence (Qnlm,cor) measuring, respectively, the number of grid
cells per MRI wavelength in the θ and φ directions and the number of grid cells per
correlation length in the radial, θ and φ directions. At late times, our simulation
has Qθ,MRI ∼ 340, Qφ,MRI ∼ 45, and rest-mass and magnetic energy densities have
Qnlm,cor ∼ 25, 20, 6 at r ∼ 8rg, indicating good r, θ resolution and marginal φ res-
olution. We also measure the thickness of the disk per unit MRI wavelength, Sd,
where Sd < 0.5 is where the MRI is suppressed. Initially Sd ∼ 0.6 at all radii, while
the time-averaged flow has Sd ∼ 0.1 out to r ∼ 60rg.
The rest-mass and internal energy densities are driven to zero near the BH
within the jet and near the axis, so we use numerical ceilings of b2/ρ0 = 300,
b2/egas = 10
9, and egas/ρ0 = 10
10. The value of b2/ρ0 is at the code’s robustness
limit for the chosen resolution.
88
4.2.4 Diagnostics
The disk’s geometric half-angular thickness (H) per radius (R) is
H
R
(r, φ) ≡ H0
R
+
(∫
θ
ρ(θ − θ0)ndAθφ
)1/n(∫
θ
ρdAθφ
)1/n , (4.7)
with n = 2, H0/R = 0, and surface differential dAθφ. One computes θ0 like H/R,
but let n = 1 and {θ0, H0/R} = pi/2.
The mass accretion rate and energy efficiency are, respectively,
M˙0 =
∣∣∣∣∫ ρ0urdAθφ∣∣∣∣ , (4.8)
η = −
∫
(T rt + ρ0u
r +Rrt )dAθφ
[M˙0]H
, (4.9)
where T is the gas stress-energy tensor and [M˙0]H is the time-averaged M˙ on the
horizon. The specific angular momentum accreted is  = (
∫
(T rφ +R
r
φ)dAθφ)/[M˙0]H .
Both η and  are composed of free particle (PAKE=kinetic+gravitational), thermal
(EN), electromagnetic (EM), and radiation (RAD) components. The jet is defined
as PAKE+EN+EM, in locations where magnetic energy density exceeds rest-mass
energy density. The wind is defined as PAKE+EN, located outside the jet. The
wind can also contain untapped EM and RAD components. The dimensionless
magnetic flux (McKinney et al., 2014) is
Υ ≈ 0.7
∫
dAθφ0.5|Br|√
[M˙0]H
, (4.10)
for radial magnetic field Br in Heaviside-Lorentz units.
To obtain the optical depth, at each instant in time we compute
τ ≈
∫
ρ0κtotdl. (4.11)
For the radial direction, dl = −fγdr, fγ ≈ ut(1− (v/c) cos θ), (v/c) ≈ 1−1/(ut)2 (as
valid at large radii), θ = 0, and the integral is from r0 = 3000 (being some radius
beyond which only transient material would contribute to the optical depth, but a
radius the disk wind has reached) to r to obtain τr(r). For the angular direction,
89
dl = fγrdθ, θ = pi/2, and the integral is from each polar axis toward the equator to
obtain τθ(θ). The flow’s “true radiative photosphere” is defined as when τr = 1, a
conservative upper limit to the radius of the photosphere for an observer, because
radiation can escape by tracking with relativistic low-density parts of the jet.
The radiative luminosity is computed at each instant as
L = −
∫
dAθφR
r
t , (4.12)
which is usually measured at r = 400rg, where we only include those angles where
the gas has τr(r) < 1.
4.3 Results and Discussion
Fig. 4.1 and Fig. 4.2 show the initial and final state of the accretion flow.
The initial disk is threaded by a weak large-scale poloidal magnetic field around
a spinning black hole with a/M = 0.8. Rotation amplifies the magnetic field via
the MRI leading to accretion of mass, energy, angular momentum, and magnetic
flux. The magnetic flux accumulates and eventually forms a quasi-steady super-
Eddington (M˙0 ∼ 400LEdd/c2 ∼ 50M˙Edd) state that is a MAD (where the MRI is
suppressed) out to r ∼ 60rg after the model is evolved for a time 31, 200rg/c. The
magnetically-compressed radiatively efficient thin disk near the BH is exposed by
the jet channel, and the magnetized wind carries a significant fraction of radiation
away from the disk.
Fig. 4.3 shows the mass accretion rate, efficiencies, disk thickness, and mag-
netic flux as time-averaged from t = 30, 000rg/c to 31, 200rg/c. The disk is in inflow
and thermal equilibrium out to r ∼ 20rg with constant fluxes of mass, energy, and
specific angular momentum vs. radius. The total (kinetic + gravitational + elec-
tromagnetic + radiative) efficiency η ≈ 50%. As in non-radiative MAD models
(Tchekhovskoy et al., 2011; McKinney et al., 2012), the efficiency is beyond the
Novikov-Thorne value of ηNT ≈ 12%.
The magnetic field threading the BH is strong, with Υ ≈ 10 at late time and
Υ ≈ 8 on average with Bz ∝ r−5/4 through the disk. The magnetic field threading
90
Figure 4.1: The initial (t = 0) state consists of a weakly magnetized radially-
extended vertically-Gaussian disk with H/R ∼ 0.2 around a spinning (a/M = 0.8)
BH. Rest-mass density is shown as color with legend, while green lines show magnetic
field lines. Density values are made dimensionless using M˙Edd, length rg, and time
rg/c. The initial disk is threaded by weak (but ordered) magnetic flux, capable of
accumulating onto the BH and leading to the MAD state.
the BH and disk leads to a magnetized wind that carries a significant amount of
radiation away from the disk and helps to avoid the classical photon trapping effect
of slim disks (Abramowicz et al., 1988). The magnetized wind at r ∼ 400rg has a
matter efficiency of about 15%, but the optically thick wind also contains about 15%
trapped radiation energy and 10% electromagnetic energy that could be tapped for
acceleration or heating. By r ∼ 400rg, the wind carries most of its radiation within
a half-opening angle of 30◦ around the polar axes.
The Blandford-Znajek (BZ) effect leads to an electromagnetic, BZ, efficiency
η ∼ 35% on the BH. Some of that energy forms a jet with η ∼ 10% at r = 50rg, but
some of that energy is absorbed by the wind, leading to a jet having only η ∼ 5%
by r = 400rg.
Radiative-loading of magnetic field lines threading the BH leads to a lower-
than-optimal BZ efficiency (relatedly, see Beskin et al. 2004; Takahashi and Ohsuga
2015) once the MAD state builds-up. The inertial loading of magnetic field lines is
due to a total energy density ρtot = ρ0 +egas + Eˆ once τ & 1 caused by ρ0 keeping up
with the MAD’s higher b2 when the numerical ceiling of b2/ρ0 = 300 is enforced. The
91
Figure 4.2: The evolved (t = 31, 200rg/c) super-Eddington MAD state, with the
figure split at x = 0 as two panels (x < 0 on left and x > 0 on right). Left
panel shows radiation-frame radiation energy density (color, with legend), radiation
velocity lines (black, fixed line width), and optical depth of unity away from each
polar axis, τθ = 1 (cyan lines). Right panel shows fluid-frame rest-mass density
(color, same legend), magnetic field lines (black, thicker lines for more magnetized
gas), where magnetic energy is equal to rest-mass energy density (red lines), and
same optical depth (cyan lines). The MAD state reaches a quasi-steady state out
to r ∼ 20rg. Radiation advects inward within the equatorial disk, outward through
the low-density jet channel, and outward along with the optically thick wind.
92
Figure 4.3: The time-averaged angle-integrated mass flux, efficiency, disk thickness,
and magnetic flux. From top to bottom, panels are: Total mass accretion rate (M˙0)
per unit M˙Edd out to where there is inflow and thermal equilibrium (r ∼ 20rg, beyond
which the line is truncated), energy efficiency η (total: solid black line, electromag-
netic: blue dashed, matter without rest-mass (i.e. kinetic+gravitational+thermal):
dark red dotted, radiation (negative within 3rg): red long-dashed, radiative luminos-
ity L: magenta dot-dashed), disk thickness H/R (evolved state: black solid, initial
state: blue dashed), and magnetic flux (total on BH and threading disk: black
solid, only threading disk: dark red dotted). The disk is magnetically compressed
near the BH, and the thin disk generates a significant radiative flux that effectively
releases starting at r ∼ 300rg. The initial disk with H/R ≈ 0.2 has puffed up to
H/R ≈ 0.3 from r ∼ 20–70rg after evolving for several thermal times. By r ∼ 400rg,
the luminosity reaches L ∼ 50LEdd.
93
optically thick radiation slows the jet magnetic field line rotation rate down by order
unity compared to the optimal BZ value, because b2/ρtot ∼ 1 in the funnel (even if
b2/ρ0  1 there) as determined by the rough condition that the jet and disk have
similar b2, yet Υ 1 implies that b2 ∼ ρ0 in the disk, and radiative energy density
Eˆ ∼ ρ0 (from a thermal estimate of H/R ∼ cs/vK ∼ 0.3). A restarted simulation
with an exponentially reduced (as b2/ρ0 approached its ceiling) opacity has a BZ
efficiency of η ∼ 100% on the BH, as consistent with Υ ≈ 8 (Tchekhovskoy et al.,
2011). A restarted simulation with a lower ceiling b2/ρ0 = 100 led to a 20% BZ
efficiency. So, mass-loading and opacity physics in the funnel are important to BZ-
driven jets and winds for super-Eddington flows where radiation crosses magnetic
field lines.
Fig. 4.4 shows the jet, wind, and “true” photosphere at large radii. The
jet-wind boundary oscillates in angle to order unity, causing the wind to absorb
some of the jet and radiation energy. We have resolved the true photosphere at
r ∼ 400rg, and the radiation escaping to an observer has a high radiative efficiency
of η ≈ 15% (comparable with the NT value of η ∼ 12%), corresponding to luminosity
L ∼ 50LEdd. Most of this radiation is within a half-opening angle of 15◦ around
the polar axes. The luminosity in this GRRMHD MAD simulation is much higher
than predicted by the slim disk solution at η ∼ 2% or seen in other viscous non-
MAD or non-MAD MHD simulations at η . 1% (Ohsuga et al., 2005; Ohsuga and
Mineshige, 2011; Sa¸dowski et al., 2014; McKinney et al., 2014). Far beyond the
photosphere, radiation might be absorbed (Sadowski and Narayan, 2015) or some
might leak out of the wind containing η ≈ 15% in trapped radiation, but our grid is
too unresolved to accurately track energy conversion for r > 600rg. The radiation
ultimately released is high in our simulation because 1) the disk near the BH is
forced to become thin; 2) a magnetized wind pulls radiation-filled material off the
disk; 3) the jet drills a channel and pushes back the opaque disk wind; and 4) the
jet becomes conical at large radii allowing free-streaming of radiation. These lead
to an effective diffusion timescale shorter than in the slim disk model.
We also performed otherwise identical simulations that are MAD with a/M =
94
0 200 400 600 800 1000
x[rg ]
0
500
1000
1500
2000
2500
3000
z[
r g
]
10−10.0
10−9.0
10−8.0
10−7.0
10−6.0
10−5.0
10−4.0
10−3.0
Figure 4.4: The evolved (t = 31, 200rg/c) super-Eddington MAD state at large radii
showing the jet, wind, and photosphere. Shows fluid-frame rest-mass density (color,
with legend), where magnetic energy equal to rest-mass energy density (red lines),
radiative photosphere at τr = 1 (black line), and the wide-angle wind component’s
magnetic field line, within which (toward the polar axis) the outflow has achieved
more than a single flow time along the line extending ∼ 1000rg (blue line). Quan-
tities are φ-averaged. The jet extends to r ∼ 30, 000rg with Lorentz factor γ ∼ 5
by r ∼ 103rg. The wind extends to r ∼ 3000rg, with most of the wind’s trapped
radiation having an opening angle of 30◦ at r ∼ 400rg. The disk at R 100rg has
not evolved much. Radiation within r . 400rg can more freely stream within the
jet’s low-density channel. By t ∼ 30, 000rg/c, at r ∼ 400rg within 15◦ around the
polar axis, the outflow become conical and optically thin.
95
0 as well as non-MAD simulations using a toroidal field in the initial disk (with
β ∼ 20 for all radii at the equator). The non-spinning MAD never produced a jet and
the radiation remains trapped by the broad disk wind. Similarly, the toroidal field
models with a/M = 0 and a/M = 0.8 show no low-density funnel region. Similar
to Sadowski and Narayan (2015), our a/M = 0 MAD model with M˙ ∼ 20M˙Edd
has a wind with NT-level efficiency of order 5% and negligible radiation flux at
large radii by r ∼ 400rg within the photosphere, because all radiation energy flux
converted to kinetic energy flux. Our toroidal field models with a/M = 0 and
a/M = 0.8 have M˙ ∼ M˙Edd with a wind efficiency of 2% and 4%, respectively, and
a radiative efficiency of about 1%. So the formation of a low-density jet channel by
a spinning BH, with enough magnetic flux to launch a jet and strong wind, helps
super-Eddington accretion become radiatively efficient.
Convergence testing was performed by restarting each previously-mentioned
model at a resolution twice lower in each dimension. The restarted simulation starts
at half-way through the higher resolution simulation, and then the two resolutions
are compared for the latter half of the simulations. We find that the total, BH,
jet and radiative efficiencies, and Υ agree to within 30% or smaller error. Our
convergence quality factors for the MRI and turbulence also suggest the simulations
are converged.
A non-relativistic MHD simulation with a/M = 0 by Jiang et al. (2014) mea-
sured a NT-level radiative efficiency of about 5% that they attributed to MRI-driven
magnetic buoyancy vertical transport of radiation. Their simulation had a small ver-
tical extent ±60rg, so τ computed from Eq. 4.11 does not include the extended wind.
Over extended distances, a significant portion of radiation energy flux can convert
into wind energy flux (Sadowski and Narayan, 2015). In our a/M = 0.8 MAD simu-
lation, we attribute the high radiative efficiency to magnetic compression of the disk
into a thin MAD, magnetized wind advection of radiation away from the disk, and
the formation of a low-density jet channel through which radiation can more freely
stream. Our simulation boundary is at r ∼ 105rg, with accurate energy conversion
out to r ∼ 600rg. Our photosphere at r ∼ 400rg accounts for wind material that
96
reached r & 3000rg. If we ignore the extended wind, for the a/M = 0.8 MAD model
we would miscalculate the radiative efficiency to be η ∼ 30% – an overestimate by
a factor of two. For the a/M = 0 MAD model, measuring at r ∼ 60rg gives η ∼ 5%
– an overestimate from η ∼ 0% (caused the wind absorbing radiation energy).
4.4 Summary
We have performed fully 3D simulations of super-Eddington accretion, includ-
ing a simulation with M˙0 ∼ 400LEdd/c2 ∼ 50M˙Edd onto a rotating black hole with
a/M = 0.8. Sufficient magnetic flux was distributed throughout the disk that the
BH and disk reached MAD levels, where magnetic forces pushing out balance gas
forces pushing in. The MAD state enables the super-Eddington accretion flow to
reach its maximum efficiency, with the total efficiency measured to be about 50%
(higher for higher BH spins). Importantly, the system has a high radiative effi-
ciency of about 15% (luminosity L ∼ 50LEdd) beyond the resolved photosphere at
r ∼ 400rg. This occurs because the magnetized wind carries radiation away from the
disk. Also, the magnetized jet creates a low-density channel for radiation to more
freely stream and pushes away the more opaque wind. These effects increase the
radiative flux escaping the disk and diminish conversion of radiation energy flux into
kinetic energy flux of the wind. This mechanism allows high radiative efficiencies
from super-Eddington accretion systems.
As applied to jetted TDEs (Tchekhovskoy et al., 2014), our super-Eddington
MAD simulations show how the jet and radiative efficiencies depend not only upon
BH spin and magnetic flux but also upon jet mass-loading physics that leads to sub-
optimal BZ efficiencies via dragging of BH field lines by optically thick radiation.
Also, while the wind absorbs some jet energy, the jet channel exposes the (otherwise
obscured) inner hot X-ray emitting disk.
97
Acknowledgments
We thank Ramesh Narayan, Alexander Tchekhovskoy, Yan-Fei Jiang, and
Aleksander Sadowski for discussions and acknowledge NASA/NSF/TCAN (NNX14AB46G),
NSF/XSEDE/TACC (TG-PHY120005), and NASA/Pleiades (SMD-14-5451).
98
Chapter 5: Role of Magnetic Field Strength and Numerical Reso-
lution in Simulations of the Heat-flux Driven Buoyancy
Instability
Abstract
The role played by magnetic fields in the intracluster medium (ICM) of
galaxy clusters is complex. The weakly collisional nature of the ICM leads to
thermal conduction that is channelled along field lines. This anisotropic heat
conduction profoundly changes the stability of the ICM atmosphere, with
convective stabilities being driven by temperature gradients of either sign.
Here, we employ the Athena magnetohydrodynamic code to investigate the
local non-linear behavior of the heat-flux driven buoyancy instability (HBI),
relevant in the cores of cooling-core clusters where the temperature increases
with radius. We study a grid of 2-d simulations that span a large range
of initial magnetic field strengths and numerical resolutions. For very weak
initial fields, we recover the previously known result that the HBI wraps the
field in the horizontal direction thereby shutting off the heat flux. However, we
find that simulations which begin with intermediate initial field strengths have
a qualitatively different behavior, forming HBI-stable filaments that resist
field-line wrapping and enable sustained vertical conductive heat flux at a level
of 10–25% of the Spitzer value. While astrophysical conclusions regarding the
role of conduction in cooling cores require detailed global models, our local
study proves that systems dominated by HBI do not necessarily quench the
conductive heat flux.
99
5.1 Introduction
The intracluster medium (ICM) is a hot, weakly-collisional plasma compris-
ing the majority of the baryonic mass in galaxy clusters. It has gained increasing
attention in the last few decades since it provides a rich environment for the study
of feedback, cosmology, low-density magnetohydrodynamics (MHD), and magneto-
genesis. One area of recent focus is the study of dynamics induced by anisotropic
conduction and anisotropic viscosity. The main goal of these studies has been to
develop a better understanding of the role of magnetic fields in moderating heat
conduction in the ICM.
Anisotropic conduction results from the fact that the electron gyroradius is
much smaller than its mean free path and renders the ICM atmosphere buoyantly
unstable to temperature gradients of either sign, relative to the direction of gravity.
The outer region of the ICM, where virialization supports a negative temperature
gradient, falls subject to the magnetothermal instability (MTI; Balbus, 2000). On
the other hand, in the inner region of cool-core clusters the plasma density is high
enough for bremsstrahlung cooling to turn over the temperature gradient and the
atmosphere is unstable to the heat-flux driven buoyancy instability (HBI; Quataert,
2008, hereafter Q08).
Previous studies of the local non-linear evolution of the HBI using numerical
simulations (Parrish and Quataert, 2008; McCourt et al., 2011) found that the in-
stability acts to quench vertical heat conduction by wrapping magnetic field lines
horizontally thereby insulating each temperature layer from one another. If this
strong insulation effect occurs in real clusters, it renders thermal conduction irrele-
vant in consideration of the thermal energy balance of the cooling core.
In the absence of significant conduction, one must invoke AGN feedback or
another type of turbulent stirring to stave off catastrophic cooling, but any such
mechanism would simultaneously reorient some field lines vertically and conduction
would again drive turbulence via the HBI. Therefore any self-consistent model of the
ICM needs to both realistically capture the physics of the HBI as well as feedback
100
processes (Balbus and Reynolds, 2008; Sharma et al., 2009; Parrish et al., 2010;
Kunz et al., 2011).
The HBI is fundamentally a weak-field instability, operating even when mag-
netic forces are completely negligible (provided that the field is strong enough to
keep the electron gyroradius smaller than its mean free path). For this reason, most
previous studies of the HBI have assumed extremely weak initial fields, β & 1011
where β ≡ P/Pmag with P being the thermal plasma pressure and Pmag the magnetic
pressure. However, the role of magnetic field/tension on the non-linear behavior of
the HBI is poorly understood. Kunz et al. (2012, hereafter K2012) performed semi-
global HBI runs with initial β = 105 and found a qualitatively different behavior
to that seen previously. Instead of the field line wrapping and quenching of con-
duction seen by previous authors, K2012 found sustained conduction via magnetic
filaments. However, it was unclear whether the semi-global nature of the simulation
or the stronger initial field was the main cause of this difference.
This leaves open several questions. Does magnetic tension play a defining role
in shaping the non-linear behavior of the HBI? Since tension forces are most relevant
for highly curved field lines, is the numerical resolution typically used in simulations
sufficient to fully capture the physics of the HBI? What role do magnetic filaments
play in the transport of heat in the ICM?
In this paper we present new 2-d simulations in the local regime which cover the
resolution-β parameter space in an attempt to answer these questions. We explore
how the non-linear dynamics of the weakly collisional ICM changes as the initial
magnetic field strength is varied from complete convective stability (HBI completely
suppressed by magnetic tension at the local scale) to the low field strength used in
previous studies. In order to determine whether current simulations have sufficient
grid resolution for all aspects of the physics of the HBI to be fully converged, we
also study trends in fundamental measures of the non-linear state of the plasma as
grid resolution is varied.
We principally find that, as suggested by inviscid simulations of K2012, the
non-linear state of the HBI is qualitatively different than in previous studies when an
101
intermediate (within a few orders of magnitude of HBI stability) initial magnetic field
strength is chosen. In particular, there is enough magnetic flux to form sustained
vertical filaments which prevent the quenching of vertical heat conduction in the
ICM, saturating at a significant fraction of the Spitzer conductivity (Spitzer, 1962).
We present the results of our parameter space exploration and describe new insights
into the physics of HBI growth as well as the formation of magnetic filaments,
including a connection between the non-linear turbulent state of the plasma and the
linear theory. Our results hold for moderate resolution 3-d simulations, validating
the applicability of a 2-d study and the physical model we present.
The paper is organized as follows. In §5.2 we briefly review the analytic de-
scription of the HBI focusing on aspects relevant for our results. In §5.3 we discuss
the numerical method, boundary conditions, and initial conditions. §5.4 presents
the results of the simulations and describes major trends we observe, which are then
placed into context of the physics of the HBI in §5.5. In §5.6 a complimentary 3-d
simulation of moderate resolution is described, and §5.7 provides a discussion of our
work in the context of existing literature.
5.2 MHD of Weakly Magnetized Plasmas
We use the standard MHD equations with a term added to the entropy equa-
tion for anisotropic thermal conduction along magnetic field lines. For easy com-
parison to recent literature we primarily follow the prescription for linear analysis
laid out in Balbus and Reynolds (2010, hereafter BR10) while excluding the term
for radiative cooling.
We choose rationalized natural units with kB = c = 0 = µ0 = 1 such that
total energy density is given by
E = + ρ
v · v
2
+
B ·B
8pi
− ρg · x
(γ − 1) (5.1)
with the internal (i.e., thermal) energy density  = P/(γ − 1), where we assume
the adiabatic index γ = 5/3 and the mean molecular weight µ = 1 (pure hydrogen,
fully ionized, plasma). The third term on the right hand side is the magnetic energy
102
density, or pressure, used in the definition of the plasma β-parameter
β =
8piP
B2
. (5.2)
The vector quantity x = xxˆ+zzˆ defines position in the 2-d simulations. zˆ points ver-
tically upward and xˆ points horizontally to the right. Then, the mass conservation,
momentum conservation, induction, and entropy equations are, respectively,
∂ρ
∂t
+∇ · (ρv) = 0, (5.3)
ρ
Dv
Dt
=
(∇×B)×B
4pi
−∇P + ρg, (5.4)
∂B
∂t
= ∇× (v ×B), (5.5)
D lnPρ−γ
Dt
= −γ − 1
P
∇ ·Q (5.6)
where ρ is the mass density, v is the fluid velocity, B is the magnetic field vector,
P is the gas pressure, g is the gravitational acceleration, and Q is the conductive
heat flux. D/Dt ≡ ∂/∂t+ v · ∇ is the Lagrangian derivative.
We define the anisotropic conductive heat flux Q, with b = B/B the unit
vector in the direction of the magnetic field, as
Q = −χb(b · ∇)T, (5.7)
where T is the kinetic gas temperature and
χ ' 6× 10−7T 5/2 erg cm−1 s−1 K−1
is the thermal conductivity (Spitzer, 1962) in physical units. As in Q08 and BR10
we define κ ≡ χT/P to be the anisotropic thermal diffusion coefficient. In natural
units, κ ≈ 10−2 is the Spitzer diffusion coefficient for a plasma of temperature 1
keV.
We ignore the ion component of the heat flux since it is smaller than the
electron contribution by a factor of (mi/me)
1/2 ≈ 42. Kunz et al. (2011) showed
that, despite this large ratio, viscosity can have a significant effect on the HBI for
conductivity at the Spitzer value, as chosen here, but that this effect is severely
103
diminished deep in cluster cores where the plasma is most collisional. For this work,
we choose not to include viscosity so that we may probe the connection between
the formation of magnetic filaments and the conduction driven physics of the HBI
alone. Regardless, the presence of anisotropic viscosity seems to only enhance the
formation of filaments (Parrish et al., 2012).
5.2.1 Background Equilibrium Conditions
We consider a vertically stratified atmosphere in which the gas is not self-
gravitating and the extent of the domain is small enough to be in the local regime.
We then approximate the gravitational acceleration to be a constant function of
position g(z) = −g0zˆ.
The initial (unperturbed) magnetic field is purely vertical for all simulations
presented. The simulations start with equilibrium initial conditions in which mag-
netic pressure is negligible compared to the thermal pressure. Moreover, a uniform
initial field provides zero magnetic pressure gradient so field strength does not affect
the equilibrium condition at all. Therefore the equilibrium state of the atmosphere
is that of hydrostatic balance with
dP
dz
= g(z)ρ = −g0ρ. (5.8)
Given an initial temperature profile T(z), in order to show distinction between
that of the HBI and that of the classical adiabatic (Schwarzschild) instability, we
employ profiles that have entropy increasing with height, i.e., that are Schwarzschild
stable. Thus, as will be shown in §5.3, we assume a simple power law representation
for equilibrium temperature and density such that ∂s/∂z > 0 where the entropy s
is
s =
P
ργ
.
5.2.2 Stability analysis
We refer the reader to BR10 for a full construction of the Wentzel-Kramers-
Brillouin (WKB) perturbation analysis using plane wave disturbances of the form
104
exp(σt+ ik · x) where x is the position vector, σ is the formal growth rate, and the
wavenumber k has Cartesian components. However, there are a few points in the
analysis which are particularly relevant to the results in this paper.
Without cooling, the linear analysis in BR10 results in the following dispersion
relation:
(σ + C) (σ2+(k · vA)2) (5.9)
+
σk2⊥N
2
k2
+ CK
g
k2
∂ lnT
∂z
= 0, (5.10)
where
C ≡
(
γ − 1
γ
)
κ(k · b)2, (5.11)
K ≡ (1− 2b2z)k2⊥ + 2bxbzkxkz, (5.12)
and the Brunt-Vaisala frequency, N, given by
N2 = − 1
γρ
∂P
∂z
∂ ln s
∂z
, (5.13)
describes the buoyant response of an adiabatic plasma. k⊥ is the component of k
perpendicular to the magnetic field direction. The square of the Alfvenic frequency
is ω2A ≡ (k · vA)2.
This cubic dispersion relation encodes three conditions for stability, one of
which encodes the HBI and MTI criterion of Balbus (2000) and Q08,
K
g
k2
∂ lnT
∂z
+ ω2A > 0, (5.14)
where it is convenient to identify the characteristic dynamical timescale with the
HBI/MTI growth rate (e-folding rate), ω2dyn ≡ g ∂ lnT∂z . Eqn. (5.14) then simplifies,
in a purely vertical magnetic field, to
− k
2
⊥
k2
g
∂ lnT
∂z
+
B2
4piρ
(k · b)2 > 0 (5.15)
This equation is essentially a quantification of the competition between the
effects of a destabilizing temperature gradient and the stabilizing influence of mag-
netic tension.
105
5.3 Numerical Setup
For our simulations we use the ATHENA MHD code (Stone et al., 2008)
which uses an unsplit Godunov method utilizing constrained transport (Gardiner
and Stone, 2005) to preserve the divergence free nature of a physical magnetic field.
To prevent unphysical transport of heat against ∇T we add a monotonic flux limiter
to our problem algorithm according to the scheme laid out in Sharma and Hammett
(2007).
Atmospheres unstable to the HBI are those with a positive temperature gra-
dient. To achieve this condition and maintain positive entropy gradient we set up
the initial equilibrium configuration with temperature, density, and pressure given
by power laws of the form
T (z) = T0
(
1 +
z
HT
)
(5.16a)
ρ(z) = ρ0
(
1 +
z
HT
)−2
(5.16b)
P (z) = P0
(
1 +
z
HT
)−1
, (5.16c)
with HT ≡ (d lnT/dz)−1 = 2 the characteristic scale height of the atmosphere with
sound speed cs ≡
√
γP/ρ ≈ 1.3. From hydrostatic equilibrium, with g = 1, we have
T0 = 2ρ0 = P0 = 2µ = 2 and we choose an initial magnetic field oriented purely
vertically (B = B0zˆ).
This equilibrium configuration has a conductive heat flux that is divergence
free,
d
dz
(
χ‖
dT
dz
)
= 0, (5.17)
since we choose a constant value for the parallel component of thermal conduction,
χ‖. Aside from ensuring that the HBI will develop with the same growth rate ev-
erywhere in the anisotropically conducting part of the atmosphere, since our models
do not account for radiative cooling and zero-divergence implies no heating, the
atmosphere is initially in thermal equilibrium.
Our simulations have dimensionless width×height of 0.1×0.3 and we limit
106
anisotropic conduction and all our measurements of the gas dynamics to the central
0.1×0.1 region. The remaining 0.1×0.1 buffer zones above and below have isotropic
conduction of the same magnitude as conduction in the active central zone. Early
testing revealed the vertical reflective momentum boundary condition to cause sup-
pression of the vertical growth by deflecting motion of the plasma before it naturally
reached a non-linear state, but addition of buffer zones of this size alleviates the ef-
fects of the boundaries.
We choose momentum reflective boundary conditions for the upper and lower
boundaries to provide mass conservation and pressure support. Periodic boundary
conditions are chosen for the left and right boundaries. The temperature is held
constant at the upper and lower boundaries where the magnetic field vector is set
identically to that in the last active cell for each column. This magnetic boundary
condition conserves total vertical magnetic flux.
In order to resolve the linear development of the HBI cleanly, we applied single-
mode perturbations to the equilibrium atmosphere in most of our simulations. For
consistency, we tested the dependence of our results on the perturbation spectrum by
running a set (the H2d128 xPS series) of simulations with identical initial conditions
but a power spectrum perturbation in momentum with |vk|2 ∝ k−α of Kolmogorov
type, α = 5/3. This produces a non-linear magnetic field topology approximately
the same as a pure-mode perturbation where two wavelengths fit in each direction
in the central square box. For all simulations exploring the β-resolution parame-
ter space we start with a single pure-mode in both velocity and magnetic field of
magnitude kx = kz = 4pi/0.1.
We performed 24 simulations to cover the β-resolution parameter space and
use the designation H2d128 5, for example, to label our 2-d HBI simulation with a
resolution (covering the active region) of 128×128 and an initial plasma β-parameter
of β = 2 × 105. To connect our work with both linear theory and previous studies
we chose five values of log10(β) ∼ 3, 5, 7, 9 and 11, and five grid resolutions R× 3R
where R = [32, 64, 128, 256, 512] 1. In addition to these 24 simulations we ran a
1We chose not to use limited computing resources on the 25th simulation, H2d512 3, since the
107
power spectrum simulation of resolution 128×384 for each β and a 3-d simulation
of resolution 128×128×384 of β = 2× 107.
5.4 Plasma Behavior as a Function of Resolution and β
In this section we show the results of our β-resolution parameter space study.
We first describe the effect of varying β (i.e., only changing magnetic field strength)
on the non-linear saturation of the HBI. We then show the dependence of these
effects, and evidence for general convergence of the physics of the HBI, with the
grid resolution.
As a precursor to our discussion of the full parameter space, we show in the first
two frames in Figures 5.1a and 5.1b the temperature and magnetic field structure of
our atmosphere under the pure-mode perturbation. These frames show the entire
height of the atmosphere at times when the HBI is in the linear growth phase and
just reaching the non-linear regime. Note that for both simulations of Fig. 5.1a and
5.1b, corresponding to values of β = 2 × 107 and 2 × 109 respectively, the linear
growth phase appears identical, in agreement with the expectations of the linear
theory.
Now we focus attention on the development of the HBI into the nonlinear
regime illustrated in the late-time frames of Figures 5.1a and 5.1b.
5.4.1 Plasma β
Central to any study of the ICM is the plasma β-parameter. In the context
of these simulations, consider the fully non-linear state of the plasma shown in the
later frames in Figure 5.1. The simulation with a moderate initial value of β ≈ 107
(H2d128 7) contrasts qualitatively to that with β ≈ 109 (H2d128 9).
The major qualitative difference is that magnetic flux bundles, or filaments,
are present near the end of the linear phase in both simulations, but persist over
many dynamical times only in the moderate β case (2 × 105 and 2 × 107). These
HBI is stable for that value of β and no dynamics is seen for any grid resolution.
108
Figure 5.1: (a) H2d128 7 at four times (in units of dynamical times, ω−1dyn). The
color bar on the right indicates temperature and magnetic field lines are overlaid.
(b) H2d128 9 at identical times as (a). Note the clear distinction between the two
after many dynamical times, shown in the last frames.
109
Figure 5.2: Horizontally averaged vertical heat flux as a function of time for 24
simulations described in the text. Each panel in the plot has identically scaled axes.
110
filaments have been seen before as the seeds of the horizontal bands of magnetic field
lines which develop in high β simulations. Horizontal motion, free from any restor-
ing force (see McCourt et al., 2011, for discussion), wraps the magnetic field lines
thereby insulating temperature layers from one another and continuously growing
in magnetic field strength. However, for the moderate values of β these filaments
persist against horizontal bulk motion and reach an entirely different saturated state.
The most important physical quantity characterizing the system, and an excel-
lent discriminator of the different non-linear states, is the vertical heat flux (VHF),
Fz = −
∫
S
χbˆ(bˆ · ∇)T · dS,
where S is a surface that cuts horizontally across the domain. Figure 5.2 shows the
horizontal average of the VHF across the center of the active domain (z = 1.5) as a
function of time for all simulations in the pure-mode study. The column containing
those simulations stable to the HBI shows constant heat flux consistent with the
magnetic field remaining vertical and conducting uniformly. The two intermediate
β columns show a nonzero heat flux in the saturated state. As shown in Fig. 5.1,
the HBI has shut off VHF in most of the plasma, but the filaments remain vertically
oriented against the turbulent motion of the plasma and support a sustained VHF.
The higher initial magnetic field strength results in a larger magnitude of sustained
VHF, which we explain via the physical model for the filaments described in §5.5.
Finally, the high β simulations lead to a non-linear state with zero VHF, duplicating
the result of previous studies. Thus, one can identify the behavior in Fig. 5.1a with
the 2nd and 3rd columns of Figure 5.2, and Fig. 5.1b with the last two.
Another interesting quantity is kinetic energy. In order to study energy equipar-
tition and turbulence in an inherently asymmetric atmosphere we break the ki-
netic energy scalar into the components contributed by motion in each direction,
KEz ≡ 12mv2z and KEx ≡ 12mv2x. Figure 5.3 shows how 〈KEz〉, the vertical kinetic
energy contribution spatially averaged across the entire active central region, varies
with β. H2d256 9 and H2d256 11 show approximately linear damping of vertical
motion, but H2d256 5 and H2d256 7 show an increase in vertical momentum with
111
Figure 5.3: The contribution to kinetic energy from vertical motion spatially aver-
aged across the entire anisotropically conducting region is shown as a function of
time for simulations H2d256 x, x∼[5 (blue-long-dash), 7 (green-dash-dot-dot-dot), 9
(black-dash-dot), 11 (red-dash)].
higher initial magnetic energy density. Both H2d256 5 and H2d256 7 saturate at a
value orders of magnitude above the weak field cases, where the vertical momentum
does not stop declining up to the end of our simulations.
However, the kinetic energy is not the only energetic discriminant for the
non-linear behavior of the plasma. Fig. 5.4 shows the vertical and horizontal contri-
butions to the volume integrated magnetic and kinetic energy within the anisotrop-
ically conducting region. For β = 2 × 109 the kinetic energy completely dominates
the dynamics of the plasma early on as horizontal motion wraps the field into layers.
It is only at late times when the field wrapping has transferred enough energy into
the magnetic field for the horizontal component B2x/8pi to overcome KEx.
The moderate β columns in Fig. 5.4 show something very different. For
sufficient resolution there is approximate equipartition between B2x/8pi ≈ KEx ≈
KEz. In fact, for initial β = 2×107, even B2z/8pi reaches equipartition with the other
energy components. In general, however, B2z/8pi is subject to a special constraint,
namely total vertical magnetic flux through the simulation domain is conserved, and
so does not necessarily come into equipartition.
While the simulations start with four different initial field strengths spanning
112
Figure 5.4: Evolution of energy components, in our rationalized units, volume inte-
grated over the active (central) region of the computational domain of dimensions
0.1x0.1x0.001. The last dimension is the arbitrary skin depth we choose for the com-
pact third dimension. Only results of HBI unstable simulations are shown, namely,
those for β = 2× 103 are omitted. KEx kinetic energy (black-line), KEz kinetic en-
ergy (black-dotted), B2x/8pi magnetic energy (red-dash), and B
2
z/8pi magnetic energy
(red-dash-dot). Each panel in the plot has identically scaled axes.
113
four orders of magnitude, the total of all the energy components (excluding B2z/8pi)
is the same for all simulations with resolution at least 64×192, whether filaments are
sustained or not. This is expected from the linear theory since the linear growth rate
of pure-mode perturbations is independent of magnetic field strength (in the weak-
field limit) and the linear growth entirely determines the amount of kinetic energy
at the time of transition into the non-linear regime. From inspection of Fig. 5.4 it
seems that the balance between volume integrated kinetic and magnetic energy soon
after transition into the non-linear regime is an indicator of whether or not filaments
are sustained. This makes sense as an excess in KE is necessary to overcome the
tension within the filaments. Production of sustained filaments is intimately tied to
the total kinetic energy, and thus the turbulent state of the plasma.
5.4.2 Grid Resolution
Now we consider the more subtle effect of grid resolution as well as its in-
terplay with β. The amplitude of VHF in the two intermediate β columns of Fig
5.2 decreases with increasing resolution. This is illustrated in Fig. 5.5 where the
conduction rate vs. time for the β = 2 × 107 column at all five grid resolutions is
smoothed and plotted in superposition. Inset is the time-averaged value of conduc-
tion during the last 18 dynamical periods for each resolution for β = 2 × 105 and
2× 107. One can see that for β ∼ 105 there is negligible change in heat flux for the
last doubling of resolution. However, for β ∼ 107 there is no such sign of conver-
gence. While the data are too sparsely sampled in resolution to conclusively make
this distinction, it is consistent with the general physical model of the filaments we
propose in the next section.
Dependence on grid resolution is also evident in Fig. 5.4. Firstly, for moderate
β, increasing resolution produces less dispersion in the values of the four volume
integrated energy components. Secondly, for high β simulations there is an increase
in the vertical energy components with increasing resolution. This second effect
could be due to the resolution of small scale turbulence.
For high β simulations it also takes less time with increasing resolution for
114
Figure 5.5: Arbitrarily scaled, horizontally averaged vertical heat flux plotted in
the β ∼ 107 column of Fig. 5.2 versus time in dynamical units. Colors, in order
of transition from black to blue, signify increasing resolution. Inset: Each data
point represents the flux averaged over the last 18ω−1dyn for each grid resolution.
Convergence is shown for simulations with initial β ∼ 105 (diamonds) but not
β ∼ 107 (stars).
KEx and B
2
x/8pi to switch dominance. This is consistent with a lower level of
numerically driven magnetic reconnection in the higher resolution simulations and
the regions of highest magnetic field curvature. Finally, in the β = 2× 105 column,
increasing resolution results in a lower saturated value of KEz. In these simulations
a significant fraction of the vertical motion is bulk flow along the filament, similar
to what was seen in the filaments of K2012. We do not elaborate on the flow in
this study, but it suffices here to simply say the width of filaments is observed to
decrease with increasing resolution. Thus, the volume integrated KEz as well as the
resolution dependence of VHF can be explained by a dependence of filament width
on resolution. In the next section we present a physical model for the filaments and
explain their dependence on the β-resolution parameter space.
5.5 Dynamics and Persistence of Filaments
Why do filaments form and what sets their field strength and width? We
hypothesize that filaments are the result of needing to pass a (conserved) vertical
115
magnetic flux through a plasma subject to an instability which attempts to remove
vertical field. We suggest that the filaments collapse until their internal field strength
renders them HBI stable, with further collapse likely prevented by the diffusive
process of turbulence.
Without viscosity, a plasma governed by ideal MHD will be driven by tur-
bulent motions toward equipartition of the total kinetic and magnetic energies via
field amplification. Therefore, after the filaments are produced, if the kinetic energy
density is sufficiently higher than the magnetic energy density, the large scale bulk
motions of the plasma (preferentially horizontal since the atmosphere is classically
Schwarzschild stable) will wrap up the magnetic field into horizontal layers. If, on
the other hand, the magnetic field dominates overall, most magnetic flux is vertical
and the filaments are sustained. While the detailed dynamics of filament stability
against the turbulent motions of the plasma could be framed in terms of horizon-
tal velocity shear and height of the filaments, the complex magnetic topology and
boundary limitations of the simulations suggest more insight may be gleaned from
an understanding of filament persistence in terms of globally averaged quantities.
Therefore, our model presents a global energetics understanding of filament per-
sistence, and a separate model for filament width and field strength in terms of
the HBI. This separation is natural since perturbations of wavelength greater than
about twice the width of the filament are better understood geometrically in terms
of bending the flux tube than in terms of the linear HBI.
A direct prediction of our model is that simulations should show that filaments
lie on the stable side of the threshold between stable and unstable regions of phase
space. Fig. 5.6 shows the result of such a measurement, obtained via the method
we now describe.
First consider the size of the filaments. Our hypothesized model predicts
filaments will have a width corresponding to a wavenumber which just satisfies HBI
stability, i.e., the filaments are stabilized when ω2A > ω
2
dyn for a region inside the
filaments of size equivalent to the scale of a random perturbation.
In our analysis for Fig. 5.6, we must be able to define and chracterize filaments
116
Figure 5.6: H2d512 5 (upper) and H2d512 7 (lower) filaments. Stability of the
filaments is demonstrated by plotting the three stability diagnostics of Eqn. (5.19)
measured for each filament. Phase space below the solid line is HBI unstable, above
is stable. The data points are evenly temporally sampled from the last 35.4ω−1dyn.
Those plotted are a subset of all extracted, for visual clarity, but the contours
represent all the data and enclose 2[4,5,...17] points.
117
Figure 5.7: Black contours (dash-dot) represent filaments from H2d128 7PS; red
(solid) represent those from H2d128 7. Contour levels are set the same as in Fig.
5.6.
in an automated manner. We consider a filament to be located wherever a series
of horizontally adjacent grid cells satisfy both |bˆz| > 0.7 and βcell < 0.1〈
√
β〉2, a
weighted average of β over the active domain. The values of magnetic energy (for
which we can use β−1 as a proxy since the average thermal pressure varies at most
by 2%), the vertical component of the temperature gradient, dT
dz
|z, and the filament
width, which we identify with wavenumber kw, are extracted for each such collection
of cells for all but the rows of our anisotropically conducting active region directly
adjacent to the buffer regions. Each set of values produces a point in Fig. 5.6.
In the plot the kw value used is calculated assuming a wavelength set at twice the
measured width of the filament. The values of vertical temperature gradient within
the filament and β−1 are unweighted horizontal averages across the filament-defining
cells.
If we assume bx is close to zero, which is generally observed to be the case
for sustained filaments, then it follows that K = −k2⊥. If we additionally assume
k⊥ ≈ k‖ for an average perturbation, then 2k2‖ ≈ k2, and Eqn. (5.15) simplifies to
− 1
2
g
T 2
dT
dz
+
1
β
k2 > 0. (5.18)
This provides a general relation giving the boundary of HBI stability in terms of our
three diagnostic measures of the filaments. In order to facilitate comparison of the
118
relation with the data in Fig. 5.6 we further simplify by setting T = 2, in our local
simulations where temperature does not significantly spatially vary. Also, since the
width of the filament is the scale at which it is stabilized to random perturbations
under our model, we identify k as the width kw. We then obtain
β−1 ∼ 1
8
dT
dz
k−2w . (5.19)
In Fig. 5.6 data points created in the method just described are extracted
from the active region for the last ∼ 35 ω−1dyn of the simulation at 200 evenly spaced
times. The magnetic filaments are clearly subject to the boundary provided by the
stability relation, thus supporting the theory that the filaments are islands of HBI
stability within an otherwise HBI susceptible plasma. Specifically, consider how the
filaments in H2d512 7 and H2d512 5 are scattered along the stable side of the curve
but generally follow it. While there is a significant vertical spread in β−1 above
the line, consider that initial values in these simulations were β−1 = 5 × 10−8 and
5×10−6 and the final spread in β−1 above the intra-filament value is consistent with
the theoretical boundary over nearly three orders of magnitude. Note that some
filament values of β−1 fall below the initial conditions value. This is consistent with
our definition of a filament having significant internal magnetic energy density over
the spatial average, and most space in the atmosphere has a very weak field in the
saturated state.
Figure 5.7 compares filaments from H2d128 7 and H2d128 7PS showing that
they cover roughly the same region in this space. This suggests the pure-mode per-
turbation recovers most of the physics of a more multi-mode initial condition, and
both lie at weaker field strengths than the simulations with higher resolution, sug-
gesting the filaments in the lower resolution simulation are prevented from shrinking
to obtain their natural peak field strength. This picture is also consistent with the
H2d128 7 data points residing at the higher values of k−2w on the plot.
This self-consistent model for the internal stability of the filaments is sufficient
to explain the width of the filaments and their dependence on magnetic field strength
and resolution, as we’ve shown. However, it also predicts the viability of magnetic
119
filaments in a plasma of arbitrarily weak magnetic field so long as there is sufficient
resolution for the filament to shrink until it becomes internally HBI stable. This
would predict that the boundary between our qualitatively different non-linear states
would shift to higher β as resolution is increased in Fig. 5.5. While it is likely we
inadequately resolve the β-resolution phase space to see such a trend to begin with,
the true dissolution of this discrepancy is evident upon closer inspection of the global
energetics in our simulations.
Recall that both the two intermediate and two high β columns of Fig. 5.2 are
well within the weak-field regime (See §5.2). The growth rate then only depends
on the pure-mode perturbation wavelength, which is identical for all simulations in
our phase-space suite. Transition to non-linearity happens at the same time in all
cases, and the exponential growth in momentum during the linear phase is driven
by a conversion of heat into kinetic energy. These facts together imply that the
total kinetic energy imprinted on the plasma at the end of the linear regime is the
integrated net heat-flux passing into the active region during that time. For these
reasons, and evident in Fig. 5.4, the peak in kinetic energy at the end of linear
growth is identical in all the simulations. Those with initial β ∼ 105 do have a
slightly higher transition value, but that is consistent with the additional magnetic
tension maintaining a pure-mode magnetic topology for slightly longer, therefore
increasing the time over which the energy conversion takes place. Clearly the peaks
in KE components for that column occur at later times than in the other columns.
In the second part of our model we suggest that the global persistence of
filaments once they’ve formed depends on how the imprinted kinetic energy just
described compares to the initial total magnetic energy. We argue that if approxi-
mately twice the magnetic energy density is greater than the kinetic energy density
imprinted at the transition to non-linearity, then sustained vertical filaments of the
form we see will be present. This is effectively a statement that either the mag-
netic field or the turbulent motion dominates and one or the other wins out early
in our simulations, leading to a distinct division in parameter space. The ‘twice’
comes from allowing an x-component of the filaments to grow to equal magnitude
120
with the vertical z-component. Beyond that, the filament would be more horizontal
than vertical, significantly quench vertical heat flux, and effectively be stretched out
beyond our definition of a filament.
This scale-invariant method for understanding filament stability is particularly
suited to a local study such as ours. Of course, it suggests that a sustained source
of turbulence may significantly alter the non-linear state of the plasma. We leave a
discussion of this and other implications for the last section.
5.6 3-D Simulation
Some studies suggest (Parrish et al., 2012) that the filaments formed in 2-
d simulations are a result of field lines not being able to slip past one another,
that is, a prevention of flux-interchange modes in the plasma. K2012 comments
on this in their H3dBrag simulation, where the filaments become more reoriented
(horizontally) in the lower part of the atmosphere where the effects of Braginskii
viscosity are less pronounced. However, this region of H3dBrag is very close to the
simulation boundary, and the effect of the boundary on the detailed field structure
is unclear.
Our model explaining the presence of the filaments, presented in the previous
section, suggests this existence is rooted in the physics of the HBI and not a feature of
dimensionality. K2012 does not provide a complimentary 3-d semi-global simulation
of the HBI without viscosity and so in extension of that work, and to confirm the
robustness of our results, we performed a 3-d simulation of 128x128x384 resolution
and moderate field strength, H3d128 7PS. The simulation was in every way identical
to H2d128 7PS save for the added dimension.
Figure 5.8 shows the magnetic structure of the non-linear state of H3d128 7PS
at late times, with coherent filaments of a qualitatively similar sort to those found
in the 2-d analog. Just as is pointed out in K2012, these filaments contain more
magnetic flux because it is easier to collect with the additional dimension during
linear growth.
121
Figure 5.8: Shown is a frame from the saturated state of H3d128 7PS. Color value
indicates magnetic field strength and field lines are integrated from grids of points
at the top and bottom of the domain. Many of the field lines thread a filament near
the center of the anisotropically conducting region.
The VHF measured in the same way as before is shown in Fig. 5.9. H3d128 7PS
clearly runs long enough for the VHF to stabilize at a value of ≈ 25% the initial
value, just as before. This value is no coincidence. The addition of a third dimen-
sion does not change the internally stabilizing β of the filaments, and so the ratio
of conducting surface area to non-conducting is still set by the initial plasma β and
total vertical flux conservation, assuming that most vertical flux is bundled into the
filaments.
5.7 Discussion and Conclusions
In this paper, we present a systematic investigation of the local non-linear
behavior of the HBI as a function of initial magnetic field strength and numerical
resolution. We identify a previously unrecognized regime of behavior at interme-
diate field strengths in which magnetic filaments are formed that can resist the
horizontal field line wrapping seen in weak-field simulations (Parrish and Quataert,
122
Figure 5.9: Vertical Heat Flux integrated over the horizontal plane slicing through
the center of the active domain in H3d128 7PS.
2008; McCourt et al., 2011) and lead to sustained conductive heat flux. We show
that the filaments themselves are HBI-stable and provide the main conduit by which
(conserved) vertical magnetic flux is passed through the domain. In addition to this
model for internal stabilization, we describe how global energetics determine the
ultimate qualitative state of the non-linear regime, that is, whether or not filaments
persist and there is a resulting significant vertical heat flux.
The robust formation of filamentary structures is tantalizing since we see this
type of structure in the cooling cores of galaxy clusters, namely in the form of
Hα filaments. For instance, HST observations of the Perseus cluster, which has
been studied closely for some time (Lynds, 1970), reveal complex braided filaments
which can be as narrow as 50 pc but may be kiloparsecs in length (Fabian et al.,
2008). More recent observations by McDonald et al. (2010, 2011) have shown that
essentially all cool-core clusters exhibit extended Hα systems which appear to extend
out to almost, but never beyond, the cooling radius of the X-ray emitting ICM. This
connection between the cooling radius and the filamentary structure suggests they
may condense out of the ICM as a result of thermal instability (Field, 1965).
However, it is the local thermal instability which produces cold filaments in
simulations, with anisotropic conduction primarily responsible for their filamentary
structure. Without a heating source to offset the background global thermal instabil-
123
ity, namely, when net cooling is purely a function of local density and temperature
of the gas, buoyancy effects stabilize the plasma to the local thermal instability
(Balbus and Soker, 1989). K2012 shows that for an atmosphere with intermediate
magnetic field strength and anisotropic viscosity, viscous heating during formation
of the filaments and a sustained vertical heat flux can significantly delay the cooling
catastrophe resulting from global thermal instability. This allows local thermal in-
stability to develop and cold filaments to naturally form. Our study supports K2012,
firmly establishing the robust formation of filaments over a range in magnetic field
strength relevant to cool-core clusters, and exploring in detail the magnitude of ver-
tical heat flux sustained by the atmosphere and its dependence on local magnetic
field strength and numerical resolution.
Though our work has not yet shown a definitive connection to the observed
Hα filaments, in part because we do not include radiative cooling in our simulations,
our findings strongly motivate further work in developing physically realistic global
simulations of cluster cores, where the ICM is likely HBI unstable. Since filaments
provide a mechanism for transporting heat (and possibly cold gas) radially in these
atmospheres, an understanding of the magnetic topology from the local to the global
scale and how these scales connect will be instrumental in understanding the role of
conduction in regulating global thermal stability. Also instrumental is the interplay
between the filaments and astrophysical turbulence, since turbulence can ultimately
destroy otherwise stable filaments. Our discussion of the presence of filaments as
a result of total energy balance provides some insight into the level of turbulence
allowed when filaments are present, but more testing, especially on the semi-global
and global scales is required.
We find, however, that to fully resolve the physics of the non-linear HBI in
the local regime requires a high grid resolution. Thus, feasible global simulations
of the ICM may require semi-analytic prescriptions for the dynamics of anisotropic
conduction as well as anisotropic viscosity. This study may aid in production of
those models since filaments have a key role in the thermodynamics of anisotropically
conducting plasma.
124
Another motivation for this work was to find, if any exist, robust measures for
convergence with resolution. Given the intimate connection between the topology of
the filaments and more global quantities such as heat conduction, it is not surprising
that we find that adequate resolution of the dynamics of the filaments is needed in
order to reach convergence in the plasma diagnostics such as VHF and equipartition.
The most stringent measure we find is to measure the position of the filaments in
HBI stability space presented in Fig. 5.6. However, a much simpler measure that
also performs well is to see how much the VHF changes with increasing resolution.
Should the HBI be driving dynamics in relatively large regions of real clusters,
numerical studies including this one suggest that it may have a key role in cluster
thermodynamics. While there is clearly an interaction between AGN and the ICM
(e.g. radio lobes, etc.), it is not known exactly how this feedback may couple
effectively to the thermodynamics. This work and K2012 suggest there may still be
a potential for conductive solutions to the cooling flow problem. Certainly at this
point, the field is still in the exploration stage of the fundamental physical processes
which work together to govern the ICM.
M.J.A., C.S.R., and T.B. acknowledge support from the NSF under ARRA
grant AST-0908212. Support for T.B. was in part provided by the National Aero-
nautics and Space Administration through Einstein Postdoctoral Fellowship Award
Number PF9-00061 issued by the Chandra X-ray Observatory Center, which is oper-
ated by the Smithsonian Astrophysical Observatory for and on behalf of the National
Aeronautics and Space Administration under contract NAS8-03060. HPC resources
were provided by the Texas Advanced Computing Center (TACC) (Ranger), as well
as the UMD OIT (Deepthought) and Dept. of Astronomy (Yorp) clusters. The
authors thank Matthew Kunz and Steve Balbus for useful conversations.
125
Chapter 6: Conclusions and a look forward
The work presented in this dissertation focusing on BH accretion flows around
a single BH has addressed several long-standing questions in the field. In the most
self-consistent context (i.e. global disk, GRMHD simulation, transport by physical
processes rather than added viscosity, etc.) to date we have demonstrated that even
a geometrically thin/optically thick disk can accumulate large-scale magnetic flux in-
ward. We have shown that strong poloidal magnetic fields can significantly alter the
radiative and total efficiency of the accretion flows. We’ve shown that magnetically
arrested disks drive powerful magneto-centrifugal winds. We’ve provided supporting
evidence that jets are not truly ‘quenched’ in the high/soft state of XRBs, but that
jet power could simply fall below detection limits. Finally, we’ve shown that the
MAD state provides a channel, literally, through which super-Eddington radiative
flux can escape for potentially arbitrarily high mass accretion rate.
Not only has our work addressed these numerous important questions using
state of the art numerical simulations, but more broadly contribute to a much needed
understanding of large-scale field effects in accretion systems which may embody a
third, primary parameter along with mass accretion rate M˙ , black hole mass, and
observer viewing angle. This could be crucial to understanding hysteresis in XRBs
and jet power in AGN.
However, this work has lead us to pose new problems that necessitate further
numerical/analytic and observational considering.
Mass loss and wind
Ongoing analysis of the simulations presented in Chapters 2 and 3 has, in part,
focused on the powerful wind. Even at radii of <50Rg the amount of mass flowing
126
outward in the wind equals the net mass flow inward to the BH. Even though the
efficiency of a part of the disk at these larger radii may be smaller, because the avail-
able gravitational potential energy is much higher it could contribute significantly
to the total radiative efficiency of the system. The efficacy of measuring efficiency
of a BH accretion flow by normalizing in terms of mass accretion rate onto the hori-
zon falls into question. While this argument is currently speculative in nature, it
strongly motivates further consideration.
The most immediate importance of better quantifying the kinetic and radiative
content of the winds is in the context of understand large-scale feedback. There is
an obvious connection between the winds in these simulations, launched from the
inner disk at 0.05−0.2c, and ultra-fast outflows (UFOs) (Tombesi et al., 2015; King
and Pounds, 2015, and references therein).
Because of the large amount of potential radiative pre-processes by the matter
in the wind it may be key to understanding the spectrum of some systems as well,
particularly in absorption.
Electro-magnetic spectrum of MAD systems
The truncation of the radiating part of the thin disk at the ISCO and the
radiative profile of the NT model is central to BH spin measurements using both
the iron line and continuum fitting methodologies. The simplicity of the iron line
approach makes it attractive for spectral modeling but its robustness is tightly cou-
pled to how closely the spin-dependence of the ISCO and disk follow the model. If,
as we have demonstrated, a thin disk can exist in the MAD state, where the rotation
of the flow near the ISCO is very sub-Keplerian and there is no clearly defined ISCO
(in fact, magnetic choking of the plunging material may enhance pressure gradients
rather than let them fall precipitously as they would free-streaming in the plunging
region), then the iron-line method may be unreliable for MAD systems.
Determining spin by fitting the continuum relies on a good understanding of
the radial distribution of radiation flux. This also remains poorly understood in
thin, MAD systems.
In tidal disruption events, it is still an open question whether the radiation
127
is primarily reprocessed emission or direct from the disk. In the prior case, both
the angular dependence of radiation and optically thick winds can be key to un-
derstanding the observations. The fact that some TDEs seem to be jetted suggests
that at least some systems could be magnetically dominated at some point in their
evolution Piran et al. (2015); Tchekhovskoy et al. (2014).
BH spin / disk misalignment
Coincident with the suit of simulations published in the paper of Chapter 2
was a complimentary set with the same initial conditions and a tilt of the BH spin
axis by ∼ 6◦. Preliminary analysis of these simulations suggests a complicated be-
havior of disk angular momentum that may go through partial alignment episodes,
if not sustained alignment or precession. Because of the clumpy and magnetically
dominated behavior of these accretion flows, a simple analysis of alignment in the
context of traditional studies of the Bardeen-Peterson effect (Bardeen and Petter-
son, 1975; Papaloizou and Pringle, 1983) is difficult. An empirical study of these
simulations is ongoing and may help to shed light on one of the, most important,
but under-studied realms of accretion astrophysics.
Connections to observations
Finally, there has been intriguing recent work connecting the MAD model
to observations (Piran et al., 2015; Zamaninasab et al., 2014; Tchekhovskoy, 2015;
Mocz and Guo, 2015). The latter shows that for a sample of AGN across the
M˙ spectrum, from very sub-Eddington LLAGN to super-Eddington quasars, the
simple relation between disk state and magnetic flux threading the horizon could be
a powerful probe into AGN population modeling. Now is the time for much more
work interweaving theory of magnetically arrested accretion flows and observations.
Building on the foundation of work laid out in this thesis, I plan to continue
addressing these open questions in the context of single-BH accretion, with a short-
term focus on how our current understanding can be applied to binary BH (BBH)
systems. The fundamental processes I’ve studied are present in BBH but with ad-
ditional complications. For instance, a so-called ‘circum-binary’ disk around two
black holes is truncated by tidal interactions with the binary at some radius. In
128
super massive BBH (SMBBH) systems it could be possible for one or both BHs
in binaries to harbor ‘mini’-disks of their own. There are fundamental aspects of
these systems that are still not understood and over the next couple years I will use
numerical simulations to model the behavior and expected observational signatures
of these systems. For SMBBHs these predictions for EM observations are especially
important since, at the time this thesis is written, space-based gravitational wave
observatories are still years away. Until then, ground-based time-domain observa-
tions hold great promise for probing SMBBH systems in conjunction with theoretical
developments.
The fields of black hole accretion and multi-messenger astrophysics promise a
rich and exciting future, of which I am honored to be a part.
129
Bibliography
M. Abramowicz, M. Jaroszynski, and M. Sikora. Title. A&A, 63:221–221, 1978.
M. A. Abramowicz, B. Czerny, J. P. Lasota, and E. Szuszkiewicz. Slim accretion
disks. ApJ, 332:646–658, September 1988. doi: 10.1086/166683.
M. A. Abramowicz, X. Chen, S. Kato, J.-P. Lasota, and O. Regev. Thermal equi-
libria of accretion disks. ApJ, 438:L37–L39, January 1995. doi: 10.1086/187709.
M. J. Avara, C. S. Reynolds, and T. Bogdanovic´. Role of Magnetic Field Strength
and Numerical Resolution in Simulations of the Heat-flux-driven Buoyancy Insta-
bility. ApJ, 773:171, August 2013. doi: 10.1088/0004-637X/773/2/171.
M. J. Avara, J. C. McKinney, and C. S. Reynolds. Efficiency of thin magnetically
arrested discs around black holes. MNRAS, 462:636–648, October 2016. doi:
10.1093/mnras/stw1643.
X.-N. Bai and J. M. Stone. Local Study of Accretion Disks with a Strong Vertical
Magnetic Field: Magnetorotational Instability and Disk Outflow. ApJ, 767:30,
April 2013. doi: 10.1088/0004-637X/767/1/30.
S. A. Balbus. Stability, Instability, and “Backward” Transport in Stratified Fluids.
ApJ, 534:420–427, May 2000. doi: 10.1086/308732.
S. A. Balbus and J. F. Hawley. A powerful local shear instability in weakly magne-
tized disks. I - Linear analysis. II - Nonlinear evolution. ApJ, 376:214–233, July
1991. doi: 10.1086/170270.
S. A. Balbus and J. F. Hawley. Instability, turbulence, and enhanced transport
in accretion disks. Reviews of Modern Physics, 70:1–53, January 1998a. doi:
10.1103/RevModPhys.70.1.
S. A. Balbus and J. F. Hawley. Instability, turbulence, and enhanced transport
in accretion disks. Rev. Mod. Phys., 70:1–53, January 1998b. doi: 10.1103/
RevModPhys.70.1.
130
S. A. Balbus and C. S. Reynolds. Regulation of Thermal Conductivity in Hot Galaxy
Clusters by MHD Turbulence. ApJ, 681:L65, July 2008. doi: 10.1086/590554.
S. A. Balbus and C. S. Reynolds. Radiative and Dynamic Stability of a Dilute
Plasma. ApJ, 720:L97–L101, September 2010. doi: 10.1088/2041-8205/720/1/
L97.
S. A. Balbus and N. Soker. Theory of local thermal instability in spherical systems.
ApJ, 341:611–630, June 1989. doi: 10.1086/167521.
J. M. Bardeen and J. A. Petterson. The Lense-Thirring Effect and Accretion Disks
around Kerr Black Holes. ApJ, 195:L65, January 1975. doi: 10.1086/181711.
K. Beckwith, J. F. Hawley, and J. H. Krolik. Estimating the Radiative Efficiency
of Magnetized Accretion Disks Around Black Holes. ArXiv/0605295, May 2006.
K. Beckwith, J. F. Hawley, and J. H. Krolik. Transport of Large Scale Poloidal Flux
in Black Hole Accretion. ArXiv e-prints, June 2009.
M. C. Begelman. Accreting Black Holes. ArXiv e-prints, October 2014.
M. C. Begelman and P. J. Armitage. A Mechanism for Hysteresis in Black Hole
Binary State Transitions. ApJ, 782:L18, February 2014. doi: 10.1088/2041-8205/
782/2/L18.
V. S. Beskin, N. L. Zakamska, and H. Sol. Radiation drag effects on magnetically
dominated outflows around compact objects. MNRAS, 347:587–600, January
2004. doi: 10.1111/j.1365-2966.2004.07229.x.
R. D. Blandford and D. G. Payne. Hydromagnetic flows from accretion discs and
the production of radio jets. MNRAS, 199:883–903, June 1982.
R. D. Blandford and R. L. Znajek. Electromagnetic extraction of energy from Kerr
black holes. MNRAS, 179:433–456, May 1977.
J. S. Bloom, D. Giannios, B. D. Metzger, S. B. Cenko, D. A. Perley, N. R. Butler,
N. R. Tanvir, A. J. Levan, P. T. O’Brien, L. E. Strubbe, F. De Colle, E. Ramirez-
Ruiz, W. H. Lee, S. Nayakshin, E. Quataert, A. R. King, A. Cucchiara, J. Guil-
lochon, G. C. Bower, A. S. Fruchter, A. N. Morgan, and A. J. van der Horst. A
Possible Relativistic Jetted Outburst from a Massive Black Hole Fed by a Tidally
Disrupted Star. Science, 333:203–, July 2011. doi: 10.1126/science.1207150.
H. Bondi. On spherically symmetrical accretion. MNRAS, 112:195, 1952.
J. R. Buchler and W. R. Yueh. Compton scattering opacities in a partially degen-
erate electron plasma at high temperatures. ApJ, 210:440–446, December 1976.
doi: 10.1086/154847.
J. K. Cannizzo. The Limit Cycle Instability in Dwarf Nova Accretion Disks, pages
6–40. World Scientific Publishing Co, 1993. doi: 10.1142/9789814350976 0002.
131
S. Chandrasekhar. Hydrodynamic and hydromagnetic stability. International Series
of Monographs on Physics, Oxford: Clarendon, 1961.
J.-P. De Villiers and J. F. Hawley. A Numerical Method for General Relativistic
Magnetohydrodynamics. ApJ, 589:458–480, May 2003. doi: 10.1086/373949.
J.-P. De Villiers, J. F. Hawley, J. H. Krolik, and S. Hirose. Magnetically Driven
Accretion in the Kerr Metric. III. Unbound Outflows. ApJ, 620:878–888, February
2005. doi: 10.1086/427142.
J. Dexter and E. Agol. Quasar Accretion Disks are Strongly Inhomogeneous. ApJ,
727:L24, January 2011. doi: 10.1088/2041-8205/727/1/L24.
J. Dexter and E. Quataert. Inhomogeneous accretion discs and the soft states of
black hole X-ray binaries. MNRAS, 426:L71–L75, October 2012. doi: 10.1111/j.
1745-3933.2012.01328.x.
J. Dexter, J. C. McKinney, S. Markoff, and A. Tchekhovskoy. Transient jet for-
mation and state transitions from large-scale magnetic reconnection in black hole
accretion discs. MNRAS, 440:2185–2190, May 2014. doi: 10.1093/mnras/stu581.
B. T. Draine. Physics of the Interstellar and Intergalactic Medium. 2011.
A. C. Fabian. Cooling Flows in Clusters of Galaxies. ARA&A, 32:277–318, 1994.
doi: 10.1146/annurev.aa.32.090194.001425.
A. C. Fabian, R. M. Johnstone, J. S. Sanders, C. J. Conselice, C. S. Crawford,
J. S. Gallagher, III, and E. Zweibel. Magnetic support of the optical emission
line filaments in NGC 1275. Nature, 454:968–970, August 2008. doi: 10.1038/
nature07169.
R. Fender and T. Belloni. GRS 1915+105 and the Disc-Jet Coupling in Ac-
creting Black Hole Systems. ARA&A, 42:317–364, September 2004. doi:
10.1146/annurev.astro.42.053102.134031.
G. B. Field. Thermal Instability. ApJ, 142:531, August 1965. doi: 10.1086/148317.
Jos A. Font. Numerical hydrodynamics and magnetohydrodynamics in general rel-
ativity. Living Reviews in Relativity, 11(7), 2008. doi: 10.1007/lrr-2008-7. URL
http://www.livingreviews.org/lrr-2008-7.
P. C. Fragile, O. M. Blaes, P. Anninos, and J. D. Salmonson. Global General
Relativistic Magnetohydrodynamic Simulation of a Tilted Black Hole Accretion
Disk. ApJ, 668:417–429, October 2007. doi: 10.1086/521092.
T. Fragos, M. Tremmel, E. Rantsiou, and K. Belczynski. Black Hole Spin-Orbit
Misalignment in Galactic X-ray Binaries. ApJ, 719:L79–L83, August 2010. doi:
10.1088/2041-8205/719/1/L79.
132
C. F. Gammie. Efficiency of Magnetized Thin Accretion Disks in the Kerr Metric.
ApJ, 522:L57–L60, September 1999. doi: 10.1086/312207.
C. F. Gammie. The Magnetorotational Instability in the Kerr Metric. ApJ, 614:
309–313, October 2004. doi: 10.1086/423443.
C. F. Gammie, J. C. McKinney, and G. To´th. HARM: A Numerical Scheme for
General Relativistic Magnetohydrodynamics. ApJ, 589:444–457, May 2003. doi:
10.1086/374594.
T. A. Gardiner and J. M. Stone. An unsplit Godunov method for ideal MHD via
constrained transport. Journal of Computational Physics, 205:509–539, May 2005.
doi: 10.1016/j.jcp.2004.11.016.
A. Harten, P.D. Lax, and B. van Leer. On upstream differencing and godunov-type
schemes for hyperbolic conservation laws. SIAM Rev., 25:35–61, 1983.
J. F. Hawley, C. F. Gammie, and S. A. Balbus. Local Three-dimensional Magneto-
hydrodynamic Simulations of Accretion Disks. ApJ, 440:742, February 1995. doi:
10.1086/175311.
J. F. Hawley, X. Guan, and J. H. Krolik. Assessing Quantitative Results in Accretion
Simulations: From Local to Global. ApJ, 738:84, September 2011. doi: 10.1088/
0004-637X/738/1/84.
J. F. Hawley, S. A. Richers, X. Guan, and J. H. Krolik. Testing Convergence for
Global Accretion Disks. ApJ, 772:102, August 2013. doi: 10.1088/0004-637X/
772/2/102.
I. V. Igumenshchev. Magnetically Arrested Disks and the Origin of Poynting Jets:
A Numerical Study. ApJ, 677:317–326, April 2008. doi: 10.1086/529025.
I. V. Igumenshchev. Magnetic Inversion as a Mechanism for the Spectral Transi-
tion of Black Hole Binaries. ApJ, 702:L72–L76, September 2009. doi: 10.1088/
0004-637X/702/1/L72.
I. V. Igumenshchev, R. Narayan, and M. A. Abramowicz. Three-dimensional Mag-
netohydrodynamic Simulations of Radiatively Inefficient Accretion Flows. ApJ,
592:1042–1059, August 2003. doi: 10.1086/375769.
K. Inayoshi, Z. Haiman, and J. P. Ostriker. Hyper-Eddington accretion flows on to
massive black holes. MNRAS, 459:3738–3755, July 2016. doi: 10.1093/mnras/
stw836.
Y.-F. Jiang, J. M. Stone, and S. W. Davis. A Godunov Method for Multidimensional
Radiation Magnetohydrodynamics Based on a Variable Eddington Tensor. ApJS,
199:14, March 2012. doi: 10.1088/0067-0049/199/1/14.
133
Y.-F. Jiang, J. M. Stone, and S. W. Davis. A Global Three-dimensional Radiation
Magneto-hydrodynamic Simulation of Super-Eddington Accretion Disks. ApJ,
796:106, December 2014. doi: 10.1088/0004-637X/796/2/106.
R. P. Kerr and A. Schild. Republication of: A new class of vacuum solutions of
the Einstein field equations. General Relativity and Gravitation, 41:2485–2499,
October 2009. doi: 10.1007/s10714-009-0857-z.
A. King and K. Pounds. Powerful Outflows and Feedback from Active
Galactic Nuclei. ARA&A, 53:115–154, August 2015. doi: 10.1146/
annurev-astro-082214-122316.
S. S. Komissarov and J. C. McKinney. The ‘Meissner effect’ and the Blandford-
Znajek mechanism in conductive black hole magnetospheres. MNRAS, 377:L49–
L53, May 2007. doi: 10.1111/j.1745-3933.2007.00301.x.
J. H. Krolik. Magnetized Accretion inside the Marginally Stable Orbit around a
Black Hole. ApJ, 515:L73–L76, April 1999. doi: 10.1086/311979.
A. K. Kulkarni, R. F. Penna, R. V. Shcherbakov, J. F. Steiner, R. Narayan,
A. Sa¸dowski, Y. Zhu, J. E. McClintock, S. W. Davis, and J. C. McKinney.
Measuring black hole spin by the continuum-fitting method: effect of deviations
from the Novikov-Thorne disc model. MNRAS, 414:1183–1194, June 2011. doi:
10.1111/j.1365-2966.2011.18446.x.
M. W. Kunz, A. A. Schekochihin, S. C. Cowley, J. J. Binney, and J. S. Sanders.
A thermally stable heating mechanism for the intracluster medium: turbulence,
magnetic fields and plasma instabilities. MNRAS, 410:2446–2457, February 2011.
doi: 10.1111/j.1365-2966.2010.17621.x.
M. W. Kunz, T. Bogdanovic´, C. S. Reynolds, and J. M. Stone. Buoyancy Instabilities
in a Weakly Collisional Intracluster Medium. ApJ, 754:122, August 2012. doi:
10.1088/0004-637X/754/2/122.
P.D. Lax and B. Wendroff. Systems of conservation laws. Commun. Pure Appl.
Math., 13:217–237, 1960.
R.J. LeVeque. Nonlinear Conservation Laws and Finite Volume Methods, pages
1–159. Springer, Berlin, Germany, 1998.
C. D. Levermore. Relating Eddington factors to flux limiters. J. Quant. Spec. Ra-
diat. Transf., 31:149–160, February 1984. doi: 10.1016/0022-4073(84)90112-2.
R. Lynds. Improved Photographs of the NGC 1275 Phenomenon. ApJ, 159, March
1970. doi: 10.1086/180500.
B. Margon, H. C. Ford, J. I. Katz, K. B. Kwitter, R. K. Ulrich, R. P. S. Stone, and
A. Klemola. The bizarre spectrum of SS 433. ApJ, 230:L41–L45, May 1979. doi:
10.1086/182958.
134
J. E. McClintock, R. Narayan, and J. F. Steiner. Black Hole Spin via Continuum
Fitting and the Role of Spin in Powering Transient Jets. Space Sci. Rev., 183:
295–322, September 2014. doi: 10.1007/s11214-013-0003-9.
J. E. McClintock et al. Measuring the spins of accreting black holes. Classical
and Quantum Gravity, 28(11):114009, June 2011. doi: 10.1088/0264-9381/28/
11/114009.
M. McCourt, I. J. Parrish, P. Sharma, and E. Quataert. Can conduction induce con-
vection? On the non-linear saturation of buoyancy instabilities in dilute plasmas.
MNRAS, 413:1295–1310, May 2011. doi: 10.1111/j.1365-2966.2011.18216.x.
M. McDonald, S. Veilleux, D. S. N. Rupke, and R. Mushotzky. On the Origin of the
Extended Hα Filaments in Cooling Flow Clusters. ApJ, 721:1262–1283, October
2010. doi: 10.1088/0004-637X/721/2/1262.
M. McDonald, S. Veilleux, D. S. N. Rupke, R. Mushotzky, and C. Reynolds. Star
Formation Efficiency in the Cool Cores of Galaxy Clusters. ApJ, 734:95, June
2011. doi: 10.1088/0004-637X/734/2/95.
J. C. McKinney. Total and Jet Blandford-Znajek Power in the Presence of an
Accretion Disk. ApJ, 630:L5–L8, September 2005. doi: 10.1086/468184.
J. C. McKinney. General relativistic magnetohydrodynamic simulations of the jet
formation and large-scale propagation from black hole accretion systems. MNRAS,
368:1561–1582, June 2006. doi: 10.1111/j.1365-2966.2006.10256.x.
J. C. McKinney and R. D. Blandford. Stability of relativistic jets from rotating, ac-
creting black holes via fully three-dimensional magnetohydrodynamic simulations.
MNRAS, 394:L126–L130, March 2009. doi: 10.1111/j.1745-3933.2009.00625.x.
J. C. McKinney and C. F. Gammie. A Measurement of the Electromagnetic Lumi-
nosity of a Kerr Black Hole. ApJ, 611:977–995, August 2004. doi: 10.1086/422244.
J. C. McKinney and R. Narayan. Disc-jet coupling in black hole accretion systems
- I. General relativistic magnetohydrodynamical models. MNRAS, 375:513–530,
February 2007. doi: 10.1111/j.1365-2966.2006.11301.x.
J. C. McKinney, A. Tchekhovskoy, and R. D. Blandford. General relativistic magne-
tohydrodynamic simulations of magnetically choked accretion flows around black
holes. MNRAS, 423:3083–3117, July 2012. doi: 10.1111/j.1365-2966.2012.21074.x.
J. C. McKinney, A. Tchekhovskoy, and R. D. Blandford. Alignment of Magnetized
Accretion Disks and Relativistic Jets with Spinning Black Holes. Science, 339:
49–, January 2013. doi: 10.1126/science.1230811.
J. C. McKinney, A. Tchekhovskoy, A. Sadowski, and R. Narayan. Three-
dimensional general relativistic radiation magnetohydrodynamical simulation of
super-Eddington accretion, using a new code HARMRAD with M1 closure. MN-
RAS, 441:3177–3208, July 2014. doi: 10.1093/mnras/stu762.
135
J. C. McKinney, L. Dai, and M. Avara. Efficiency of Super-Eddington Magnetically-
Arrested Accretion. ArXiv/1508.02433, August 2015a.
J. C. McKinney, L. Dai, and M. J. Avara. Efficiency of super-Eddington
magnetically-arrested accretion. MNRAS, 454:L6–L10, November 2015b. doi:
10.1093/mnrasl/slv115.
D. L. Meier. The Association of Jet Production with Geometrically Thick Accretion
Flows and Black Hole Rotation. ApJ, 548:L9–L12, February 2001. doi: 10.1086/
318921.
M. C. Miller and E. J. M. Colbert. Title. 2004. Int. J. Mod. Phys. and D13, 1.
B. Mishra, P. C. Fragile, L. C. Johnson, and W. Kluz´niak. Three-dimensional,
global, radiative GRMHD simulations of a thermally unstable disc. MNRAS,
463:3437–3448, December 2016. doi: 10.1093/mnras/stw2245.
P. Mocz and X. Guo. Interpreting MAD within multiple accretion regimes. MNRAS,
447:1498–1503, February 2015. doi: 10.1093/mnras/stu2555.
D. Morales Teixeira, P. C. Fragile, V. V. Zhuravlev, and P. B. Ivanov. Conservative
GRMHD Simulations of Moderately Thin, Tilted Accretion Disks. ApJ, 796:103,
December 2014. doi: 10.1088/0004-637X/796/2/103.
R. Narayan and I. Yi. Advection-dominated accretion: A self-similar solution. ApJ,
428:L13–L16, June 1994. doi: 10.1086/187381.
R. Narayan and I. Yi. Advection-dominated Accretion: Underfed Black Holes and
Neutron Stars. ApJ, 452:710, October 1995. doi: 10.1086/176343.
R. Narayan, I. V. Igumenshchev, and M. A. Abramowicz. Magnetically Arrested
Disk: an Energetically Efficient Accretion Flow. PASJ, 55:L69–L72, December
2003. URL http://adsabs.harvard.edu/abs/2003PASJ...55L..69N.
R. Narayan, A. SA¨ dowski, R. F. Penna, and A. K. Kulkarni. GRMHD simulations
of magnetized advection-dominated accretion on a non-spinning black hole: role
of outflows. MNRAS, 426:3241–3259, November 2012. doi: 10.1111/j.1365-2966.
2012.22002.x.
R. Narayan, Y. Zhu, D. Psaltis, and A. Sad¸owski. HEROIC: 3D general relativis-
tic radiative post-processor with comptonization for black hole accretion discs.
MNRAS, 457:608–628, March 2016. doi: 10.1093/mnras/stv2979.
S. C. Noble, J. H. Krolik, and J. F. Hawley. Direct Calculation of the Radiative
Efficiency of an Accretion Disk Around a Black Hole. ApJ, 692:411–421, February
2009. doi: 10.1088/0004-637X/692/1/411.
S. C. Noble, J. H. Krolik, and J. F. Hawley. Dependence of Inner Accretion Disk
Stress on Parameters: The Schwarzschild Case. ApJ, 711:959–973, March 2010.
doi: 10.1088/0004-637X/711/2/959.
136
S. C. Noble, J. H. Krolik, J. D. Schnittman, and J. F. Hawley. Radiative Efficiency
and Thermal Spectrum of Accretion onto Schwarzschild Black Holes. ApJ, 743:
115, December 2011. doi: 10.1088/0004-637X/743/2/115.
M.L. Norman and K.-H.A. Winkler. Why ultrarelativistic hydrodynamics is dif-
ficult. In M.L. Norman and K.-H.A. Winkler, editors, Astrophysical Radiation
Hydrodynamics: Proceedings of the NATO Advanced Workshop held in Garching,
August 2-13, 1982, volume 188 of NATO ASI Series C, pages 449–476, Dordrecht,
Netherlands, 1986. Reidel.
I. D. Novikov and K. S. Thorne. Astrophysics of black holes. In C. Dewitt and B. S.
Dewitt, editors, Black Holes (Les Astres Occlus), pages 343–450, 1973a.
I. D. Novikov and K. S. Thorne. In Black Holes-Les Astres Occlus, ed. C. De Witt
& B. S. De Witt. New York: Gordon & Breach, 1973b. URL http://adsabs.
harvard.edu/abs/1973blho.conf..343N.
K. Ohsuga and S. Mineshige. Global Structure of Three Distinct Accretion
Flows and Outflows around Black Holes from Two-dimensional Radiation-
magnetohydrodynamic Simulations. ApJ, 736:2, July 2011. doi: 10.1088/
0004-637X/736/1/2.
K. Ohsuga, M. Mori, T. Nakamoto, and S. Mineshige. Supercritical Accretion Flows
around Black Holes: Two-dimensional, Radiation Pressure-dominated Disks with
Photon Trapping. ApJ, 628:368–381, July 2005. doi: 10.1086/430728.
J. C. B. Papaloizou and J. E. Pringle. The time-dependence of non-planar accretion
discs. MNRAS, 202:1181–1194, March 1983. doi: 10.1093/mnras/202.4.1181.
I. J. Parrish and E. Quataert. Nonlinear Simulations of the Heat-Flux-driven Buoy-
ancy Instability and Its Implications for Galaxy Clusters. ApJ, 677:L9, April
2008. doi: 10.1086/587937.
I. J. Parrish, E. Quataert, and P. Sharma. Turbulence in Galaxy Cluster Cores:
A Key To Cluster Bimodality? ApJ, 712:L194–L198, April 2010. doi: 10.1088/
2041-8205/712/2/L194.
I. J. Parrish, M. McCourt, E. Quataert, and P. Sharma. The effects of anisotropic
viscosity on turbulence and heat transport in the intracluster medium. MNRAS,
422:704–718, May 2012. doi: 10.1111/j.1365-2966.2012.20650.x.
R. F. Penna, J. C. McKinney, R. Narayan, A. Tchekhovskoy, R. Shafee, and J. E.
McClintock. Simulations of magnetized discs around black holes: effects of black
hole spin, disc thickness and magnetic field geometry. MNRAS, 408:752–782,
October 2010a. doi: 10.1111/j.1365-2966.2010.17170.x.
R. F. Penna, J. C. McKinney, R. Narayan, A. Tchekhovskoy, R. Shafee, and J. E.
McClintock. Simulations of magnetized discs around black holes: effects of black
137
hole spin, disc thickness and magnetic field geometry. MNRAS, 408:752–782,
October 2010b. doi: 10.1111/j.1365-2966.2010.17170.x.
R. F. Penna, A. Sa¸dowski, and J. C. McKinney. Thin-disc theory with a non-
zero-torque boundary condition and comparisons with simulations. MNRAS, 420:
684–698, February 2012. doi: 10.1111/j.1365-2966.2011.20084.x.
M. E. Pessah, C.-k. Chan, and D. Psaltis. Angular Momentum Transport in Accre-
tion Disks: Scaling Laws in MRI-driven Turbulence. ApJ, 668:L51–L54, October
2007. doi: 10.1086/522585.
T. Piran, A. Sa¸dowski, and A. Tchekhovskoy. Jet and disc luminosities in tidal
disruption events. MNRAS, 453:157–165, October 2015. doi: 10.1093/mnras/
stv1547.
P. Polko and J. C. McKinney. Electromagnetic vs. Lense-Thirring alignment of
black hole accretion discs. ArXiv e-prints, December 2015.
E. Quataert. Buoyancy Instabilities in Weakly Magnetized Low-Collisionality Plas-
mas. ApJ, 673:758-762, February 2008. doi: 10.1086/525248.
R. A. Remillard and J. E. McClintock. X-Ray Properties of Black-Hole Bina-
ries. ARA&A, 44:49–92, September 2006. doi: 10.1146/annurev.astro.44.051905.
092532.
C. S. Reynolds et al. Trapping of Magnetic Flux by the Plunge Region of a Black
Hole Accretion Disk. ApJ, 651:1023–1030, November 2006. doi: 10.1086/507691.
D. M. Rothstein and R. V. E. Lovelace. Advection of Magnetic Fields in Accretion
Disks: Not So Difficult After All. ApJ, 677:1221–1232, April 2008. doi: 10.1086/
529128.
D. M. Russell, J. C. A. Miller-Jones, T. J. Maccarone, et al. Testing the Jet Quench-
ing Paradigm with an Ultradeep Observation of a Steadily Soft State Black Hole.
ApJ, 739:L19, September 2011. doi: 10.1088/2041-8205/739/1/L19.
M. Ruszkowski, H.-Y. K. Yang, and E. Zweibel. Global simulations of galactic winds
including cosmic ray streaming. ArXiv e-prints, February 2016.
G. B. Rybicki and A. P. Lightman. Radiative Processes in Astrophysics. June 1986.
A. Sadowski and R. Narayan. Powerful radiative jets in super-critical accretion disks
around non-spinning black holes. arxiv: 1503.00654, March 2015.
Y. Sakurai, K. Inayoshi, and Z. Haiman. Hyper-Eddington mass accretion on to a
black hole with super-Eddington luminosity. MNRAS, 461:4496–4504, October
2016. doi: 10.1093/mnras/stw1652.
138
G. Salvesen, J. B. Simon, P. J. Armitage, and M. C. Begelman. Accretion disc dy-
namo activity in local simulations spanning weak-to-strong net vertical magnetic
flux regimes. MNRAS, 457:857–874, March 2016. doi: 10.1093/mnras/stw029.
A. Sa¸dowski, R. Narayan, J. C. McKinney, and A. Tchekhovskoy. Numerical simu-
lations of super-critical black hole accretion flows in general relativity. MNRAS,
439:503–520, March 2014. doi: 10.1093/mnras/stt2479.
A. Sa¸dowski, R. Narayan, A. Tchekhovskoy, D. Abarca, Y. Zhu, and J. C. McK-
inney. Global simulations of axisymmetric radiative black hole accretion discs
in general relativity with a mean-field magnetic dynamo. MNRAS, 447:49–71,
February 2015. doi: 10.1093/mnras/stu2387.
A. Sa¸dowski, J.-P. Lasota, and R. Abramowicz, M. A. andNarayan. Energy flows
in thick accretion discs and their consequences for black hole feedback. MNRAS,
456:3915–3928, March 2016. doi: 10.1093/mnras/stv2854.
K. Schwarzschild. ber das Gravitationsfeld eines Massenpunktes nach der Einstein-
schen Theorie. Sitzungsberichte der Kniglich Preussischen Akademie der Wis-
senschaften, 1:189–196, 1916.
C. K. Seyfert. Nuclear Emission in Spiral Nebulae. ApJ, 97:28, January 1943. doi:
10.1086/144488.
R. Shafee, J. C. McKinney, R. Narayan, A. Tchekhovskoy, C. F. Gammie, and J. E.
McClintock. Three-Dimensional Simulations of Magnetized Thin Accretion Disks
around Black Holes: Stress in the Plunging Region. ApJ, 687:L25–L28, November
2008. doi: 10.1086/593148.
N. I. Shakura and R. A. Sunyaev. Black holes in binary systems. Observational
appearance. A&A, 24:337–355, 1973. URL http://adsabs.harvard.edu/abs/
1973A%26A....24..337S.
S. L. Shapiro. Spin, Accretion, and the Cosmological Growth of Supermassive Black
Holes. ApJ, 620:59–68, February 2005. doi: 10.1086/427065.
P. Sharma and G. W. Hammett. Preserving monotonicity in anisotropic diffusion.
Journal of Computational Physics, 227:123–142, November 2007. doi: 10.1016/j.
jcp.2007.07.026.
P. Sharma, B. D. G. Chandran, E. Quataert, and I. J. Parrish. Buoyancy Instabilities
in Galaxy Clusters: Convection Due to Adiabatic Cosmic Rays and Anisotropic
Thermal Conduction. ApJ, 699:348–361, July 2009. doi: 10.1088/0004-637X/
699/1/348.
K. Shibata and Y. Uchida. A magnetodynamic mechanism for the formation of
astrophysical jets. I - Dynamical effects of the relaxation of nonlinear magnetic
twists. PASJ, 37:31–46, 1985.
139
K. Shibata and Y. Uchida. A magnetodynamic mechanism for the formation of
astrophysical jets. II - Dynamical processes in the accretion of magnetized mass
in rotation. PASJ, 38:631–660, 1986.
H. Shiokawa, J. C. Dolence, C. F. Gammie, and S. C. Noble. Global General
Relativistic Magnetohydrodynamic Simulations of Black Hole Accretion Flows:
A Convergence Study. ApJ, 744:187, January 2012. doi: 10.1088/0004-637X/
744/2/187.
F. H. Shu. The physics of astrophysics. Volume 1: Radiation. 1991.
A. Soltan. Masses of quasars. MNRAS, 200:115–122, July 1982.
K. A. Sorathia, C. S. Reynolds, and P. J. Armitage. Connections Between Local
and Global Turbulence in Accretion Disks. ApJ, 712:1241–1247, April 2010. doi:
10.1088/0004-637X/712/2/1241.
K. A. Sorathia, C. S. Reynolds, J. M. Stone, and K. Beckwith. Global Simulations
of Accretion Disks. I. Convergence and Comparisons with Local Models. ApJ,
749:189, April 2012. doi: 10.1088/0004-637X/749/2/189.
L. Spitzer. Physics of Fully Ionized Gases. 1962.
J. M. Stone and M. L. Norman. ZEUS-2D: A radiation magnetohydrodynamics code
for astrophysical flows in two space dimensions. I - The hydrodynamic algorithms
and tests. ApJS, 80:753–790, June 1992.
J. M. Stone, J. F. Hawley, C. F. Gammie, and S. A. Balbus. Three-dimensional Mag-
netohydrodynamical Simulations of Vertically Stratified Accretion Disks. ApJ,
463:656, June 1996. doi: 10.1086/177280.
J. M. Stone, T. A. Gardiner, P. Teuben, J. F. Hawley, and J. B. Simon. Athena:
A New Code for Astrophysical MHD. ApJS, 178:137–177, September 2008. doi:
10.1086/588755.
H. R. Takahashi and K. Ohsuga. Radiation drag effects in black hole outflows from
super-critical accretion disks via special relativistic radiation magnetohydrody-
namics simulations. PASJ, February 2015. doi: 10.1093/pasj/psu145.
S. Takeuchi, K. Ohsuga, and S. Mineshige. A Novel Jet Model: Magnetically Colli-
mated, Radiation-Pressure Driven Jet. PASJ, 62:L43, October 2010.
A. Tchekhovskoy. Launching of Active Galactic Nuclei Jets. In I. Contopoulos,
D. Gabuzda, and N. Kylafis, editors, The Formation and Disruption of Black
Hole Jets, volume 414 of Astrophysics and Space Science Library, page 45, 2015.
doi: 10.1007/978-3-319-10356-3 3.
A. Tchekhovskoy and J. C. McKinney. Prograde and retrograde black holes: whose
jet is more powerful? MNRAS, 423:L55–L59, June 2012. doi: 10.1111/j.
1745-3933.2012.01256.x.
140
A. Tchekhovskoy, R. Narayan, and J. C. McKinney. Efficient generation of jets
from magnetically arrested accretion on a rapidly spinning black hole. MNRAS,
418:L79–L83, November 2011. doi: 10.1111/j.1745-3933.2011.01147.x.
A. Tchekhovskoy, J. C. McKinney, and R. Narayan. General Relativistic Modeling
of Magnetized Jets from Accreting Black Holes. Journal of Physics Conference
Series, 372(1):012040, July 2012. doi: 10.1088/1742-6596/372/1/012040.
A. Tchekhovskoy, B. D. Metzger, D. Giannios, and L. Z. Kelley. Swift J1644+57
gone MAD: the case for dynamically important magnetic flux threading the black
hole in a jetted tidal disruption event. MNRAS, 437:2744–2760, January 2014.
doi: 10.1093/mnras/stt2085.
R. Teyssier. Grid-Based Hydrodynamics in Astrophysical Fluid Flows. ARA&A,
53:325–364, August 2015. doi: 10.1146/annurev-astro-082214-122309.
F. Tombesi, M. Mele´ndez, S. Veilleux, J. N. Reeves, E. Gonza´lez-Alfonso, and C. S.
Reynolds. Wind from the black-hole accretion disk driving a molecular outflow in
an active galaxy. Nature, 519:436–438, March 2015. doi: 10.1038/nature14261.
E. P. Velikhov. Stability of an ideally conducting liquid flowing between rotating
cylinders in a magnetic field. J. Exper. Theor. Phys., 36:1398–1404, 1959.
K. Watarai, T. Mizuno, and S. Mineshige. Title. ApJ, 549, 2001. L77-L77.
J. R. Wilson. Magnetohydrodynamics near a black hole. In R. Ruffini, editor, 1st
Marcel Grossmann Meeting on General Relativity, pages 393–413, 1977.
L. Woltjer. Emission Nuclei in Galaxies. ApJ, 130:38, July 1959. doi: 10.1086/
146694.
H.-Y. K. Yang and C. S. Reynolds. How AGN Jets Heat the Intracluster Medium
- Insights from Hydrodynamic Simulations. ApJ, 829:90, October 2016. doi:
10.3847/0004-637X/829/2/90.
N. L. Zakamska and R. Narayan. Models of Galaxy Clusters with Thermal Conduc-
tion. ApJ, 582:162–169, January 2003. doi: 10.1086/344641.
M. Zamaninasab, E. Clausen-Brown, T. Savolainen, and A. Tchekhovskoy. Dynam-
ically important magnetic fields near accreting supermassive black holes. Nature,
510:126–128, June 2014. doi: 10.1038/nature13399.
Ya. B. Zeldovich. The Fate of a Star and the Evolution of Gravitional Energy upon
Accretion. Soviet Physics Doklady, 9:195, 1964.
Y. Zhu, S. W. Davis, R. Narayan, A. K. Kulkarni, R. F. Penna, and J. E. McClintock.
The eye of the storm: light from the inner plunging region of black hole accretion
discs. MNRAS, 424:2504–2521, August 2012. doi: 10.1111/j.1365-2966.2012.
21181.x.
141
Y. Zhu, R. Narayan, A. Sadowski, and D. Psaltis. HERO - A 3D general relativistic
radiative post-processor for accretion discs around black holes. MNRAS, 451:
1661–1681, August 2015. doi: 10.1093/mnras/stv1046.
142