Mon. Not. R. Astron. Soc. 408, 752–782 (2010)

doi:10.1111/j.1365-2966.2010.17170.x

Simulations of magnetized discs around black holes: effects of black hole spin, disc thickness and magnetic ﬁeld geometry
Robert F. Penna,1 Jonathan C. McKinney,2 † Ramesh Narayan,1 Alexander Tchekhovskoy,1 Rebecca Shafee3 and Jeffrey E. McClintock1
1Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA 2Department of Physics and Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Stanford, CA 94305-4060, USA 3Center for Brain Science, Harvard University, 52 Oxford Street, Cambridge, MA 02138, USA
Accepted 2010 June 9. Received 2010 May 18; in original form 2010 March 4

ABSTRACT The standard general relativistic model of a razor-thin accretion disc around a black hole, developed by Novikov & Thorne (NT) in 1973, assumes the shear stress vanishes at the radius of the innermost stable circular orbit (ISCO) and that, outside the ISCO, the shear stress is produced by an effective turbulent viscosity. However, astrophysical accretion discs are not razor thin; it is uncertain whether the shear stress necessarily vanishes at the ISCO, and the magnetic ﬁeld, which is thought to drive turbulence in discs, may contain large-scale structures that do not behave like a simple local scalar viscosity. We describe 3D general relativistic magnetohydrodynamic simulations of accretion discs around black holes with a range of spin parameters, and we use the simulations to assess the validity of the NT model. Our ﬁducial initial magnetic ﬁeld consists of multiple (alternating polarity) poloidal ﬁeld loops whose shape is roughly isotropic in the disc in order to match the isotropic turbulence expected in the poloidal plane. For a thin disc with an aspect ratio |h/r| ∼ 0.07 around a non-spinning black hole, we ﬁnd a decrease in the accreted speciﬁc angular momentum of 2.9 per cent relative to the NT model and an excess luminosity from inside the ISCO of 3.5 per cent. The deviations in the case of spinning black holes are also of the same order. In addition, the deviations decrease with decreasing |h/r|. We therefore conclude that magnetized thin accretion discs in X-ray binaries in the thermal/high-soft spectral state ought to be well described by the NT model, especially at luminosities below 30 per cent of Eddington where we expect a very small disc thickness |h/r| 0.05. We use our results to determine the spin equilibrium of black hole accretion discs with a range of thicknesses and to determine how electromagnetic stresses within the ISCO depend upon black hole spin and disc thickness. We ﬁnd that the electromagnetic stress and the luminosity inside the ISCO depend on the assumed initial magnetic ﬁeld geometry. We consider a second geometry with ﬁeld lines following density contours, which for thin discs leads to highly radially elongated magnetic ﬁeld lines. This gives roughly twice larger deviations from NT for both the accreted speciﬁc angular momentum and the luminosity inside the ISCO. Lastly, we ﬁnd that the disc’s corona (including any wind or jet) introduces deviations from NT in the speciﬁc angular momentum that are comparable to those contributed by the disc component, while the excess luminosity of bound gas from within the ISCO is dominated by only the disc component. Based on these indications, we suggest that differences in results between our work and other similar work are due to differences in the assumed initial magnetic ﬁeld geometry as well as the inclusion of disc

E-mail: rpenna@cfa.harvard.edu (RFP); jmckinne@stanford.edu (JCM); rnarayan@cfa.harvard.edu (RN); atchekho@cfa.harvard.edu (AT); shafee@fas.harvard.edu (RS); jmcclintock@cfa.harvard.edu (JEM) † Chandra Fellow
C 2010 The Authors. Journal compilation C 2010 RAS

Simulations of magnetized discs 753
gas versus all the gas when comparing the speciﬁc angular momentum from the simulations with the NT model. Key words: accretion, accretion discs – black hole physics – gravitation – hydrodynamics – MHD – methods: numerical.

1 INTRODUCTION
Accreting black holes (BHs) are among the most powerful astrophysical objects in the Universe. Although they have been the target of intense research for a long time, many aspects of BH accretion theory remain uncertain to this day. Pioneering work by Bardeen (1970), Shakura & Sunyaev (1973), Novikov & Thorne (1973) (hereafter NT), Page & Thorne (1974) and others indicated that BH accretion through a razor-thin disc can be highly efﬁcient, with up to 42 per cent of the accreted rest-mass-energy being converted into radiation. These authors postulated the existence of a turbulent viscosity in the disc, parametrized via the famous α-prescription (Shakura & Sunyaev 1973). This viscosity causes outwards transport of angular momentum; in the process, it dissipates energy and produces the radiation. The authors also assumed that, within the innermost stable circular orbit (ISCO) of the BH, the viscous torque vanishes and material plunges into the BH with constant energy and angular momentum ﬂux per unit rest-mass ﬂux. This is the so-called ‘zero-torque’ boundary condition.
Modern viscous hydrodynamical calculations of discs with arbitrary thicknesses suggest that the zero-torque condition is a good approximation when the height (h) to radius (r) ratio of the accreting gas is small: |h/r| 0.1 (Paczyn´ski 2000; Afshordi & Paczyn´ski 2003; Shafee, Narayan & McClintock 2008b; Sa¸dowski 2009; Abramowicz et al. 2010). Radiatively efﬁcient discs in active galactic nuclei (AGN) and X-ray binaries are expected to have disc thickness |h/r| < 0.1 whenever the luminosity is limited to less than about 30 per cent of the Eddington luminosity (McClintock et al. 2006). The above hydrodynamical studies thus suggest that systems in this limit should be described well by the standard relativistic thin disc theory as originally developed by NT.
In parallel with the above work, it has for long been recognized that the magnetic ﬁeld could be a complicating factor that may signiﬁcantly modify accretion dynamics near and inside the ISCO (Thorne 1974). This issue has become increasingly important with the realization that angular momentum transport in discs is entirely due to turbulence generated via the magnetorotational instability (MRI) (Balbus & Hawley 1991, 1998). However, the magnetic ﬁeld does not necessarily behave like a local viscous hydrodynamical stress. Near the BH, the magnetic ﬁeld may have large-scale structures (MacDonald 1984), which can induce stresses across the ISCO (Krolik 1999; Gammie 1999; Agol & Krolik 2000), leading to changes in, e.g., the density, velocity, and amount of dissipation and emission. Unlike turbulence, the magnetic ﬁeld can transport angular momentum without dissipation (e.g. Li 2002), or it can dissipate in current sheets without transporting angular momentum. In Agol & Krolik (2000), the additional electromagnetic stresses are treated simply as a freely tunable model parameter on top of an otherwise hydrodynamical model. A more complete magnetohydrodynamical (MHD) model of a magnetized thin disc has been developed by Gammie (1999). In this model, the controlling free parameter is the speciﬁc magnetic ﬂux, i.e. magnetic ﬂux per unit rest-mass ﬂux. Larger values of this parameter lead to larger deviations from NT due to electromagnetic stresses, but the exact value of the parameter

for a given accretion disc is unknown. For instance, it is entirely possible that electromagnetic stresses become negligible in the limit when the disc thickness |h/r| → 0. The value of the speciﬁc magnetic ﬂux is determined by the non-linear turbulent saturation of the magnetic ﬁeld, so accretion disc simulations are the best way to establish its magnitude.
The coupling via the magnetic ﬁeld between a spinning BH and an external disc, or between the hole and the corona, wind and jet (hereafter corona-wind-jet), might also play an important role in modifying the accretion ﬂow near the BH. The wind or jet (hereafter wind-jet) can transport angular momentum and energy away from the accretion disc and BH (Blandford 1976; Blandford & Znajek 1977; McKinney & Gammie 2004; McKinney 2006b; McKinney & Narayan 2007b; Komissarov & McKinney 2007). The wind-jet power depends upon factors such as the BH spin (McKinney 2005; Hawley & Krolik 2006; Komissarov et al. 2007), disc thickness (Meier 2001; Tchekhovskoy, McKinney & Narayan 2008, 2009a; Tchekhovskoy, Narayan & McKinney 2009b, 2010), and the strength and large-scale behaviour of the magnetic ﬁeld (McKinney & Gammie 2004; Beckwith, Hawley & Krolik 2008a; McKinney & Blandford 2009), and these can affect the angular momentum transport through an accretion disc. In this context, we note that understanding how such factors affect disc structure may be key in interpreting the distinct states and complex behaviours observed for BH X-ray binaries (Fender, Belloni & Gallo 2004; Remillard & McClintock 2006). These factors also affect the BH spin history (Gammie, Shapiro & McKinney 2004), and so must be taken into account when considering the effect of accretion on the cosmological evolution of BH spin (Hughes & Blandford 2003; Gammie et al. 2004; Berti & Volonteri 2008).
Global simulations of accretion discs using general relativistic magnetohydrodynamical (GRMHD) codes (e.g. Gammie, McKinney & To´th 2003; De Villiers, Hawley & Krolik 2003) currently provide the most complete understanding of how turbulent magnetized accretion ﬂows around BHs work. Most simulations have studied thick (|h/r| 0.15) discs without radiative cooling. Such global simulations of the inner accretion ﬂow have shown that ﬂuid crosses the ISCO without any clear evidence that the torque vanishes at the ISCO, i.e. there is no apparent ‘stress edge’ (McKinney & Gammie 2004; Krolik, Hawley & Hirose 2005; Beckwith, Hawley & Krolik 2008b). Similar results were previously found with a pseudo-Newtonian potential for the BH (Krolik & Hawley 2002). In these studies, a plot of the radial proﬁle of the normalized stress within the ISCO appears to indicate a signiﬁcant deviation from the NT thin disc theory (Krolik et al. 2005; Beckwith et al. 2008b), and it was thus expected that much thinner discs might also deviate signiﬁcantly from NT. A complicating factor in drawing ﬁrm conclusions from such studies is that the assumed initial global magnetic ﬁeld geometry and strength can signiﬁcantly change the magnitude of electromagnetic stresses and associated angular momentum transport inside the ISCO (McKinney & Gammie 2004; Beckwith et al. 2008a).
The implications of the above studies for truly thin discs (|h/r| 0.1) remain uncertain. Thin discs are difﬁcult to resolve numerically,

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

754 R. F. Penna et al.
and simulations have been attempted only recently. Simulations of thin discs using a pseudo-Newtonian potential for the BH reveal good agreement with standard thin disc theory (Reynolds & Fabian 2008). The ﬁrst simulation of a thin (|h/r| ≈ 0.05) disc using a full GRMHD model was by Shafee et al. (2008a, hereafter S08). They considered a non-spinning (a/M = 0) BH and an initial ﬁeld geometry consisting of multiple opposite-polarity bundles of poloidal loops within the disc. They found that, although the stress proﬁle appears to indicate signiﬁcant torques inside the ISCO, the actual angular momentum ﬂux per unit rest-mass ﬂux through the disc component deviates from the NT prediction by only 2 per cent, corresponding to an estimated deviation in the luminosity of only about 4 per cent. The study by S08 was complemented by Noble, Krolik & Hawley (2009, hereafter N09), who considered a thin (|h/r| ≈ 0.1) disc around an a/M = 0.9 BH and an initial ﬁeld geometry consisting of a single highly elongated poloidal loop bundle whose ﬁeld lines follow the density contours of the thin disc. They found 6 per cent more luminosity than predicted by NT. More recently, Noble, Krolik & Hawley (2010, hereafter N10), considered a thin (|h/r| ≈ 0.07) disc around a non-spinning (a/M = 0) BH and reported up to 10 per cent deviations from NT in the speciﬁc angular momentum accreted through the accretion ﬂow.
In this paper, we extend the work of S08 by considering a range of BH spins, disc thicknesses, ﬁeld geometries, box sizes, numerical resolutions, etc. Our primary result is that we conﬁrm S08, viz., geometrically thin discs are well-described by the NT model. We show that there are important differences between the dynamics of the gas in the disc and in the corona-wind-jet. In addition, we ﬁnd that the torque and luminosity within the ISCO can be signiﬁcantly affected by the geometry and strength of the initial magnetic ﬁeld, a result that should be considered when comparing simulation results to thin disc theory. In this context, we discuss likely reasons for the apparently different conclusions reached by N09 and N10.
The equations we solve are given in Section 2, diagnostics are described in Section 3 and our numerical setup is described in Section 4. Results for our ﬁducial thin disc model for a non-rotating BH are given in Section 5, and we summarize convergence studies in Section 6. Results for a variety of BH spins and disc thicknesses are presented in Section 7 and for thin discs with different magnetic ﬁeld geometries and strengths in Section 8. We compare our results with previous studies in Section 9, discuss the implications of our results in Section 10 and conclude with a summary of the salient points in Section 11.

2 GOVERNING EQUATIONS

The system of interest to us is a magnetized accretion disc around a rotating BH. We write the BH Kerr metric in Kerr–Schild (KS; horizon-penetrating) coordinates (Font, Iba´n˜ez & Papadopoulos 1998; Papadopoulos & Font 1998), which can be mapped to Boyer– Lindquist (BL) coordinates or an orthonormal basis in any frame (McKinney & Gammie 2004). We work with Heaviside–Lorentz units, set the speed of light and gravitational constant to unity (c = G = 1), and let M be the BH mass. We solve the GRMHD equations of motion for rotating BHs (Gammie et al. 2003) with an additional cooling term designed to keep the simulated accretion disc at a desired thickness.
Mass conservation gives

∇μ(ρ0uμ) = 0,

(1)

where ρ0 is the rest-mass density, corresponding to the mass density in the ﬂuid frame, and uμ is the contravariant four-velocity. Note

that we write the orthonormal three-velocity as vi (the covariant three-velocity is never used below). Energy-momentum conservation gives

∇μTνμ = Sν ,

(2)

where

the

stress

energy

tensor

T

μ ν

includes

both

matter

and

electro-

magnetic terms,

Tνμ = (ρ0 + ug + pg + b2)uμuν + (pg + b2/2)δνμ − bμbν ,

(3)

where ug is the internal energy density and pg = ( − 1)ug is the ideal gas pressure with = 4/3.1 The contravariant ﬂuid-frame magnetic

four-ﬁeld is given by bμ, and is related to the lab-frame three-ﬁeld

via

bμ

=

Bν hμν /ut

where

hμν

=

uμuν

+

δ

μ ν

is

a

projection

tensor,

and δμν is the Kronecker delta function. We write the orthonormal

three-ﬁeld as Bi (the covariant three-ﬁeld is never used below). The

magnetic energy density (ub) and magnetic pressure (pb) are then

given by ub = pb = bμbμ/2 = b2/2. Note that the angular velocity

of the gas is = uφ/ut . Equation (2) has a source term

Sν =

dU dτ

uν ,

(4)

which is a radiation four-force corresponding to a simple isotropic comoving cooling term given by dU/dτ . We ignore radiative transport effects such as heat currents, viscous stresses or other effects that would enter as additional momentum sources in the comoving frame. In order to keep the accretion disc thin, we employ the same ad hoc cooling function as in S08:

dU dτ

=

−ug

log

(K /Kc ) τcool

S [θ ],

(5)

where τ is the ﬂuid proper time, K = pg/ρ0 is the entropy constant, Kc = 0.000 69 is set to be the same entropy constant as the torus
atmosphere and is the entropy constant we cool the disc towards, and K0 Kc is the entropy constant of the initial torus2. The gas cooling time is set to τ cool = 2π / K, where K = (1/M)/[(a/M) + (R/M)3/2] is the Keplerian angular frequency and R = r sin θ is the

cylindrical radius (We consider variations in the cooling time-scale

in section 5.7.). We use a shaping function given by the quantity S[θ ] = exp[−(θ − π /2)2/(2(θ nocool)2)], where we set θ nocool = {0.1, 0.3, 0.45, 0.45} for our sequence of models with target thickness of |h/r| = {0.07, 0.1, 0.2, 0.3}, although we note that the thickest model with target |h/r| = 0.3 has no cooling turned on.

The shaping function S[θ ] is used to avoid cooling in the low-

density funnel-jet region where the thermodynamics is not accurately evolved and where the gas is mostly unbound (see Fig. 13 in Section 5.7). In addition, we set the cooling function dU/dτ =

0 if (1) the time-step, dt, obeys dt > τ cool, which ensures that the cooling does not create negative entropy gas; or (2) log(K/Kc) <

1Models with = 5/3 show some minor differences compared to models with = 4/3 (McKinney & Gammie 2004; Mignone & McKinney 2007). 2 We intended to have a constant K everywhere at t = 0, but a normalization issue led to Kc K0. Because of this condition, the disc cools towards a slightly thinner equilibrium at the start of the simulation, after which the
cooling proceeds as originally desired by cooling towards the ﬁducial value K = Kc. Our models with |h/r| ≈ 0.07 are least affected by this issue. Also, since we do not make use of the cooling-based luminosity near t = 0, this issue does not affect any of our results. We conﬁrmed that this change leads
to no signiﬁcant issues for either the magnitude or scaling of quantities with thickness by repeating some simulations with the intended Kc = K0. The otherwise similar simulations have thicker discs as expected (very minor
change for our thin disc model as expected), and we ﬁnd consistent results
for a given measured thickness in the saturated state.

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

Simulations of magnetized discs 755

0, which ensures the gas is only cooled, never heated. Photon cap-

ture by the BH is not included, so the luminosity based upon this

cooling function is an upper limit for radiation from the disc. The

above cooling function drives the speciﬁc entropy of the gas towards

the reference speciﬁc entropy Kc. Since speciﬁc entropy always in-

creases due to dissipation, this cooling function explicitly tracks

dissipation. Hence, the luminosity generated from the cooling func-

tion should not be considered as the true luminosity, but instead

should be considered as representing the emission rate in the limit

that all dissipated energy is lost as radiation. Any other arbitrary

cooling function that does not track dissipation would require full

radiative transfer to obtain the true luminosity.

Magnetic ﬂux conservation is given by the induction equation

∂t

√ ( −gB

i

)

=

−∂j

√ [ −g(B

i

vj

−

Bj vi )],

(6)

where vi = ui/ut , and g = Det (gμν) is the determinant of the metric. In steady state, the cooling is balanced by heating from shocks, gridscale reconnection and grid-scale viscosity. No explicit resistivity or viscosity is included.

3 DIAGNOSTICS

In this section, we describe several important diagnostics that we have found useful in this study. First, since we regulate the disc height via an ad hoc cooling function, we check the scaleheight of the simulated disc as a function of time and radius to ensure that our cooling function operates properly. Secondly, the equations we solve consist of seven time-dependent ideal MHD equations, corresponding to four relevant conserved quantities.3 Using these quantities we construct three dimensionless ﬂux ratios corresponding to the accreted speciﬁc energy, speciﬁc angular momentum and speciﬁc magnetic ﬂux. Thirdly, we check what the duration of the simulations should be in order to reach a quasi-steady-state (‘inﬂow equilibrium’) at any given radius. Finally, we describe how we compute the luminosity.
When the speciﬁc ﬂuxes are computed as a spatial or temporal average/integral, we always take the ratio of averages/integrals of ﬂuxes (i.e. dxF1/ dxF2) rather than the average/integral of the ratio of ﬂuxes (i.e. dx(F1/F2)). The former is more appropriate for capturing the mean behaviour, while the latter can be more appropriate when investigating ﬂuxes with signiﬁcant phase-shifted correlations between each other. As relevant for this study, the accretion disc has signiﬁcant vertical stratiﬁcation and the local value of the ratio of ﬂuxes can vary considerably without any effect on the bulk accretion ﬂow. Similarly, potentially one ﬂux can (e.g.) nearly vanish over short periods, while the other ﬂux does not, which leads to unbounded values for the ratio of ﬂuxes. However, the time-averaged behaviour of the ﬂow is not greatly affected by such short periods of time. These cases demonstrate why the ratio of averages/integrals is always performed for both spatial and temporal averages/integrals.
When comparing the ﬂux ratios or luminosities from the simulations against the NT model, we evaluate the per cent relative difference D[f ] between a quantity f and its NT value as follows:

D[f ]

≡

f 100

−

f [NT] .

f [NT]

(7)

3 The energy momentum of the ﬂuid is not strictly conserved because of radiative cooling; however, the ﬂuid component of the energy-momentum equations still proves to be useful. Only energy conservation of the ﬂuid is strongly affected for our types of models.

For a density-weighted time-averaged value of f , we compute

f ρ0 ≡

f ρ0(r, θ, φ)dAθφdt , ρ0(r, θ, φ)dAθφdt

(8)

√ where dAθφ ≡ −gdθ dφ is an area element in the θ − φ plane, and

the integral over dt is a time average over the duration of interest,

which corresponds to the period when the disc is in steady state.

For a surface-averaged value of f , we compute

f≡

f dAθφ . dAθ φ

(9)

3.1 Disc thickness measurement

We deﬁne the dimensionless disc thickness per unit radius, |h/r|, as the density-weighted mean angular deviation of the gas from the midplane,

hπ

≡ θ−

.

r 2 ρ0

(10)

(This quantity was called θ abs in S08.) Notice that we assume the accretion disc plane is on the equator (i.e. we assume θ ρ0 = π/2). As deﬁned above, |h/r| is a function of r. When we wish to characterize the disc by means of a single estimate of its thickness, we use the value of |h/r| at r = 2rISCO, where rISCO is the ISCO radius (rISCO = 6M for a non-spinning BH and rISCO = M for a maximally spinning BH; Shapiro & Teukolsky 1983). As we show in
Section 5.4, this choice is quite reasonable. An alternative thickness
measure, given by the root-mean-square thickness (h/r)rms, allows us to estimate how accurate we can be about our deﬁnition of
thickness. This quantity is deﬁned by

h ≡
r rms

π 2 1/2 θ− .
2 ρ0

(11)

The range of θ for the disc thickness integrals in the above equations

is from 0 to π .

3.2 Fluxes of mass, energy and angular momentum
The mass, energy and angular momentum conservation equations give the following ﬂuxes:

M˙ (r, t) = −

ρ0ur dAθφ ,

θφ

(12)

e

≡

E˙ (r, t) M˙ (r, t)

=

θ

φ
M˙

Ttr dAθ (r, t)

φ

,

(13)

j≡

J˙(r, t) M˙ (r, t)

=−

θ

φ Tφr dAθφ M˙ (r, t)

.

(14)

The above relations deﬁne the rest-mass accretion rate (sometimes just referred to as the mass accretion rate), M˙ ; the accreted energy ﬂux per unit rest-mass ﬂux, or speciﬁc energy, e; and the accreted angular momentum ﬂux per unit rest-mass ﬂux, or speciﬁc angular momentum, j . Positive values of these quantities correspond to an inward ﬂux through the BH horizon.
The BH spin evolves due to the accretion of mass, energy and angular momentum, which can be described by the dimensionless spin-up parameter s,

s

≡

d(a /M ) dt

M M˙

=

j

− 2 a e, M

(15)

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

756 R. F. Penna et al.

