# Spin and Accretion Rate Dependence of Black Hole X-Ray Spectra

Brooks E. Kinch CCS-2: Computational Physics and Methods
Los Alamos National Laboratory
NM 87545, USA  Jeremy D. Schnittman Gravitational Astrophysics Laboratory,
Code 663
NASA Goddard Space Flight Center
Greenbelt, MD 20771, USA  Scott C. Noble Gravitational Astrophysics
Laboratory, Code 663
NASA Goddard Space Flight Center
Greenbelt, MD 20771, USA  Timothy R. Kallman X-Ray Astrophysics Laboratory,
Code 662
NASA Goddard Space Flight Center
Greenbelt, MD 20771, USA  Julian H. Krolik Department of Physics and Astronomy
Johns Hopkins University
Baltimore, MD 21218, USA

###### Abstract

We present a survey of how the spectral features of black hole X-ray binary
systems depend on spin, accretion rate, viewing angle, and Fe abundance when
predicted on the basis of first principles physical calculations. The power
law component hardens with increasing spin. The thermal component strengthens
with increasing accretion rate. The Compton bump is enhanced by higher
accretion rate and lower spin. The Fe K$\alpha$ equivalent width grows sub-
linearly with Fe abundance. Strikingly, the K$\alpha$ profile is more
sensitive to accretion rate than to spin because its radial surface brightness
profile is relatively flat, and higher accretion rate extends the production
region to smaller radii. The overall radiative efficiency is at least 30–-100%
greater than as predicted by the Novikov-Thorne model.

Magnetohydrodynamical simulations (1966) — General relativity (641) —
Accretion (14) — X-ray binary stars (1811)

††journal: AJ

##  1 Introduction

The X-ray photons which escape from the innermost regions of a black hole
system convey information about the mass and spin of the black hole itself as
well as the properties of the accretion flow in its immediate environment.
X-ray spectra observed from stellar-mass black holes in X-ray binary systems
typically show the following common features: a thermal peak superimposed on a
power law which extends to high energy and eventually rolls over to an
exponential decay, an extremely broad Fe K$\alpha$ fluorescent emission line
centered at 6–7 keV, and a wide “bump” superimposed on the power law above the
K$\alpha$ feature. In recent years, much has been learned from General
Relativistic magnetohydrodynamic (GRMHD) simulations about the dynamics of the
accretion flow: how it is regulated, in part, by the black hole’s mass and
spin, and how, at a qualitative level, it could create these spectral
features.

Despite these advances, the standard picture of the accretion flow geometry
often used for data analysis is little more than a cartoon: a razor thin but
optically thick disk [with inner cutoff at exactly the innermost stable
circular orbit (ISCO) radius] emits blackbody radiation according to a
prescribed radial temperature profile [usually that of Shakura & Sunyaev
(1973)], while it is illuminated by hard X-rays generated in a stationary
“lamppost” located at a single point on the rotation axis some height above
the black hole. In reality, it is virtually impossible for plasma to remain
stationary at such a location, and both physical and observational arguments
point strongly to an extended (Wilkins et al., 2016; Chainakun et al., 2019;
Zoghbi et al., 2020, 2021; Zdziarski et al., 2021) and dynamic (Kara et al.,
2019) corona. Moreover, these features are by no means independent: the hard
X-rays are created by Compton upscattering of thermal disk photons in the
corona, and are “reprocessed” if they strike the disk to create both the
K$\alpha$ line and the Compton bump. The profile of the K$\alpha$ line is
therefore determined by how the hard X-ray flux varies over the surface of the
disk in concert with how the effectiveness of reprocessing varies over that
same surface, as well as Compton scattering in the corona as they escape.

All this physics is well-understood, but the many interactions make a complex
problem. In order to calculate black hole accretion self-consistently from
first principles, one must contend with the relevant physics governing the
processes just described: General Relativity, magnetohydrodynamics, atomic
physics, relativistic Compton scattering, and radiation transport. Our method
for doing so in as fundamental a way as possible is documented in Kinch et al.
(2016, 2019, 2020), and is reviewed briefly in the next section. In short, we
post-process snapshots taken from state of the art three-dimensional GRMHD
simulations using codes which self-consistently solve the equilibrium
radiation transport and thermal plus photoionization balance problems,
accounting for all relativistic effects.

In this paper, we use this established method to investigate how the features
of the synthetic spectra predicted on the basis of this first principles
approach depend on the properties of the modeled black hole system: the mass
and spin of the black hole, the rate at which the inflowing material accretes
onto the black hole, the elemental abundance of Fe, and the angle at which the
disk is inclined relative to the observer’s line of sight.

By beginning with self-consistent physics, this method provides a number of
advantages with respect to the usual methods of spectral analysis, which
involve fitting to data the model spectra generated by codes like reflionx
(Ross & Fabian, 2005), xillver (García & Kallman, 2010; García et al., 2011,
2013), and relxill (García et al., 2014). While the models generated by these
codes do produce good fits to real spectral data, they rely on parameterized
descriptions of the accretion flow geometry, illumination pattern, temperature
structure, etc.—not surprisingly, they therefore require many more parameters
in order to specify a spectrum. By post-processing simulation data, we not
only replace an ad hoc description of the black hole environment with a
physically-motivated one, we also directly connect our smaller set of
parameters to the properties and physics of the real systems we seek to
understand.

##  2 Method

Our starting point is the extremely well-resolved ThinHR series of simulations
(Noble et al., 2010) performed by the 3D general relativistic MHD code harm3d
(Noble et al., 2009). These simulations studied the MHD evolution of an
accretion disk in a Kerr spacetime. The matter was assumed to have an ideal
gas equation of state with adiabatic index $c_{P}/c_{V}=5/3$ [see Schnittman
et al. (2013) for a discussion concerning the choice of equation of state],
cooling according to a target-temperature cooling function in the disk body
(Noble et al., 2011) and an inverse Compton (IC) cooling function in the
corona (Kinch et al., 2020).

