{'fraction_non_alphanumeric': 0.029080932784636488, 'fraction_numerical': 0.01196159122085048, 'mean_word_length': 4.588307440719542, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dependence of inner accretion disk stress on parameters: the Schwarzschild case

We explore the parameter dependence of inner disk stress in black hole accretion by contrasting the results of a number of simulations, all employing 3-d general relativistic MHD in a Schwarzschild spacetime. Five of these simulations were performed with the intrinsically conservative code HARM3D, which allows careful regulation of the disk aspect ratio, H/R; our simulations span a range in H/R from 0.06 to 0.17. We contrast these simulations with two previously reported simulations in a Schwarzschild spacetime in order to investigate possible dependence of the inner disk stress on magnetic topology. In all cases, much care was devoted to technical issues: ensuring adequate resolution and azimuthal extent, and averaging only over those time-periods when the accretion flow is in approximate inflow equilibrium. We find that the time-averaged radial-dependence of fluid-frame electromagnetic stress is almost completely independent of both disk thickness and poloidal magnetic topology. It rises smoothly inward at all radii (exhibiting no feature associated with the ISCO) until just outside the event horizon, where the stress plummets to zero. Reynolds stress can also be significant near the ISCO and in the plunging region; the magnitude of this stress, however, depends on both disk thickness and magnetic topology. The two stresses combine to make the net angular momentum accreted per unit rest-mass 7-15% less than the angular momentum of the ISCO.


Introduction
At the very beginning of accretion disk studies, their overall properties were analyzed by applying the constraints of energy and angular momentum conservation to the simplest reasonable approximation to their structure: they were assumed to be time-steady and axisymmetric, and any internal vertical structure was integrated over (Novikov & Thorne 1973;Shakura & Sunyaev 1973;Page & Thorne 1974). The equation of energy conservation can be closed by counting the energy carried by photons to infinity, but no such ready closure exists for the angular momentum equation; angular momentum can be conserved for any rate of angular momentum transport through the disk provided it does not vary with radius. Consequently, it was necessary to guess the angular momentum flux in order to complete the solution. A convenient way to parametrize this guess is in terms of the net angular momentum accreted onto the black hole per accreted rest-mass, j net . The choice made by the original papers, and still widely-used today, is to suppose that no stresses act on the flow from the innermost stable circular orbit (ISCO) inward to the event horizon; if so, j net = u φ (ISCO), the orbital angular momentum of a test-particle at the last stable orbit (here u µ is the covariant four-velocity).
However, this was never more than a heuristic guess. As remarked by Thorne (1974), although the zero-stress boundary condition is plausibly motivated by hydrodynamic reasoning-the inertia of matter inside the ISCO should always be much less than that in the stable-orbit portion of the disk outside the ISCO-it might well be invalid if magnetic fields are important. In fact, one of the crucial things we have learned in the years since the 1970s is that magnetic fields are, in fact, essential to accretion due to the presence of the magnetorotational instability (MRI) (Balbus & Hawley 1998).
On that basis, this traditional boundary condition has been questioned (Krolik 1999b;Gammie 1999), and a number of numerical simulations of global disks (Hawley & Krolik initial field was entirely vertical. Before presenting the results of these simulations, we will also discuss the importance of a number of technical considerations-particularly having to do with spatial resolution and the establishment of inflow equilibrium-to obtaining meaningful results.

Simulation Details
The new simulations reported here were made using the code HARM3D, an intrinsically conservative method to solve the equations of 3-d MHD in an arbitrary metric. This new code is described in Noble et al. (2009); see also Gammie et al. (2003) for a description of the earlier axisymmetric version HARM, and Noble et al. (2006) for additional details on the primitive variable solver. We employ the same methodology as before, with only a few exceptions. In the following summary of this code's techniques, we emphasize those points that are different from the previous description or particular to the simulations discussed in this paper.
One of the principal aims of the present work is to study the influence of disk thickness H on the stress at the ISCO. We define it as the density moment in the coordinate frame: where g µν is the metric, g is the determinant of the metric, and ρ is the rest-mass density. When the density profile follows a Gaussian distribution with standard deviation H G , H = 2/πH G = 0.798H G . As in Noble et al. (2009), we regulate the thickness by cooling bound portions of the disk when the local temperature is greater than some target temperature T * (r). In terms of intensive quantities, bound matter satisfies (ρ + u + P ) u t > −ρ and gas has temperature above the target when (Γ − 1) u/ρ > T * ; here, P is the gas pressure, u is the internal energy density, u µ is the fluid's 4-velocity, and Γ is the adiabatic index of the equation of state: P = (Γ − 1) u. The relativistic enthalpy h ≡ 1 + (u + P )/ρ. We set Γ = 5/3 throughout. The optically-thin cooling function is implemented by modifying the stress-energy conservation equation to include a sink term: The fluid-frame emissivity, L, and T * are designed so as to keep the density scale height at the desired value. The emissivity is the same as before but we slightly modified T * to include a neglected relativistic correction. The new target temperature is where R z is the relativistic correction to the vertical component of gravity (Abramowicz et al. 1997;Krolik 1999a) 2 : Here, l k and ǫ k are the specific angular momentum (u φ ) and energy (u t ) of circular time-like geodesics in the equator of a black hole with spin parameter a. For r < r ISCO , l k and ǫ k remain at their ISCO values.
All of the new simulations were performed in a Schwarzschild spacetime (a = 0) described in terms of Kerr-Schild coordinates and run for durations of 12000M to 15000M (in our units, G = c = 1, so that both time and distance have units M, the mass of the central black hole). In all cases, the initial condition was a hydrostatic torus, but we examined two varieties of this state: in one the radial coordinate of the pressure maximum r pmax = 35M and the inner edge r in = 20M; in the other, r pmax = 25M and r in = 15M. The former set of parameters were chosen to match those of Shafee et al. (2008a), the latter to match those of De . In the remainder of this paper, we will refer to the former set as "HR" and to the latter set as "LR" because the grid schemes used for the former were in general higher resolution than for the latter. The q parameter determines the angular velocity profile in the initial torus (Ω ∝ r −q ). The choice of q along with the choice of pressure maximum and inner torus edge determines the characteristic vertical thickness of the initial torus. For LR simulations, q was set to 1.68, while for HR, q was determined by requiring the initial disk's thickness at r = r pmax to be equal to the run's target thickness. To study the effect of different disk thicknesses, we ran HR simulations for three different target temperatures, chosen to make H/R ≃ 0.05, 0.08, and ≃ 0.16; these were designated "Thin", "Medium", and "Thick", respectively. We ran only LR simulations for the Thin and Medium cases. Finally, in all five simulations, the initial magnetic field consisted of dipole poloidal loops, with field lines following density contours, and with amplitude set such that the mean initial plasma β = 100. Turbulence was seeded by adding random perturbations to u at the 1% level.
Boundary conditions were imposed through assignment of primitive variables in ghost zones; the primitive variables are ρ, u, v i (spatial velocity components), and B i (magnetic field components). The v i are spatial components of the four-velocity as seen by observers in the ZAMO frame; the B i are the spatial components of the magnetic field as represented in the Maxwell field tensor, i.e., B i ≡ * F it / √ 4π. Outflow boundary conditions were taken both at r = r max and r = r min : all primitive variables are extrapolated at 0 th -order into the ghost zones, but u r is set to zero-and v i recalculated-whenever it points into the domain. In order to prevent numerical boundary effects from propagating outside the trapped surface, we chose r min so that the numerical domain extended 5 to 25 cells inside the event horizon.
When there was a cutout around the polar axis (see grid details below), reflective boundary conditions are imposed on the perpendicular vector components while all other quantities are extrapolated beyond the cutout with 0 th -order accuracy. If the cutout size is negligible, ghost zone values are set so as to make the variables continuous across the pole. In order to gain a factor of four in computer resources, we simulated only a quarter of the azimuthal domain, employing periodic boundary conditions linking φ = 0 and φ = π/2.
Adequate resolution throughout the accretion flow and throughout the duration of the simulation is vital to ensure quantitative accuracy, particularly for magnetic effects. To explore the consequences of resolution effects, we used two different grid schemes, one for HR, the other for LR. Both were designed with an eye toward satisfying several conflicting criteria. On the one hand, it is always desirable to have the finest feasible spatial resolution.
Toward that end, the grid scheme should provide at least several dozen cells per scale height on either side of the equatorial plane, and the poloidal cell aspect ratio should never be too large. In addition, the number of cells within a wavelength of the fastest-growing MRI mode should only occasionally fall below ≃ 6, the minimum resolution level at which the mode grows at the correct rate (Sano et al. 2004). On the other hand, to fit within existing computational resources, the total cell count cannot be too great nor the time step too small.
Each run used a unique discretization tailored to its particular target thickness.
HARM3D accomplishes nonuniform discretization through continuous coordinate transformations from an underlying uniform mesh. The finite volume equations are discretized with respect to points uniformly distributed in coordinates, x µ , which are placed nonuniformly in the spherical coordinate system r, θ, and φ. The center of cell C ijk is so that ∆r/r is the same everywhere. For N 1 grid cells along the x 1 axis and minimum and maximum radii r min and r max , ∆x 1 = 1 N 1 log rmax r min .
The relationship between x 2 and θ is determined differently in the HR and LR simulations. In the LR group, we follow previous work (e.g., Gammie et al. 2003;Noble et al. 2007Noble et al. , 2009), but introduce a "cutout" or excised region around the polar axis: where θ c is the approximate size of the cutout, θ 0 and s control the nonlinearity of the transformation, and ξ is the amplitude of the nonlinear part. One drawback to this transformation is that if one tries to place a majority of the points within the first two scale heights from the equator, the minimum ∆θ, which is always found at θ = π/2, is so small that the time step is prohibitively small.
In the HR simulations, we follow the method of Shafee et al. (2008a), adapted here to include a cutout. In this method, the sinusoidal nonlinear term is replaced with a polynomial: where n is a positive odd integer, θ c is the size of the excised region, and ξ is still the relative amplitude of the nonlinear term. Note that near the equator, where most of the points are located, the linear term dominates and ∆θ(x 2 ) is nearly uniform. Periodicity of θ ∈ [0, 2π] is ensured by making x 2 a periodic triangle function over [0,2]. We set n = 9 whenever equation (6) is used for the runs presented here.
Because gridscale dissipation scales with the ratio of cell dimension to the length scale on which physical quantities vary, cells that are far from cubical may have effective dissipation properties that are anisotropic. This might produce unphysical results because physical dissipation mechanisms are unlikely to have this property. For this reason, we strive to limit the degree of anisotropy in our cell shapes, although cells longer by factors of a few in the φ-coordinate than in the others are acceptable because orbital shear tends to draw out features in the azimuthal direction. At the same time, it is important, especially for small H/R simulations, to make the θ-direction cell thickness small enough to put an adequate number of cells within a vertical scale height. Achieving these conflicting goals is easier with the HR grid scheme than with the LR scheme (as illustrated by the data in Table 1. Our HR runs use grids with ∆r : r∆θ ≃ 2 : 1 while achieving at least 60 cells per scale height. In LR simulations, although the aspect ratios are acceptable at an altitude ∼ H away from the equatorial plane, they are rather extreme close to the plane, and-of these runs-only ThinLR has a comparable number of cells per vertical scale height to HR runs (see Table 2).
The time step ∆t = ∆x 0 is set equal to 0.8 times the shortest cell crossing time for the fastest MHD characteristic from the previous time step.
The defining parameters of the simulations are collected in Table 1. In the first column we state the names of the runs. The body of the name corresponds to the disk thickness, the suffix refers to the grid resolution. The target aspect ratio H/R is shown in the second column. Note that the "radius" to which the scale height is compared is the radial coordinate, which is not identical to a length in any frame of reference. Since we regulate the temperature, not H directly, the actual value will not in general exactly coincide with the target. The observed value of H/R is given in Table 2 and is obtained by averaging over the time interval shown in its last column (whose origin is discussed in § 3.2) and radially-averaging from the ISCO (r = 6M) to the run's r pmax . The third column, cell-count, is self-explanatory. Columns 4 to 8 define the poloidal discretization. Columns 9 and 10 state the typical cell aspect ratio at the equator and at two scale heights from the equator, and Column 11 lists the size of the smallest poloidal extent-∆θ-of a cell, which always occurs at θ = π/2.
We will also analyze two older simulations, both done with the general relativistic 3-d MHD code GRMHD (originally described in De ). Both began with the gas in a hydrostatic torus. The first, called "KD0c," was described in Hawley & Krolik (2006) and had initial conditions similar to the LR models presented here.
The other, called "VD0" and analyzed in Beckwith et al. (2009), began with a constant intensity vertical magnetic field filling the annulus 35M ≤ r ≤ 55M and running through an initial torus with a pressure maximum at r = 40M. Unlike the HARM3D simulations there was no explicit cooling function included. Rather, the internal energy equation was evolved and the only heating included was shock heating captured by an artificial viscosity term. The resulting H/R values are given in Table 2 and, in terms of the descriptive terms used here, qualify these simulations as "thick."

Quality and Consistency Checks
We are investigating the stress at the ISCO in a steady state accretion disk with a focus on the role of vertical scale height. To model the accretion disk system, we use the equations of conservative ideal MHD along with an ad hoc cooling function to control the disk's thickness. We can, however, run only a discrete set of simulations, evolved for a limited time from particular initial conditions, using a grid with restricted spatial extent and modest resolution. The conclusions we obtain must be assessed within the context of the limitations of those simulations.
In this section we consider the effects of some of these numerical limitations by developing several quantitative diagnostics to measure their possible significance. We begin by considering the adequacy of the grid resolution, and then look at how closely the inner disk approximates a statistical steady state. Finally, since this study focuses on the effect of disk scale height for the ISCO stress, we check to see how well we are able to control that variable.

Resolution
Inadequate resolution can cause a number of numerical artifacts. For example, the growth rate of the underlying MRI can be suppressed if there are fewer than ≃ 6 zones within a wavelength of the fastest growing mode (Sano et al. 2004). The MRI produces turbulence, but only a small range in wavenumber space can be captured, possibly distorting the properties of that turbulence. The rate at which nonlinear mode-mode couplings transfer energy from large-scale motions to small may be altered. And resolution that is too coarse may drive an artificially large rate of magnetic numerical dissipation.
Unfortunately, there is no a priori standard by which we can measure whether a given grid scheme either exhibits excessive magnetic reconnection or improperly evaluated nonlinear mode couplings. Only through a numerical convergence test, in which increasingly better-resolved simulations give quantitatively consistent results, can one demonstrate that resolution artifacts are not influencing the outcomes. As a practical matter, however, one carries out simulations such as these at the highest feasible resolution. Lower resolution simulations might provide some information, but with the computational resources at our disposal it has not been possible to improve sufficiently upon the resolution used in the production simulations to establish formal convergence.
It is possible, however, to check whether the grid resolution satisfies certain physical criteria, such as having a sizable number of cells per vertical scale height and sufficient cells across the fastest growing MRI wavelength that the linear growth of these modes is correctly described. Both criteria may be met in the initial conditions, but the data must be examined throughout the relevant volume (the main disk body) and for the duration of the simulation to ensure that they continued to be satisfied.
Data on the mean number of cells per scale height are displayed in To quantify the quality of resolution of the MRI, we evaluated the parameter Q ≡ λ MRI /∆z, where both λ MRI and ∆z are computed in the fluid frame, and whereê µ andê µ are contra-and co-variant tetrad systems in the local fluid-frame, respectively. To assure ourselves that the simulations always satisfied this criterion, we created animations of Q with frames every 20M in time. In Figure 1, we show sample stills from the simulations ThinHR, ThinLR, MediumHR, and ThickHR, each exhibiting a vertical slice at fixed azimuthal angle. All of these simulations exhibited better than adequate resolution at all times: the typical number of cells across the fastest-growing mode in the disk body was ∼ > 20. Azimuthally-averaged versions of the data shown in Figure 1 display ratios of fastest-growing wavelength to cell-size greater than 12 throughout the entire region within two scale heights of the equatorial plane. It is important to recognize, however, that even in a superbly-resolved simulation there will always be the occasional region in which the local poloidal field strength is small simply because this is a chaotic system in which fields are free to have either sign.
The importance of Q comes from the requirement that the MRI maintain turbulence in the face of continual dissipation. If the MRI growth rate is artificially reduced due to poor resolution, 3 and field amplification by MHD turbulence is too weak, the field strength can become increasingly anemic. Thus it is possible for a simulation that began with ample cells per MRI wavelength to suddenly find itself under-resolved. The end-result is a field that dies away and a cessation of accretion due entirely to inadequate resolution.
This appears to happen in the MediumLR simulation beginning at t ≃ 9000M.
MediumLR demonstrates the utility of the resolution diagnostic. Figure 2 shows two stills, at times 6000M and 12000M, from the resolution animation for this simulation. Initially well-resolved, MediumLR becomes more poorly resolved as the field strength diminishes.
As a result, the accretion rate also decays.
The two GRMHD simulations can be analyzed in similar fashion. By the standards of the HARM3D simulations, both have relatively few cells per scale height, only ≃ 12-16, not quite as many as in MediumLR. We also show sample snapshots of Q in Figure 3, in each case showing the final time-step of the simulation. KD0c was performed several years ago, and we did not save enough snapshots to create an animation; of all those available (every 80M from t = 8000M to t = 10000M), the one shown in Figure 3 appears to be the least well-resolved. In fact, there is an indication of a secular worsening of resolution beginning at t ≃ 9000M. In the case of VD0, we possess more data and have confirmed that at no time was the resolution substantially poorer than shown here. As can be seen, in neither case was the resolution as good as in the best of the HARM3D simulations, although VD0 at its end-point was better than KD0c at its.

Inflow equilibrium
In the end, we hope to use these evolving simulations to describe time-steady accretion flows. Starting from our initial condition, this state can never be achieved at all radii because the angular momentum removed from accreting material will be transferred to matter at larger radius, causing that matter to move outward. Moreover, because only a fraction of the initial torus mass is accreted in the duration of the simulations, the radial surface density profile can relax to the one associated with inflow equilibrium only within a short distance outside the initial inner radius of the torus (20M for HR, 15M for LR).
Nonetheless, it is possible to identify a period of time for which an approximate state of inflow equilibrium really does obtain over a reasonable dynamic range in radius, subject, of course, to the sorts of fluctuations that occur in statistically stationary turbulent systems.
To identify that region in both time and space, we impose several tests relying on the conservation of mass and angular momentum. Whenever the disk is in inflow equilibrium, the radial fluxes of mass and angular momentum should be constant as a function of radius, but because of turbulent fluctuations, they are constant only in a time-averaged sense. The radial range of the equilibrium region is determined by the range over which the time-averaged values of these quantities are nearly constant. To identify the time periods in which equilibrium obtains, we consider the mass interior to several specific radii M(< r; t) and the time-dependence of the specific angular momentum accreted onto the black hole, j net . In inflow equilibrium, these should be roughly constant over time.
The mass interior to radius r is defined as the integral of the mass density over the computational volume from the horizon to radial coordinate r, or Thus, once inflow equilibrium is established, the mass within a given radius should stay the same, as the amount of mass entering from outside is matched by the mass accreted onto the black hole at the center.
Figures 4-6 display the history of mass inside r = 10M, 15M, and 20M for each of the three disk aspect ratio categories studied. In all cases, at early times the mass in the disk grows as accretion from the initial torus fills in rings at smaller radii. Eventually, the mass-interior curves level off, signaling the approach to equilibrium with respect to this criterion. In the case of MediumLR, the mass in the inner disk declines at late times, a symptom of the diminution in the accretion rate which we attribute to the artificial decay of the magnetic field.
Similar data for the two GRMHD simulations is given in Figure 7. Both simulations reach inflow equilibrium with respect to this criterion, after ≃ 3000M in the case of KD0c, after ≃ 10000M in the case of VD0.
Another test of inflow equilibrium is provided by the history of the specific angular momentum accreted into the black hole, j net . This is determined by dividing the total angular momentum flux by the accretion rate: where the brackets represent the shell integration of the bracketed quantity: The quantity j net (r, t) will be constant in r where there is inflow equilibrium. If the accretion disk is truly in a statistically steady state in regard to its angular momentum flow, this quantity should exhibit no trends in time, varying only modestly due to fluctuations intrinsic to the turbulence. Figure  The time-dependence of j net in MediumLR is instructive. Overall, especially after 10000M, j net rises toward u φ (ISCO). As discussed in §3.1, the average number of grid zones per most unstable wavelength dropped in this simulation as the field weakened. Inadequate resolution leads to ever-weaker magnetic field which, in turn, results in a rise of the specific accreted angular momentum toward the ISCO value as the magnetic stress is reduced. Combining what we have seen in the mass-interior plots with those in the j net histories, we define the averaging periods, ∆t ave , for these simulations as the time when both criteria for inflow equilibrium are met. The results of this analysis are shown in the last column of is the contrast in Ṁ = ρu r inside r = 20M as much as 30%; in one (MediumLR) it is ≃ 10%; in all the others (ThinLR, MediumHR, ThickHR, KD0c, VD0), it is no more than a few percent.

Scale-height regulation
Finally we examine the actual time-averaged scale heights achieved in the various simulations; these are shown in Figure 10. The scale-height regulation employed in the HARM3D simulations is quite successful at enforcing a fixed ratio H/R (except in the plunging region in ThickHR), but, as shown in Table 2, the actual value obtained can be different from the target, ≃ 20% greater in the cases of ThinHR, MediumHR, and MediumLR, but ≃ 70% greater in the case of ThinLR. Part of this consistent offset can be attributed to magnetic support, and part to the fact that the temperature is typically slightly greater than the target. The large offset in ThinLR is a consequence of its initial condition, in which the gas was given a thickness almost twice as great as the target.
Although its initial mean plasma β was 100, when cooling compresses the disk by a sizable factor, the magnetic field strengthens while the gas pressure falls. As a result, much of the disk mass of ThinLR outside a thin midplane layer was supported magnetically, and its scale height was substantially increased beyond what its gas pressure alone could support.
This fact emphasizes the importance of choosing initial conditions relatively close to the expected time-averaged state. Inside the ISCO, the aspect ratio decreases slightly because the inflow time is so short that the disk cannot maintain vertical hydrostatic equilibrium.
The time-averaged scale heights of the two GRMHD simulations are shown in Figure 11.
This code evolves the internal energy equation, rather than the total energy equation, and its simulations did not employ an explicit cooling function. Rather, entropy is conserved except where local shock heating is captured by an artificial viscosity term. Thus the scale height is not controlled directly and is determined primarily by the initial condition. In the end, the two simulations have very similar scale-height profiles, with H/R ≃ 0.15 for r ≥ 10M, but declining inward, reaching ≃ 0.12 at the ISCO (r = 6M) and ≃ 0.06 just outside the horizon. The mean H/R for both is ≃ 0.14.

Results: Stress in the Inner Disk
With this background in mind, we can now discuss the results for time-averaged stress in the inner disk. We will present them in two ways: in terms of the radial profile of the spherical shell-integrated fluid-frame electromagnetic stress, and in terms of the angular momentum flux and the quantity j net defined in the introduction.

Fluid-frame electromagnetic stress profile
We begin with the radial profile of the electromagnetic stress. Both for the purpose of highlighting the physics and for the purpose of contrasting with the Novikov-Thorne model, it is best to compute it in the fluid frame: where b µ is the magnetic 4-vector. Each component of the vector dx (µ) = e (µ) ν dx ν represents the extent of a cell's dimension as measured in the fluid element's rest frame, and e (µ) ν is the orthonormal tetrad that transforms vectors in the Boyer-Lindquist coordinate frame to the local fluid frame (see Beckwith et al. (2008b) for explicit expressions for the tetrad).
The physical significance of the electromagnetic fluid-frame stress profile is that it describes the rate at which angular momentum is carried outward by electromagnetic fields.
The net angular momentum accreted by the black hole is diminished to the degree that stresses like these convey angular momentum outward even while inflowing matter carries its orbital angular momentum inward. The Novikov-Thorne model assumes that the stress begins to decline outside of the ISCO, reaching zero at that point. Agol & Krolik (2000) showed how changing that boundary condition to account for non-zero stress at the ISCO, i.e., j net < u φ (ISCO), could significantly alter the shape of the stress profile. Even in the disk body, well outside the ISCO, a smaller j net can lift the time-averaged stress above the Novikov-Thorne curve. applies, here r ≤ 20M. Outside the ISCO, but inside the domain of inflow equilibrium, the fluid frame stress usually (but not always) lies slightly above the Novikov-Thorne prediction. Although the match is not perfect, it is somewhat better described by the Agol-Krolik model. As the flow approaches the ISCO, where the Novikov-Thorne model would predict that the stress begins to fall, the electromagnetic stress rises steadily. Inside the ISCO, in the plunging region, the stress continues to rise inward, with a slope that is similar to or slightly steeper than outside the ISCO. Just outside the event horizon, the stress falls sharply to zero: because this stress component is nothing more than the radial flux of angular momentum of rotation in the disk plane, when the black hole has no angular momentum (i.e., does not rotate), it cannot act as a source of angular momentum and the stress immediately outside it must go to zero.
The corresponding profiles for the two GRMHD simulations are shown in Figure 13.
They are very similar to one another, and qualitatively similar to, but quantitatively different from, the HARM3D profiles. Like the HARM3D profiles, the radial slope of the stress is nearly constant in the disk outside the ISCO; unlike the HARM3D profiles, in these two the stress rises somewhat more steeply inside the ISCO. As a result, the peak stress in these two simulations is ≃ 2 times greater than seen in the HARM3D cases.
In conclusion, despite the range of temperatures and scale heights covered by these models, there is a great similarity between the stress profiles of all the HARM3D models as well as KD0c and VD0. In other words, with respect to this measure of the stress, there appears to be no dependence on H/R whatsoever.
Furthermore, the similarity between KD0c and VD0 suggests that the change of magnetic topology from closed dipolar loops on the scale of the disk thickness to net vertical field also makes little difference to the radial variation of accretion stress. One possible explanation for this insensitivity is the fact that reconnection events in the corona largely decoupled the magnetic field in the inner disk of VD0 from the large-scale flux (Beckwith et al. 2009). The remaining field in the disk then has a topology not so different from that in the other simulations. In §5 we review the results from previously published simulations to explore further the possible role of field topology on ISCO stress.

Specific accreted angular momentum
The value of j net , the mean angular momentum accreted per unit rest-mass, summarizes the net angular momentum flow in the system. It is determined by several effects. In the accretion disk body, the orbital angular momentum, u φ , is close to the value associated with a circular test-particle orbit at that radius, but can be altered by an amount ∼ (H/R) 2 by radial pressure gradients, both gas and magnetic. Stresses, both electromagnetic (Maxwell) and fluid (Reynolds) move angular momentum through the accretion flow; to the degree that they have a net divergence, they can either add or remove angular momentum from the fluid. In the classical Novikov-Thorne model, the fluid's angular momentum is assumed to match the local circular orbit angular momentum at all locations outside the ISCO, but is fixed at the ISCO angular momentum at all smaller radii. The stresses, whether Maxwell or Reynolds, are constrained to be exactly what they need to be to produce this result: finite in the disk body, zero in the plunging region. As a result, j net is predicted to be exactly u φ (ISCO).
As can be seen from the data listed in Table 2, this is not the case in the simulations. We find that, over the range of thicknesses and magnetic geometries studied, j net ranges from ≃ 2.93 to ≃ 3.21. Thus, in all these cases, the accreted angular momentum per unit accreted rest-mass is 7-15% below u φ (ISCO) = 3.464. Interestingly, the largest j net by far was seen in KD0c, a simulation performed on a comparatively coarse grid; if it were discounted, the depression of j net below u φ would be 10-15%.
The separate elements contributing to this departure from the Novikov-Thorne prediction are shown by the curves in Figures 14 and 15. The solid lines correspond to the net specific angular momentum flux; in inflow equilibrium this should be constant with radius. The dashed curves show the time-averaged specific angular momentum carried by the accreting matter, i.e., ρhu r u φ / ρu r . The difference between these curves and the net angular momentum flux is due to the electromagnetic stress. The dot-dash curves show the time-averaged mass-weighted mean angular momentum at each radius, ρhu φ / ρ . These last two quantities can be compared with the u φ corresponding to a circular orbit, shown by the dotted line (which is held constant at the ISCO value inside of that radius).
For all the models, the mass-weighted mean angular momentum generally follows the circular orbit value outside of the ISCO, but continues to decline inside of that radius rather than holding steady. In the thin disk models the offset from the circular orbit value is small, while in the thicker cases this offset is somewhat larger. This is precisely what one would expect; the gas in the disk is partially supported by the outward radial decline in pressure, primarily magnetic, in disks with a larger H/R.
In every case, the curve of mean angular momentum accreted by the matter lies well below the curve of the local mass-weighted mean angular momentum. One way to view the origin of this offset stems from the fact that the material in the disk is turbulent; local fluid elements have angular momenta that fluctuate from time to time and from place to place.
It should be no surprise that fluid elements with angular momentum slightly smaller than the value that would support a local circular orbit tend to move inward faster than those with larger angular momentum. In other words, the mean accreted angular momentum is systematically biased toward lower values by orbital mechanics that sorts the fluid elements according to their place in the local angular momentum distribution.
This effect can also be identified with turbulent Reynolds stress. To see this identification directly, consider the equation of angular momentum equilibrium integrated over spherical shells and averaged in time: where M r φ is the Maxwell stress. The term ∝ u r u φ can be broken into two pieces, one reflecting the advection of the mean angular momentum (weighted by enthalpy), the other reflecting departures from the mean. Dividing through by the mass accretion rate, the previous equation becomes where δX ≡ X − X . That is, the net rate at which angular momentum is carried inward per unit rest-mass accreted is the local mean angular momentum reduced by the ratio of the total stress, Maxwell plus turbulent Reynolds, to accretion rate. Comparing this formalism to Figures 14 and 15, we see that the dot-dashed curves show ρhu φ / ρ ≃ hu φ , while the dashed curves show ρhu r u φ /(− ρu r ). Their offset can then be attributed to the turbulent Reynolds stress normalized to the accretion rate, (δρu r )δ(hu φ ) /(− ρu r ). This turbulent Reynolds stress can be quantitatively significant, particularly in VD0 and to a lesser degree in ThickHR.
The separation between the curve of the accretion-weighted mean angular momentum and j net is the electromagnetic angular momentum flux, and in every simulation but VD0 it clearly makes the largest contribution to the outward angular momentum flux. In all the other cases, the only place where the Maxwell stress does not outweigh the Reynolds stress is in the immediate vicinity of the event horizon. There, the matter's angular momentum flux becomes almost exactly the total because a non-spinning black hole has no angular momentum to lose. For the same reason, just outside the horizon ||b|| 2 u r u φ comes to exceed b r b φ in magnitude, and the net electromagnetic angular momentum flux turns (weakly) negative.
To conclude, then, all these models show values of j net that are reduced below the ISCO value, regardless of H/R. The electromagnetic stress in the fluid frame hardly varies at all from one simulation to the next; consequently, its contribution to the net angular momentum flux is likewise nearly the same in all cases. Scale height does seem to have an effect on the run of u φ through the disk: u φ is reduced below the circular orbit value in proportion to the magnitude of radial pressure support. It also appears that the Reynolds stress levels may be partially controlled by the scale height in the disk and partially by the magnetic topology.

Review of Previous Simulations
In this section we summarize some results from previous simulations. Some of these, although not designed specifically to study the influence of scale height, can nevertheless provide additional information about how the inner disk stress depends on other parameters, notably magnetic field topology. Others are more directly comparable to our HR series of simulations.
We have already discussed models KD0c (Krolik et al. 2005) which began with simple dipole loops in the initial torus, and VD0 (Beckwith et al. 2009), a model that began with a vertical field piercing the initial torus. Beckwith et al. (2008a) feature three simulations all computed in a Kerr (a/M = 0.9) spacetime with identical fluid initial states, but differing initial magnetic fields: in KDPg, the field configuration was nested dipolar loops similar to those in KD0c; in QDP, there was a pair of quadrupolar field loops above and below the equator whose associated currents had opposite signs; and in TDPa, the initial condition held only toroidal field. Figure 2 of Beckwith et al. (2008a) provides data for investigating how these different initial field configurations affect inner disk stress. Panel f shows the coordinate-frame electromagnetic angular momentum flux (shell-and time-averaged T r φ (EM)) as a function of radius. The curves for KDPg and QDPa, the large dipolar field and quadrupolar simulations, are nearly identical. Both curves are generally a factor of several higher than that of TDPa, the one whose initial field was purely toroidal. Figure 2c of Beckwith et al. (2009) shows the radial run of net specific angular momentum (what we call j net here, but is labeled L in that figure). Both KDPg and QDPa are slightly below (by ≃ 0.05) the ISCO value of 2.1. TDPa, on the other hand, has an average accreted specific angular momentum of 2.13, very nearly the exact ISCO angular momentum, when averaged over the last 9000M in time. The similarity between quadrupole and dipole initial fields carries over to simulations with a Schwarzschild hole. Another GRMHD simulation, QD0 (described in (Beckwith et al. 2008b)), also has quadrupolar loops in the initial torus. It was run for 10 4 M in time, and averaging over the last 2000M in time gives a value of j net = 3.21 for the specific angular momentum accreted into the hole, similar to KD0c. As a whole, these results suggest that the crucial distinction may be between a field with significant poloidal character (either with or without net flux) and one that is only toroidal.
There is another case where we can clearly see the significance of a particular field configuration. The vertical field model VD0 was evolved in both two-and three-dimensions (Beckwith et al. 2009). In the present context the contrast between the two is interesting (see Fig. 7 of Beckwith et al. 2009). In the axisymmetric simulation the value of j net shows strong fluctuations. Between the times of 10 4 M and 1.5 × 10 4 M (the end of this simulation), the mean of j net is 2.85, but with a standard deviation of 0.43. Values as low as ∼ 1.0 are reached at particular moments. The mean value of j net over the same interval in VD0 is similar, 2.89, but the standard deviation is only 0.18 and the minimum value reached is 2.43. Axisymmetric simulations with vertical fields typically show strong MRI "channel modes" characterized by extended radial flows accompanied by radial magnetic field. It is the presence of those extended radial fields through the plunging region that provide the strong torques. Thus, the two-dimensional simulation illustrates the basic principle obtaining in three-dimensions, but in exaggerated form. In this particular case, the presence of a net vertical field (which cannot be reconnected away within the disk) also prevents the antidynamo effect from dissipating the turbulence, allowing the stress in the plunging region to remain over the full evolution. used a temperature regulation scheme based on maintaining constant entropy rather than a target temperature, but-compared to our ThinHR simulation-the resulting aspect ratio was only slightly thicker in the mean and somewhat less constant as a function of radius.
In their initial hydrodynamic conditions, our HR simulations were almost identical to those of Shafee et al., differing only in the q parameter (theirs was chosen to give an initial state with H/R ≃ 0.1, ours had variously H/R ≃ 0.05, 0.08, and 0.16). However, they did differ in the initial state of the magnetic field. Whereas our initial magnetic field was a set of nested dipole loops centered on the pressure maximum, they imposed two sets of loops, centered on r = 28M and r = 38M, which they then perturbed randomly with ≃ 50% fractional amplitude.
The spatial grid used by Shafee et al. was also very similar to the one we employed for ThinHR. They used a grid with 512 × 128 × 32 cells, whereas ThinHR used 912 × 160 × 64 cells. Both radial grids were logarithmic, but their grid extended to a slightly smaller radius than ours (50M as opposed to 70M); their radial cell size was therefore about 1.6 times larger than ours. The polar-angle grid scheme for ThinHR was finer than that of Shafee et al. near the equator, resulting in ≃ 30% more cells per scale height. The two simulations had identical azimuthal resolutions because, even though we used twice as many cells, our azimuthal extent was also twice as great (π/2 as opposed to π/4). Shafee et al. reported that in their initial condition, there were ∼ 10 cells per fastest-growing MRI wavelength, but do not say how that number changed during the simulation.
Despite this very close similarity, Shafee et al. arrived at a result that was quite different from ours. They found that the fluid-frame electromagnetic stress followed the Novikov-Thorne prediction very closely all the way to r = 9M, and then maintained more or less that amplitude all the way to r ≃ 2.5M. In contrast, as Figure 12 shows, in ThinHR (and all our other simulations), the fluid-frame electromagnetic stress is tens of percent above the Novikov-Thorne prediction in the disk body, and rises steeply inward inside r = 10M, reaching a level ≃ 5 times greater than the stress found in the Shafee et al. calculation. Similarly, although we found u φ (ISCO) − j net ≃ 0.33, Shafee et al. found a value less than half as large, only ≃ 0.14.
The origin of this contrast is uncertain. Although we have not yet performed a simulation with several dipolar loops, the very minor contrast between Kerr simulations with dipolar and quadrupolar initial field suggests that different forms of poloidal field are not, by themselves, significant. However, different field structures can place different demands on spatial resolution; as clearly demonstrated by MediumLR (and perhaps by KD0), inadequate resolution can lead to substantial artificial suppression of magnetic field strength. The figures illustrating results from MediumLR that we have shown in this paper all draw on data from the period during that simulation when it remained well-resolved; at later times, as its resolution quality failed, its electromagnetic stresses steadily weakened and u φ (ISCO) − j net diminished. Because it entails more small scale structure, a pair of dipole loops, as in the Shafee et al. initial condition, may create a turbulent magnetic field more vulnerable to reconnection than our single dipolar loop, particularly when perturbed by 50% and studied with a grid having larger radial cells.
It is possible that two other considerations may also play a role in creating these contrasting conclusions. First, their simulation ran only to a time of 10000M, and they presented no data demonstrating how well, and over what range of radii, it reached a state of inflow equilibrium. Because we found that a disk of this thickness takes ≃ 10000M to reach equilibrium, their shorter duration may be problematic. Second, the Shafee et al.
simulation had an azimuthal range of only π/4. Schnittman et al. (2006) showed that the characteristic azimuthal coherence length of features in full 2π global disk simulations was ≃ 1 radian. This result suggests that an azimuthal extent of only π/4 might misrepresent the MHD turbulence.

Discussion
It is almost twenty years since the first recognition (Balbus & Hawley 1991) that magnetic fields can produce significant stress with disks; it has been ten years since the first suggestions (Krolik 1999b;Gammie 1999) that this stress could continue within the ISCO, with implications for the overall efficiency and luminosity of accretion disk. However, there remains controversy about how large these effects may be and how they depend on disk parameters. Because the first general relativistic MHD accretion simulations demonstrated that these stresses are significant in disks with thickness H/R ≃ 0.15 (Krolik et al. 2005), the discussion in the last several years has centered on whether they might diminish with decreasing disk thickness (Reynolds & Fabian 2008;Shafee et al. 2008b). In this paper we have carried out a series of simulations carefully designed to isolate the effect of varying H/R; we find in all of our simulations that the electromagnetic fluid-frame stress increases steeply inward almost all the way to the event horizon. Indeed, even for a contrast of a factor of 3 from the thickest to the thinnest disk, we find an almost imperceptibly small change in the fluid-frame electromagnetic stress profile. The most natural interpretation of our results is that the radial distribution of electromagnetic stress depends at most only weakly on disk thickness. In our simulations, these stresses diminish j net , the net angular momentum per unit rest-mass accreted, by ∼ 10%. This result supports quantitatively the very crude qualitative argument given in Krolik (1999b) that the Alfvén speed in the plunging region would always become marginally relativistic, more or less independent of the accretion rate, so that magnetic stresses there would always be significant, but would always be dominated by gravitational forces.
It is also in keeping with other lines of qualitative reasoning. One might begin by asking, "If stirring of MHD turbulence by the MRI leads to significant magnetic stresses in the disk body, what might change near the ISCO?" As shown by Gammie & Popham (1998), the orbital shear in the plunging region in a Kerr spacetime differs from the Newtonian value by only a number of order unity. Consequently, if the magnetic pressure continues to be smaller than the gas pressure in the disk's equatorial plane, one would expect linear growth of the MRI to behave in very much the same way as in the disk body. If the gas pressure falls relative to the magnetic pressure so that the plasma β drops below unity, thereby quenching linear growth of the MRI, then magnetic stresses are surely important. However, nonlinear development of the turbulence can be expected to change as the infall time becomes as short or shorter than the eddy turnover time, the time for the energy of turbulent motions to move from long length scales to short. Indeed, earlier work (Hirose et al. 2004) has shown that the magnetic field becomes both much smoother and somewhat less tightly-wound (i.e., the radial component grows somewhat relative to the toroidal component) just inside the ISCO. Such a change in structure could alter both the rate of field amplification and the rate of field dissipation. On the basis of simple arguments like these, however, it is difficult to say whether these changes should lead to a larger or smaller mean field intensity and therefore stress. In the limit that both amplification and dissipation become weaker, flux-freezing results in growing magnetic stresses (Krolik 1999b;Gammie 1999). Moreover, with the principal dynamics-orbital shear and a growing radial velocity-all acting in the equatorial plane, there is no obvious place for a dependence on disk thickness. This physical argument is bolstered by the fact that the different simulation versions differ only slightly. Significantly different grid schemes, contrasting initial conditions, and even wholly different codes make only slight differences in the outcome. Even the topological contrast of substantial net vertical flux versus none at all seems to change the stress by only a modest amount.
On the other hand, we have also found that fluid effects in the inner disk can also contribute to a diminution in j net . Pressure support of the matter in the disk is, by definition, proportional to (H/R) 2 . Consequently, the mass-weighted mean angular momentum at any location in the disk is smaller than the angular momentum of a test-particle at that radius by a comparable amount. More significantly, Reynolds stresses can reduce it further, by an amount that increases both with increasing disk thickness and with net magnetic flux trapped on the horizon.
A few notes of caution should be injected into this discussion, however: Previous work studying accretion in a Kerr (a/M = 0.9) geometry (Beckwith et al. 2008a) suggests that a disk magnetic field that has no net poloidal content might produce weaker electromagnetic stress both throughout the disk and in the ISCO region. This effect should be explored more thoroughly in future work. We have also shown that before interpreting simulation results in terms of their implications for steady-state accretion, it is important to check carefully both that the simulation approximates inflow equilibrium and that the simulation's spatial grid provides adequate resolution throughout the period studied.
We conclude, then, that there appears to be little evidence for a strong dependence of near-ISCO electromagnetic stress on either disk thickness or the net magnetic flux. Because disk thickness is really a function of accretion rate, substantial near-ISCO electromagnetic stress should be seen in black hole accretion systems whether they are accreting at a rate near Eddington or far below. Its weak dependence on net magnetic flux suggests that the impact of electromagnetic stresses should be similarly independent of external magnetic boundary conditions. At the same time, we also find a supplemental reduction of the net accreted angular momentum due to Reynolds stresses, and this depends on both H/R and magnetic geometry. When the Reynolds stress is weakest, so only the near-universal electromagnetic stress acts, j net is reduced below u φ (ISCO) by ≃ 7-10%; when the Reynolds stress is strongest, the reduction is as large as 15%.
Fully quantitative conclusions, however, await several extensions of this work: to rotating black holes, to disks with more complex magnetic topologies, and to disks whose scale heights are constant, rather than proportional to radius. In the radiation-dominated regime (which should apply to the inner regions of accretion disks around black holes whenever the accretion rate is near Eddington, particularly when the central mass is large: Shakura & Sunyaev (1973)), H is constant as a function of radius. Therefore, a flat-topped disk is a more realistic model for disks as we are likely to find them in Nature. All of these extensions should be feasible in the near-term. This work was supported by NSF grant AST-0908869 and NASA grant NNX09AD14G (JFH), and by NSF grant AST-0507455 (JHK). Some of the simulations described were carried out on the Teragrid Ranger system at TACC and the Teragrid Abe system at NCSA, both supported by the National Science Foundation. Other simulations were run on the Johns Hopkins Homewood High-Performance Computing Center cluster.               a Average of H(r, t)/R in time (over each run's ∆t ave ) and in radius (weighted uniformly with respect to log r) from r ISCO to r pmax .
b Averaging period used here was 8000M − 10000M instead of ∆t ave because only these data were saved.