where the angular integrals used to compute j and e include all θ and φ angles (Gammie et al. 2004). For s = 0 the BH is in so-called ‘spin equilibrium’, corresponding to when the dimensionless BH spin, a/M, does not change in time.
The ‘nominal’ efﬁciency, corresponding to the total loss of speciﬁc energy from the ﬂuid, is obtained by removing the rest-mass term from the accreted speciﬁc energy:

e˜ ≡ 1 − e.

(16)

The time-averaged value of e˜ at the horizon (r = rH) gives the total nominal efﬁciency: e˜(rH) , which is an upper bound on the total photon radiative efﬁciency.
The range of θ over which the ﬂux density integrals in the above equations are computed depends on the situation. In S08, we limited the θ range to δθ = ±0.2 corresponding to 2–3 density scaleheights in order to focus on the disc and to avoid including the disc wind or BH jet. In this paper, we are interested in studying how the contributions to the ﬂuxes vary as a function of height above the equatorial plane. Our expectation is that the disc and corona-windjet contribute differently to these ﬂuxes. Thus, we consider different ranges of θ in the integrals, e.g. from (π /2) − 2|h/r| to (π /2) + 2|h/r|, (π /2) − 4|h/r| to (π /2) + 4|h/r|, or 0 to π . The ﬁrst and third are most often used in later sections.

3.3 Splitting angular momentum ﬂux into ingoing and outgoing components

For a more detailed comparison of the simulation results with the NT model, we decompose the ﬂux of angular momentum into an ingoing (‘in’) term which is related to the advection of mass-energy into the BH and an outgoing (‘out’) term which is related to the forces and stresses that transport angular momentum radially outward. These ingoing and outgoing components of the speciﬁc angular momentum are deﬁned by

jin(r, t) ≡

(ρ0 + ug + b2/2)ur −ρ0ur

uφ ,

(17)

jout(r, t) ≡ j − jin(r, t).

(18)

By this deﬁnition, the ‘in’ quantities correspond to inward transport of the comoving mass-energy density of the disc, uμuνTμν = ρ0 + ug + b2/2. Note that ‘in’ quantities are products of the mean velocity ﬁelds ur and uμ and not the combination ur uμ ; the latter includes a contribution from correlated ﬂuctuations in ur and uμ, which corresponds to the Reynolds stress. The residual of the total
ﬂux minus the ‘in’ ﬂux gives the outward, mechanical transport
by Reynolds stresses and electromagnetic stresses. One could also
consider a similar splitting for the speciﬁc energy. The above de-
composition most closely matches our expectation that the inward ﬂux should agree with the NT result as |h/r| → 0. Note, however, that our conclusions in this paper do not require any particular de-
composition. This decomposition is different from S08 and N10 where the entire magnetic term (b2ur uφ − br bφ) is designated as the ‘out’ term. Their choice overestimates the effect of electromag-
netic stresses, since some magnetic energy is simply advected into
the BH. Also, the splitting used in S08 gives non-monotonic j in versus radius for some BH spins, while the splitting we use gives
monotonic values for all BH spins.

3.4 The magnetic ﬂux
The no-monopoles constraint implies that the total magnetic ﬂux ( = S B · d A) vanishes through any closed surface or any open

surface penetrating a bounded ﬂux bundle. The magnetic ﬂux conservation equations require that magnetic ﬂux can only be transported to the BH or through the equatorial plane by advection. The absolute magnetic ﬂux ( S |B · d A|) has no such restrictions and can grow arbitrarily due to the MRI. However, the absolute magnetic ﬂux can saturate when the electromagnetic ﬁeld comes into force balance with the matter. We are interested in such a saturated state of the magnetic ﬁeld within the accretion ﬂow and threading the BH.
We consider the absolute magnetic ﬂux penetrating a spherical surface and an equatorial surface given, respectively, by

r(r, θ, t) =

|Br |dAθ φ ,

θφ

(19)

r =r

θ (r, θ, t) =

|Bθ |dAr φ .

r =rH φ

(20)

Nominally, r has an integration range of θ = 0 to θ = θ when measured on the BH horizon, while when computing quantities

around the equatorial plane θ has the range θ ± θ . One useful

normalization of the magnetic ﬂuxes is by the total ﬂux through one

hemisphere of the BH plus through the equator

tot(r, t) ≡ r(r = rH, θ = 0 . . . π/2, t) + θ (r, θ = π/2, t),

(21)

which gives the normalized absolute radial magnetic ﬂux

˜ r(r, θ, t) ≡

r(r, θ, t)

,

tot(r = Rout, t = 0)

(22)

where Rout is the outer radius of the computational box. The normalized absolute magnetic ﬂux measures the absolute magnetic ﬂux

on the BH horizon or radially through the equatorial disc per unit

absolute ﬂux that is initially available.

The Gammie (1999) model of a magnetized thin accretion ﬂow

suggests another useful normalization of the magnetic ﬂux is by the

absolute mass accretion rate

M˙ G(r, t) ≡

ρ0|ur |dAθφ ,

θφ

(23)

which gives the normalized speciﬁc absolute magnetic ﬂuxes

(r, t)

=

M˙

r(r, t) G(r, t)

,

(24)

√ ϒ(r, t) ≡ 2

(r, t)

M

M˙ G(r = rH, t) , SAH

(25)

where SA = (1/r2) θ φ dAθφ is the local solid angle, SAH = SA(r = rH) is the local solid angle on the horizon, (r, t) is the radial magnetic ﬂux per unit rest-mass ﬂux (usually speciﬁc magnetic ﬂux) and ϒ(r, t)c3/2/G is a particular dimensionless normalization of the speciﬁc magnetic ﬂux that appears in the MHD accretion model developed by Gammie (1999). Since the units used for the magnetic ﬁeld are arbitrary, any constant factor can be applied to and one would still identify the quantity as the speciﬁc magnetic ﬂux. So to simplify the discussion we henceforth call ϒ the speciﬁc magnetic ﬂux. To obtain equation (25), all involved integrals should have a common θ range around the equator. These quantities all have absolute magnitudes because a sign change does not change the physical effect. The quantities j, e, e˜, and ϒ are each conserved along poloidal ﬁeld-ﬂow lines for stationary ideal MHD solutions (Bekenstein & Oron 1978; Takahashi et al. 1990).

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

Gammie’s (1999) model of a magnetized accretion ﬂow within the ISCO assumes: (1) a thin equatorial ﬂow; (2) a radial poloidal ﬁeld geometry (i.e. |Bθ | |Br|); (3) a boundary condition at the ISCO corresponding to zero radial velocity and (4) no thermal contribution. The model reduces to the NT solution within the ISCO for ϒ → 0, and deviations from NT’s solution are typically small (less than 12 per cent for j across all BH spins; see Appendix A)√for ϒ 1. We have deﬁned the quantity ϒ in equation (24) with the 2 factor, the square root of the total mass accretion rate through the horizon per unit solid angle and Heaviside–Lorentz units for Br so that the numerical value of ϒ at the horizon is identically equal to the numerical value of the free parameter in Gammie (1999), i.e. their Fθφ normalized by FM = −1. As shown in that paper, ϒ directly controls deviations of the speciﬁc angular momentum and speciﬁc energy away from the non-magnetized thin disc theory values of the NT model. Even for discs of ﬁnite thickness, the parameter shows how electromagnetic stresses control deviations between the horizon and the ISCO. Note that the ﬂow can deviate from NT at the ISCO simply due to ﬁnite thermal pressure (McKinney & Gammie 2004). In Appendix (Table A1), we list numerical values of j and e˜ for Gammie’s (1999) model, and show how these quantities deviate from NT for a given BH spin and ϒ.
We ﬁnd ϒ to be more useful as a measure of the importance of the magnetic ﬁeld within the ISCO than our previous measurement in S08 of the α-viscosity parameter,

T φˆ rˆ α= ,
pg + pb

(26)

where T φˆ rˆ = eφˆ μerˆ ν T μν is the orthonormal stress-energy tensor components in the comoving frame, and eμνˆ is the orthonormal contravariant tetrad system in the local ﬂuid-frame. This is related
to the normalized stress by

W M˙

=

M˙

T φˆ rˆ dAθ φ dLφ

φ

,

(27)

where dAθφ = eθˆ μeφˆ ν dθ μdφν is the comoving area element, dLφ = eφˆ νdφν evaluated at θ = π /2 is the comoving φ length element, θ μ = {0, 0, 1, 0} and φν = {0, 0, 0, 1}. This form for W is a simple generalization of equation (5.6.1b) in NT73, and note that the NT solution for W /M˙ is given by equation (5.6.14a) in NT73. In S08, W was integrated over ﬂuid satisfying −ut(ρ0 + ug + pg + b2)/ρ0 < 1 (i.e. only approximately gravitationally bound ﬂuid and no wind-jet). We use the same deﬁnition of bound in this paper. As shown in S08, a plot of the radial proﬁle of W /M˙ or α within the
ISCO does not necessarily quantify how much the magnetic ﬁeld
affects the accretion ﬂow properties, since even apparently large
values of this quantity within the ISCO do not cause a signiﬁcant
deviation from NT in the speciﬁc angular momentum accreted. On
the other hand, the Gammie (1999) parameter ϒ does directly re-
late to the electromagnetic stresses within the ISCO and is an ideal
MHD invariant (so constant versus radius) for a stationary radial
ﬂow. One expects that appropriately time-averaged simulation data
could be framed in the context of this stationary model in order to
measure the effects of electromagnetic stresses.

3.5 Inﬂow equilibrium
When the accretion ﬂow has achieved steady state inside a given radius, the quantity M˙ (r, t) will (apart from normal ﬂuctuations due to turbulence) be independent of time, and if it is integrated over

Simulations of magnetized discs 757

all θ angles will be constant within the given radius.4 The energy and angular momentum ﬂuxes have a non-conservative contribution due to the cooling function and therefore are not strictly constant. However, the cooling is generally a minor contribution (especially in the case of the angular momentum ﬂux), and so we may still measure the non-radiative terms to verify inﬂow equilibrium.
The radius out to which inﬂow equilibrium can be achieved in a given time can be estimated by calculating the mean radial velocity vr and then deriving from it a viscous time-scale −r/vr. From standard accretion disc theory and using the deﬁnition of α given in equation (27), the mean radial velocity is given by

vr ∼ −α

h2 r vK,

(28)

where vK ≈ (r/M)−1/2 is the Keplerian speed at radius r and α is the standard viscosity parameter given by equation (26) (Frank, King & Raine 1992). Although the viscous time-scale is the nominal time needed to achieve steady state, in practice it takes several viscous times before the ﬂow really settles down, e.g. see the calculations reported in Shapiro (2010). In the present paper, we assume that inﬂow equilibrium is reached after two viscous times, and hence we estimate the inﬂow equilibrium time, tie, to be

r r 3/2

tie

∼ −2 vr

∼2

M

1 α|h/r |2

r 3/2

∼ 5000

,

M

(29)

where, in the right-most relation, we have taken a typical value of α ∼ 0.1 for the gas in the disc proper (i.e. outside the ISCO) and we have set |h/r| ≈ 0.064, as appropriate for our thinnest disc models.
A simulation must run until t ∼ tie before we can expect inﬂow equilibrium at radius r. According to the above Newtonian estimate, a thin disc simulation with |h/r| ∼ 0.064 that has been run for a time of 30000M will achieve steady state out to a radius of only ∼3 M. However, this estimate is inaccurate since it does not allow for the boundary condition on the ﬂow at the ISCO. Both the boundary condition and the effects of GR are included in the formula for the radial velocity given in equation (5.9.8) of NT, which we present for completeness in Appendix B. That more accurate result, which is what we use for all our plots and numerical estimates, shows that our thin disc simulations should reach inﬂow equilibrium within r/M = 9, 7, 5.5, 5, respectively, for a/M = 0, 0.7, 0.9, 0.98. These estimates are roughly consistent with the radii out to which we have a constant j versus radius in the simulations discussed in Section 7.

3.6 Luminosity measurement

We measure the radiative luminosity of the accreting gas directly from the cooling function dU/dτ . At a given radius, r, in the steady region of the ﬂow, the luminosity per unit rest-mass accretion rate interior to that radius is given by

L(< r) M˙ (r, t)

=

1 M˙ (r, t)(tf

−

ti)

tf t =ti

r r =rH

π θ =0

φ

dU dτ

ut dVtr θφ ,

(30) where dVtr θφ = √−gdt dr dθ dφ and the 4D integral goes from
the initial time ti to the ﬁnal time tf over which the simulation

4 If we integrate over a restricted range of θ , then there could be additional mass ﬂow through the boundaries in the θ direction and M˙ (r, t) will no
longer be independent of r, though it would still be independent of t.

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

758 R. F. Penna et al.

results are time-averaged, from the radius rH of the horizon to the radius r of interest, and usually over all θ and φ. We ﬁnd it useful to compare the simulations with thin disc theory by computing the ratio of the luminosity emitted inside the ISCO (per unit rest-mass accretion rate) to the total radiative efﬁciency of the NT model:

L˜ in

≡

L(< rISCO) M˙ e˜[NT]

.

(31)

This ratio measures the excess radiative luminosity from inside the ISCO in the simulation relative to the total luminosity in the NT model (which predicts zero luminosity here). We also consider the excess luminosity over the entire inﬂow equilibrium region

L˜ eq

≡

L(r

<

req) − L(r < M˙ e˜[NT]

req)[NT] ,

(32)

which corresponds to the luminosity (per unit mass accretion rate)
inside the inﬂow equilibrium region (i.e. r < req, where req is the radius out to which inﬂow equilibrium has been established) sub-
tracted by the NT luminosity all divided by the total NT efﬁciency. Large per cent values of L˜ in or L˜ eq would indicate large per cent deviations from NT.

4 PHYSICAL MODELS AND NUMERICAL METHODS
This section describes our physical models and numerical methods. Table 1 provides a list of all our simulations and shows the physical and numerical parameters that we vary. Our primary models are labelled by names of the form AxHRy, where x is the value of the BH spin parameter and y is approximately equal to the disc thickness |h/r|. For instance, our ﬁducial model A0HR07 has a non-spinning BH (a/M = 0) and a geometrically thin disc with

|h/r| ∼ 0.07. We discuss this particular model in detail in Section 5. Table 1 also shows the time span (from Ti/M to Tf/M) used to perform the time-averaging, and the last column shows the actual value of |h/r| in the simulated model as measured during inﬂow equilibrium, e.g. |h/r| = 0.064 for model A0HR07.
4.1 Physical models
This study considers BH accretion disc systems with a range of BH spins: a/M = 0, 0.7, 0.9, 0.98 and a range of disc thicknesses: |h/r| = 0.07, 0.13, 0.2, 0.3. The initial mass distribution is an isentropic equilibrium torus (Chakrabarti 1985a,b; De Villiers et al. 2003). All models have an initial inner torus edge at rin = 20M, while the torus pressure maximum for each model is located between Rmax = 35M and Rmax = 65M. We chose this relatively large radius for the initial torus because S08 found that placing the torus at smaller radii caused the results to be sensitive to the initial mass distribution. We initialize the solution so that ρ0 = 1 is the maximum rest-mass density. In S08, we set q = 1.65 ( ∝ r−q in non-relativistic limit) and K = 0.00034 with = 4/3, while in this paper we adjust the initial angular momentum proﬁle such that the initial torus has the target value of |h/r| at the pressure maximum. For models with |h/r| = 0.07, 0.13, 0.2, 0.3, we ﬁx the speciﬁc entropy of the torus by setting, respectively, K = K0 ≡ {0.00034, 0.0035, 0.009, 0.009} in the initial polytropic equation of state given by p = K0ρ0 . The initial atmosphere surrounding the torus has the same polytropic equation of state with nearly the same entropy constant of K = 0.000 69, but with an initial restmass density of ρ0 = 10−6(r/M)−3/2, corresponding to a Bondi-like atmosphere.
Recent GRMHD simulations of thick discs indicate that the results for the disc (but not the wind-jet) are roughly independent of the initial ﬁeld geometry (McKinney & Narayan 2007a,b; Beckwith

Table 1. Simulation parameters.

Model name Ti/M– Tf /M

a M

Nr Nθ Nφ

Rin rH

Y

φ

Rmax M

Rout M

q Cooling

λﬁeld 2π r

h r

A0HR07 12500–27350

0 256 64

A7HR07 12500–20950 0.7 256 64

A9HR07 14000–23050 0.9 256 64

A98HR07 14000–19450 0.98 256 64

A0HR1 5050–14150

0 256 64

A7HR1 5050–12550 0.7 256 64

A9HR1 5050–10000 0.9 256 64

A98HR1 12000–13600 0.98 256 64

A0HR2 6000–13750

0 256 64

A7HR2 12000–17900 0.7 256 64

A9HR2 12000–15200 0.9 256 64

A98HR2

6100–7100 0.98 256 64

A0HR3

4700–7900

0 256 64

A7HR3 10000–11900 0.7 256 64

A9HR3

4700–7900 0.9 256 64

A98HR3

4700–7900 0.98 256 64

C0 6000–10000

0 512 128

C1 12500–18900

0 256 64

C2 12500–22500

0 256 64

C3 12500–19500

0 256 64

C4 12500–21700

0 256 64

C5 12500–20000

0 256 32

C6 12500–20800

0 256 128

A0HR07LOOP1 17000–22100

0 256 64

A0HR3LOOP1

3000–8000

0 256 64

32 0.9 0.13 π /2

32 0.92 0.13 π /2

32 0.92 0.13 π /2

32 0.92 0.13 π /2

32 0.9 0.37 3π /4

32 0.9 0.37 3π /4

32 0.9 0.37 3π /4

32 0.9 0.37 3π /4

32 0.9 0.65

π

32 0.9 0.65

π

32 0.88 0.65

π

32 0.9 0.65

π

32 0.9 0.65

π

32 0.9 0.65

π

32 0.88 0.65

π

32 0.9 0.65

π

32 0.9 0.15 π /4

16 0.9 0.13 π /4

64 0.9 0.13

π

16 0.9 0.13 π /2

64 0.9 0.13 π /2

32 0.9 0.13 π /2

32 0.9 0.13 π /2

32 0.9 0.13 π /2

32 0.9 0.65

π

35 50 1.65 35 50 1.65 35 50 1.65 35 50 1.65 45 120 1.94 45 120 1.92 45 120 1.91 45 120 1.91 65 200 1.97 65 200 1.97 65 200 1.97 65 200 1.97 65 200 1.97 65 200 1.97 65 200 1.97 65 200 1.97 35 50 1.65 35 50 1.65 35 50 1.65 35 50 1.65 35 50 1.65 35 50 1.65 35 50 1.65 35 50 1.65 65 200 1.97

Yes 0.065 0.064 Yes 0.065 0.065 Yes 0.065 0.054 Yes 0.065 0.059 Yes 0.25 0.12 Yes 0.25 0.09 Yes 0.25 0.13 Yes 0.25 0.099 Yes 0.28 0.18 Yes 0.28 0.16 Yes 0.28 0.21 Yes 0.28 0.18 No 0.28 0.35 No 0.28 0.34 No 0.28 0.341 No 0.28 0.307 Yes 0.16 0.052 Yes 0.065 0.063 Yes 0.065 0.062 Yes 0.065 0.061 Yes 0.065 0.061 Yes 0.065 0.052 Yes 0.065 0.065 Yes 0.25 0.048 No 0.5 0.377

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

et al. 2008a). However, a detailed study for thin discs is yet to be performed.
We consider a range of magnetic ﬁeld geometries described by the vector potential Aμ which is related to the Faraday tensor by Fμν = Aν,μ − Aμ,ν. As in S08, we consider a general multiple-loop ﬁeld geometry corresponding to N separate loop bundles stacked radially within the initial disc. The vector potential we use is given by

Aφ,N ∝ Q2 sin

log(r /S ) λﬁeld/(2π r)

[1 + w(ranc − 0.5)] ,

(33)

where ranc is a random number generator for the domain 0 to 1 (see

below for a discussion of perturbations.). All other Aμ are initially zero. All our multi-loop and one-loop simulations have S = 22M,

and the values of λﬁeld/(2π r) are listed in Table 1. For multi-loop models, each additional ﬁeld loop bundle has opposite polarity. We use Q = (ug/ug,max − 0.2) (r/M)3/4, and set Q = 0 if either r < S or Q < 0, and ug,max is the maximum value of the internal energy density in the torus. By comparison, in S08, we set S = 1.1rin, rin = 20M, λﬁeld/(2π r) = 0.16, such that there were two loops centred at r = 28M and 38M. The intention of introducing multiple-loop

bundles is to keep the aspect ratio of the bundles roughly 1:1

in the poloidal plane, rather than loop(s) that are highly elon-

gated in the radial direction. For each disc thickness, we tune

λﬁeld/(2π r) in order to obtain initial poloidal loops that are roughly isotropic.

As in S08, the magnetic ﬁeld strength is set such that the plasma

β parameter satisﬁes βmaxes ≡ pg,max/pb,max = 100, where pg,max is the maximum thermal pressure and pb,max is the maximum magnetic pressure in the entire torus. Since the two maxima never occur at

the same location, β = pg/pb varies over a wide range of values within the disc. This approach is similar to how the magnetic ﬁeld

was normalized in other studies (Gammie et al. 2003; McKinney

& Gammie 2004; McKinney 2006b; McKinney & Narayan 2007b;

Komissarov & McKinney 2007). It ensures that the magnetic ﬁeld

is weak throughout the disc. Care must be taken with how one

normalizes any given initial magnetic ﬁeld geometry. For example,

for the one-loop ﬁeld geometry used by McKinney & Gammie

(2004), if one initializes the ﬁeld with a mean (volume-averaged) β¯ = 100, then the inner edge of the initial torus has β ∼ 1 and the

initial disc is not weakly magnetized.

For most models, the vector potential at all grid points was ran-

domly perturbed by 2 per cent (w in equation 33) and the inter-

nal energy density at all grid points was randomly perturbed by 10 per cent.5 If the simulation starts with perturbations of the vector potential, then we compute tot (used to obtain ˜ r) using the
pre-perturbed magnetic ﬂux in order to gauge how much ﬂux is

dissipated due to the perturbations. Perturbations should be large

enough to excite the non-axisymmetric MRI in order to avoid the

axisymmetric channel solution, while they should not be so large as

to induce signiﬁcant dissipation of the magnetic energy due to grid-

scale magnetic dissipation just after the evolution begins. For some

models, we studied different amplitudes for the initial perturbation

in order to ensure that the amplitude does not signiﬁcantly affect

our results. For a model with |h/r| ∼ 0.07, a/M = 0, and a single

polarity ﬁeld loop, one simulation was initialized with 2 per cent

vector potential perturbations and 10 per cent internal energy per-

turbations, while another otherwise similar simulation was given no