The disk body radiates a roughly thermal spectrum into the hot corona, which
is energized by magnetic dissipation. This thermal spectrum is Comptonized in
the corona. Some of the upscattered photons strike the disk, where their
energy is reprocessed, altering its quasi-thermal spectrum.

To predict the shape of the spectrum emerging from this complex feedback
system, we post-processed snapshot output data using the codes pandurata and
ptransx as described in Kinch et al. (2019): the geodesic-integrating, ray-
tracing code pandurata (Schnittman & Krolik, 2013) solves the time-independent
radiation transport problem in the optically thin, hot corona; for the denser,
cooler disk, the plane-parallel radiative transfer solver ptransx (Kinch et
al., 2016), together with xstar (Kallman & Bautista, 2001), finds the thermal-
plus photoionization-balancing seed photon spectrum and disk albedo. Applied
iteratively, these codes build up a globally self-consistent solution for the
radiation field, temperature structure, and ionization balance of the
computational volume, accounting for all relativistic effects, and including
the spectrum as would be seen by a distant observer. A schematic summary of
the simulation-to-spectrum pipeline is shown in Figure 1.

This technique is, at present, the most realistic, first principles-based
means by which to generate theoretical predictions for the X-ray spectra of
black hole systems. Furthermore it is the _only_ method which does so starting
from a description of the accretion flow environment as it develops within a
magnetorotational instability (MRI)-resolving (Balbus & Hawley, 1991; Hawley &
Balbus, 1991) simulation—widely agreed to be the actual mechanism for angular
momentum transport in real disks—while also reproducing all the hallmark
spectral features observed in real moderately-accreting black holes. While
Narayan et al. (2016) also generated synthetic distant observer spectra
starting from simulation output, they did so without the machinery required to
incorporate disk reprocessing effects (e.g., K$\alpha$ lines), and their base
simulations had too little coronal luminosity to produce substantial power law
components in their final spectra.

Each synthetic spectrum our method generates is characterized entirely by only
five parameters: the central black hole mass ($M$), the black hole’s
dimensionless spin parameter ($a$; sometimes written $a^{*}$ or $a/M$) the
mass accretion rate expressed in Eddington units
($\dot{m}=\dot{M}/\dot{M}_{\mathrm{Edd}}$), the Fe abundance relative to its
solar value ($[\mathrm{Fe}]/[\mathrm{Fe}]_{\odot}$; in principle the entire
elemental abundance list is a free parameter, but Fe contributes the most
easily observable line emission in X-ray binary spectra), and the viewer
inclination angle ($i$, where $i=0^{\circ}$ is a face-on view of the disk,
while $i=90^{\circ}$ is an edge-on view). These parameters directly quantify
independent physical or geometric properties of real black hole systems, and
therefore cannot be reduced to an even smaller set. The process by which scale
free simulation results are translated into physical (cgs) units for post-
processing is thoroughly detailed in Kinch et al. (2019, 2020): for the
purposes of this paper, it suffices to explain that the ThinHR simulations are
run choosing only a value for the dimensionless spin $a$ until they achieve
approximate inflow equilibrium. At that point, the accretion rate must be
chosen in order to employ the IC cooling function, and the simulations are run
for an additional $1000M$ to achieve a new quasi-steady state. A choice for
$M$ at this stage allows the full translation of these last $1000M$ of
simulation data into cgs units for post-processing, which also requires
specification of the Fe abundance. Finally, because pandurata is an emitter-
to-observer ray-tracing code, spectra for all viewing angles $i$ are computed
simultaneously.

An important technical caveat we briefly discuss here [again see Kinch et al.
(2020) for details] are harm3d’s “zero cooling” conditions. Each time step,
harm3d must recover a series of familiar “primitive” variables (density,
internal energy, pressure, etc.) from a more natural (in GR) set of
“conserved” variables—those quantities which are actually evolved by harm3d.
The conserved-to-primitive recovery procedure is performed via a nonlinear
solver, and is subject to failure in extremely magnetically-dominated cells.
In these regions, a separate entropy-conserving integration is performed. The
thresholds beyond which we choose to invoke the entropy equation are (in
dimensionless code units) $B^{2}/u>10^{4}$ (i.e, below a minimum
plasma-$\beta=(c_{P}/c_{V}-1)8\pi 10^{-4}$, where $c_{P}/c_{V}$ is the
adiabatic index) or $B^{2}/\rho>1$, where $B^{2}$ is the magnetic energy
density, $u$ is the internal energy density, and $\rho$ is the (code units)
density. In practice, these conditions occur only near the evacuated polar
“jet” region of the simulations. In addition, while these simulations are
spatially well-resolved in the disk and corona, the rapid changes in density
and temperature as a function of polar angle approaching the $z$-axis are
difficult to resolve with any standard grid [see Kinch et al. (2020) for
details]. The net result of these limitations are patches of uncooled and
often unphysically hot material near the $z$-axis. However, as we will
demonstrate below, the unreliable thermodynamic treatment of the jet gas _does
not_ affect the synthetic X-ray spectrum we generate from a given set of
simulation data. These choices are made for the stability of the GRMHD
integrator, and these thresholds are not part of our parameter set because
they do not affect the spectral region of interest. A future study, focusing
on emission from the jet region, would require a more physically robust
treatment thereof: pair production must be included to correctly treat the
thermodynamics of the extremely hot, diffuse plasma [see Wong et al. (2021)
for an implementation in a harm-like code]; additionally, a scheme like
patchwork (Shiokawa et al., 2018) might be used to better resolve the
jet/corona boundary.

Figure 1: A schematic representation of our pandurata+ptransx pipeline, which
generates synthetic model spectra out of harm3d simulation data.

##  3 Results

###  3.1 Simulation Output