5 In S08, we had a typo saying we perturbed the ﬁeld by 50 per cent, while it was actually perturbed the same as these models, i.e.: 2 per cent vector potential perturbations and 10 per cent internal energy perturbations.

Simulations of magnetized discs 759
seed perturbations. Both become turbulent at about the same time t ∼ 1500M. The magnetic ﬁeld energy at that time is negligibly different, and there is no evidence for signiﬁcant differences in any quantities during inﬂow equilibrium.
4.2 Numerical methods
We perform simulations using the GRMHD code HARM that is based upon a conservative shock-capturing Godunov scheme. One key feature of our code is that we use horizon-penetrating KS coordinates for the Kerr metric (Gammie, McKinney & To´th 2003; McKinney & Gammie 2004; McKinney 2006a; Noble et al. 2006; Mignone & McKinney 2007; Tchekhovskoy, McKinney & Narayan 2007), which avoids any issues with the coordinate singularity in BL coordinates. Even with KS coordinates, one must ensure that the inner-radial boundary of the computational domain is outside the so-called inner horizon (at r/M ≡ 1 − 1 − (a/M)2) so that the equations remain hyperbolic, and one must ensure that there are plenty of grid cells spanning the region near the horizon in order to avoid numerical diffusion out of the horizon.
Another key feature of our code is the use of a third-order accurate (fourth-order error) PPM scheme for the interpolation of primitive quantities (i.e. rest-mass density, four-velocity relative to a ZAMO observer and lab-frame three-magnetic ﬁeld) (McKinney 2006b). The interpolation is similar to that described in Colella & Woodward (1984), but we modiﬁed it to be consistent with interpolating through point values of primitives rather than average values. We do not use the PPM steepener, but we do use the PPM ﬂattener that only activates in strong shocks (e.g. in the initial bow shock off the torus surface, but rarely elsewhere). The PPM scheme attempts to ﬁt a monotonic third-order polynomial directly through the grid face where the dissipative ﬂux enters in the Godunov scheme. Only if the polynomial is non-monotonic does the interpolation reduce order and create discontinuities at the cell face, and so only then does it introduce dissipative ﬂuxes. It therefore leads to extremely small dissipation compared to the original schemes used in HARM, such as the ﬁrst-order accurate (second-order error) minmod or monotonized central (MC) limiter type schemes that always create discontinuities (and so dissipative ﬂuxes) at the cell face regardless of the monotonicity for any primitive quantity that is not linear in space.
Simulations of fully 3D models of accreting BHs producing jets using our 3D GRMHD code show that this PPM scheme leads to an improvement in effective resolution by at least factors of roughly 2 per dimension as compared to the original HARM MC limiter scheme for models with resolution 256 × 128 × 32 (McKinney & Blandford 2009). The PPM method is particularly well-suited for resolving turbulent ﬂows since they rarely have strong discontinuities and have most of the turbulent power in long wavelength modes. Even moving discontinuities are much more accurately resolved by PPM than minmod or MC. For example, even without a steepener, a simple moving contact or moving magnetic rotational discontinuity is sharply resolved within about four cells using the PPM scheme as compared to being diffusively resolved within about 8–15 cells by the MC limiter scheme.
A second-order Runge–Kutta method-of-lines scheme is used to step forward in time, and the time-step is set by using the fast magnetosonic wavespeed with a Courant factor of 0.8. We found that a fourth-order Runge–Kutta scheme does not signiﬁcantly improve accuracy, since most of the ﬂow is far beyond the grid cells inside the horizon that determine the time-step. The standard HARM HLL

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

760 R. F. Penna et al.
scheme is used for the dissipative ﬂuxes, and the standard HARM To´th scheme is used for the magnetic ﬁeld evolution.

4.3 Numerical model setup

The code uses uniform internal coordinates (t, x(1), x(2), x(3)) mapped to the physical coordinates (t, r, θ , φ). The radial grid mapping is

r(x(1)) = R0 + exp (x(1)),

(34)

which spans from Rin to Rout. The parameter R0 = 0.3M controls the resolution near the horizon. Absorbing (outﬂow, no inﬂow allowed) boundary conditions are used. The θ-grid mapping is

θ (x(2)) = [Y (2x(2) − 1) + (1 − Y )(2x(2) − 1)7 + 1](π/2), (35)

where x(2) ranges from 0 to 1 (i.e. no cut-out at the poles) and Y is an adjustable parameter that can be used to concentrate grid zones toward the equator as Y is decreased from 1 to 0. Roughly half of the θ resolution is concentrated in the disc region within ±2|h/r| of the midplane. The HR07 and HR2 models listed in Table 1 have 11 cells per |h/r|, while the HR1 and HR3 models have seven cells per |h/r|. The high-resolution run, C6, has 22 cells per |h/r|, while the low-resolution model, C5, has ﬁve cells per |h/r|. For Y = 0.15 this grid gives roughly six times more angular resolution compared to the grid used in McKinney & Gammie (2004) given by equation (8) with h = 0.3. Reﬂecting boundary conditions are used at the polar axes.
The φ-grid mapping is given by φ(x(3)) = 2π x(3), such that x(3) varies from 0 to 1/8, 1/4, 3/8, 1/2 for boxes with φ = π /4, π /2, 3π /4, π , respectively. Periodic boundary conditions are used in the φ-direction. In all cases, the spatial integrals are renormalized to refer to the full 2π range in φ, even if our computational box size is limited in the φ-direction. We consider various φ in order to check whether this changes our results. Previous GRMHD simulations with the full φ = 2π extent suggest that φ = π /2 is sufﬁcient since coherent structures only extend for about one radian (see ﬁg. 12 in Schnittman, Krolik & Hawley 2006). However, in other GRMHD studies with φ = 2π , the m = 1 mode was found to be dominant, so this requires further consideration (McKinney & Blandford 2009). Note that S08 used φ = π /4, while both N09 and N10 used φ = π /2.
The duration of our simulations with the thinnest discs varies from approximately 20000M to 30000M in order to reach inﬂow equilibrium and to minimize ﬂuctuations in time-averaged quantities. We ensure that each simulation runs for a couple of viscous times in order to reach inﬂow equilibrium over a reasonable range of radius. Note that the simulations cannot be run for a duration longer than tacc ∼ Mdisc(t = 0)/M˙ ∼ 105M, corresponding to the time-scale for accreting a signiﬁcant fraction of the initial torus. We are always well below this limit.
Given ﬁnite computational resources, there is a competition between duration and resolution of a simulation. Our simulations run for relatively long durations, and we use a numerical resolution of Nr × Nθ × Nφ = 256 × 64 × 32 for all models (except those used for convergence testing). In S08 we found this resolution to be sufﬁcient to obtain convergence compared to a similar 512 × 128 × 32 model with φ = π /4. In this paper, we explicitly conﬁrm that our resolution is sufﬁcient by convergence testing our results (see section 6). Near the equatorial plane at the ISCO, the grid aspect ratio in dr : rdθ : r sin θ dφ is 2:1:7, 1:1:4, 1:1:3 and 1:1:3, respectively, for our HR07, HR1, HR2 and HR3 models. The 2:1:7 grid aspect ratio for the HR07 model was found to be sufﬁcient in S08. A grid aspect ratio of 1:1:1 would be preferable in order to ensure the dissipation

is isotropic in Cartesian coordinates, since in Nature one would not expect highly anisotropic dissipation on the scale resolved by our grid cells. However, ﬁnite computational resources require a balance between a minimum required resolution, grid aspect ratio and duration of the simulation.
As described below, we ensure that the MRI is resolved in each simulation both as a function of space and as a function of time by measuring the number of grid cells per fastest growing MRI mode:

QMRI ≡ λMRI ≈ 2π |vAθˆ |/|

(r, θ)| ,

θˆ θˆ

(36)

where θˆ ≡ |eθˆ μdxμ| is the comoving orthonormal θ -directed

grid cell length, eμνˆ is the orthonormal contravariant tetrad system

iAnlftvhee´nloscpaeleﬂd,ubidθˆ-f≡rameθˆeμ,b|μvAθˆis|

= the

bθˆ bθˆ /(b2 comoving

+ ρ0 + ug + pg) is the orthonormal θ-directed

four-ﬁeld and | (r, θ)| is the temporally and azimuthally averaged

absolute value of the orbital frequency.

During the simulation, the rest-mass density and internal en-

ergy densities can become quite low beyond the corona, but the

code only remains accurate and stable for a ﬁnite value of b2/ρ0, b2/ug, and ug/ρ0 for any given resolution. We enforce b2/ρ0 104, b2/ug 104 and ug/ρ0 104 by injecting a sufﬁcient amount

of mass or internal energy into a ﬁxed zero angular momentum

observer (ZAM√O) frame with four-velocity uμ = {−α, 0, 0, 0}, where α = 1/ −gtt is the lapse. In some simulations, we have to

use stronger limits given by b2/ρ0 10, b2/ug 102 and ug/ρ0 10, in order to maintain stability and accuracy. Compared to our

older method of injecting mass-energy into the comoving frame,

the new method avoids run-away injection of energy momentum

in the low-density regions. We have conﬁrmed that this procedure

of injecting mass-energy does not contaminate our results for the

accretion rates and other diagnostics.

5 FIDUCIAL MODEL OF A THIN DISC AROUND A NON-ROTATING BLACK HOLE
Our ﬁducial model, A0HR07, consists of a magnetized thin accretion disc around a non-rotating (a/M = 0) BH. This is similar to the model described in S08; however, here we consider a larger suite of diagnostics, a resolution of 256 × 64 × 32, and a computational box with φ = π /2. As mentioned in Section 4.1, the initial torus parameters are set so that the inner edge is at r = 20M, the pressure maximum is at r = 35M and |h/r| 0.1 at the pressure maximum (see Fig. 1).
The initial torus is threaded with magnetic ﬁeld in the multi-loop geometry as described in Section 4.1. For this model, we use four loops in order to ensure that the loops are roughly circular in the poloidal plane. Once the simulation begins, the MRI leads to MHD turbulence which causes angular momentum transport and drives the accretion ﬂow to a quasi-steady-state.
The ﬁducial model is evolved for a total time of 27350M. We consider the period of steady state to be from Ti = 12500M to Tf = 27350M and of duration T = 14850M. All the steady state results described below are obtained by time-averaging quantities over this steady-state period, which corresponds to about 160 orbital periods at the ISCO, 26 orbits at the inner edge of the initial torus (r = 20M) and 11 orbits at the pressure maximum of the initial torus (r = 35M).

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

Simulations of magnetized discs 761

Figure 1. The initial state of the ﬁducial model (A0HR07) consists of weakly magnetized gas in a geometrically thin torus around a non-spinning (a/M = 0) BH. Colour maps have red as the highest values and blue as the lowest values. Panel (a): linear colour map of rest-mass density, with solid lines showing the thickness |h/r| of the initial torus. Note that the BH horizon is at r = 2M, far to the left of the plot, so the torus is clearly geometrically thin. Near the pressure maximum |h/r| 0.1, and elsewhere |h/r| is even smaller. Panel (b): contour plot of b2 overlaid on linear colour map of rest-mass density shows that the initial ﬁeld consists of four poloidal loops centred at r/M = 29, 34, 39, 45. The wiggles in b2 are due to the initial perturbations. Panel (c): linear colour map of the plasma β shows that the disc is weakly magnetized throughout the initial torus. Panel (d): linear colour map of the number of grid cells per fastest growing MRI wavelength, QMRI, shows that the MRI is properly resolved for the primary two loops at the centre of the disc.
5.1 Initial and evolved disc structure
Fig. 1 shows contour plots of various quantities in the initial solution projected on the (R, z) = (r sin θ , r cos θ )-plane. Notice the relatively small vertical extent of the torus. The disc has a thickness of |h/r| ∼ 0.06–0.09 over the radius range containing the bulk of the mass. The four magnetic loops are clearly delineated. The plot of QMRI indicates that the MRI is well-resolved within the two primary loops. The left-most and right-most loops are marginally underresolved, so a slightly slower-growing MRI mode is expected to control the

Figure 2. The evolved state of the ﬁducial model (A0HR07) consists of a weakly magnetized thin disc surrounded by a strongly magnetized corona. All plots show quantities that have been time-averaged over the period 12500M to 27350M. Colour maps have red as highest values and blue as lowest values. Panel (a): linear colour map of rest-mass density, with solid lines showing the disc thickness |h/r|. Note that the rest-mass density drops off rapidly inside the ISCO. Panel (b): linear colour map of b2 shows that a strong magnetic ﬁeld is present in the corona above the equatorial disc. Panel (c): linear colour map of plasma β shows that the β values are much lower than in the initial torus. This indicates that considerable ﬁeld ampliﬁcation has occurred via the MRI. The gas near the equatorial plane has β ∼ 10 far outside the ISCO and approaches β ∼ 1 near the BH. Panel (d): linear colour map of the number of grid cells per fastest growing MRI wavelength, QMRI, shows that the MRI is properly resolved within most of the accretion ﬂow. Note that QMRI (determined by the vertical magnetic ﬁeld strength) is not expected to be large inside the plunging region where the ﬁeld is forced to become mostly radial or above the disc within the corona where the ﬁeld is mostly toroidal.
dynamics in this region. However, the two primary loops tend to dominate the overall evolution of the gas.
Fig. 2 shows the time-averaged solution during the quasi-steadystate period from Ti = 12500M to Tf = 27350M. We refer to the disc during this period as being ‘evolved’ or ‘saturated’. The evolved disc is in steady state up to r ∼ 9M, as expected for the duration of our simulation. The rest-mass density is concentrated in the disc midplane within ±2|h/r|, while the magnetic energy

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

762 R. F. Penna et al.

Figure 3. Magnetic ﬁeld lines (red vectors) and magnetic energy density (grey-scale map) are shown for the ﬁducial model (A0HR07). Panel (a): snapshot of the magnetic ﬁeld structure at time 27200M shows that the disc is highly turbulent for r > rISCO = 6M and laminar for r < rISCO. Panel (b): time-averaged magnetic ﬁeld in the saturated state shows that for r 9M, viz., the region of the ﬂow that we expect to have achieved inﬂow equilibrium, the geometry of the time-averaged magnetic ﬁeld closely resembles that of a split-monopole. The dashed, vertical line marks the position of the ISCO.

Figure 4. Flow stream lines (red vectors) and rest-mass density (grey-scale map) are shown for the ﬁducial model (A0HR07). Panel (a): snapshot of the velocity structure and rest-mass density at time 27200M clearly show MRIdriven turbulence in the interior of the disc. The rest-mass density appears more diffusively distributed than the magnetic energy density shown in Fig. 3(a). Panel (b): time-averaged streamlines and rest-mass density show that for r 9M the velocity ﬁeld is mostly radial with no indication of a steady outﬂow. Time-averaging smooths out the turbulent ﬂuctuations in the velocity. The dashed, vertical line marks the position of the ISCO.

density is concentrated above the disc in a corona. The MRI is properly resolved with QMRI ≈ 6 in the disc midplane.6 The gas in the midplane has plasma β ∼ 10 outside the ISCO and β ∼ 1 near the BH, indicating that the magnetic ﬁeld has been ampliﬁed beyond the initial minimum of β ∼ 100.
Fig. 3 shows the time-averaged structure of the magnetic ﬁeld during the quasi-steady-state period. The ﬁeld has a smooth splitmonopole structure near and inside the ISCO. Beyond r ∼ 9M, however, the ﬁeld becomes irregular, reversing direction more than once. At these radii, the simulation has not reached inﬂow equilibrium.

5.2 Velocities and the viscous time-scale
Fig. 4 shows the velocity structure in the evolved model. The snapshot indicates well-developed turbulence in the interior of the disc at radii beyond the ISCO (r > 6M), but laminar ﬂow inside the ISCO and over most of the corona. The sudden transition from turbulent to laminar behaviour at the ISCO, which is seen also in the magnetic ﬁeld (Fig. 3a), is a clear sign that the ﬂow dynamics are quite
6 Sano et al. (2004) found that having about six grid cells per wavelength of the fastest growing MRI mode during saturation leads to convergent behaviour for the electromagnetic stresses, although their determination of six cells was based upon a second-order van Leer scheme that is signiﬁcantly more diffusive than our PPM scheme. Also, the (time-averaged or single time value of) vertical ﬁeld is already (at any random spatial position) partially sheared by the axisymmetric MRI, and so may be less relevant than the (e.g.) maximum vertical ﬁeld per unit orbital time at any given point that is not yet sheared and so represents the vertical component one must resolve. These issues imply we may only need about four cells per wavelength of the fastest growing mode (as deﬁned by using the time-averaged absolute vertical ﬁeld strength).

Figure 5. Flow stream lines are shown for the ﬁducial model (A0HR07). Panel (a): snapshot of the velocity structure at time 27200M clearly shows MRI-driven turbulence in the interior of the disc. Panel (b): time-averaged streamlines show that for r 9M the velocity ﬁeld is mostly radial. The dashed, vertical line marks the position of the ISCO.
different in the two regions. Thus the ISCO clearly has an effect on the accreting gas. The time-averaged ﬂow shows that turbulent ﬂuctuations are smoothed out within r ∼ 9M. Fig. 5 shows the velocity stream lines using the line integral convolution method to illustrate vector ﬁelds. This ﬁgure again conﬁrms that the accretion ﬂow is turbulent at radii larger than rISCO but it becomes laminar inside the ISCO, and it again shows that time-averaging smooths out turbulent ﬂuctuations out to r ∼ 9M.

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

Simulations of magnetized discs 763

Figure 6. The time-averaged, angle-averaged, rest-mass density-weighted three-velocities and viscous time-scale in the ﬁducial model (A0HR07) are compared with the NT model. Angle-averaging is performed over the disc gas lying within ±2|h/r| of the midplane. Top panel: the orthonormal radial three-velocity (solid line), and the analytical GR estimate given in equation (B7) of Appendix B (dashed line). Agreement for r > rISCO between the simulation and NT model is found when we set α|h/r|2 ≈ 0.000 33. At smaller radii, the gas dynamics is no longer determined by viscosity and hence the two curves deviate. Middle panel: the orthonormal azimuthal three-velocity vφ (solid line) and the corresponding Keplerian three-velocity (dashed line). Bottom panel: the inﬂow equilibrium time-scale tie ∼ −2r/vr (solid line) of the disc gas is compared to the analytical GR thin disc estimate (dashed line). At r ∼ 9M, we see that tie ∼ 2 × 104M. Therefore, the simulation needs to be run for this time period (which we do) before we can reach inﬂow equilibrium at this radius.
Fig. 6 shows components of the time-averaged velocity that are angle-averaged over ±2|h/r| around the midplane (thick dashed lines in Fig. 8). By limiting the range of the θ integral, we focus on the gas in the disc, leaving out the corona-wind-jet. Outside the ISCO, the radial velocity from the simulation agrees well with the analytical GR estimate (equation B7 in Appendix B). By making this comparison, we found α|h/r|2 ≈ 0.00033. For our disc thickness |h/r| = 0.064, this corresponds to α ≈ 0.08, which is slightly smaller than the nominal estimate α ∼ 0.1 we assumed in Section 3.5. As the gas approaches the ISCO, it accelerates rapidly in the radial direction and ﬁnally free falls into the BH. This region of the ﬂow is not driven by viscosity and hence the dynamics here are not captured by the analytical formula.
Fig. 6 also shows the inﬂow equilibrium time tie, which we take to be twice the GR version of the viscous time: tie = −2r/vr . This is our estimate of the time it will take for the gas at a given radius to reach the steady state. We see that, in a time of ∼27350M, the total duration of our simulation, the solution can be in steady state only inside a radius of ∼9M. Therefore, in the time-averaged results described below, we consider the results to be reliable only over this range of radius.
5.3 Fluxes versus time
Fig. 7 shows various ﬂuxes versus time that should be roughly constant once inﬂow equilibrium has been reached. The ﬁgure shows

Figure 7. The ﬁgure shows for the ﬁducial model (A0HR07) the timedependence at the horizon of the mass accretion rate, M˙ (top panel); nominal efﬁciency, e˜, with dashed line showing the NT value (next panel); accreted speciﬁc angular momentum, j , with dashed line showing the NT value (next panel); absolute magnetic ﬂux relative to the initial absolute magnetic ﬂux, ˜ r (next panel); and dimensionless speciﬁc magnetic ﬂux, ϒ (bottom panel). All quantities have been integrated over all angles. The mass accretion rate varies by factors of up to 4 during the quasi-steady-state phase. The nominal efﬁciency is close to, but on average slightly lower than, the NT value. This means that the net energy loss through photons, winds and jets is below the radiative efﬁciency of the NT model. The speciﬁc angular momentum is clearly lower than the NT value, which implies that some stresses are present inside the ISCO. The absolute magnetic ﬂux at the BH horizon grows until it saturates due to local force balance. The speciﬁc magnetic ﬂux ϒ 1, indicating that electromagnetic stresses inside the ISCO are weak and cause less than 7 per cent deviations from NT in j .
the mass ﬂux, M˙ (rH, t), nominal efﬁciency, e˜(rH, t), speciﬁc angular momentum, j (rH, t), normalized absolute magnetic ﬂux, ˜ r(rH, t) (normalized using the unperturbed initial total ﬂux) and speciﬁc magnetic ﬂux, ϒ(rH, t), all measured at the event horizon (r = rH). These ﬂuxes have been integrated over the entire range of θ from 0 to π . The quantities M˙ , e˜ and j appear to saturate already at t ∼ 7000M. However, the magnetic ﬁeld parameters saturate only at ∼12500M. We consider the steady-state period of the disc to begin only after all these quantities reach their saturated values.
The mass accretion rate is quite variable, with root-mean-square (rms) ﬂuctuations of the order of 2. The nominal efﬁciency e˜ is fairly close to the NT efﬁciency, while the speciﬁc angular momentum j is clearly below the NT value. The results indicate that torques are present within the ISCO, but do not dissipate much energy or cause signiﬁcant energy to be transported out of the ISCO. The absolute magnetic ﬂux per unit initial absolute ﬂux, ˜ r, threading the BH grows to about 1 per cent, which indicates that the magnetic ﬁeld strength near the BH is not just set by the amount of magnetic ﬂux in the initial torus. This suggests our results are insensitive to the total absolute magnetic ﬂux in the initial torus. The speciﬁc magnetic ﬂux is ϒ ≈ 0.86 on average. Magnetic stresses are relatively weak since ϒ 1, which implies the magnetic ﬁeld contributes no more than 7 per cent to deviations from NT in j (Gammie 1999) (see Appendix A). During the quasi-steady-state period, the small

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

764 R. F. Penna et al.
deviations from NT in j are correlated in time with the magnitude of ϒ. This is consistent with the fact that the speciﬁc magnetic ﬂux controls these deviations. Also note that ˜ r is roughly constant in time while ϒ varies in time. This is clearly because M˙ is varying in time and also consistent with the fact that ϒ and M˙ are anti-correlated in time.