In Figures 2–6, we show various azimuthally-averaged simulation quantities
$1000M$ after switching on the IC coronal cooling function; also shown is the
$\phi$-averaged location of the upper and lower photosphere surfaces. Note
that as spin increases, the photosphere surfaces move closer to the midplane.
This is due, in part, to how the density is scaled from dimensionless “code”
units to physical (cgs) units: $\rho\propto\dot{m}/\eta_{\mathrm{NT}}$, where
$\dot{m}$ is the accretion rate in Eddington units and $\eta_{\mathrm{NT}}$ is
the Novikov & Thorne (1973) radiative efficiency ($L=\eta\dot{M}c^{2}$)
calculated using analytic accretion disk theory, assuming a zero torque
boundary condition at the ISCO. See Table 1 for the values of
$\eta_{\mathrm{NT}}$ for the spins considered here.

As the spin increases, the density in the jet cone decreases (Figure 2),
leading to an increasing number of coronal cells subject to evolution via the
entropy equation and therefore not cooled (Figure 3). This is due to an
increase in the magnetization $B^{2}/\rho$ as shown in Figure 4; the
plasma-$\beta$ in fact increases with increasing spin, shown (as its
reciprocal quantity $B^{2}/u$) in Figure 5. The net result is that the jet
cone becomes more unphysically hot as the spin increases, as shown in Figure
6.

Figure 2: The $\phi$-averaged density structure of each simulation, using a
snapshot $1000M$ after switching on the IC cooling function. The white lines
indicate the upper and lower photosphere surfaces. The simulations are run
with the accretion rate and spin indicated, all with a central black hole mass
of $10M_{\odot}$. The radial extent of these figures is $30M$, and each are
plotted using the same color scale.  Figure 3: The $\phi$-averaged cooling
rate for each simulation at $t=+1000M$. The cooling rate used is the target-
temperature rate in the disk (between the photosphere surfaces) and the IC
cooling rate in the corona (outside the photosphere surfaces). The white
regions indicate zero cooling _in all $\phi$ cells_ for this snapshot.  Figure
4: The $\phi$-averaged magnetization (magnetic energy density divided by mass-
energy density) for each simulation at $t=+1000M$. On this scale, blue
indicates $\langle B^{2}/\rho\rangle>1$, where we expect many cells to have
zeroed-out cooling rates.  Figure 5: The $\phi$-averaged ratio of magnetic
energy density to internal energy density (equal to $(c_{P}/c_{V}-1)8\pi$
divided by the plasma-$\beta$) for each simulation at $t=+1000M$.  Figure 6:
The $\phi$-averaged electron temperature for each simulation at $t=+1000M$.
For cells within the disk body, the effective temperature at the corresponding
photosphere is used for $T_{e}$, adjusted according to our detailed
photoionization calculation just inside the photosphere.

The unphysical thermodynamics within the jet region, however, _do not_ affect
the X-ray spectra we calculate from these simulations. In Figure 7 we show how
the corona’s total cooling budget is distributed among gas of different
temperature; likewise, we show in Figure 8 how the corona’s mass varies in
temperature. The general trend is a greater portion of the coronal cooling and
material at higher temperatures for the higher spin runs—though the
$\dot{m}=0.1$ simulation results display a more complex structure relative to
the $\dot{m}=0.01$ results. Note that while a fair share of the total coronal
cooling is due to gas above 1 MeV, this same gas accounts for essentially none
of the coronal mass budget. As explained in Kinch et al. (2020),
$L_{\mathrm{IC}}\propto\rho\Theta_{e}^{2}$ for $\Theta_{e}\gg 1$, where
$\Theta_{e}$ is the dimensionless electron temperature
($\Theta_{e}=k_{B}T_{e}/m_{e}c^{2}$; $m_{e}c^{2}=511$ keV); thus for the super
hot gas which resides in the evacuated polar regions, the inverse Compton
cooling rate (when nonzero) is often extreme, even when the density is very
low.

Figure 7: For each of the six simulations as indicated, the curves represent
the relative fraction of the total coronal cooling at a given electron
temperature. These harm3d-calculated cooling rate values are integrated over
the coronal volume and averaged over the entire $1000M$ run. All six curves
are normalized so that the integrals with respect to $\log T_{e}$ are unity,
for ease of comparison.

Figure 8: Like Figure 7, but the curves indicate the relative fraction of the
total coronal mass at a given electron temperature. Note that while a
substantial fraction of the total cooling budget is due to gas hotter than 1
MeV, it accounts for essentially none of the corona’s _mass_.

In the next section we will discuss X-ray spectra at length—here we briefly
explore how the anomalously hot polar gas affects spectral predictions. In
Figure 9 we show, for two of the six simulations, the results of post-
processing the $t=+1000M$ snapshots with the pandurata+ptransx pipeline
described above. After its coronal $T_{e}$ map is determined, a value
$\Theta_{e}^{\mathrm{cut}}$ is imposed before the final ray-tracing operation:
inverse Compton scattering events are not simulated in cells with
$k_{B}T_{e}/m_{e}c^{2}>\Theta_{e}^{\mathrm{cut}}$. As we see, increasing
$\Theta_{e}^{\mathrm{cut}}$—i.e., allowing pandurata to probe hotter
gas—affects only the very high energy portion of the spectrum. Because this
hot gas is so diffuse, pandurata _under_ -samples it even if it more than
sufficiently samples the cooler, nearer to the disk regions of the corona.
Increasing the number of photon packets by a factor of 20 leads to an observed
spectrum that is indistinguishable from one constructed with the standard
number all the way from 0.1 keV to 1 MeV; even at 2 MeV, the departure is only
10%.

In practice, extraordinarily rare but large amplification scattering events
can prevent pandurata from achieving a global solution—this is because, with
the current code structure, it is not computationally feasible to simulate
enough scattering events in the polar regions to reliably reconstruct a
meaningful (Monte Carlo) cooling rate to then match to the (analytic) harm3d
cooling rate. A biasing/splitting procedure could be employed to better sample
these rare events—a method specialized for relativistic Compton scattering is
discussed in Kinch & Cleveland (in press). However, because we are confident
that the spectrum below 200 keV is unaffected by these rare events, we proceed
by simply applying the no-scattering rule for any cell with $\Theta_{e}>1$.

Figure 9: The post-processed spectral luminosity (measured at infinity) for
the two simulations indicated. For each, the spectra are shown after applying
three different values of $\Theta_{e}^{\mathrm{cut}}$ in pandurata’s ray-
tracing calculation, as described in the text. The value chosen for
$\Theta_{e}^{\mathrm{cut}}$ affects the spectra only above 200 keV.

In Figure 10, we show the ray-traced luminosity and volume-integrated cooling
rates as functions of time for each of the six simulations. Note that the
corona is relatively more luminous than the disk as the spin increases, and
relatively less so at higher accretion rate. This trend is summarized by the
time-averaged $L_{\mathrm{disk}}/L_{\mathrm{tot}}$ values listed in Table 1.
The higher spin runs tend to have more variable cooling rates, especially
along the jet cone boundary. The prominent spikes in the $\dot{m}=0.1$,
$a=0.5$ run, as well as the more numerous smaller spikes in the $a=0.9$ runs,
are due to very hot material falling below the magnetization threshold for
evolution via the entropy equation and therefore rapidly cooling down
following a sustained period of zero cooling. Figure 10 also shows the ray-
traced value of the total luminosity, i.e., at what luminosity distant
observers would measure the system. The ray-traced luminosity is always less
than the volume-integrated value because the former accounts for gravitational
redshift, capture of photons by the black hole, and—especially so for the
higher spin runs—an increasing fraction of the total cooling occurring in the
(deliberately or otherwise) inadequately sampled hot, diffuse polar region.
The ray-tracing is performed only every $100M$ (whereas the three other curves
in Figure 10 are sampled every $1M$) starting at $t=+100M$, so these curves
appear artificially less variable. The “spiking” behavior typically occurs in
gas at sufficiently high temperature that it does not affect the ray-traced
X-ray luminosity.

Figure 10: The ray-traced luminosity, as well as the volume-integrated cooling
rates, expressed in ratio to the Eddington luminosity, as functions of time
after switching on the IC cooling function. The dotted lines indicate the
nominal luminosity expected for the choice of accretion rate.

Figure 11 shows the time-averaged, net mass inflow rate as a function of
radius (in units of the ISCO radius) for each of the six simulations. Time-
averaged mass inflow equilibrium is established out to at least a few
multiples of the ISCO radius for all runs.

Figure 11: The time-averaged net mass inflow rate, expressed in ratio to the
Eddington accretion rate (assuming $\eta=\eta_{\mathrm{NT}}$). The thin lines
above and below indicate plus-or-minus one standard deviation from the time-
averaged value. The dotted lines indicate the nominal net mass inflow rate
expected for this choice of accretion rate scaling (again, assuming
$\eta=\eta_{\mathrm{NT}}$). The $x$-axes are expressed in multiples of the
ISCO radius: for $a=0$, 0.5, and 0.9, $r_{\mathrm{ISCO}}/M=6$, 4.23, and 2.32,
respectively.

The time-averaged values for ray-traced luminosity, the total volume-
integrated cooling rate, and the mass inflow rate (measured at the event
horizon) allow us to calculate post hoc values for the radiative efficiency,
listed in Table 1. The volume-integrated cooling rate value is reported down
the “$\eta$” column; the values computed using the ray-traced luminosity are
labeled “$\eta_{\mathrm{rt}}$,” and the classical geodesic thin disk
efficiency given by Novikov & Thorne (1973) is “$\eta_{\mathrm{NT}}$.” Note
that the $\eta$ values tend to be a factor $\sim$ 1.5–3 times the
$\eta_{\mathrm{NT}}$ values, while the $\eta_{\mathrm{rt}}$ values are not
quite so large (but still always greater than $\eta_{\mathrm{NT}}$). The
values for $\eta_{\mathrm{rt}}$ are lower than the direct volume-integrated
$\eta$ values because, as discussed above, there is a substantial fraction of
the total cooling budget which is not reflected in the ray-tracing procedure,
especially for high spin, as the low density in the jet cone inhibits adequate
Monte Carlo sampling of the inverse Compton power. A secondary but much
smaller effect is the capture of photons by the black hole, which is accounted
for in determining $\eta_{\mathrm{rt}}$, but not in $\eta$ or
$\eta_{\mathrm{NT}}$ [though see Schnittman et al. (2016)]. Thus
$\eta_{\mathrm{rt}}$ is in this sense a lower bound for the predicted
radiative efficiency.

In addition, the inferred radiative efficiency values are a strong increasing
function of the spin (like $\eta_{\mathrm{NT}}$) but vary little with the
nominal accretion rate, as is expected since the target-temperature cooling
function, the IC cooling function, and the mass inflow rate all depend
linearly on density. The distribution of cooling between corona and disk,
however, does depend on the nominal accretion rate choice—as we see below,
this division significantly impacts the resulting X-ray continuum spectra.

Table 1: Summary of Key Simulation Results:

$\dot{m}$ | $a$ | $L_{\mathrm{disk}}/L_{\mathrm{tot}}$ | $\eta$ | $\eta_{\mathrm{rt}}$ | $\eta_{\mathrm{NT}}$
---|---|---|---|---|---
0.01 | 0 | 0.357 | 0.1284 | 0.1086 | 0.0572
0.01 | 0.5 | 0.319 | 0.2460 | 0.2096 | 0.0821
0.01 | 0.9 | 0.131 | 0.4164 | 0.2170 | 0.1558
0.1 | 0 | 0.612 | 0.0917 | 0.0910 | 0.0572
0.1 | 0.5 | 0.407 | 0.2739 | 0.2094 | 0.0821
0.1 | 0.9 | 0.299 | 0.4010 | 0.2292 | 0.1558

Note. — For the six simulations performed, characterized by their accretion
rate ($\dot{m}$) and spin ($a$), we report time-averaged values for the
fraction of total luminosity in the disk
($L_{\mathrm{disk}}/L_{\mathrm{tot}}$), the inferred radiative efficiency
using the volume-integrated cooling rate ($\eta$), and the inferred radiative
efficiency using the ray-traced luminosity ($\eta_{\mathrm{rt}}$). For
comparison, the Novikov-Thorne radiative efficiency for each spin value
($\eta_{\mathrm{NT}}$) is also listed. These values are calculated for a
$10M_{\odot}$ black hole.