5.4 Disc thickness and ﬂuxes versus radius Fig. 8 shows the time-averaged disc thickness of the ﬁducial model as a function of radius. Both measures of thickness deﬁned in Section 3.1 are shown; they track each other. As expected, our primary thickness measure, |h/r|, is the smaller of the two. This thickness measure varies by a small amount across the disc, but it is generally consistent with the following ﬁducial value, viz., the value |h/r| = 0.064 at r = 2rISCO = 12M.
Fig. 9 shows the behaviour of various ﬂuxes versus radius for the full θ integration range (0 to π ). We see that the mass accretion rate, M˙ , and the speciﬁc angular momentum ﬂux, j , are constant up to a radius r ∼ 9M. This is exactly the distance out to which we expect inﬂow equilibrium to have been established, given the inﬂow velocity and viscous time-scale results discussed in Section 5.2. The consistency of these two measurements gives us conﬁdence that the simulation has truly achieved steady-state conditions inside r = 9M. Equally clearly, and as also expected, the simulation is not in steady state at larger radii.
The second panel in Fig. 9 shows that the inward angular momentum ﬂux, j in, agrees reasonably well with the NT prediction. It falls below the NT curve at large radii, i.e. the gas there is sub-Keplerian. This is not surprising since we have included the contribution of the corona-wind-jet gas which, being at high latitude, does not rotate at the Keplerian rate. Other quantities, described below, show a sim-
Figure 8. The time-averaged scaleheight, |h/r|, versus radius in the ﬁducial model (A0HR07) is shown by the solid lines. The above-equator and belowequator values of the disc thickness are |h/r| ∼ 0.04–0.06 in the inﬂow equilibrium region r < 9M. We use the speciﬁc value of |h/r| = 0.064 as measured at r = 2rISCO (light dashed lines) as a representative thickness for the entire ﬂow. Twice this representative thickness (thick dashed lines) is used to ﬁx the θ range of integration for averaging when we wish to focus only on the gas in the disc instead of the gas in the corona-wind-jet. The root mean square thickness (h/r)rms ∼ 0.07–0.13 is shown by the dotted lines.

Figure 9. Mass accretion rate and speciﬁc ﬂuxes are shown as a function of radius for the ﬁducial model (A0HR07). From top to bottom the panels show the following. Top panel: mass accretion rate. Second panel: the accreted speciﬁc angular momentum, j (dotted line), j in (solid line) and the NT proﬁle (dashed line). Third panel: the nominal efﬁciency e˜ (solid line) and the NT proﬁle (dashed line). Bottom panel: the speciﬁc magnetic ﬂux ϒ. For all quantities the integration range includes all θ . The mass accretion rate and j are roughly constant out to r ∼ 9M, as we would expect for inﬂow equilibrium. The proﬁle of j in lies below the NT value at large radii because we include gas in the slowly rotating corona. At the horizon, j and e˜ are modestly below the corresponding NT values. The quantity ϒ ∼ 0.86 and is roughly constant out to r ∼ 6M, indicating that electromagnetic stresses are weak inside the ISCO.
ilar effect due to the corona. At the horizon, j in = 3.286, which is 5 per cent lower than the NT value. This deviation is larger than that found by S08. Once again, it is because we have included the gas in the corona-wind-jet, whereas S08 did not.
The third panel in Fig. 9 shows that the nominal efﬁciency e˜ at the horizon lies below the NT prediction. This implies that the full accretion ﬂow (disc+corona+wind+jet) is radiatively less efﬁcient than the NT model. However, the overall shape of the curve as a function of r is similar to the NT curve. The ﬁnal panel in Fig. 9 shows the value of ϒ versus radius. We see that ϒ ≈ 0.86 is constant out to r ∼ 6M. A value of ϒ ∼ 1 would have led to 7 per cent deviations from NT in j , and only for ϒ ∼ 6.0 would deviations become 50 per cent (see Appendix A). The fact that ϒ ∼ 0.86 1 indicates that electromagnetic stresses are weak and cause less than 7 per cent deviations from NT in j . Note that one does not expect ϒ to be constant7 outside the ISCO where the magnetic ﬁeld is dissipating due to MHD turbulence and the gas is forced to be nearly Keplerian despite a sheared magnetic ﬁeld.
As we have hinted above, we expect large differences between the properties of the gas that accretes in the disc proper, close to the

7 We also ﬁnd that tationdθldaφw√’ −ogf |Bﬁerl|d,

the ideal MHD lines, F(r) ≡

invardiθandtφ√re−lagte|dvr

to Bφ

the ‘isoro− vφBr| /

is nearly Keplerian outside the ISCO and is (as pre-

dicted by the Gammie 1999 model) roughly constant from the ISCO to

the horizon (see also McKinney & Gammie 2004; McKinney & Narayan

2007a).

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

Simulations of magnetized discs 765
In summary, a comparison of Figs 9 and 10 shows that all aspects of the accretion ﬂow in the ﬁducial simulation agree much better with the NT prediction when we restrict our attention to regions close to the midplane. In other words, the gas in the disc proper, deﬁned here as the region lying within ±2|h/r| of the midplane, is well described by the NT model. The deviation of the angular momentum ﬂux j in or j at the horizon relative to NT is 3 per cent, similar to the deviation found by S088, while the nominal efﬁciency e˜ agrees to within ∼1 per cent.

Figure 10. Similar to Fig. 9, but here the integration range only includes angles within ±2|h/r| = ±0.128 radians of the midplane. This allows us to focus on the disc gas. The mass accretion rate is no longer constant because streamlines are not precisely radial. The quantities shown in the second and third panels are not affected by the non-constancy of M˙ because they are ratios of time-averaged ﬂuxes within the equatorial region and are related to ideal MHD invariants. As compared to Fig. 9, here we ﬁnd that j , j in and e˜ closely follow the NT model. For example, j (rH) = 3.363 is only 2.9 per cent less than NT. This indicates that the disc and coronal regions behave quite differently. As one might expect, the disc region behaves like the NT model, while the corona-wind-jet does not. The speciﬁc magnetic ﬂux is even smaller than in Fig. 9 and is ϒ ∼ 0.45, which indicates that electromagnetic stresses are quite weak inside the disc near the midplane.
midplane, and that which ﬂows in the corona-wind-jet region. To focus just on the disc gas, we show in Fig. 10 the same ﬂuxes as in Fig. 9, except that we have restricted the θ range to π /2 ± 2|h/r|. The mass accretion rate is no longer perfectly constant for r < 9M. This is simply a consequence of the fact that the ﬂow streamlines do not perfectly follow the particular constant 2|h/r| disc boundary we have chosen. The non-constancy of M˙ does not signiﬁcantly affect the other quantities plotted in this ﬁgure since they are all normalized by the local M˙ .
The speciﬁc angular momentum, speciﬁc energy and speciﬁc magnetic ﬂux are clearly shifted closer to the NT values when we restrict the angular integration range. Compared to the NT value, viz., j NT(rH) = 3.464, the ﬁducial model gives j (rH) = 3.363 (2.9 per cent less than NT) when integrating over ±2|h/r| around the midplane (i.e. only over the disc gas) and gives j (rH) = 3.266 (5.7 per cent less than NT) when integrating over all θ (i.e. including the corona-wind-jet). Even though the mass accretion rate through the corona-wind-jet is much lower than in the disc, still this gas contributes essentially as much to the deviation of the speciﬁc angular momentum as the disc gas does. In the case of the speciﬁc magnetic ﬂux, integrating over ±2|h/r| around the midplane we ﬁnd ϒ ≈ 0.45, while when we integrate over all angles ϒ ≈ 0.86. The Gammie (1999) model of an equatorial (thin) magnetized ﬂow within the ISCO shows that deviations in the speciﬁc angular momentum are determined by the value of ϒ. We ﬁnd that the measured values of ϒ are able to roughly predict the measured deviations from NT in j .

5.5 Comparison with Gammie (1999) model
Fig. 11 shows a comparison between the ﬁducial model and the Gammie (1999) model of a magnetized thin accretion ﬂow within the ISCO (see also Appendix A). Quantities have been integrated within ±2|h/r| of the midplane and time-averaged over a short period from t = 17400M to t = 18400M. Note that time-averaging b2, ρ0, etc. over long periods can lead to no consistent comparable solution if the value of ϒ varies considerably during the period used for averaging. Also, note that the presence of vertical stratiﬁcation, seen in Figs 9 and 10 showing that ϒ depends upon height, means the vertical-averaging used to obtain ϒ can sometimes make it difﬁcult to compare the simulations with the Gammie (1999) model which has no vertical stratiﬁcation. In particular, using equation (24) over this time period, we ﬁnd that ϒ ≈ 0.2 0.3, 0.44, 0.7, 0.8 for integrations around the midplane of, respectively, ±0.01, ±0.05, ±2|h/r|, ±π /4, ±π /2, with best matches to the Gammie model (i.e. b2/2 and other quantities match) using an actual value of ϒ = 0.2, 0.33, 0.47, 0.8, 0.92. This indicates that stratiﬁcation likely causes our diagnostic to underestimate the best match with the Gammie model once the integration is performed over highly stratiﬁed regions. However, the consistency is fairly good considering how much ϒ varies with height.
Overall, Fig. 11 shows how electromagnetic stresses control the deviations from NT within the ISCO. The panels with D[j ] and D[e] show how the electromagnetic ﬂux starts out large at the ISCO and drops to nearly zero on the horizon. This indicates the electromagnetic ﬂux has been converted into particle ﬂux within the ISCO by ideal (non-dissipative) electromagnetic stresses.9 The simulated magnetized thin disc agrees quite well with the Gammie solution, in contrast to the relatively poor agreement found for thick discs (McKinney & Gammie 2004). Only the single parameter ϒ determines the Gammie solution, so the agreement with the value and radial dependence among multiple independent terms is a strong validation that the Gammie model is working well. Nevertheless, there are some residual deviations near the ISCO where the thermal pressure dominates the magnetic pressure. Even if deviations from NT are present right at the ISCO, the total deviation of the particle ﬂux between the ISCO and horizon equals the deviation predicted by the Gammie (1999) model, as also found in McKinney & Gammie (2004) for thick discs. This indicates that the Gammie (1999) model accurately predicts the effects of electromagnetic stresses inward of the ISCO.
8 The quantities j in and j are nearly equal at the horizon in the calculations reported here whereas they were different in S08. This is because S08 used an alternate deﬁnition of j in. If we had used that deﬁnition here, we would have found a deviation of ∼2 per cent in j in, just as in S08. 9 This behaviour is just like that seen in ideal MHD jet solutions, but inverted with radius.

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

766 R. F. Penna et al.

Figure 11. Comparison between the accretion ﬂow (within ±2|h/r| around the midplane) in the ﬁducial model (A0HR07), shown by solid lines, and
the model of a magnetized thin accretion disc (inﬂow solution) within the
ISCO by Gammie (1999), shown by dashed lines. In all cases the red vertical
line shows the location of the ISCO. Top-left panel shows the radial fourvelocity, where the Gammie solution assumes ur = 0 at the ISCO. Finite thermal effects lead to non-zero ur at the ISCO for the simulated disc.
Bottom-left panel shows the rest-mass density (ρ0, black line), the internal energy density (ug, magenta line) and magnetic energy density (b2/2, green line). Top-right and bottom-right panels show the per cent deviations from
NT for the simulations and Gammie solution for the speciﬁc particle kinetic ﬂux (uμ, black line), speciﬁc enthalpy ﬂux ((ug + pg) uμ/ρ0, magenta line) and speciﬁc electromagnetic ﬂux [(b2ur uμ − br bμ)/(ρ0 ur ), green line], where for j we use μ = φ and for e we use μ = t. As usual, the simulation result for the speciﬁc ﬂuxes is obtained by a ratio of ﬂux integrals instead
of the direct ratio of ﬂux densities. The total speciﬁc ﬂux is constant versus
radius and is a sum of the particle, enthalpy and electromagnetic terms. This ﬁgure is comparable to ﬁg. 10 for a thick (|h/r| ∼ 0.2–0.25) disc in McKinney & Gammie (2004). Finite thermal pressure effects cause the
ﬁducial model to deviate from the inﬂow solution near the ISCO, but the
solutions rapidly converge inside the ISCO and the differences between
the simulation result and the Gammie model (relative to the total speciﬁc
angular momentum or energy) are less than 0.5 per cent.

Figure 12. Luminosity per unit rest-mass accretion rate versus radius (top panel) and the logarithmic derivative of this quantity (bottom panel) are shown for the ﬁducial model (A0HR07). The integration includes all θ angles. The simulation result (solid lines, truncated into dotted lines outside the radius of inﬂow equilibrium) shows that the accretion ﬂow emits more radiation than the NT prediction (dashed lines) at small radii. However, the excess luminosity within the ISCO is only L˜ in ≈ 3.5 per cent, where e˜[NT] is the NT efﬁciency at the horizon (or equivalently at the ISCO).
deviation of the simulation results from the NT model, we estimate what fraction of the disc luminosity is emitted inside the ISCO; recall that the NT model predicts zero luminosity here. The ﬁducial simulation gives L(< rISCO)/M˙ = 0.0021, which is 3.5 per cent of the nominal efﬁciency e˜[NT] = 0.058 of a thin NT disc around a non-spinning BH. This shows that the excess luminosity radiated within the ISCO is quite small. The relative luminosity within the ISCO is L˜ in = 3.5 per cent and the relative luminosity within the inﬂow equilibrium region is L˜ out = 8.0 per cent. Hence, we conclude that, for accretion discs which are as thin as our ﬁducial model, viz., |h/r| ∼ 0.07, the NT model provides a good description of the luminosity proﬁle.

Finally, note that the electromagnetic stresses within the ISCO are ideal and non-dissipative in the Gammie model. Since the ﬂow within the ISCO in the simulation is mostly laminar leading to weak non-ideal (resistive or viscous) effects, the dissipative contribution (which could lead to radiation) can be quite small. An exception to this is the presence of extended current sheets, present near the equator within the ISCO in the simulations, whose dissipation requires a model of the (as of yet, poorly understood) rate of relativistic reconnection.
5.6 Luminosity versus radius
Fig. 12 shows radial proﬁles of two measures of the disc luminosity: L(< r)/M˙ , which is the cumulative luminosity inside radius r, and d(L/M˙ )/d ln r, which gives the local luminosity at r. We see that the proﬁles from the simulation are quite close to the NT prediction, especially in the steady-state region. As a way of measuring the

5.7 Luminosity from disc versus corona-wind-jet
The ﬁducial model described so far includes a tapering of the cooling rate as a function of height above the midplane, given by the function S[θ ] (see equation 5). We introduced this taper in order to only cool bound [−ut(ρ0 + ug + pg + b2)/ρ0 < 1] gas and to avoid including the emission from the part of the corona-wind-jet that is prone to excessive numerical dissipation due to the low resolution used high above the accretion disc. This is a common approach that others have also taken when performing GRMHD simulations of thin discs (N09; N10). However, since our tapering function does not explicitly refer to how bound the gas is, we need to check that it is consistent with cooling only bound gas. We have explored this question by re-running the ﬁducial model with all parameters the same except that we turned off the tapering function altogether, i.e. we set S[θ ] = 1. This is the only model for which the tapering function is turned off.

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

Figure 13. Enclosed luminosity versus radius for models with different cooling prescriptions and θ integration ranges. The black dashed line corresponds to the NT model. The luminosity for the ﬁducial model A0HR07, which includes a tapering of the cooling with disc height as described in Section 2, is shown integrated over ±2|h/r| from the midplane (black dotted line), integrated over all bound gas (black long dashed line), and integrated over all ﬂuid (black solid line). Essentially all the gas is bound and so the black solid and long dashed lines are indistinguishable. The red lines are for a model that is identical to the ﬁducial run, except that no tapering is applied to the cooling. For this model the lines are: red solid line: all angles, all ﬂuid; red dotted line: ±2|h/r| around the midplane; red long dashed line: all bound gas. The main result is that the luminosity from bound gas is nearly the same (especially at the ISCO) whether or not we include tapering (compare the red long dashed line and the black long dashed line).
Fig. 13 shows a number of luminosity proﬁles for the ﬁducial model and the no-tapering model. This comparison shows that, whether or not we include a taper, the results for the luminosity from all the bound gas is nearly the same. Without a tapering, there is some luminosity at high latitudes above ±8|h/r| corresponding to emission from the low-density jet region (black solid line). This region is unbound and numerically inaccurate, and it is properly excluded when we use the tapering function. Another conclusion from the above test is that, as far as the luminosity is concerned, it does not matter much whether we focus on the midplane gas ((π /2) ± 2|h/r|) or include all the bound gas. The deviations of the luminosity from NT in the two cases are similar – changes in the deviation are less than 1 per cent.
An important question to ask is whether the excess luminosity from within the ISCO is correlated with, e.g. deviations from NT in j , since D[j ] could then be used as a proxy for the excess luminosity. We investigate this in the context of the simulation with no tapering. For an integration over ±2|h/r| around the midplane (which we identify with the disc component), or over all bound gas, or over all the gas (bound and unbound), the excess luminosity inside the ISCO is L˜ in = 3.3 per cent, 4.4 per cent, 5.4 per cent, and the deviation from NT in j is D[j ] = −3.6 per cent, −6.7 per cent, −6.7 per cent, respectively. We ignore the luminosity from unbound gas since this is mostly due to material in a very low density region of the simulation where thermodynamics is not evolved accurately. Considering the rest of the results, we see that

Simulations of magnetized discs 767
D[j ] is 100 per cent larger when we include bound gas outside the disc compared to when we consider only the disc gas, whereas the excess luminosity increases by only 32 per cent. Therefore, when we compute j by integrating over all bound gas and then assess the deviation of the simulated accretion ﬂow from the NT model, we strongly overestimate the excess luminosity of the bound gas relative to NT. A better proxy for the latter is the deviations from NT in j integrated only over the disc component (i.e. over ±2|h/r| around the midplane).
Furthermore, we note that the gas that lies beyond ±2|h/r| from the disc midplane consists of coronal gas, which is expected to be optically thin and to emit a power-law spectrum of photons. For many applications, we are not interested in this component but rather care only about the thermal blackbody-like emission from the optically thick region of the disc. For such studies, the most appropriate diagnostic from the simulations is the radiation emitted within ±2|h/r| of the midplane. According to this diagnostic, the excess emission inside the ISCO is only L˜ in = 3.4 per cent in the model without tapering, and 3.5 per cent in the ﬁducial model that includes tapering.
Lastly, we consider variations in the cooling time-scale, τ cool, which is another free parameter of our cooling model that we generally set to 2π / K. However, we consider one model that is otherwise identical to the ﬁducial model except we set τ cool to be ﬁve times shorter so that the cooling rate is ﬁve times faster. We ﬁnd that L˜ in = 4.2 per cent, which is slightly larger than the ﬁducial model with L˜ in = 3.5 per cent. Even though the cooling rate is ﬁve times faster than an orbital rate, there is only 20 per cent more luminosity from within the ISCO. This is likely due to the ﬂow within the ISCO being mostly laminar with little remaining turbulence to drive dissipation and radiation.
6 CONVERGENCE WITH RESOLUTION AND BOX SIZE
The ﬁducial model described earlier was computed with a numerical resolution of 256 × 64 × 32, using an azimuthal wedge of π /2. This is to be compared with the simulation described in S08, which made use of a 512 × 128 × 32 grid and used a π /4 wedge. These two runs give very similar results, suggesting that the details of the resolution and wedge size are not very important. To conﬁrm this, we have run a number of simulations with different resolutions and wedge angles. The complete list of runs is: 256 × 64 × 32 with
φ = π /2 (ﬁducial run, model A0HR07)), 512 × 128 × 32 with φ = π /4 (S08, model C0), 256 × 128 × 32 with φ = π /2 (model C6), 256 × 32 × 32 with φ = π /2 (model C5), 256 × 64 × 64 with φ = π /2 (model C4), 256 × 64 × 64 with φ = π (model C2), 256 × 64 × 16 with φ = π /2 (model C3) and 256 × 64 × 16 with φ = π /4 (model C1). Fig. 14 shows the accreted speciﬁc angular momentum, j , ingoing component of the speciﬁc angular momentum, j in, and the nominal efﬁciency e˜ as functions of radius for all the models used for convergence testing. Fig. 15 similarly shows the cumulative luminosity L(< r)/M˙ and differential luminosity d(L/M˙ )/d ln r as functions of radius. The overwhelming impression from these plots is that the sequence of convergence simulations agrees with one another quite well. Also, the average of all the runs matches the NT model very well; this is especially true for the steady-state region of the ﬂow, r < 9M. Thus, qualitatively, we conclude that our results are well-converged. For more quantitative comparison, Fig. 16 shows the proﬁle of j versus r for the various models, this time with each model separately

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

768 R. F. Penna et al.

Figure 14. This plot shows j , j in and e˜ for a sequence of simulations that are similar to the ﬁducial run (A0HR07), viz., |h/r| ≈ 0.07, a/M = 0, but use different radial resolutions, or θ resolutions, or box sizes. The integration range in θ is over ±2|h/r| around the midplane. Only the region of the ﬂow in inﬂow equilibrium, 2M < r < 9M, is shown in the case of j . The different lines are as follows: black dashed line: NT model; black solid line: ﬁducial model A0HR07; blue solid line: model C0 (S08); magenta dotted line: model C1; magenta solid line: model C2; red dotted line: model C3; red solid line: model C4; green dotted line: model C5; green solid line: model C6. Note that changes in the numerical resolution or other computational parameters lead to negligible changes in the values of j , j in and e˜ in the region of the ﬂow that is in inﬂow equilibrium, r < 9M. For r 9M, the ﬂow has not achieved steady state, which explains the large deviations in e˜. Only the lowest resolution models are outliers.

Figure 15. Similar to Fig. 14, but for the normalized luminosity, L(< r)/M˙ , and its logarithmic derivative, d(L/M˙ )/d ln r, both shown versus radius. We
see that all the models used to test convergence show consistent luminosity
proﬁles over the region that is in inﬂow equilibrium, r < 9M. The wellconverged models have L˜ in 4 per cent, which indicates only a low level of luminosity inside the ISCO.