###  3.2 Predicted Spectra

In Figure 12, we show the spectral intensity for a series of observers at
infinity (characterized by their viewing angle $i$), as calculated by our
pandurata+ptransx post-processing pipeline for snapshots of each of the six
simulations $1000M$ after switching on the IC coronal cooling function.
Superimposed on each spectrum is a power law fit to the NuSTAR X-ray band,
3–79 keV, with associated photon index $\Gamma$. The overall trend is
immediately apparent: the spectrum _hardens_ with increasing spin, and
_softens_ with increased accretion rate. The hardening with spin is due to the
increased coronal temperature at higher spins (see Figure 6), while the
softening with accretion rate is due to the larger share of total cooling
occurring in the disk (and therefore radiating thermally) at higher accretion
rate (see Table 1).

Not surprisingly, the smaller solid angle the disk subtends for the more
inclined viewers results in a diminished overall flux as $i$ increases from
$0^{\circ}$ to $90^{\circ}$. Also, the peak hardness typically occurs for
viewers inclined at $\sim 60^{\circ}$; this is due to the larger scattering
optical depth traversed by photons emerging more nearly parallel to the plane
of the disk (relative to those that escape through, e.g., the extremely
optically thin jet cone). This trend is less well-established for the
$\dot{m}=0.1$ runs, in part because these spectra more substantially deviate
from power laws in the NuSTAR band. This deviation is depicted in Figure 13 by
dividing the distant observer spectra by their power law fits.

The prominent Fe K$\alpha$ emission feature at $\sim 6$ keV we discuss further
below. Note that the “Compton bump” between 20–50 keV is enhanced for the
higher accretion rate runs. This feature arises from high energy (already
Compton upscattered) corona photons undergoing Compton downscattering when
they reflect off the cool disk, and is enhanced further by the Fe K-edge
photoionization absorption opacity above $\sim 10$ keV—the absorption of
Compton downscattering photons by Fe forms the boundary between the Fe
K$\alpha$ emission feature and the Compton bump.

The elevated photosphere surfaces of the higher accretion rate runs (recall
Figure 2) increase the likelihood of coronal photons colliding with the disk
surface while traversing the corona; in addition, the softer incident spectra
striking the cooler outer disk layers results in a greater fraction of the
disk-incident photons undergoing K-edge absorption. Both of these effects
contribute to the enhanced Compton bumps present in the $\dot{m}=0.1$ spectra.
Note that the $\dot{m}=0.1$, $a=0.9$ spectra are so poorly-described by a
power law in the NuSTAR band that the bump is not as well-defined in this
view.

Figure 12: The spectral intensity as viewed by distant observers at a sampling
of inclination angles, for the $t=+1000M$ snapshots from each of the six
simulations as indicated. The dashed red lines are power law fits to the
NuSTAR band for each spectrum, with associated photon index reported.  Figure
13: The distant observer spectra divided by their power law fits at a sampling
of viewing angles. The feature at $\sim$ 6 keV is the Fe K$\alpha$ emission
line, while the broader deviation at 20–50 keV is the “Compton bump.”

Figure 14 shows the “true” Fe K$\alpha$ emission line feature. We define
“true” Fe K$\alpha$ emission as those line photons originating from a bound-
bound transition with photon energies in the range 6.3–7.0 keV; by tagging
these photons at creation in ptransx and tracking them through pandurata, we
are able to identify the “line” and “continuum” components of a predicted
spectrum by their physical origin, without resorting to a fitting procedure
for isolating the continuum. The corresponding “true” equivalent width of
these K$\alpha$ emission features all lie in the range 100–500 eV, which are
in line with typical observed strengths.

The line profiles shown in Figure 14 are analogous to a typical continuum-
subtracted representation, but with the continuum known exactly (Figure 13 is
more akin to what is usually done in practice, i.e., a fit to the underlying
continuum). Though the line profiles we predict are relativistically-
broadened, the occasionally-observed “double-horned” feature (due to Doppler
shifting of advancing and receding material in the disk) is essentially
completely smoothed over by Compton scattering in the corona.

Figure 14: The “true” continuum-subtracted Fe K$\alpha$ line profiles, for a
sampling of viewer inclination angles. Only those photons whose physical
origin is a bound-bound transition in the range 6.3–7.0 keV are included.

Unlike the continuum, the Fe K$\alpha$ emission profiles we calculate do not
show strong, immediately apparent trends. The red wing in particular—often the
primary driver for spin measurements in X-ray reflection spectroscopy
techniques—is nearly unchanged as a function of spin when normalized for total
K$\alpha$ emission. Figure 15 shows how the disk surface emission of Fe
K$\alpha$ varies with multiples of the ISCO radius: a typical emission model
assumes a power law-decaying surface brightness profile with a sharp inner
cutoff at the ISCO radius. While our predicted surface brightness profiles do
in fact terminate near (but not exactly at) the ISCO radius (and the
termination is more abrupt at $\dot{m}=0.01$ than for $\dot{m}=0.1$), they do
not decline rapidly with increasing radius—rather, they vary haphazardly and
only gently decline with radius (except for the $a=0.9$ curves, which are
roughly flat), reflecting the turbulent accretion flow and the spatially-
extended corona. The most substantial trends concerning the Fe K$\alpha$
emission are: its equivalent width and peak contrast tend to increase with
accretion rate, and the _blue_ wing tends to strengthen at higher spin.

Figure 15: The $\phi$-, top- and bottom-averaged Fe K$\alpha$ surface
brightness for each of the six simulations at $t=+1000M$. The $x$-axis is
expressed in multiples of the ISCO radius: for $a=0$, 0.5, and 0.9,
$r_{\mathrm{ISCO}}/M=6$, 4.23, and 2.32, respectively. While K$\alpha$
production terminates _near_ the ISCO radius, it does not do so precisely _at_
the ISCO, nor does it decline steeply with increasing radius as is usually
assumed with a “lamppost”-style corona model. Moreover, the location and
sharpness of the inner K$\alpha$ cut-off depend on accretion rate and spin.