identiﬁed. It is clear that j has converged, since there are very minor deviations from our highest resolution/largest box size to our next highest resolution/next largest box size. All other quantities, including e˜, jin and ϒ are similarly converged. The model with Nφ = 64 shows slightly less deviations from NT in j than our other models. However, it also shows slightly higher luminosity than our other models. This behaviour is likely due to the stochastic temporal behaviour of all quantities versus time, but this could also be due to the higher φ-resolution causing a weaker ordered magnetic ﬁeld to be present leading to weaker ideal electromagnetic stresses, smaller deviations from NT in j within the ISCO, but with the remaining turbulent ﬁeld being dissipated giving a higher luminosity. The Nφ resolution appears to be the limit on our accuracy.
Further quantitative details are given in Table 2, where we list numerical results for all the convergence test models, with the θ integration performed over both ±2|h/r| around the midplane and over all angles. We see that there are some trends as a function of resolution and/or φ. Having only 32 cells in θ or 16 cells in φ gives somewhat poor results, so these runs are underresolved. However, even for these runs, the differences are not large. Note that ϒ reaches a steady state much later than all other quantities, and our C? (where ? is 0 through 6) models did not run as long as the ﬁducial model. This explains why ϒ is a bit lower for the C? models. Overall, we conclude that our choice of resolution 256 × 64 × 32 for the ﬁducial run (A0HR07) is adequate to reach convergence.

Figure 16. This is a more detailed version of Fig. 14, showing j versus r for individually labelled models. The models correspond to the ﬁducial resolution (solid lines), a higher resolution run (dot–dashed lines) and a lower resolution run (dotted lines). Generally, there are only minor differences between the ﬁducial and higher resolution models.
7 DEPENDENCE ON BLACK HOLE SPIN AND DISC THICKNESS
In addition to the ﬁducial model and the convergence runs described in previous sections, we have run a number of other simulations to explore the effect of the BH spin parameter a/M and the disc thickness |h/r| on our various diagnostics: j, jin, e˜, the luminosity

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

Simulations of magnetized discs 769

Table 2. Convergence study.

Model name |h/r|

M˙

± 2|h/r| A0HR07
C0 C1 C2 C3 C4 C5 C6 All θ A0HR07 C0 C1 C2 C3 C4 C5 C6

0.064 0.052 0.063 0.062 0.061 0.061 0.052 0.065
0.064 0.052 0.063 0.062 0.061 0.061 0.052 0.065

0.066 0.064 0.041 0.063 0.026 0.054 0.008 0.088
0.074 0.071 0.042 0.069 0.029 0.057 0.009 0.092

e˜
0.058 0.059 0.056 0.058 0.058 0.058 0.055 0.059
0.054 0.057 0.055 0.056 0.057 0.057 0.053 0.058

D[e˜]
−0.829 −2.708
2.204 −0.699 −1.969 −1.406
3.955 −3.355
4.723 0.738 3.680 1.664 0.893 0.811 7.687 −1.813

j
3.363 3.363 3.415 3.360 3.358 3.385 3.417 3.355
3.266 3.312 3.398 3.307 3.325 3.358 3.338 3.334

D[j ]
−2.913 −2.905 −1.406 −3.016 −3.060 −2.286 −1.364 −3.155
−5.717 −4.392 −1.894 −4.541 −4.008 −3.075 −3.636 −3.761

j in
3.355 3.351 3.408 3.333 3.339 3.378 3.427 3.333
3.281 3.291 3.389 3.262 3.305 3.351 3.353 3.306

D[j in]
−3.153 −3.259 −1.621 −3.789 −3.610 −2.489 −1.067 −3.778
−5.275 −5.002 −2.178 −5.846 −4.580 −3.255 −3.218 −4.560

s
3.363 3.363 3.415 3.360 3.358 3.385 3.417 3.355
3.266 3.312 3.398 3.307 3.325 3.358 3.338 3.334

L˜ in
0.035 0.019 0.034 0.047 0.039 0.019 0.034 0.054
0.035 0.032 0.037 0.053 0.042 0.020 0.039 0.057

L˜ eq
0.080 0.042 0.072 0.109 0.087 0.018 0.070 0.103
0.053 0.049 0.062 0.080 0.064 0.009 0.059 0.086

100 ˜ r
1.355 0.399 0.584 0.925 0.727 0.714 0.315 0.933
6.677 2.18
0.940 5.757 1.401 2.690 0.726 5.091

ϒ
0.450 0.359 0.223 0.323 0.449 0.296 0.322 0.256
0.863 0.480 0.142 0.710 0.309 0.359 0.243 0.534

and ϒ. We consider four values of the BH spin parameter, viz., a/M = 0, 0.7, 0.9, 0.98, and four disc thicknesses, viz., |h/r| = 0.07, 0.1, 0.2, 0.3. We summarize here our results for this 4 × 4 grid of models.
Geometrically thick discs are expected on quite general grounds to deviate from the standard thin disc model. The inner edge of the disc, as measured for instance by the location of the sonic point, is expected to deviate from the ISCO, the shift scaling roughly as |rin − rISCO| ∝ (cs/vK)2 [cs is sound speed, where c2s = pg/ (ρ0 + ug + pg)]. This effect is seen in hydrodynamic models of thick discs, e.g. Narayan, Kato & Honma (1997) and Abramowicz et al. (2010), where it is shown that rin can move either inside or outside the ISCO; it moves inside when α is small and outside when α is large. In either case, these hydrodynamic models clearly show that, as |h/r| → 0, i.e. as cs/vK → 0, the solution always tends towards the NT model (Shafee et al. 2008b).
While the hydrodynamic studies mentioned above have driven much of our intuition on the behaviour of thick and thin discs, it is an open question whether or not the magnetic ﬁeld plays a signiﬁcant role. In principle, magnetic effects may cause the solution to deviate signiﬁcantly from the NT model even in the limit |h/r| → 0 (Krolik 1999; Gammie 1999). One of the major goals of the present paper is to investigate this question. We show in this section that, as |h/r| → 0, magnetized discs do tend towards the NT model. This statement appears to be true for a range of BH spins. We also show that the speciﬁc magnetic ﬂux ϒ inside the ISCO decreases with decreasing |h/r| and remains quite small. This explains why the magnetic ﬁeld does not cause signiﬁcant deviations from NT in thin discs.
Fig. 17 shows the speciﬁc angular momentum, j , and the ingoing component of this quantity, j in, versus radius for the 4 × 4 grid of models. The θ integral has been taken over ±2|h/r| around the midplane in order to focus on the equatorial disc properties. The value of j is roughly constant out to a radius well outside the ISCO, indicating that we have converged solutions in inﬂow equilibrium extending over a useful range of radius. As discussed in Section 3.5, inﬂow equilibrium is expected within r/M = 9, 7, 5.5, 5, respectively, for a/M = 0, 0.7, 0.9, 0.98. This is roughly consistent with the radius out to which the quantity j (integrated over all angles) is constant, and this motivates why in all such plots we only show j over the region where the ﬂow is in inﬂow equilibrium. The

Figure 17. The net accreted speciﬁc angular momentum, j (the nearly horizontal dotted lines), and the ingoing component of this quantity, j in (the sloping curved lines), as a function of radius for the 4 × 4 grid of models. Each panel corresponds to a single BH spin, a/M = 0, 0.7, 0.9 or 0.98, and shows models with four disc thicknesses, |h/r| ≈ 0.07, 0.1, 0.2, 0.3 (see legend). The θ integral has been taken over ±2|h/r| around the midplane. In each panel, the thin dashed black line, marked by two circles which indicate the location of the horizon and the ISCO, shows the NT solution for j in. As expected, we see that thicker discs exhibit larger deviations from NT. However, as a function of spin, there is no indication that deviations from NT become any larger for larger spins. In the case of the thinnest models with |h/r| ≈ 0.07, the NT model works well for gas close to the midplane for all spins.
four panels in Fig. 17 show a clear trend, viz., deviations from NT are larger for thicker discs, as expected. Interestingly, for higher BH spins, the relative deviations from NT actually decrease.
Fig. 18 shows the nominal efﬁciency, e˜, as a function of radius for the 4 × 4 grid of models. Our thickest disc models (|h/r| ≈

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

770 R. F. Penna et al.

Figure 18. Similar to Fig. 17, but for the nominal efﬁciency, e˜. For thin (|h/r| 0.1) discs, the results are close to NT for all BH spins. As expected, the thicker models deviate signiﬁcantly from NT. In part this is because the ad hoc cooling function we use in the simulations is less accurate for thick discs, and in part because the models with |h/r| ≈ 0.3 have no cooling and start with marginally bound/unbound gas that implies e˜ ∼ 0. The a/M = 0.98 models show erratic behaviour at large radii where the ﬂow has not achieved inﬂow equilibrium.
0.3) do not include cooling, so the efﬁciency shown is only due to losses by a wind-jet. We see that the efﬁciency is fairly close to the NT value for all four thin disc simulations with |h/r| ∼ 0.07; even in the worst case, viz., a/M = 0.98, the deviation from NT is only ∼5 per cent. In the case of thicker discs, the efﬁciency shows larger deviations from NT and the proﬁle as a function of radius also looks different. For models with |h/r| ≈ 0.3, there is no cooling so large deviations are expected.
Fig. 19 shows the luminosity, L(< r)/M˙ , versus radius for our 4 × 4 grid of models, focusing just on the region that has reached inﬂow equilibrium. The luminosity is estimated by integrating over all θ angles. Our thickest disc models (|h/r| ≈ 0.3) do not include cooling and so are not plotted. The various panels show that, as |h/r| → 0, the luminosity becomes progressively closer to the NT result in the steady-state region of the ﬂow near and inside the ISCO. Thus, once again, we conclude that the NT luminosity proﬁle is valid for geometrically thin discs even when the accreting gas is magnetized.
A ﬁgure (not shown) that is similar to Fig. 17 but for the speciﬁc magnetic ﬂux indicates that ϒ ≤ 1 within ±2|h/r| near the ISCO for all BH spins and disc thicknesses. For our thinnest models, ϒ ≤ 0.45, for which the model of Gammie (1999) predicts that the speciﬁc angular momentum will deviate from NT by less than 1.9 per cent, 3.0 per cent, 3.8 per cent, 4.2 per cent for BH spins a/M = 0, 0.7, 0.9, 0.98, respectively (see Appendix A). The numerical results from the simulations show deviations from NT that are similar to these values. Thus, overall, our results indicate that electromagnetic stresses are weak inside the ISCO for geometrically thin discs.
Finally, for all models we look at plots (not shown) of M(< r) (mass enclosed within radius), M˙ (r) (total mass accretion

Figure 19. Similar to Fig. 17, but for the normalized luminosity, L(< r)/M˙ . For thin discs, the luminosity deviates only slightly from NT near and inside the ISCO. There is no strong evidence for any dependence on the BH spin. The region at large radii has not reached inﬂow equilibrium and is not shown.
rate versus radius) and [h/r](r) (disc scaleheight versus radius). We ﬁnd that these are consistently ﬂat to the same degree and to the same radius as the quantity j (r) is constant as shown in Fig. 17. This further indicates that our models are in inﬂow equilibrium out to the expected radius.
7.1 Scaling laws versus a/M and |h/r|
We now consider how the magnitude of j, e˜, L(< rISCO) and ϒ scales with disc thickness and BH spin. Table 3 lists numerical results corresponding to θ integrations over ±2|h/r| around the midplane and over all angles.10 Fig. 20 shows selected results corresponding to models with a non-rotating BH for quantities integrated over ±2|h/r|. We see that the deviations of various diagnostics from the NT values scale roughly as |h/r|. In general, the deviations are quite small for the thinnest model with |h/r| ≈ 0.07.
Next, we consider ﬁts of our simulation data as a function of BH spin and disc thickness to reveal if, at all, these two parameters control how much the ﬂow deviates from NT. In some cases we directly ﬁt the simulation results instead of their deviations from NT, since for thick discs the actual measurement values can saturate independent of thickness leading to large non-linear deviations from NT. Before making the ﬁts, we ask how quantities might scale with a/M and |h/r|. With no disc present, the rotational symmetry forces any scaling to be an even power of BH spin (McKinney & Gammie 2004). However, the presence of a rotating disc breaks this symmetry, and any accretion ﬂow properties, such as deviations from NT’s model, could depend linearly upon a/M (at least for small spins). This motivates performing a linear ﬁt in a/M. Similarly, the
10 Some thicker disc models without cooling show small or slightly negative efﬁciencies, e˜, which signiﬁes the accretion of weakly unbound gas. This can occur when a magnetic ﬁeld is inserted into a weakly bound gas in hydrostatic equilibrium.

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

Simulations of magnetized discs 771

Table 3. BH spin and disc thickness study.

Model name
±2|h/r| A0HR07 A7HR07 A9HR07 A98HR07 A0HR1 A7HR1 A9HR1 A98HR1 A0HR2 A7HR2 A9HR2 A98HR2 A0HR3 A7HR3 A9HR3 A98HR3
All θ A0HR07 A7HR07 A9HR07 A98HR07 A0HR1 A7HR1 A9HR1 A98HR1 A0HR2 A7HR2 A9HR2 A98HR2 A0HR3 A7HR3 A9HR3 A98HR3

|h/r|
0.064 0.065 0.054 0.059 0.12 0.09 0.13 0.099 0.18 0.16 0.21 0.18 0.350 0.34 0.341 0.307
0.064 0.065 0.054 0.059 0.12 0.09 0.13 0.099 0.18 0.16 0.21 0.18 0.350 0.34 0.341 0.307

M˙
0.066 0.050 0.045 0.013 4.973 2.443 2.133 2.372 48.286 31.665 40.603 29.410 44.066 41.045 35.852 24.486
0.074 0.065 0.060 0.021 6.036 2.907 2.777 3.235 59.025 41.327 53.746 43.815 49.207 47.146 42.733 30.293

e˜
0.058 0.107 0.156 0.225 0.056 0.099 0.142 0.213 0.055 0.049 0.090 0.191 −0.003 −0.007 −0.001 −0.015
0.054 0.093 0.135 0.171 0.053 0.093 0.128 0.197 0.050 0.046 0.085 0.154 −0.004 −0.007 −0.004 −0.014

D[e˜]
−0.829 −2.899 −0.157
3.897 1.470 4.808 9.014 8.810 4.134 52.330 41.922 18.496 104.582 106.384 100.582 106.206
4.723 10.132 13.294 26.735 6.579 10.343 17.577 15.950 12.631 55.244 45.564 34.174 106.180 107.191 102.370 105.873

j
3.363 2.471 2.042 1.643 3.138 2.446 1.947 1.626 2.774 2.412 1.946 1.588 2.309 1.967 1.722 1.783
3.266 2.282 1.860 1.559 2.908 2.344 1.823 1.425 2.465 2.186 1.739 1.411 2.128 1.788 1.530 1.597

D[j ]
−2.913 −4.465 −2.762 −2.335 −9.424 −5.447 −7.261 −3.393 −19.916 −6.736 −7.315 −5.650 −33.331 −23.970 −17.999
5.987
−5.717 −11.776 −11.404
−7.348 −16.046
−9.376 −13.164 −15.291 −28.830 −15.499 −17.164 −16.120 −38.572 −30.878 −27.117
−5.105

j in
3.355 2.527 2.074 1.679 3.162 2.524 2.124 1.703 2.872 2.576 2.155 1.870 2.408 2.236 1.998 1.886
3.281 2.511 2.209 1.799 2.980 2.457 2.069 1.880 2.596 2.394 2.058 1.876 2.220 2.065 1.869 1.655

D[j in]
−3.153 −2.294 −1.213 −0.199 −8.724 −2.409
1.157 1.200 −17.098 −0.425 2.624 11.117 −30.473 −13.549 −4.832 12.102
−5.275 −2.933
5.222 6.905 −13.961 −4.988 −1.449 11.744 −25.067 −7.453 −2.001 11.497 −35.919 −20.166 −10.977 −1.624

s
3.363 1.220 0.523 0.124 3.138 1.184 0.402 0.084 2.774 1.081 0.309 0.001 2.309 0.557 −0.080 −0.205
3.266 1.012 0.303 −0.065 2.908 1.074 0.254 −0.149 2.465 0.851 0.092 −0.247 2.128 0.377 −0.276 −0.390

L˜ in
0.035 0.048 0.041 0.069 0.084 0.060 0.068 0.062 0.167 0.050 0.011 0.052 0.000 0.000 0.000 0.000
0.035 0.040 0.035 0.048 0.087 0.068 0.066 0.078 0.164 0.045 0.012 0.045 0.000 0.000 0.000 0.000

L˜ eq
0.080 0.084 0.082 0.127 0.134 0.108 0.107 0.112 0.235 0.034 −0.026 0.068 −0.049 −0.060 −0.073 −0.104
0.053 0.040 0.031 0.039 0.110 0.093 0.075 0.094 0.197 0.015 −0.031 0.033 −0.049 −0.060 −0.073 −0.104

100 ˜ r
1.355 0.919 0.455 0.276 2.976 0.909 1.064 0.451 2.518 0.919 0.795 0.651 3.039 1.956 1.437 0.369
6.677 8.496 8.945 2.460 8.880 2.295 4.256 6.599 5.798 2.679 3.799 2.072 4.724 3.433 2.652 0.668

ϒ
0.450 0.393 0.218 0.228 0.871 0.406 0.466 0.254 1.235 0.631 0.557 0.459 1.182 0.746 0.543 0.246
0.863 1.156 1.299 0.626 1.247 0.525 0.735 1.349 1.771 0.923 1.093 0.887 1.331 0.952 0.914 0.273

thickness relates to a dimensionless speed: cs/vK ∼ |h/r|, while there are several different speeds in the accretion problem that could force quantities to have an arbitrary dependence on |h/r|. Although, in principle, deviations might scale as some power of |h/r|, we assume here a linear scaling ∝ |h/r|. This choice is driven partly by simplicity and partly by Fig. 20 which shows that the simulation results agree well with this scaling.
These rough arguments motivate obtaining explicit scaling laws for a quantity’s deviations from NT as a function of a/M and |h/r|. For all quantities we use the full 4 × 4 set of models, except for the luminosity and efﬁciency we exclude the two thickest models in order to focus on the luminosity for thin discs with cooling. We perform a linear least-squares ﬁt in both a/M and |h/r|, and we report the absolute per cent difference between the upper 95 per cent conﬁdence limit (C+) and the best-ﬁtting parameter value (f ) given by E = 100|C+ − f |/|f |. Note that if E > 100 per cent, then the best-ﬁtting value is no different from zero to 95 per cent conﬁdence (such parameter values are not reported). After the linear ﬁt is provided, the value of E is given for each parameter in order of appearance in the ﬁt. Only the statistically signiﬁcant digits are shown.
First, we consider how electromagnetic stresses depend upon a/M and |h/r|. Gammie (1999) has shown that the effects of electromagnetic stresses are tied to the speciﬁc magnetic ﬂux, ϒ, and that for ϒ 1 there are weak electromagnetic stresses causing only

minor deviations (less than 12 per cent for j across all BH spins)

from N√T. Let us consider how ϒ should scale with |h/r|, where ϒ = 2(r/M)2Br /( −(r/M)2ρ0ur ) in the equatorial plane and is assumed to be constant from the ISCO to the horizon. For sim-

plicity, let us study the case of a rapidly rotating BH. First, con-

sider the boundary conditions near the ISCO provided by the disc,

where cs/vK ∼ |h/r| and the Keplerian rotation speed reaches vK ∼ 0.5. This implies cs ∼ 0.5|h/r|. Secondly, consider the ﬂow that connects the ISCO and the horizon. The gas in the disc beyond

the ISCO has β ∼ (cs/vA)2 ∼ 10, but reaches β ∼ 1 inside the

ISCO in any GRMHD simulations of turbulent magnetized discs,

which that ϒ

gives that ∼ 1.4Br

/c√s ρ∼0

vA. Thus, vA ∼ 0.5|h/r|. Finally, notice at the horizon where ur ∼ −1 and r =

M. The Keplerian rotation at the ISCO leads to a magnetic ﬁeld

with orthonormal radial (|Br | ∼ |Br |) and toroidal (|Bφ|) com-

ponents with similar values near |Br | ∼ |Br | ∼ |Bφ| ∼ |b| and so
Alfve´n three-speed is vA = |b|/

the ϒ b2

+∼ISρC10.O4+|bau|n/gd√+hρpo0.griFz∼ounr|t,bhe|g/ri,√vitnhρge0

in any massive disc, so that ϒ ∼ 1.4 vA ∼ 0.7|h/r| for a rapidly

rotating BH. Extending these rough arguments to all BH spins at a

ﬁxed disc thickness also gives that ϒ ∝ −0.8(a/M) for a/M 0.7.

These arguments demonstrate three points: (1) ϒ 1 gives b2/ρ0

1, implying a force-free magnetosphere instead of a massive

accretion disc; (2) ϒ ∝ |h/r| and (3) ϒ ∝ −(a/M). Since the local

condition for the magnetic ﬁeld ejecting mass is b2/ρ0 1 (see

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

772 R. F. Penna et al.

Figure 20. The relative difference between j in the simulation and in the NT model (top panel), the relative difference between the nominal efﬁciency, e˜, and the NT value (middle panel), and the luminosity inside the ISCO normalized by the net radiative efﬁciency of the NT model, L˜ in (bottom panel), where e˜[NT] has been evaluated at the horizon (equivalently at the ISCO). There is a rough linear dependence on |h/r| for all quantities, where a linear ﬁt is shown as a dotted line in each panel. Note that the thicker disc models are not expected to behave like NT, and actually have j roughly similar across all spins. For |h/r| ≈ 0.07, the excess luminosity from within the ISCO is less than 4 per cent of the total NT efﬁciency.

e.g. Komissarov & Barkov 2009), this shows that ϒ ∼ 1 deﬁnes a boundary that the disc component of the ﬂow cannot signiﬁcantly pass beyond without eventually incurring disruption by the strong magnetic ﬁeld within the disc.
We now obtain the actual ﬁt, which for an integration over ±2|h/r| gives

ϒ ≈ 0.7 + h − 0.6 a , rM

(37)

with E = 33 per cent, 70 per cent, 40 per cent, indicating a reasonable ﬁt. There is essentially 100 per cent conﬁdence in the sign of the ﬁrst and third parameters and 98 per cent conﬁdence in the sign of the second parameter. This ﬁt is consistent with our basic analytical estimate for the scaling. Since most likely ϒ ≤ 0.9 in the limit that |h/r| → 0 across all BH spins, the electromagnetic stresses are weak and cause less than 12 per cent deviation from NT in j . This means that the NT solution is essentially recovered for magnetized thin discs. For an integration over all angles, ϒ ≈ 1 with E = 35 per cent, and there is no statistically signiﬁcant trend with disc thickness or BH spin. The value of ϒ ∼ 1 is consistent with the presence of the highly magnetized corona-wind-jet above the disc component (McKinney & Gammie 2004).
Next, we consider whether our simulations can determine the equilibrium value of the BH spin as a function of |h/r|. The spin evolves as the BH accretes mass, energy and angular momentum, and it can stop evolving when these come into a certain balance leading to d(a/M)/dt = 0 (see equation 15). In spin equilibrium, the spin-up parameter s = j − 2(a/M)e has s = 0 and solving for a gives the equilibrium spin aeq/M = j /(2e). For the NT solution, s is fairly linear for a/M > 0 and aeq/M = 1. In Appendix A, we