Thus far we have used a solar elemental composition. In Figure 16, we show the
K$\alpha$ line profiles which result from post-processing the $t=+1000M$
snapshot of the $\dot{m}=0.01$, $a=0$ simulation, but with the Fe abundance at
half, equal to, and three times the solar value. As one might expect, the Fe
K$\alpha$ line strength increases with Fe abundance, though sub-linearly.
Above some abundance the K$\alpha$ strength would saturate, as the
availability of K-edge photons (not Fe atoms) would become the limiting factor
for K$\alpha$ production. The underlying continuum (not shown) does change
slightly with Fe abundance, as the K-edge opacity and corresponding
photoionization heating affects the temperature structure of the disk.

Figure 16: The “true” (continuum-subtracted) Fe K$\alpha$ line profiles for
the same simulation ($\dot{m}=0.01$, $a=0$), at the same viewing angle
($i=30^{\circ}$), for three different values of the Fe abundance. Note that
the strength of the feature increases sub-linearly with the Fe abundance;
scaling to the overall K$\alpha$ flux, the line profile _shapes_ are
essentially identical.

So far we have always set $M=10M_{\odot}$. Figure 17 shows the post-processed
spectrum of the $\dot{m}=0.01$, $a=0$ run, with central black hole masses
equal to 3, 10, and 30 $M_{\odot}$. On the scale shown (relative to the
Eddington luminosity, which also scales with mass), the change is minor: as
the mass increases, the thermal peak recedes—as expected, as
$T_{\mathrm{disk}}\propto M^{-1/4}$—while the power law component hardens
slightly due to the lower and redder thermal peak. The effect of central black
hole mass on the Fe K$\alpha$ feature is very small: the Fe K$\alpha$
equivalent width changes by only $\lesssim 10\%$ over the given mass range,
and not in a uniform way; additionally, the continuum-subtracted line shape
remains nearly unchanged after controlling for the overall linear scaling of
flux with mass.

Figure 17: The spectral luminosity for a distant observer, for the $t=+1000M$
snapshot of the $\dot{m}=0.01$, $a=0$ simulation, scaled to the three
indicated values of the central black hole mass.

##  4 Discussion

As emphasized in the introduction and also in our preceding methods and “proof
of concept” papers, our 3D GRMHD simulation-derived spectral predictions are
the most directly connected to the underlying physics of any X-ray binary
spectral model available. The disk (thermal) and coronal (power law)
contributions to the continuum, the strength and shape of the Fe K$\alpha$
emission line, the magnitude of the Compton bump reflection feature, and
location of the high energy rollover—all these features emerge naturally from
the self-consistent treatment of the physical processes responsible, against a
background accretion flow geometry provided by high resolution, MRI-resolving
simulations. Because the parameters of our model all specify the actual
properties of the black hole system—and the spectral properties follow suit
without additional parameterization—we have made _forward_ predictions for how
we expect these underlying properties to affect observations. Below we review
the key results of this survey:

  1. 1.

_The Fe K $\alpha$ line profile shape is much more strongly affected by the
accretion rate than by the spin._

Contrary to the usual assumption in X-ray spectral modeling, the inner radial
cutoff of K$\alpha$ emission we predict (for the range of accretion rates and
spins we consider) is always outside the ISCO, though it does move inward (and
terminates less sharply) for the higher accretion rate runs. The spatially-
extended coronae our simulations develop illuminate the disks uniformly: the
result is a radial K$\alpha$ surface brightness profile which, outside the
cutoff, is fairly flat, in contrast to the usually-assumed steep profiles
(Wilkins & Fabian, 2012; Dauser et al., 2013) associated with lamppost
geometries. Flat radial profiles result in relatively unchanging K$\alpha$ red
wings with increasing spin as the fraction of the total K$\alpha$ photons
which originate deep in the potential increases only a very small amount.

The K$\alpha$ profiles do, however, broaden and strengthen for the higher
accretion rate runs. Consider the line profiles in Figure 14: the effect of
raising the accretion rate of a non-spinning black hole is to enhance the red
wing by a greater degree than does increasing the spin to 0.9. This is because
the location of the inner cutoff is in fact determined not _only_ by the spin
(through the ISCO radius), but also by the overall density scale as set by the
accretion rate. An overall greater disk column density means more radiating Fe
atoms at smaller radii, thereby pushing the inner cutoff farther inward.

As expected, the blue wing of the Fe K$\alpha$ line profile tends to
strengthen (relative to the underlying continuum) at more edge-on views, due
to relativistic beaming of K$\alpha$ photons from the rapidly rotating inner
disk (Fabian et al., 1989) and the corresponding enhanced electron scattering
optical depth along more inclined lines of sight.

The Fe K$\alpha$ line strength (as measured by its equivalent width) increases
sub-linearly with Fe abundance. For the example shown, raising the Fe
abundance by a factor of 6 results in an equivalent width which grows only by
a factor of 2.5. We note also that the K$\alpha$ line profile _shape_ is
nearly identical after controlling for total K$\alpha$ photons seen by a given
distant observer—the relative distribution of K$\alpha$ production over the
disk surface is not substantially affected by the Fe abundance.

  2. 2.

_The overall radiative efficiency is at least 30–100% greater than predicted
by the Novikov-Thorne model._

From analytic accretion disk theory, we expect an overall increase in the
radiative efficiency at higher spin—we find the same trend in our simulations,
though the radiative efficiency values we calculate are substantially greater
than the analytic results (compare the $\eta_{\mathrm{rt}}$ values we infer
from our calculations to the $\eta_{\mathrm{NT}}$ column in Table 1).