note that for ϒ ∼ 0.2–1 that the deviations from NT roughly scale as ϒ. Since ϒ ∝ |h/r|, one expects s to also roughly scale with |h/r|. This implies that deviations from NT in the spin equilibrium should scale as |h/r|. Hence, one should have 1 − aeq/M ∝ |h/r|.
Now we obtain the actual ﬁt. We consider two types of ﬁts. In
one case, we ﬁt s (with ﬂuxes integrated over all angles) and solve s = 0 for aeq/M. This gives

s ≈ 3.2 − 2.5 h − 2.9 a , rM

(38)

with E = 8 per cent, 36 per cent, 8 per cent, indicating quite a good ﬁt. There is an essentially 100 per cent conﬁdence for the sign of all parameters, indicating the presence of well-deﬁned trends. Solving the equation s = 0 for a/M shows that the spin equilibrium value, aeq/M, is given by

1 − aeq ≈ −0.08 + 0.8 h . Mr

(39)

In the other case, we ﬁt j /(2e) and re-solve for aeq/M, which gives directly

1 − aeq ≈ −0.10 + 0.9 h , Mr

(40)

with E = 9 per cent, 38 per cent with a 99.99 per cent conﬁdence in the sign of the |h/r| term. Both of these procedures give a similar ﬁt (the ﬁrst ﬁt is statistically better) and agree within statistical errors, which indicates a linear ﬁt is reasonable. For either ﬁt, one should set aeq/M = 1 when the above formula gives aeq/M > 1 to be consistent with our statistical errors and the correct physics. Note that the overshoot aeq/M > 1 in the ﬁt is consistent with a linear extrapolation of the NT dependence of s for a/M > 0, which also overshoots in the same way due to the progressively non-linear behaviour of s above a/M ≈ 0.95.
These spin equilibrium ﬁts imply that, within our statistical errors, the spin can reach aeq/M → 1 as |h/r| → 0. Thus, our results are consistent with NT by allowing maximal BH spin for thin discs.11 Our results are also roughly consistent with the thick disc one-loop ﬁeld geometry study by Gammie et al. (2004). Using our deﬁnition of disc thickness, their model had |h/r| ∼ 0.2–0.25 and they found aeq/M ∼ 0.9, which is roughly consistent with our scaling law. The ﬁt is also consistent with results for even thicker discs (|h/r| ∼ 0.4 near the horizon) with aeq/M ∼ 0.8 (Abramowicz, Jaroszynski & Sikora 1978; Popham & Gammie 1998).
Overall, the precise scaling relations given for ϒ and aeq should be considered as suggestive and preliminary. More work is required to test the convergence and generality of the actual coefﬁcients. While we explicitly tested convergence for the a/M = 0 ﬁducial model, the other a/M were not tested as rigorously. A potential issue is that we ﬁnd the saturated state has fewer cells per (vertical magnetic ﬁeld) fastest growing mode for the axisymmetric MRI in models with a/M = 0.9, 0.98 than in models with a/M = 0, 0.7 due to the relative weakness of the vertical ﬁeld in the saturated state for the high spin models. However, both the rough analytical arguments and the numerical solutions imply that electromagnetic stresses scale somewhat linearly with BH spin. This consistency suggests that many measurements for the simulations, such as ϒ and aeq, may be independent of smallness of the vertical ﬁeld. This fact could be due to these quantities being only directly related to the radial and toroidal magnetic ﬁeld strengths rather than the vertical

11 Here, we do not include BH spin changes by photon capture, which gives a limit of aeq/M = 0.998 (Thorne 1974).

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

magnetic ﬁeld strength. Further, our thick disc models resolve the axisymmetric MRI better than the thinnest disc model. This suggests that the scaling of ϒ and aeq with disc thickness may be a robust result.
Lastly, we consider how the speciﬁc angular momentum, nominal efﬁciency and luminosity from within the ISCO deviate from NT as functions of spin and thickness. Overall, ﬁtting these quantities does not give very strong constraints on the actual parameter values, but we can state the conﬁdence level of any trends. For each of e˜, j, jin and L˜ in, the deviation from NT as |h/r| → 0 is less than 5 per cent with a conﬁdence of 95 per cent. For j integrated over ±2|h/r|, D[j ] decreases with |h/r| and increases with a/M both with 99 per cent conﬁdence. When integrating j over all angles, D[j ] only decreases with |h/r| to 99 per cent conﬁdence. For j in integrated over ±2|h/r|, D[j in] only increases with a/M with 99.8 per cent conﬁdence and only decreases with |h/r| with 97 per cent conﬁdence. When integrating j in over all angles, D[j in] only increases with a/M to essentially 100 per cent conﬁdence and only decreases with |h/r| to 99.8 per cent conﬁdence. For e˜ integrated over ±2|h/r|, D[e˜] only increases with |h/r| with 98 per cent conﬁdence with no signiﬁcant trend with a/M. When integrating e˜ over all angles, D[e˜] only increases with a/M with 95 per cent conﬁdence with no signiﬁcant trend with |h/r|. For L˜ in, there is a 98 per cent conﬁdence for this to increase with |h/r| with no signiﬁcant trend with a/M. Overall, the most certain statement that can be made is that our results are strongly consistent with all deviations from NT becoming less than a few per cent in the limit that |h/r| → 0 across the full range of BH spins.
8 T H I N D I S C S W I T H VA RY I N G M AG N E T I C FIELD GEOMETRY
We now consider the effects of varying the initial ﬁeld geometry. Since the magnetic ﬁeld can develop large-scale structures that do not act like a local scalar viscosity, there could in principle be longlasting effects on the accretion ﬂow properties as a result of the initial ﬁeld geometry. This is especially a concern for geometrically thin discs, where the one-loop ﬁeld geometry corresponds to a severely squashed and highly organized ﬁeld loop bundle with long-range coherence in the radial direction, whereas our ﬁducial four-loop model corresponds to nearly circular loops which impose much less radial order on the MRI-driven turbulence. To investigate this question we have simulated a model similar to our ﬁducial run except that we initialized the gas torus with a one-loop type ﬁeld geometry instead of our usual multi-loop geometry.
Fig. 21 shows the radial dependence of j, jin, e˜ and ϒ for the two ﬁeld geometries under consideration, and Table 4 reports numerical estimates of various quantities at the horizon. Consider ﬁrst the solid lines (four-loop ﬁducial run) and dotted lines (oneloop run) in Fig. 21, both of which correspond to integrations in θ over ±2|h/r| around the midplane. The simulation with four-loops is clearly more consistent with NT than the one-loop simulation. The value of j at the horizon in the four-loop model deviates from NT by −2.9 per cent. Between the times of 12900M and 17300M, the one-loop model deviates by −5.6 per cent, while at late time over the saturated period the one-loop model deviates by −7.2 per cent. The long-dashed lines show the effect of integrating over all θ for the one-loop model. This introduces yet another systematic deviation from NT (as already noted in Section 5.7); now the net deviation of j becomes −10.7 per cent for times 12900M to 17300M and becomes −15.8 per cent for the saturated state. Overall, this implies that the assumed initial ﬁeld geometry has a considerable

Simulations of magnetized discs 773
Figure 21. Radial proﬁles of j and j in (top panel), e˜ (middle panel) and ϒ (bottom panel) are shown for two different initial ﬁeld geometries. Results for the ﬁducial four-loop ﬁeld geometry (model A0HR07) integrated over ±2|h/r| around the midplane are shown by solid lines, for the one-loop ﬁeld geometry (model A0HR07LOOP1) integrated over ±2|h/r| around the midplane by dotted lines and the one-loop model integrated over all angles by long-dashed lines. The short-dashed lines in the top two panels show the NT result. We see that the one-loop ﬁeld geometry shows larger deviations from NT in j and ϒ compared to the four-loop geometry. The panels also reemphasize the point that including all θ angles in the angular integration leads to considerable changes in j and ϒ due to the presence of magnetic ﬁeld in the corona-wind-jet.
impact on the speciﬁc angular momentum proﬁle and the stress inside the ISCO. This also indicates that the saturated state is only reached after approximately 17000M, and it is possible that the oneloop model may never properly converge because magnetic ﬂux of the same sign (how much ﬂux is initially available is arbitrary due to the arbitrary position of the initial gas pressure maximum) may continue to accrete on to the BH and lead to a qualitatively different accretion state [as seen in Igumenshchev, Narayan & Abramowicz (2003a) and McKinney & Gammie (2004) for their vertical ﬁeld model]. At early times, the nominal efﬁciency, e˜, shows no signiﬁcant dependence on the ﬁeld geometry, and remains near the NT value for both models. At late time in the one-loop model, e˜ rises somewhat, which may indicate the start of the formation of a qualitatively different accretion regime.
Fig. 22 shows the normalized luminosity. We see that the oneloop model produces more luminosity inside the ISCO. For times 12900M to 17300M, L˜ in = 5.4 per cent (integrated over all θ ) compared to 3.5 per cent for the four-loop ﬁeld geometry. Thus there is 50 per cent more radiation coming from inside the ISCO in this model. At late time during the saturated state, L˜ in = 4.6 per cent (integrated over all θ ). Thus there is approximately 31 per cent more radiation coming from inside the ISCO in this model during the late phase of accretion.
Table 4 also reports the results for thick (|h/r| ≈ 0.3) disc models initialized with the multi-loop and one-loop ﬁeld geometries. This again shows that the deviations from NT are inﬂuenced by the initial magnetic ﬁeld geometry and scale with |h/r| in a way expected by our scaling laws. The one-loop models show deviations from NT in

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

774 R. F. Penna et al.
Figure 22. Similar to Fig. 21 for the initial four-loop and one-loop ﬁeld geometries, but here we show the luminosity (top panel) and log-derivative of the luminosity (bottom panel). The luminosity is slightly higher for the one-loop model compared to the four-loop model. j are larger as related to the larger value of ϒ. The deviations from NT are less affected by the initial magnetic ﬁeld geometry for thicker discs, because the deviations from NT are also driven by thermal effects and Reynolds stresses rather than primarily electromagnetic stresses as for thin discs.
These effects can be partially understood by looking at the speciﬁc electromagnetic stress, ϒ, shown in Fig. 21. We ﬁnd ϒ ≈ 0.45 for the four-loop ﬁeld geometry. For times 12900M to 17300M, ϒ ≈ 0.71 in the one-loop ﬁeld geometry, and during the saturated state ϒ ≈ 1.28. For times 12900M to 17300M, the 50 per cent larger ϒ appears to be the reason for the 50 per cent extra luminosity inside the ISCO in the one-loop model. The magnetized thin disc model of Gammie (1999) predicts that, for a/M = 0, speciﬁc magnetic ﬂuxes of ϒ = 0.45, 0.71 should give deviations from NT of −D[j ] = 1.9, 3.9, respectively. These are close to the deviations seen in the simulations, but they are not a perfect match for reasons we can explain. First, the details of how one spatially averages quantities (e.g. average of ratio versus ratio of averages) when computing ϒ lead to moderate changes in its value, and, for integrations outside the midplane, comparisons to the Gammie model can require slightly

higher ϒ than our diagnostic reports. Secondly, the ﬁnite thermal pressure at the ISCO leads to (on average over time) a deviation already at the ISCO that is non-negligible compared to the deviation introduced by electromagnetic stresses between the ISCO and horizon. This thermal component is not always important, e.g. see the comparison in Fig. 11. Still, as found in McKinney & Gammie (2004) for thick discs at least, the deviations from NT contributed by the thermal pressure are of the same order as the deviations predicted by the Gammie model. These results motivate extending the Gammie (1999) model to include a ﬁnite (but still small) thermal pressure such that the boundary conditions at the ISCO lead to a non-zero radial velocity.
Within the ISCO, we ﬁnd that the time-averaged and volumeaveraged comoving ﬁeld strength for the four-loop geometry roughly follows |b| ∝ r−0.7 within ±2|h/r| of the disc midplane, while at higher latitudes we have a slightly steeper scaling. For times 12900M to 17300M, the one-loop geometry has |b| ∝ r−1.1 within ±2|h/r| of the disc midplane, and again a slightly steeper scaling in the corona. Other than this scaling, there are no qualitative differences in the distribution of any comoving ﬁeld component with height above the disc. While the Gammie (1999) solution does not predict a power-law dependence for |b|, for a range between ϒ = 0.4 and 0.8, the variation near the horizon is approximately |b| ∝ r−0.7 − r−0.9, which is roughly consistent with the simulation results. The slightly steeper slope we obtain for the one-loop geometry is consistent with a higher speciﬁc magnetic ﬂux, although the variations in ϒ for integration over different ranges of angle imply stratiﬁcation and a non-radial ﬂow which the Gammie (1999) model does not account for. This fact and the rise in ϒ with decreasing radius seen in Fig. 21 indicate a non-trivial degree of angular compression as the ﬂow moves towards the horizon. Overall, our results suggest that deviations from NT depend on the assumed ﬁeld geometry, and that the Gammie (1999) model roughly ﬁts the simulations.
Fig. 23 shows the same type of plot as in Fig. 7, but here we compare the ﬁducial four-loop model with the one-loop model. As mentioned above, the one-loop geometry has a larger deviation in j from the NT value, corresponding to larger stresses inside the ISCO. The absolute magnetic ﬂux (per unit initial total absolute magnetic ﬂux) on the BH ˜ r is of the order of 1/2, suggesting that the inner half of the initial ﬁeld bundle accreted on to the BH, while the other half was advected to larger radii. This is consistent with what is seen in simulations of thick tori (McKinney & Gammie 2004; Beckwith et al. 2008a). This suggests that using the one-loop geometry leads to results that are sensitive to the initial absolute magnetic ﬂux, while the multiple-loop geometry leads to results that are insensitive to the initial absolute magnetic ﬂux. Such dependence of the

Table 4. Field geometry study.

Model name
±2|h/r| A0HR07 A0HR07LOOP1 A0HR3 A0HR3LOOP1
All θ A0HR07 A0HR07LOOP1 A0HR3 A0HR3LOOP1

|h/r|
0.064 0.048 0.350 0.377
0.064 0.048 0.350 0.377

M˙
0.066 0.036 44.066 32.577
0.074 0.040 49.207 39.382

e˜
0.058 0.066 −0.003 0.001
0.054 0.075 −0.004 −0.001

D[e˜]
−0.829 −14.846 104.582
98.892
4.723 −31.857 106.180 102.499

j
3.363 3.215 2.309 1.823
3.266 2.915 2.128 1.575

D[j ]
−2.913 −7.193 −33.331 −47.389
−5.717 −15.847 −38.572 −54.526

j in
3.355 3.234 2.408 1.984
3.281 2.928 2.220 1.734

D[j in]
−3.153 −6.637 −30.473 −42.717
−5.275 −15.478 −35.919 −49.932

s
3.363 3.215 2.309 1.823
3.266 2.915 2.128 1.575

L˜ in
0.035 0.049 0.000 0.000
0.035 0.046 0.000 0.000

L˜ eq
0.080 0.059 −0.049 −0.049
0.053 0.048 −0.049 −0.049

100 ˜ r
1.355 6.198 3.039 7.599
6.677 43.935
4.724 11.444

ϒ
0.450 1.281 1.182 2.246
0.863 3.464 1.331 2.523

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

Simulations of magnetized discs 775

Figure 23. Similar to Fig. 7, but here we compare the initial four-loop ﬁducial model (black solid lines) and the one-loop model (dashed magenta lines). The horizontal black dashed lines in the second and third panels show the predictions of the NT model. The mass accretion rate, M˙ , has larger root-mean-squared ﬂuctuations in the one-loop model, which is indicative of more vigorous turbulence. The nominal efﬁciency, e˜, shows no clear difference. The speciﬁc angular momentum, j , is lower in the one-loop model compared to the four-loop model. This indicates that the one-loop ﬁeld leads to larger stress within the ISCO. The absolute magnetic ﬂux (per unit initial total absolute ﬂux) on the BH is larger in the one-loop model than in the four-loop model. Since ˜ r ∼ 1/2 for the one-loop model, essentially half of the initial loop was advected on to the BH, while the other half gained angular momentum and has been advected away. This may indicate that the one-loop geometry is a poor choice for the initial ﬁeld geometry, since the magnetic ﬂux that ends up on the BH is determined by the initial conditions. The value of ϒ is about twice higher in the one-loop model, which implies about twice greater electromagnetic stresses within the ISCO.
electromagnetic stress on initial magnetic ﬁeld geometry has also been reported in 3D pseudo-Newtonian simulations by Reynolds & Armitage (2001) and in 3D GRMHD simulations by Beckwith et al. (2008a).
Fig. 24 shows the electromagnetic stress as computed by equation (27) for the multiple-loop ﬁducial model (A0HR07) and the otherwise identical one-loop model (A0HR07LOOP1). We only show the electromagnetic part of the stress, and within the ISCO this is to within a few per cent the same as the total stress obtained by including all terms in the stress-energy tensor. Outside the ISCO, the total stress agrees more with the NT model. The ﬁgure shows the full-angle integrated electromagnetic stress, the electromagnetic stress integrated over only ±2|h/r|, the NT stress and the Gammie (1999) electromagnetic stress for ϒ = 0.60, 0.89, 0.90, 1.21 (we choose ϒ, the only free parameter of the model, such that the peak magnitude of the stress agrees with the simulation). The chosen ϒ values are close to our diagnostic’s value of ϒ for these models, which demonstrates that the Gammie (1999) model is consistently predicting the simulation’s results with a single free parameter. The stress is normalized by the radially dependent M˙ (r) that is computed over the same θ integration range. We do not restrict the integration to bound material as done in S08 (in S08, the stress is integrated over ±2|h/r| and only for bound material, while in N10 the

Figure 24. Normalized electromagnetic stress, W /M˙ , as a function of radius for the ﬁducial model (black lines) and the otherwise identical one-loop model (magenta lines). The solid lines correspond to a θ integration over all angles, while the dotted lines correspond to a θ integration over ±2|h/r|. The dashed black line shows the NT result, while the dashed green lines show the Gammie (1999) result for ϒ = 0.60, 0.89, 0.90, 1.21 for lines from the bottom to the top. The one-loop model shows about 50 per cent larger peak normalized stress (integrated over all angles) compared to the multi-loop model (integrated over all angles), which is consistent with the one-loop model leading to larger deviations from NT (about 50 per cent larger luminosity over the period used for time averaging). The large differences between the solid and dotted lines again highlight the fact that the stress within the disc is much smaller than the stress over all θ that includes the corona+wind+jet. As pointed out in S08, even though such a plot of the electromagnetic stress appears to indicate large deviations from NT within the ISCO, this is misleading because one has not speciﬁed the quantitative effect of the non-zero value of W /M˙ on physical quantities within the ISCO. Apparently, large values of W /M˙ do not necessarily correspond to large deviations from NT. For example, quantities such as j , e and the luminosity only deviate by a few per cent from NT for the multi-loop model. The Gammie (1999) model gives a reasonable ﬁt to the simulation’s stress proﬁle within the ISCO.
stress12 is only over bound material). The stress for the ﬁducial model was time-averaged over the saturated state, while the oneloop model was time-averaged from time 12900M to 17300M in order to best compare with the early phase of accretion for the one-loop model studied in N10.
Fig. 24 shows that the simulation and NT stress do not agree well, and it suggests there is an apparently large stress within the ISCO. However, as ﬁrst pointed out by S08 and discussed in Section 3.4, this stress does not actually correspond to a large deviation from NT in physically relevant quantities such as the speciﬁc angular momentum, speciﬁc energy and luminosity. This point is clariﬁed by making a comparison to the Gammie (1999) model’s stress, which agrees reasonably well with the simulation stress inside the ISCO. Even though the stress may appear large inside the ISCO, the stress corresponding to the Gammie model with (e.g.) ϒ =
12 N10’s ﬁgs 12 and 13 show stress versus radius, but some of the integrals they computed were not re-normalized to the full 2π when using their simulation φ-box size of π /2, so their stress curves are all a constant factor of four times larger than the actual stress (Noble, private communication).

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

776 R. F. Penna et al.
0.60 only translates into a few per cent deviations from NT. This ﬁgure also demonstrates that the initial magnetic ﬁeld geometry affects the amplitude of the stress in the same direction as it affects other quantities and is reasonably well predicted by the Gammie (1999) model. The initial magnetic ﬁeld sets the saturated value of ϒ, which is directly related to the electromagnetic stresses within the ISCO. The one-loop model leads to a peak stress (integrated over all angles) within the ISCO that is about 50 per cent larger than the multi-loop model (integrated over all angles), which is likely related to the extra 50 per cent luminosity in the one-loop model compared to the multi-loop model. The fact that the stress normalization changes with initial ﬁeld geometry is consistent with other 3D GRMHD simulations of thick discs by Beckwith et al. (2008a). This ﬁgure again shows how the stress within the disc (±2|h/r|) is much smaller than the total disc+corona+wind+jet (all θ ).
Finally, we discuss previous results obtained for other ﬁeld geometries using an energy-conserving 2D GRMHD code (McKinney & Gammie 2004). While such 2D simulations are unable to sustain turbulence, the period over which the simulations do show turbulence agrees quite well with the corresponding period in 3D simulations. This implies that the turbulent period in the 2D simulations may be qualitatively correct. The ﬁducial model of McKinney & Gammie (2004) was of a thick (|h/r| ∼ 0.2–0.25) disc with a one-loop initial ﬁeld geometry around an a/M = 0.9375 BH. This model had ϒ ∼ 1 near the midplane within the ISCO and ϒ ∼ 2 when integrated over all θ angles. Their measured value of j ≈ 1.46 integrated over all angles, |b| ∝ r−1.3 within the ISCO within the disc midplane (McKinney & Narayan 2007b), along with ϒ ∼ 1–2 are roughly consistent with the Gammie (1999) model prediction of j ≈ 1.5. Similarly, the strong vertical ﬁeld geometry model they studied had ϒ ∼ 2 near the midplane within the ISCO and ϒ ∼ 6 integrated over all θ angles. Their measurement of j ≈ −1 integrated over all angles is again roughly consistent with the model prediction of j ≈ −1.2 for ϒ ∼ 6. Note that in this model, ϒ rises (as usual to reach saturation) with time, but soon after ϒ 2 in the midplane, the disc is pushed away by the BH and then ϒ is forced to be even larger. Evidently, the accumulated magnetic ﬂux near the BH pushes the system into a force-free magnetosphere state – not an accretion state. This shows the potential danger of using strongﬁeld initial conditions (like the one-loop ﬁeld geometry), since the results are sensitive to the assumed initial ﬂux that is placed on (or rapidly drops on to) the BH. Even while the disc is present, this particular model exhibits net angular momentum extraction from the BH. This interesting result needs to be conﬁrmed using 3D simulations of both thick and thin discs.
9 COMPARISONS WITH OTHER RESULTS
The results we have obtained in the present work are consistent with those of Armitage, Reynolds & Chiang (2001) and Reynolds & Fabian (2008), who carried out pseudo-Newtonian studies, and with the results of S08, who did a full GRMHD simulation. These studies found only minor deviations from NT for thin accretion discs with a multi-loop initial ﬁeld geometry. However, more recently, N09 and N10 report apparently inconsistent results, including factors of up to 5 larger deviations from NT in the speciﬁc angular momentum (2 per cent in S08 versus 10 per cent in N10) for the same disc thickness of |h/r| ∼ 0.07. They also found a 50 per cent larger deviation from NT in the luminosity (4 per cent in S08 versus 6 per cent in N09). Furthermore, in N10 they concluded that the electromagnetic stresses have no dependence on disc thick-

ness or initial magnetic ﬁeld geometry, whereas we ﬁnd that the electromagnetic stresses have a statistically signiﬁcant dependence on both disc thickness and magnetic ﬁeld geometry.
We have considered several possible explanations for these differences, as we now describe. We attempt to be exhaustive in our comparison with the setup and results by N09 and N10, because our works and their works seek accuracies much better than order two in measuring deviations from NT. Thus, any deviations between our results by factors of 2 or more must be investigated further in order to ensure a properly understood and accurate result.
First, we brieﬂy mention some explanations that N10 propose as possible reasons for the discrepant results, viz., differences in (1) numerical algorithm or resolution; (2) box size in φ-direction:
φ; (3) amplitude of initial perturbations; (4) accuracy of inﬂow equilibrium and (5) duration of the simulations. Our algorithms are similar except that their PPM interpolation scheme assumes primitive quantities are cell averages (standard PPM), while ours assumes they are point values (as required to be applied in a higherorder scheme). They used LAXF dissipative ﬂuxes, while we used HLL ﬂuxes that are about twice more accurate for shocks and may be more accurate in general. On the other hand, they used parabolic interpolation for the Toth electric ﬁeld, while we use the standard Toth scheme. Given these facts, we expect that the accuracy of our algorithms is similar. Overall, our convergence testing and other diagnostics (see Section 6) conﬁrm that none of their proposed issues can be the cause of differences between S08 and N10.
We have shown that inﬂow equilibrium must include saturation of the speciﬁc magnetic ﬂux, ϒ, which generally saturates later in time than other quantities. By running our ﬁducial model A0HR07 to a time of nearly 30000M, we ensure that we have a long period of steady-state conditions to compute our diagnostic quantities. The fact that we need to run our ﬁducial thin disc simulation for such a long time to reach inﬂow equilibrium up to a radius r ∼ 9M is completely consistent with our analytical estimate of the time-scale calculated using equation (B7) of Appendix B (see the earlier discussion in Section 3.5 and Fig. 6). In the comparison between the numerical and analytical results shown in Fig. 6, we found agreement by setting α|h/r|2 ≈ 0.00033 which, for our disc with |h/r| ≈ 0.064, corresponds to α ≈ 0.08. With this value of α|h/r|2, we would have to run the simulation until t ∼ 83000M, 160000M, to reach inﬂow equilibrium out to 15M, 20M, respectively, corresponding to a couple viscous time-scales at each radius. N10 stated that they reach inﬂow equilibrium within r ∼ 15M– 20M in a time of only t ∼ 10000M. Since their disc thickness is |h/r| ≈ 0.06, even a single viscous time-scale would require their simulations to have α ∼ 0.38 to reach an inﬂow equilibrium up to r ∼ 15M, and an even larger value of α for r ∼ 20M. This seems unlikely. We can partially account for their result by considering our one-loop model, which up to t ∼ 17000M has α|h/r|2 twice as large and α about 70 per cent larger than in the ﬁducial four-loop run. However, this still falls far short by a factor of roughly 3 of what N10 would require for inﬂow equilibrium up to 15M– 20M. Further, our A0HR07LOOP1 model, which is similar to their model, only reaches a saturated state by 17000M, and only the ϒ quantity indicates that saturation has been reached. If we were to measure quantities from 10000M to 15000M as in N10, we would have underestimated the importance of magnetic ﬁeld geometry on the electromagnetic stresses.
Since all these simulations are attempting to obtain accuracies better than the factors of 2 in the results, this inﬂow equilibrium issue should be explored further. A few possible resolutions include: (1) N10’s higher resolution leads to a much larger α; (2) their disc

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

has a larger ‘effective’ thickness, e.g. |h/r| ∼ 0.13, according to equation (5.9.8) in NT (see equation B7 of Appendix B); (3) some aspects of their solution have not yet reached inﬂow equilibrium within a radius much less than r ∼ 15M, such as the value of ϒ versus time that saturates much later than other quantities; or (4) they achieve constant ﬂuxes versus radius due to transient nonviscous effects – although one should be concerned that the actual value of such ﬂuxes continues to secularly evolve in time and one still requires evolution on the longer viscous (turbulent diffusion) time-scale to reach true inﬂow equilibrium.
Secondly, we considered various physical setup issues, including differences in: (1) range of BH spins considered; (2) range of disc thicknesses studied; (3) ad hoc cooling function and (4) equation of state. We span the range of BH spins and disc thicknesses studied by N10, so this is unlikely to explain any of the differences. Some differences could be due to the disc thickness versus radius proﬁles established by the ad hoc cooling functions in the two studies. N10’s cooling function is temperature-based and allows cooling even in the absence of any dissipation, while ours is based upon the speciﬁc entropy and cools the gas only when there is dissipation. Both models avoid cooling unbound gas. In S08 and in the present paper, we use an ideal gas equation of state with = 4/3, while N09 and N10 used = 5/3. The properties of the turbulence do appear to depend on the equation of state (Mignone & McKinney 2007), so it is important to investigate further the role of in thin discs.
Thirdly, the assumed initial ﬁeld geometry may introduce critical differences in the results. Issues with the initial ﬁeld geometry include how many ﬁeld reversals are present, how isotropic the ﬁeld loops are in the initial disc, how the electromagnetic ﬁeld energy is distributed in the initial disc and how the magnetic ﬁeld is normalized. In S08 and here, we have used a multi-loop geometry in the initial torus consisting of alternating polarity poloidal ﬁeld bundles stacked radially. We ensure that the ﬁeld loops are roughly isotropic within the initial torus. We set the ratio of maximum gas pressure to maximum magnetic pressure to βmaxes = 100, which gives us a volume-averaged mean β within the dense part of the torus (ρ0/ρ0,max ≥ 0.2) of β¯ ∼ 800. Our procedure ensures that all initial local values of β within the disc are much larger than the values in the evolved disc, i.e. there is plenty of room for the magnetic ﬁeld to be ampliﬁed by the MRI.
We have also studied a one-loop geometry that is similar to the one-loop geometry used in N09 and N10. Their initial φ-component of the vector potential is Aφ ∝ MAX(ρ0/ρ0,max − 0.25, 0) (Noble, private communication). They initialize the magnetic ﬁeld geometry by ensuring that the volume-averaged gas pressure divided by the volume-averaged magnetic pressure is βaverages = 100 (Noble, private communication). (They stated that the mean initial plasma β is β¯ = 100.) For their thin disc torus model parameters, this normalization procedure leads to a portion of the inner radial region of the torus to have a local value of β∼ 3–8, which may be a source of differences in our results since such a small β is lower than present in the saturated disc. N10 make use of an older set of simulations from a different non-energy-conserving code (Hawley & Krolik 2006; Beckwith et al. 2008a) to investigate the effect of other ﬁeld geometries. The results from this other code have strong outliers, e.g. the KD0c model, and so we are unsure if these other simulations can be used for such a study.
N10 state that they ﬁnd no clear differences in the electromagnetic stresses for different initial ﬁeld geometries. As shown in their ﬁgs 12 and 13, the Agol & Krolik (2000) model captures the smoothing of the stress outside the ISCO, but it is not a model for the

Simulations of magnetized discs 777
behaviour of the stress inside the ISCO. We ﬁnd that electromagnetic stresses have a clear dependence on both disc thickness and the initial magnetic ﬁeld geometry, with a trend that agrees with the Gammie (1999) model of a magnetized thin disc. Our Fig. 24 shows that the stress within the ISCO is reasonably well modelled by the Gammie (1999) model. Our one-loop thin disc model gives a peak normalized stress (integrated over all angles) of about 3.2 × 10−3 for times 12900M to 17300M, which is comparable to the one-loop thin disc model in N10 with peak normalized stress (integrated over all angles) of about 2.5 × 10−3 (after correcting for their φ-box size). Hence, we are able to account for the results of their one-loop model.
In addition, we used the speciﬁc magnetic ﬂux, ϒ, an ideal MHD invariant that is conserved within the ISCO, to identify how electromagnetic stresses scale with disc thickness and magnetic ﬁeld geometry. In the saturated state, the value of ϒ, which controls the electromagnetic stresses, is different for different initial magnetic ﬁeld geometries. The larger the value of ϒ, the larger the deviations from NT in j . We ﬁnd that j within the disc (±2|h/r| from the midplane) deviates from NT by −3 per cent in our four-loop model and −6 per cent in our one-loop model for times 12900M to 17300M. Integrating over all angles, j deviates by −6 per cent for the four-loop model and −11 per cent for the one-loop model for times 12900M to 17300M. Thus, we ﬁnd a clear factor of 2 change, depending on the assumed initial ﬁeld geometry and the range of integration. The excess luminosity is 3.5 per cent for the four-loop model and 5.4 per cent for the one-loop model for times 12900M to 17300M. Recalling that N10 ﬁnd a deviation from NT of about −10 per cent in j (integrated over all angles) and a luminosity excess beyond NT of about 6%, this shows we can completely account for the apparent inconsistencies mentioned by N10 by invoking dependence of the results on the initial ﬁeld geometry and the presence of extra stress beyond the disk component of the accretion ﬂow.
Fourth, let us consider measurement and interpretation differences. Our ultimate goal is to test how well the NT model describes a magnetized thin accretion disc. The primary quantity that is used to measure this effect in S08 and N10 is the speciﬁc angular momentum j . However, the measurements are done differently in the two studies. In S08 as well as in this paper, we focus on the disc gas by limiting the range of θ over which we compute the averaging integrals (±2|h/r| from the midplane). In contrast, N10 compute their integrals over a much wider range of θ which includes both the disc and the corona-wind-jet (Noble, 2010, private communications). We have shown in Section 5.7 that the disc and corona-wind-jet contribute roughly equally to deviations of j from the NT value. In principle, the luminosity from the corona-wind-jet could be important, but we have shown that the excess luminosity of bound gas within the ISCO is dominated by the disc. This means that the measure used by N10, consisting of integrating over all gas to obtain j , cannot be used to infer the excess luminosity of bound gas within the ISCO. Further, the corona would largely emit non-thermal radiation, so for applications in which one is primarily interested in the thermal component of the emitted radiation, one should evaluate the accuracy of the NT model by restricting the angular integration range to the disc component within ±2|h/r|.
Fifth, let us consider how the results from N10 scale with disc thickness for the speciﬁc case of a non-spinning (a/M = 0) BH. We have performed a linear least-squares ﬁt of their simulation results, omitting model KD0c which is a strong outlier. For j integrated over all θ , their relative difference follows D[j ] ≈ −7 − 45|h/r| with conﬁdence of 95 per cent that these coefﬁcients, respectively, only deviate by ±67 per cent and ±89 per cent. These ﬁts imply that, as

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

778 R. F. Penna et al.
|h/r| → 0, the relative deviation of j from the NT value is about −7 per cent, but they could easily be as low as −2 per cent. Their results do not indicate a statistically signiﬁcant large deviation from NT as |h/r| → 0. Since the total deviation in j from NT includes the effects of electromagnetic (and all other) stresses, this implies that their models are consistent with weak electromagnetic stresses as |h/r| → 0.
Further, we have already established that the one-loop geometry gives (at least) twice the deviation from NT compared to the four-loop geometry, plus there is another factor of 2 arising from including the corona-wind-jet versus not including it. This net factor of 4 applied to N10’s results implies that j would deviate by about −2 per cent or even as low as −0.5 per cent from NT in the limit |h/r| → 0 if they were to consider a four-loop ﬁeld geometry and focus only on the disc gas. Thus, their models show no statistically signiﬁcant large deviations from NT. In addition, our results in Section 7.1 show that, whether we consider an integral over all angles or only over the disc, there is no statistically signiﬁcant large deviation from NT as |h/r| → 0.
In summary, we conclude that the apparent differences between the results obtained in S08 and the present paper on the one hand, and those reported in N09 and N10 on the other, are due to (1) dependence on initial magnetic ﬁeld geometry (multi-loop versus one-loop); (2) dependence upon the initial magnetic ﬁeld distribution and normalization and (3) measurement and interpretation differences (disc versus corona-wind-jet). Note in particular that the one-loop initial ﬁeld geometry is severely squashed in the vertical direction and elongated in the radial direction for thin discs, and it is not clear that such a geometry would ever arise naturally. There are also indications from our simulation that the one-loop geometry may actually never reach a converged state due to the arbitrary amount of magnetic ﬂux accreted on to the BH due to the single polarity of the initial magnetic ﬁeld. Finally, if one is trying to test how well-simulated thin accretion discs compare with NT, then it is important to restrict the comparison to disc gas near the midplane. One should not expect the gas in the corona-wind-jet to agree with the NT model.
10 DISCUSSION
We now discuss some important consequences of our results and also consider issues to be addressed by future calculations. First, we discuss the relevance to BH spin measurements.
In recent years, BH spin parameters have been measured in several BH X-ray binaries by ﬁtting their X-ray continuum spectra in the thermal (or high-soft) spectral state (Zhang, Cui & Chen 1997; Shafee et al. 2006; McClintock et al. 2006; Davis, Done & Blaes 2006; Liu et al. 2008; Gou et al. 2009, 2010). This method is based on several assumptions that require testing (Narayan, McClintock & Shafee 2008), the most critical being the assumption that an accretion disc in the radiatively efﬁcient thermal state is well-described by the Novikov–Thorne model of a thin disc. More speciﬁcally, in analysing and ﬁtting the spectral data, it is assumed that the radial proﬁle of the radiative ﬂux, or equivalently the effective temperature, in the accretion disc closely follows the prediction of the NT model.
Practitioners of the continuum-ﬁtting method generally restrict their attention to relatively low-luminosity systems below 30 per cent of the Eddington luminosity. At these luminosities, the maximum height of the disc photosphere above the midplane is less than 10 per cent of the radius, i.e. (h/r)photosphere ≤ 0.1

(McClintock et al. 2006). For a typical disc, the photospheric disc thickness is approximately twice the mean absolute thickness |h/r| that we consider in this paper. Therefore, the discs that observers consider for spin measurement have |h/r| 0.05, i.e. they are thinner than the thinnest disc (|h/r|min ∼ 0.06) that we (S08, this paper) and others (N09, N10) have simulated.
The critical question then is the following. Do the ﬂux proﬁles of very thin discs match the NT prediction? At large radii the two will obviously match very well since the ﬂux proﬁle is determined simply by energy conservation.13 However, in the region near and inside the ISCO, analytic models have to apply a boundary condition, and the calculated ﬂux proﬁle in the inner region of the disc depends on this choice. The conventional choice is a ‘zero-torque’ boundary condition at the ISCO. Unfortunately, there is disagreement on the validity of this assumption. Some authors have argued that the magnetic ﬁeld strongly modiﬁes the zero-torque condition and that, therefore, real discs might behave very differently from the NT model near the ISCO (Krolik 1999; Gammie 1999). Other authors, based either on heuristic arguments or on hydrodynamic calculations, ﬁnd that the NT model is accurate even near the ISCO so long as the disc is geometrically thin (Paczyn´ski 2000; Afshordi & Paczyn´ski 2003; Abramowicz et al. 2010; S08). Investigating this question was the primary motivation behind the present study.
We described in this paper GRMHD simulations of geometrically thin (|h/r| ∼ 0.07) accretion discs around BHs with a range of spins: a/M = 0, 0.7, 0.9, 0.98. In all cases, we ﬁnd that the speciﬁc angular momentum j of the accreted gas as measured at the horizon (this quantity provides information on the dynamical importance of torques at the ISCO) shows only minor deviations at the level of ∼2–4 per cent from the NT model. Similarly, the luminosity emitted inside the ISCO is only ∼3–7 per cent of the total disc luminosity. When we allow for the fact that a large fraction of this radiation will probably be lost into the BH because of relativistic beaming as the gas plunges inward (an effect ignored in our luminosity estimates), we conclude that the region inside the ISCO is likely to be quite unimportant. Furthermore, our investigations indicate that deviations from the NT model decrease with decreasing |h/r|. Therefore, since the discs of interest to observers are generally thinner than the thinnest discs we have simulated, the NT model appears to be an excellent approximation for modelling the spectra of BH discs in the thermal state.
One caveat needs to be mentioned. Whether or not the total luminosity of the disc agrees with the NT model is not important since, in spectral modelling of data, one invariably ﬁts a normalization (e.g. the accretion rate M˙ in the model KERRBB; Li et al. 2005) which absorbs any deviations in this quantity. What is important is the shape of the ﬂux proﬁle versus radius. In particular, one is interested in the radius at which the ﬂux or effective temperature is maximum (Narayan et al. 2008; McClintock et al. 2009). Qualitatively, one imagines that the fractional shift in this radius will be of the order of the fractional torque at the ISCO, which is likely to be of the order of the fractional error in j . We thus guess that, in the systems of interest, the shift is nearly always below 10 per cent. We plan to explore this question quantitatively in a future study.
Another issue is the role of the initial magnetic ﬁeld topology. We ﬁnd that, for a/M = 0, starting with a one-loop ﬁeld geometry gives an absolute relative deviation in j of 7.1 per cent, and
13 This is why the formula for the ﬂux as a function of radius in the standard thin disc model does not depend on details like the viscosity parameter α (Frank et al. 1992).

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

an excess luminosity inside the ISCO of 4.9 per cent, compared to 2.9 per cent and 3.5 per cent for our standard four-loop geometry. Thus, having a magnetic ﬁeld distribution with long-range correlation in the radial direction seems to increase deviations from the NT model, though even the larger effects we ﬁnd in this case are probably not a serious concern for BH spin measurement. Two comments are in order on this issue. First, the four-loop geometry is more consistent with nearly isotropic turbulence in the poloidal plane and, therefore, in our view a more natural initial condition. Secondly, the one-loop model develops a stronger ﬁeld inside the ISCO and around the BH and might therefore be expected to produce a relativistic jet with measurable radio emission. However, it is well-known that BH X-ray binaries in the thermal state have no detectable radio emission. This suggests that the magnetic ﬁeld is probably weak, i.e. more consistent with our four-loop geometry.
Next, we discuss the role of electromagnetic stresses on the dynamics of the gas in the plunging region inside the ISCO. In order to better understand this issue, we have extracted for each of our simulations the radial proﬁle of the speciﬁc magnetic ﬂux, ϒ. This quantity appears as a dimensionless free parameter (called Fθφ) in the simple MHD model of the plunging region developed by Gammie (1999). The virtues of the speciﬁc magnetic ﬂux are its well-deﬁned normalization and its constancy with radius for stationary ﬂows (Takahashi et al. 1990). In contrast, quantities like the stress W or the viscosity parameter α have no well-deﬁned normalization; W can be normalized by any quantity that has an energy scale, such as ρ0, M˙ or b2, while α could be deﬁned with respect to the total pressure, the gas pressure or the magnetic pressure. The numerical values of W or α inside the ISCO can thus vary widely, depending on which deﬁnition one chooses. For instance, although S08 found α ∼ 1 inside the ISCO, the speciﬁc angular momentum ﬂux, j , deviated from NT by no more than a few per cent. Further, Fig. 24 shows that (even for the multi-loop model) the stress appears quite large within the ISCO, but this is misleading because the effects of the stress are manifested in the speciﬁc angular momentum, speciﬁc energy and luminosity – all of which agree with NT to within less than 10 per cent for the multi-loop model. Since W and α do not have a single value within the ISCO or a unique normalization, we conclude that they are not useful for readily quantifying the effects of the electromagnetic stresses within the ISCO.
Gammie’s (1999) model shows how the value of ϒ is directly related to the electromagnetic stresses within the ISCO. Unfortunately, the actual value of ϒ is a free parameter which cannot be easily determined from ﬁrst principles. It is possible that accretion discs might have ϒ 1, in which case, the model predicts large deviations from NT. For example, if ϒ = 6, then for an a/M = 0 BH j is lowered by 56 per cent relative to the NT model. We have used our 3D GRMHD simulations which include self-consistent MRI-driven turbulence to determine the value of ϒ for various BH spins, disc thicknesses and ﬁeld geometries. For the multiple-loop ﬁeld geometry, we ﬁnd that the speciﬁc magnetic ﬂux varies with disc thickness and spin as

ϒ ≈ 0.7 + h − 0.6 a , rM

(41)

within the disc component, which indicates that electromagnetic stresses are weak and cause less than 8 per cent deviations in j in the limit |h/r| → 0 for all BH spins. Our rough analytical arguments for how ϒ should scale with |h/r| and a/M are consistent with the above formula. Even for the one-loop ﬁeld geometry, ϒ 1 for thin discs, so electromagnetic stresses cause only minor deviations from NT for all BH spins (for ϒ 1, less than 12 per cent in j ). Not all

Simulations of magnetized discs 779
aspects of the Gammie (1999) model agree with our simulations. As found in McKinney & Gammie (2004), the nominal efﬁciency, e˜, does not match well and for thin discs is quite close to NT. Since the true radiative efﬁciency is limited to no more than e˜, this predicts only weak deviations from NT in the total luminosity even if j has non-negligible deviations from NT. Also, this highlights that the deviations from NT in j are due to non-dissipated electromagnetic stresses and cannot be used to directly predict the excess luminosity within the ISCO. The assumption of a radial ﬂow in a split-monopole ﬁeld is approximately valid, but the simulations do show some nonradial ﬂow and vertical stratiﬁcation, a non-zero radial velocity at the ISCO and thermal energy densities comparable to magnetic energy densities. Inclusion of these effects is required for better consistency with simulation results inside the ISCO.
Next, we consider how our results lend some insight into the spin evolution of BHs. Standard thin disc theory with photon capture predicts that an accreting BH spins up until it reaches spin equilibrium at aeq/M ≈ 0.998 (Thorne 1974). On the other hand, thick nonradiative accretion ﬂows deviate signiﬁcantly from NT and reach equilibrium at aeq/M ∼ 0.8 for a model with α ∼ 0.3 and |h/r| ∼ 0.4 near the horizon (Popham & Gammie 1998). GRMHD simulations of moderately thick (|h/r| ∼ 0.2–0.25) magnetized accretion ﬂows give aeq/M ≈ 0.9 (Gammie et al. 2004). In this paper, we ﬁnd from our multi-loop ﬁeld geometry models that spin equilibrium scales as