The classical radiative efficiency of the Novikov-Thorne model is based upon
the assumption of zero stress at the ISCO radius, which automatically
identifies the radiative efficiency with the binding energy of a particle on
an ISCO orbit. The justification for this assumption has been challenged in
many papers (Thorne, 1974; Krolik, 1999; Gammie, 1999; Noble et al., 2009,
2010; Avara et al., 2016) on the grounds that magnetic stresses are likely to
be continuous across the ISCO. Based solely on these grounds, the actual
efficiency has previously been estimated to be modestly elevated above NT
(Noble et al., 2009, 2011; Schnittman et al., 2016) when the magnetic flux
near the horizon is not too large, and possibly elevated rather more if it is
very large [Avara et al. (2016), but see Morales Teixeira et al. (2018) for a
contrary view]. The enhancement we find here is much larger than in this
previous work, even though our methods are very close to those of Noble et al.
(2009). We therefore attribute this change to the one new element relative to
our previous work: the physical cooling function applied in the corona.

A similar result was reported for the non-spinning case in Kinch et al.
(2020)—as explained there, the more physically realistic inverse Compton
coronal cooling function turns out to be much more efficient at radiating away
heat than was the old target-temperature formulation. The result is that more
magnetic heating in the corona is radiated away as upscattered photons instead
of being advected through the event horizon or carried out of the simulation
volume. Importantly, this is a substantial increase in radiative efficiency
above NT in the presence of only moderately-strong magnetic fields, in
contrast to the “magnetically arrested disk” state typically associated with
similar increases above NT (Avara et al., 2016).

  3. 3.

_The power law component of the X-ray spectrum hardens with increasing spin,
while the thermal component strengthens with increasing accretion rate. The
Compton bump is enhanced by higher accretion rate and lower spin._

These are the primary ways the X-ray continuum depends on the parameters
varied. In our simulation-to-spectrum pipeline, the code units density
converts to cgs units according to $\rho\propto\dot{m}/\eta_{\mathrm{NT}}$,
where $\eta_{\mathrm{NT}}$ is an increasing function of the dimensionless spin
(listed in Table 1). Even if we scaled the density using our inferred $\eta$
values instead, the trend with spin is similar. This scaling relationship of
density with accretion rate is a central assumption in our calculation, though
a roughly linear relationship is supported by recent simulations which include
radiation forces (Jiang et al., 2019). For a fixed $\dot{m}$ value, the
physical density scale decreases with increasing spin, resulting in
photosphere surfaces which are closer to the midplane, and therefore contain a
smaller portion of the (overall more efficiently cooling) simulation volume.
For increasing $\dot{m}$ values at a fixed spin, the photospheres move farther
from the midplane, and a greater portion of the (again increased) total
cooling is captured by the disk. In this way, values for $\dot{m}$ and $a$
dictate the division between disk and corona and, thereby, the broad trends of
the X-ray continuum spectum.

The deviation from a pure power law component, i.e., the magnitude of the
Compton bump, is determined by the interaction between the disk’s outer layers
and the disk-incident coronal X-ray flux. It is a complex calculation handled
by the iterations between ptransx and pandurata. Broadly speaking, the
elevated photosphere surfaces at higher $\dot{m}$ result in more coronal
photon trajectories intercepted by the disk, and therefore a greater portion
of the photons which reach infinity having experienced some Compton
downscattering when reflecting off the cooler disk surface.

  4. 4.

_The power law slope varies non-monotonically with viewer inclination angle._

The overall flux declines with increasing viewer inclination (i.e., toward
more edge-on views), with peak hardness occurring for viewers at $i\sim
60^{\circ}$, largely independent of the other parameter values. This is due to
a larger Compton scattering optical depth through the extended corona for
viewers so inclined, relative to the face-on or edge-on observers (who are
more likely to see only those photons which escape not having scattered even
once). The structure of the corona in our 3D GRMHD simulations develops
naturally as a consequence of the buoyancy of magnetic fields. Despite the
persistence of point source lamppost-style modeling in the analysis of time-
integrated spectra, the spatially-extended coronal structure we find enjoys
considerable observational support in the form of variability studies (Wilkins
et al., 2016; Chainakun et al., 2019; Zoghbi et al., 2020, 2021).

##  5 Conclusion

To sum up: the continuum features are more clearly and dramatically affected
by varying the black hole spin and accretion rate than are the K$\alpha$
features. Importantly, the red wing of the Fe K$\alpha$ profile—long used as a
proxy for the black hole spin in X-ray reflection spectroscopy studies (Miller
& Miller, 2015)—actually depends more on the accretion rate than it does on
the spin. In fact, we conclude from our parameter study that an observer
interested in _measuring_ values for these underlying black hole properties
using an X-ray spectrum would be best off relying on the relative strengths of
the thermal to power law components, and the power law’s photon index, in
order to constrain the accretion rate and the spin, while using the overall
magnitude and the relative strength of the blue wing of the Fe K$\alpha$ line
profile to infer the Fe abundance and the viewer inclination angle,
respectively.

In this paper, we sought to understand and explain how the _physics_ of black
hole accretion is affected by the properties of the system and how these
physical changes translate into spectral properties—thus the discussion has
focused on general trends. However, the predictions we make are quantitative
and specific. Moving forward, the template spectra shown in this survey will
form the core of a simulation- and physics-based spectral modeling library to
be used for the analysis of real X-ray telescope data with xspec (Arnaud,
1996).

BEK was supported by the U.S. Department of Energy Advanced Simulation and
Computing Program’s Metropolis Fellowship, through the Los Alamos National
Laboratory, and used resources provided by the Los Alamos National Laboratory
Institutional Computing Program. Los Alamos National Laboratory is operated by
Triad National Security, LLC, for the National Nuclear Security Administration
of the U.S. Department of Energy (Contract No. 89233218CNA000001). JHK was
partially supported by NSF grants PHY-1707826, AST-1715032, and AST-2009260.