aeq ≈ 1.1 − 0.8 h , Mr

(42)

where one should set aeq/M = 1 if the above formula gives aeq/M > 1. This gives a result consistent with the above-mentioned studies of thick discs, and it is also consistent with our rough analytical prediction based upon our scaling of ϒ and using the Gammie model prediction for the spin equilibrium. This result also agrees with the NT result in the limit |h/r| → 0 within our statistical errors, and shows that magnetized thin discs can approach the theoretical limit of aeq/M ≈ 1, at least in the multi-loop case. In the singleloop ﬁeld geometry, because of the presence of a more radially elongated initial poloidal ﬁeld, we ﬁnd a slightly stronger torque on the BH. However, before a time of order 17000M, the deviations in the equilibrium spin parameter, aeq/M, between the four-loop and one-loop ﬁeld geometries appear to be minor, so during that time the scaling given above roughly holds. Of course, it is possible (even likely) that radically different ﬁeld geometries or anomalously large initial ﬁeld strengths will lead to a different scaling.
Lastly, we mention a number of issues which we have neglected but are potentially important. A tilt between the angular momentum vector of the disc and the BH rotation axis might signiﬁcantly affect the accretion ﬂow (Fragile et al. 2007). We have not accounted for any radiative transfer physics, nor have we attempted to trace photon trajectories (see e.g. N09 and Noble & Krolik 2009). In principle, one may require the simulation to be evolved for hundreds of orbital times at a given radius in order to completely erase the initial conditions (Sorathia, Reynolds & Armitage 2010), whereas we only run the model for a couple of viscous time-scales. New pseudo-Newtonian simulations show that convergence may require resolving several disc scaleheights with high resolution (Sorathia et al. 2010), and a similar result has been found also for shearing box calculations with no net ﬂux and no stratiﬁcation (Stone 2010, private communication). In contrast, we resolve only a couple of scaleheights. Also, we have only studied two different types of initial

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

780 R. F. Penna et al.
ﬁeld geometries. Future studies should consider whether alternative ﬁeld geometries change our results.
11 CONCLUSIONS
We set out in this study to test the standard model of thin accretion discs around rotating BHs as developed by NT. We studied a range of disc thicknesses and BH spins and found that magnetized discs are consistent with NT to within a few per cent when the disc thickness |h/r| 0.07. In addition, we noted that deviations from the NT model decrease as |h/r| goes down. These results suggest that BH spin measurements via the X-ray continuum-ﬁtting method (Zhang et al. 1997; Shafee et al. 2006; McClintock et al. 2006; Davis et al. 2006; Liu et al. 2008; Gou et al. 2009, 2010), which are based on the NT model, are robust to model uncertainties so long as |h/r| 0.07. At luminosities below 30 per cent of Eddington, we estimate disc thicknesses to be |h/r| 0.05, so the NT model is perfectly adequate.
These results were obtained by performing global 3D GRMHD simulations of accreting BHs with a variety of disc thicknesses, BH spins and initial magnetic ﬁeld geometries in order to test how these affect the accretion disc structure, angular momentum transport, luminosity and the saturated magnetic ﬁeld. We explicitly tested the convergence of our numerical models by considering a range of resolutions, box sizes and amplitude of initial perturbations. As with all numerical studies, future calculations should continue to clarify what aspects of such simulations are converged by performing more parameter space studies and running the simulations at much higher resolutions. For example, it is possible that models with different BH spins require more or less resolution than the a = 0 models, while we ﬁxed the resolution for all models and only tested convergence for the a = 0 models.
We conﬁrmed previous results by S08 for a non-spinning (a/M = 0) BH, which showed that thin (|h/r| 0.07) discs initialized with multiple poloidal ﬁeld loops agree well with the NT solution once they reach steady state. For the ﬁducial model described in the present paper, which has similar parameters as the S08 model, we ﬁnd 2.9 per cent relative deviation in the speciﬁc angular momentum accreted through the disc, and 3.5 per cent excess luminosity from inside the ISCO. Across all BH spins that we have considered, viz., a/M = 0, 0.7, 0.9, 0.98, the relative deviation from NT in the speciﬁc angular momentum is less than 4.5 per cent, and the luminosity from inside the ISCO is less than 7 per cent (typically smaller, and much of it is likely lost to the hole). In addition, all deviations from NT appear to be roughly proportional to |h/r|.
We found that the assumed initial ﬁeld geometry modiﬁes the accretion ﬂow. We investigated this effect by considering two different ﬁeld geometries and quantiﬁed it by measuring the speciﬁc magnetic ﬂux, ϒ, which is an ideal MHD invariant (like the speciﬁc angular momentum or speciﬁc energy). The speciﬁc magnetic ﬂux can be written as a dimensionless free parameter that enters the magnetized thin disc model of Gammie (1999). This parameter determines how much the ﬂow deviates from NT as a result of electromagnetic stresses. We found that ϒ allows a quantitative understanding of the ﬂow within the ISCO, while the electromagnetic stress (W) has no well-deﬁned normalization and varies widely within the ISCO. While a plot of the stress may appear to show large stresses within the ISCO, the actual deviations from NT can be small. This demonstrates that simply plotting W is not a useful diagnostic for measuring deviations from NT. We found that the

speciﬁc magnetic ﬂux of the gas inside the ISCO was substantially larger when we used a single poloidal magnetic loop to initialize the simulation compared to our ﬁducial four-loop run. For a/M = 0 and |h/r| 0.07, the early saturated phase (times 12900M to 17300M) of the evolution for the one-loop geometry gave 5.6 per cent relative deviation in the speciﬁc angular momentum and 5.8 per cent excess luminosity inside the ISCO. These deviations are approximately twice as large as the ones we found for the four-loop simulation. At late times, the one-loop model generates signiﬁcant deviations from NT, which is a result similar to that found in a vertical ﬁeld model in McKinney & Gammie (2004). However, we argued that the multiple-loop geometry we used is more natural than the single-loop geometry, since for a geometrically thin disc the magnetic ﬁeld in the one-loop model is severely squashed vertically and highly elongated radially. The one-loop model is also likely to produce a strong radio jet.
More signiﬁcant deviations from NT probably occur for discs with strong ordered magnetic ﬁeld, as found in 2D GRMHD simulations by McKinney & Gammie (2004). Of course, in the limit that the magnetic ﬁeld energy density near the BH exceeds the local rest-mass density, a force-free magnetosphere will develop and deviations from the NT model will become extreme. We argued that this corresponds to when the speciﬁc magnetic ﬂux ϒ 1 near the disc midplane. Our one-loop model appears to be entering such a phase at late time after accumulation of a signiﬁcant amount of magnetic ﬂux. Such situations likely produce powerful jets that are not observed in BH X-ray binaries in the thermal state. However, transition between the thermal state and other states with a strong power-law component (Fender et al. 2004; Remillard & McClintock 2006) may be partially controlled by the accumulation of magnetic ﬂux causing the disc midplane (or perhaps just the corona) to breach the ϒ ∼ 1 barrier. Such a behaviour has been studied in the non-relativistic regime (Narayan, Igumenshchev & Abramowicz 2003; Igumenshchev, Narayan & Abramowicz 2003b; Igumenshchev 2009), but more work using GRMHD simulations is required to validate the behaviour.
We also found that the apparently different results obtained by N10 were mostly due to measurement and interpretation differences. We found that both the disc and the corona-wind-jet contribute nearly equally to deviations in the total speciﬁc angular momentum relative to the NT model. However, the corona-wind-jet contributes much less to the luminosity than the disc component. Therefore, if one is interested in comparing the luminous portion of the disc in the simulations against the NT model, the only fair procedure is to consider only the disc gas, i.e. gas within a couple of scaleheights of the midplane. This is the approach we took in this study (also in S08). N10 on the other hand included the corona-wind-jet gas in their calculation of the speciﬁc angular momentum. The dynamics of the coronal gas differs considerably from the NT model. Therefore, while it does not contribute to the luminosity of bound gas, it doubles the deviation of the speciﬁc angular momentum from the NT model. In addition, N10 used a one-loop initial ﬁeld geometry for their work which, as discussed above, further enhanced deviations.
ACKNOWLEDGMENTS
We thank Phil Armitage, the referee, for comments that greatly improved the paper’s presentation, and we thank Scott Noble for detailed discussions about his work. We thank Charles Gammie, Chris Reynolds, Jim Stone, Kris Beckwith, John Hawley, Julian Krolik, Chris Done, Chris Fragile, Martin Pessah and Niayesh Afshordi

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

for useful discussions. This work was supported in part by NASA grant NNX08AH32G (AT and RN), NSF grant AST-0805832 (AT and RN), NASA Chandra Fellowship PF7-80048 (JCM), an NSF Graduate Research Fellowship (RFP) and the NSF through TeraGrid resources provided by NCSA (Abe), LONI (QueenBee), NICS (Kraken) under grant numbers TG-AST080025N and TGAST080026N.
REFERENCES
Abramowicz M., Jaroszynski M., Sikora M., 1978, A&A, 63, 221
Abramowicz M. A., Jaroszynski M., Kato S., Lasota J.-P., Rozanska A., Sadowski A., 2010, A&A, in press (arXiv:1003.3887)
Afshordi N., Paczyn´ski B., 2003, ApJ, 592, 354 Agol E., Krolik J. H., 2000, ApJ, 528, 161 Armitage P. J., Reynolds C. S., Chiang J., 2001, ApJ, 548, 868 Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214 Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 1 Bardeen J. M., 1970, ApJ, 161, 103 Beckwith K., Hawley J. F., Krolik J. H., 2008a, ApJ, 678, 1180 Beckwith K., Hawley J. F., Krolik J. H., 2008b, MNRAS, 390, 21 Bekenstein J. D., Oron E., 1978, Phys. Rev. D, 18, 1809 Berti E., Volonteri M., 2008, ApJ, 684, 822 Blandford R. D., 1976, MNRAS, 176, 465 Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433 Chakrabarti S. K., 1985a, ApJ, 288, 1 Chakrabarti S. K., 1985b, ApJ, 294, 383 Colella P., Woodward P. R., 1984, J. Comp. Phys., 54, 174 Davis S. W., Done C., Blaes O. M., 2006, ApJ, 647, 525 De Villiers J.-P., Hawley J. F., Krolik J. H., 2003, ApJ, 599, 1238 Fender R. P., Belloni T. M., Gallo E., 2004, MNRAS, 355, 1105 Font J. A., Iba´n˜ez J. M., Papadopoulos P., 1998, ApJ, 507, L67 Fragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007, ApJ, 668,
417 Frank J., King A., Raine D., 1992, Accretion power in astrophysics. Cam-
bridge Univ. Press, Cambridge, UK Gammie C. F., 1999, ApJ, 522, L57 Gammie C. F., McKinney J. C., To´th G., 2003, ApJ, 589, 444 Gammie C. F., Shapiro S. L., McKinney J. C., 2004, ApJ, 602, 312 Gou L. et al., 2009, ApJ, 701, 1076 Gou L., McClintock J. E., Steiner J. F., Narayan R., Cantrell A. G., Bailyn
C. D., Orosz J. A., 2010, ApJ, 718, L122 Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103 Hughes S. A., Blandford R. D., 2003, ApJ, 585, L101 Igumenshchev I. V., 2009, ApJ, 702, L72 Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003a, ApJ, 592, 1042 Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003b, ApJ, 592, 1042 Komissarov S. S., Barkov M. V., 2009, MNRAS, 397, 1153 Komissarov S. S., McKinney J. C., 2007, MNRAS, 377, L49 Komissarov S. S., Barkov M. V., Vlahakis N., Ko¨nigl A., 2007, MNRAS,
380, 51 Krolik J. H., 1999, ApJ, 515, L73 Krolik J. H., Hawley J. F., 2002, ApJ, 573, 754 Krolik J. H., Hawley J. F., Hirose S., 2005, ApJ, 622, 1008 Li L., 2002, Phys. Rev. D, 65, 084047 Li L., Zimmerman E. R., Narayan R., McClintock J. E., 2005, ApJS, 157,
335 Liu J., McClintock J. E., Narayan R., Davis S. W., Orosz J. A., 2008, ApJ,
679, L37 MacDonald D. A., 1984, MNRAS, 211, 313 McClintock J. E., Shafee R., Narayan R., Remillard R. A., Davis S. W., Li
L., 2006, ApJ, 652, 518 McClintock J. E., Narayan R., Gou L., Liu J., Penna R. F., Steiner J. F., 2009,
in Comastri A., Angelini L., Cappi M., eds, AIP Conf. Proc. Vol. 1248, X-ray Astronomy 2009: Present Status, Multi-Wavelength Approach and Future Perspectives. Am. Inst. Phys., New York, p. 101

Simulations of magnetized discs 781
McKinney J. C., 2005, ApJ, 630, L5 McKinney J. C., 2006a, MNRAS, 367, 1797 McKinney J. C., 2006b, MNRAS, 368, 1561 McKinney J. C., Blandford R. D., 2009, MNRAS, 394, L126 McKinney J. C., Gammie C. F., 2004, ApJ, 611, 977 McKinney J. C., Narayan R., 2007a, MNRAS, 375, 513 McKinney J. C., Narayan R., 2007b, MNRAS, 375, 531 Meier D. L., 2001, ApJ, 548, L9 Mignone A., McKinney J. C., 2007, MNRAS, 378, 1118 Narayan R., Kato S., Honma F., 1997, ApJ, 476, 49 Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ, 55,
L69 Narayan R., McClintock J. E., Shafee R., 2008, in Yuan Y.-F., Li X.-D., Lai
D., eds, AIP Conf. Ser. Vol. 968, Astrophysics of Compact Objects. Am. Inst. Phys., New York, p. 265 Noble S. C., Krolik J. H., 2009, ApJ, 703, 964 Noble S. C., Gammie C. F., McKinney J. C., Del Zanna L., 2006, ApJ, 641, 626 Noble S. C., Krolik J. H., Hawley J. F., 2009, ApJ, 692, 411 (N09) Noble S. C., Krolik J. H., Hawley J. F., 2010, ApJ, 711, 959 (N10) Novikov I. D., Thorne K. S., 1973, in DeWitt C., DeWitt B. S., eds, Black Holes. Les astres occlus. Gordon and Breach, New York, p. 343 (NT) Paczyn´ski B., 2000, preprint (astro-ph/0004129) Page D. N., Thorne K. S., 1974, ApJ, 191, 499 Papadopoulos P., Font J. A., 1998, Phys. Rev. D, 58, 024005 Popham R., Gammie C. F., 1998, ApJ, 504, 419 Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49 Reynolds C. S., Armitage P. J., 2001, ApJ, 561, L81 Reynolds C. S., Fabian A. C., 2008, ApJ, 675, 1048 Sano T., Inutsuka S., Turner N. J., Stone J. M., 2004, ApJ, 605, 321 Sa¸dowski A., 2009, ApJS, 183, 171 Schnittman J. D., Krolik J. H., Hawley J. F., 2006, ApJ, 651, 1031 Shafee R., McClintock J. E., Narayan R., Davis S. W., Li L.-X., Remillard R. A., 2006, ApJ, 636, L113 Shafee R., McKinney J. C., Narayan R., Tchekhovskoy A., Gammie C. F., McClintock J. E., 2008a, ApJ, 687, L25 (S08) Shafee R., Narayan R., McClintock J. E., 2008b, ApJ, 676, 549 Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337 Shapiro S. L., 2010, Phys. Rev. D, 81, 024019 Shapiro S. L., Teukolsky S. A., 1983, Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects. Research supported by the National Science Foundation. Wiley-Interscience, New York, p. 663 Sorathia K. A., Reynolds C. S., Armitage P. J., 2010, ApJ, 712, 1241 Takahashi M., Nitta S., Tatematsu Y., Tomimatsu A., 1990, ApJ, 363, 206 Tchekhovskoy A., McKinney J. C., Narayan R., 2007, MNRAS, 379, 469 Tchekhovskoy A., McKinney J. C., Narayan R., 2008, MNRAS, 388, 551 Tchekhovskoy A., McKinney J. C., Narayan R., 2009a, ApJ, 699, 1789 Tchekhovskoy A., Narayan R., McKinney J. C., 2009b, New Astron., 15, 749 Tchekhovskoy A., Narayan R., McKinney J. C., 2010, ApJ, 711, 50 Thorne K. S., 1974, ApJ, 191, 507 Zhang S. N., Cui W., Chen W., 1997, ApJ, 482, L155
APPENDIX A: EXAMPLE SOLUTIONS AND SCALINGS FOR THE GAMMIE (1999) MODEL
Table A1 gives representative solutions for the Gammie (1999) model of a magnetized thin accretion ﬂow. The columns correspond to the BH spin, a; the speciﬁc magnetic ﬂux, ϒ; the nominal efﬁciency, e˜; per cent deviation of e˜ from the NT value; the speciﬁc angular momentum, j ; per cent deviation of j from NT; and the normalized rate of change of the dimensionless BH spin, s (see equation 15). For ϒ 0.5 and across all BH spins, the relative change in the speciﬁc angular momentum is less than 5 per cent and the relative change in the efﬁciency is less than 9 per cent. For small values of ϒ 1, the deviations of j and e˜ from NT behave

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782

782 R. F. Penna et al.

Table A1. Thin magnetized inﬂow solutions.

a M

ϒ e˜

D[e˜] j

D[j ] s

0 0.1 0.0576 0.709 3.46

−0.172 3.46

0 0.2 0.0584 2.14 3.45

−0.52 3.45

0 0.3 0.0595 4.08 3.43

−0.991 3.43

0 0.5 0.0624 9.17 3.39

−2.23 3.39

0 1 0.0727 27.1 3.24

−6.57 3.24

0 1.5 0.0859 50.2 3.04

−12.2 3.04

0 6 0.19 232 1.51

−56.4 1.51

0.7 0.1 0.105 1.03 2.58

−0.286 1.33

0.7 0.2 0.107 3.07 2.56

−0.853 1.31

0.7 0.3 0.11 5.8

2.54

−1.61 1.3

0.7 0.5 0.117 12.8 2.49

−3.57 1.26

0.7 1 0.142 36.7 2.32

−10.2 1.12

0.7 1.5 0.172 66.3 2.11

−18.5 0.95

0.7 6 0.477 360 −0.00714 −100 −0.74

0.9 0.1 0.157 1.17 2.09

−0.386 0.576

0.9 0.2 0.161 3.37 2.08

−1.11 0.567

0.9 0.3 0.165 6.29 2.06

−2.07 0.555

0.9 0.5 0.177 13.7 2.01

−4.5

0.524

0.9 1 0.215 38.3 1.84

−12.6 0.423

0.9 1.5 0.262 68.3 1.63

−22.5 0.3

0.9 6 0.845 443 −0.958 −146 −1.24

0.98 0.1 0.236 0.949 1.68

−0.397 0.179

0.98 0.2 0.241 2.86 1.66

−1.2

0.174

0.98 0.3 0.246 5.36 1.64

−2.25 0.168

0.98 0.5 0.261 11.6 1.6

−4.9

0.152

0.98 1 0.309 32.2 1.45

−13.6 0.1

0.98 1.5 0.368 57.1 1.28

−24.1 0.0379

0.98 6 1.21 416 −1.27

−175 −0.862

0.998 0.1 0.319 −0.63 1.4

0.344 0.0374

0.998 0.2 0.327 2.02 1.38

−1.11 0.0342

0.998 0.3 0.332 3.66 1.36

−2 0.0322

0.998 0.5 0.345 7.73 1.33

−4.22 0.0273

0.998 1 0.388 20.9 1.23

−11.4 0.0113

0.998 1.5 0.439 37

1.11

−20.2 −0.00819

0.998 6 1.19 272 −0.675 −148 −0.292

systematically and one can derive simple ﬁtting functions. For j we ﬁnd

log10 [−D[j ]] ≈ 0.79 + 0.37(a/M) + 1.60 log10 ϒ

(A1)

∼ (4/5) + (1/3)(a/M) + (8/5) log10 ϒ,

(A2)

with an L2 error norm of 0.7 per cent, 0.7 per cent, respectively, for the ﬁrst and second relations, while for e˜ we ﬁnd

log10 [D[e˜]] ≈ 1.44 + 0.12(a/M) + 1.60 log10 ϒ

(A3)

∼ (3/2) + (1/10)(a/M) + (8/5) log10 ϒ,

(A4)

with an L2 error norm of 0.9 per cent, 1 per cent, respectively, for the ﬁrst and second relations. These results indicate that the deviations

from NT scale as ϒ8/5 for ϒ 1. For ϒ 1, the index on ϒ depends on the spin parameter. In the span from ϒ ∼ 0.2 to ϒ ∼ 1, a linear ﬁt across all BH spins gives −D[j ] ∼ −1 + 11ϒ and D[e˜] ∼ −4 + 33ϒ, which are rough, though reasonable looking, ﬁts.

APPENDIX B: INFLOW EQUILIBRIUM TIME-SCALE IN THE NOVIKOV–THORNE MODEL

The radius out to which inﬂow equilibrium is achieved in a given
time may be estimated by calculating the mean radial velocity vr and then deriving from it a viscous time-scale −r/vr.
When the ﬂow has achieved a steady state, the accretion rate,

M˙ = −2π r vrD1/2,

(B1)

is a constant independent of time and position. Here, we derive an expression for vr corresponding to the general relativistic NT thin disc model. In what follows, capital script letters denote standard functions of r and a (cf. equations 14 and 35 in Page & Thorne 1974) which appear as relativistic corrections in otherwise Newtonian expressions. They reduce to unity in the limit r/M → ∞.
The vertically integrated surface density may be deﬁned as = 2hρ, where h is the disc scaleheight and ρ is the rest-mass density at the midplane. In equilibrium, density is related to pressure by

dp = ρ × (‘acceleration of gravity’) dz

(B2)

=

ρ

Mz r3

B2DE A2C

,

the vertically integrated solution of which is

h = (p/ρ)1/2/| |AB−1C1/2D−1/2E −1/2.

(B3) (B4)

The pressure may be parametrized in terms of the viscous stress, |trˆφˆ | = αp, which is a known function of r and a:

M˙ C1/2Q W = 2htrˆφˆ = 2π BD . The surface density is then

(B5)

= 1 M˙ A2B−3C3/2D−2E −1Q. 2π αh2| |

(B6)

Substituting this in equation (B1), the radial velocity is

vr = −α|h/r|2| |rA−2B3C−3/2D3/2E Q−1.

(B7)

This result is independent of the exact form of the pressure and opacity and so is valid in all regions of the disc. The inﬂow equilibrium time may be estimated as tie ∼ −2r/vr.

This paper has been typeset from a TEX/LATEX ﬁle prepared by the author.

C 2010 The Authors. Journal compilation C 2010 RAS, MNRAS 408, 752–782