## References

  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Avara et al. (2016) Avara, M. J., McKinney, J. C., & Reynolds, C. S. 2016, MNRAS, 462, 636, doi: 10.1093/mnras/stw1643
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Chainakun et al. (2019) Chainakun, P., Watcharangkool, A., Young, A. J., & Hancock, S. 2019, MNRAS, 487, 667, doi: 10.1093/mnras/stz1319
  * Dauser et al. (2013) Dauser, T., Garcia, J., Wilms, J., et al. 2013, MNRAS, 430, 1694, doi: 10.1093/mnras/sts710
  * Fabian et al. (1989) Fabian, A. C., Rees, M. J., Stella, L., & White, N. E. 1989, MNRAS, 238, 729, doi: 10.1093/mnras/238.3.729
  * Gammie (1999) Gammie, C. F. 1999, ApJ, 522, L57, doi: 10.1086/312207
  * García et al. (2013) García, J., Dauser, T., Reynolds, C. S., et al. 2013, ApJ, 768, 146, doi: 10.1088/0004-637X/768/2/146
  * García & Kallman (2010) García, J., & Kallman, T. R. 2010, ApJ, 718, 695, doi: 10.1088/0004-637X/718/2/695
  * García et al. (2011) García, J., Kallman, T. R., & Mushotzky, R. F. 2011, ApJ, 731, 131, doi: 10.1088/0004-637X/731/2/131
  * García et al. (2014) García, J., Dauser, T., Lohfink, A., et al. 2014, ApJ, 782, 76, doi: 10.1088/0004-637X/782/2/76
  * Hawley & Balbus (1991) Hawley, J. F., & Balbus, S. A. 1991, ApJ, 376, 223, doi: 10.1086/170271
  * Jiang et al. (2019) Jiang, Y.-F., Blaes, O., Stone, J. M., & Davis, S. W. 2019, ApJ, 885, 144, doi: 10.3847/1538-4357/ab4a00
  * Kallman & Bautista (2001) Kallman, T., & Bautista, M. 2001, ApJS, 133, 221, doi: 10.1086/319184
  * Kara et al. (2019) Kara, E., Steiner, J. F., Fabian, A. C., et al. 2019, Nature, 565, 198, doi: 10.1038/s41586-018-0803-x
  * Kinch & Cleveland (in press) Kinch, B. E., & Cleveland, M. A. in press, ANS M&C 2021
  * Kinch et al. (2020) Kinch, B. E., Noble, S. C., Schnittman, J. D., & Krolik, J. H. 2020, ApJ, 904, 117, doi: 10.3847/1538-4357/abc176
  * Kinch et al. (2016) Kinch, B. E., Schnittman, J. D., Kallman, T. R., & Krolik, J. H. 2016, ApJ, 826, 52, doi: 10.3847/0004-637X/826/1/52
  * Kinch et al. (2019) —. 2019, ApJ, 873, 71, doi: 10.3847/1538-4357/ab05d5
  * Krolik (1999) Krolik, J. H. 1999, ApJ, 515, L73, doi: 10.1086/311979
  * Miller & Miller (2015) Miller, M. C., & Miller, J. M. 2015, Phys. Rep., 548, 1, doi: 10.1016/j.physrep.2014.09.003
  * Morales Teixeira et al. (2018) Morales Teixeira, D., Avara, M. J., & McKinney, J. C. 2018, MNRAS, 480, 3547, doi: 10.1093/mnras/sty2044
  * Narayan et al. (2016) Narayan, R., Zhu, Y., Psaltis, D., & Saḑowski, A. 2016, MNRAS, 457, 608, doi: 10.1093/mnras/stv2979
  * Noble et al. (2009) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, ApJ, 692, 411, doi: 10.1088/0004-637X/692/1/411
  * Noble et al. (2010) —. 2010, ApJ, 711, 959, doi: 10.1088/0004-637X/711/2/959
  * Noble et al. (2011) Noble, S. C., Krolik, J. H., Schnittman, J. D., & Hawley, J. F. 2011, ApJ, 743, 115, doi: 10.1088/0004-637X/743/2/115
  * Novikov & Thorne (1973) Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus), ed. C. Dewitt & B. S. Dewitt, 343–450
  * Ross & Fabian (2005) Ross, R. R., & Fabian, A. C. 2005, MNRAS, 358, 211, doi: 10.1111/j.1365-2966.2005.08797.x
  * Schnittman & Krolik (2013) Schnittman, J. D., & Krolik, J. H. 2013, ApJ, 777, 11, doi: 10.1088/0004-637X/777/1/11
  * Schnittman et al. (2013) Schnittman, J. D., Krolik, J. H., & Noble, S. C. 2013, ApJ, 769, 156, doi: 10.1088/0004-637X/769/2/156
  * Schnittman et al. (2016) —. 2016, ApJ, 819, 48, doi: 10.3847/0004-637X/819/1/48
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Shiokawa et al. (2018) Shiokawa, H., Cheng, R. M., Noble, S. C., & Krolik, J. H. 2018, ApJ, 861, 15, doi: 10.3847/1538-4357/aac2dd
  * Thorne (1974) Thorne, K. S. 1974, ApJ, 191, 507, doi: 10.1086/152991
  * Wilkins et al. (2016) Wilkins, D. R., Cackett, E. M., Fabian, A. C., & Reynolds, C. S. 2016, MNRAS, 458, 200, doi: 10.1093/mnras/stw276
  * Wilkins & Fabian (2012) Wilkins, D. R., & Fabian, A. C. 2012, MNRAS, 424, 1284, doi: 10.1111/j.1365-2966.2012.21308.x
  * Wong et al. (2021) Wong, G. N., Ryan, B. R., & Gammie, C. F. 2021, ApJ, 907, 73, doi: 10.3847/1538-4357/abd0f9
  * Zdziarski et al. (2021) Zdziarski, A. A., Dziełak, M. A., Marco, B. D., Szanecki, M., & Niedźwiecki, A. 2021, The Astrophysical Journal Letters, 909, L9, doi: 10.3847/2041-8213/abe7ef
  * Zoghbi et al. (2020) Zoghbi, A., Kalli, S., Miller, J. M., & Mizumoto, M. 2020, ApJ, 893, 97, doi: 10.3847/1538-4357/ab7dc8
  * Zoghbi et al. (2021) Zoghbi, A., Miller, J. M., & Cackett, E. 2021, arXiv e-prints, arXiv:2103.04994.  https://arxiv.org/abs/2103.04994