The initial torus was in hydrodynamical equilibrium, but it evolves rapidly as the simulation proceeds. The overall evolution is very similar to that of the GT4 of Hawley 2000 (see Fig. 13 of that paper). The radial field is sheared by the differential rotation, creating toroidal field. Turbulence develops within the disk with the onset of the magnetorotational instability. The resulting Maxwell stresses drive angular momentum transfer; the disk begins to accrete into the central hole and expands radially outward. Relatively rapidly the inner portion of the disk achieves an approximate steady-state. By , the mass accretion rate reaches roughly the long-time mean; after , the shape of the disk's average radial density distribution no longer changes appreciably. However, it must be borne in mind that there are always significant fluctuations around all these mean states. Thus, we observe the disk in a quasi-stationary state for several dozen orbital periods at the marginally stable orbit. In much of what follows, we will present a detailed analysis of the state of the disk during this time frame, focusing particularly on the end of the simulation, using that ``snapshot" as a fairly typical sample of conditions in the quasi-stationary disk.

A contour plot of the azimuthally-averaged gas density at the initial
and final times is presented in Figure 1. As might be expected,
the
density diminishes inside the marginally stable orbit and outside the
initial outer radial boundary of the disk, while also decreasing with
increasing altitude away from the equatorial plane. Although the
vertical height of the highest contour levels at the final time are
roughly independent of radius within the disk proper, the
azimuthally-averaged exponential scale-height of the pressure,
*H*,
increases very nearly linearly with radius from
out to ;
at the marginally stable orbit (*R* = 3), ,
while at *R* = 10, .
This particular disk shape is
almost certainly a consequence of the equation of state, in which all
fluid elements retain the same specific entropy except for shock
heating. In a real disk, one in which turbulent dissipation
introduces
heat and radiation vents it, the shape might be rather different.

**Figure 1:**
The azimuthally-averaged gas density distribution at the beginning
and the end (*t*=1500) of the simulation. There are 14 evenly spaced
contours. The maximum density at *t*=0 is 34.2; at *t*=1500,
it is 14.3.

Density gradients in radius and altitude are determined primarily
by the overall dynamical effects of approximate hydrostatic equilibrium.
No such requirements exist for azimuthal fluctuations, whose
amplitudes characterize the turbulence. We
quantify the azimuthal fluctuation level by defining

where
refers to the azimuthal average at
fixed *R* and *z* of the quantity *X*. Measured in
this way,
- 0.6 over most of the problem volume,
with smaller values typically nearer the equatorial plane. In a few
places,
rises to be greater than unity. Thus, quite
sizable density fluctuations are generic.

**Figure 2:**
Density in the equatorial plane for *R* < 10
at *t*=1500.

These fluctuations are primarily sheared strips whose radial extent is
about 0.1*R*, but extend ~1 radian in azimuthal angle
(Figure 2).
This pattern (tightly wrapped spiral features) is characteristic of
fluctuations in almost all quantities. Part of their nature is
revealed by constructing spacetime diagrams in which the surface
density at a specific radius is plotted as a function of *t* and
.
From this visualization we are able to measure the angular speed of
these features; it is generally quite close to the local orbital
frequency
.
We can also measure their persistence at any particular radius;
it is usually less than an orbital period. In this sense, they can be
regarded as fluctuations that merely ``ride" with the local orbital
motion, but see further discussion in §3.5.

The magnetic field distribution is rather more complicated than the
density distribution. A poloidal projection (i.e., averaged
azimuthally) shows the magnetic pressure to have its largest value
outside the dense regions of the torus, with large variations inside
the torus (Figure 3). In the shearing-box simulations of Miller &
Stone (2000), the field strength likewise peaked away from the midplane,
but we find somewhat greater contrast. Field-strength
fluctuations are also larger than
those for density: the *rms* fractional fluctuations in the
azimuthal direction are about 1 to 1.5
in the disk body, falling to ~ 0.25
several gas scale-heights immediately above and below the
disk. Like the density, the characteristic shape of fluctuations is
sheared arcs that run for ~1 radian before losing coherence.
However, the peaks in magnetic field strength tend to lie adjacent to,
rather than coincident with, peaks in gas density.

**Figure 3:**
The azimuthally-averaged magnetic pressure
plotted on a logarithmic scale from 10^{-5}
to 10^{-2}. The mean value of magnetic pressure
within the region where the density is greater than 0.01
is 7.1 × 10^{-4}. Using the mean gas pressure this
corresponds to an average plasma beta value of 13.

The structure of the field in the inner part of the disk is illustrated by Figure 4, the azimuthally-averaged poloidal magnetic field vectors. There is a complicated field pattern within the disk, which is surrounded by a stronger and more organized field of roughly dipolar character. Figure 4 illustrates the poloidal field, but in general, the toroidal field is several times larger; ratios anywhere from ~1 to ~10 are common.

As one expects from magnetic fields in shearing plasma, there is a
mean
correlation between *B*_{R} and
in the sense that the ratio of
the Maxwell stress to the magnetic pressure,
and has (approximately) constant magnitude
- 0.3 in the body of the disk (see Figure 12). Inside the
marginally stable orbit, this ratio increases slowly in magnitude as
radial accretion ``combs'' the field out into radial lengths that are
strongly sheared, increasing the correlation between
*B*_{R}
and
.
In fact, the Maxwell stress is so correlated with the
average magnetic pressure
,
that a figure
showing the stress is almost indistinguishable from one showing the
magnetic pressure.

**Figure 4:**
Azimuthally-averaged poloidal magnetic field vectors
in the inner region of the simulation overlaid on density contours. The
length of each arrow is proportional to the magnitude of the poloidal
field.

Perhaps the single most basic quantity of interest is the accretion
rate, defined in equation (10). The disk is never perfectly
time-steady; v_{R} varies substantially from place to place
and from time to time. In fact, disk quantities tend to be quite
nonaxisymmetric, and much of the final accretion inside
*r*_{ms} takes place in a spiral flow. However,
time-averaging can remove much of this variability, and Figure 5 shows
the time-averaged accretion rate as a function of *R* along with
the instantaneous values of M-dot
at *t*= 1000 and 1500. The time-averaged accretion rate
is nearly constant as a function of radius inside *R*=10,
and equal to 5 in code units. We shall adopt
this as the fiducial steady state accretion rate for the inner disk.
Between *R*=9 and 14 there is inflow at a decreasing rate; outside
*R*=14 the net flux is outward. In this simulation, the initial disk
is entirely contained within the grid, and the outer part of the disk
must move outward in order to ``soak up" the angular momentum received
from the inner disk material. Over the course of the simulation the
total disk mass decreases by 14%; about 90% of this has gone into the
black hole.

**Figure 5: **
The accretion rate (in code units) as a function of
radius at times *t*=1000
and *t*=1500 (solid lines). The dashed curve is the time-average of the
accretion rate between those two times. Note that
negative accretion rate means outward flow.

Figure 6 shows a detailed time-history of the accretion rate
through
the inner radial boundary. As can be seen from this figure, the
accretion rate shows significant time variations after its initial
growth phase. The nature of the time-variability can be studied using
the Fourier power spectrum of the accretion rate
during the latter portion of the simulation (beyond *t*=600,
i.e., after theinitial growth phase is over). This is shown in
Figure 7. It is
most simply described as a ``red" power-law; roughly speaking,
.
However, there is also some curvature, in the
sense that it steepens with increasing frequency. Because the
accretion
rate does appear to have a well-defined mean value,
must
be > -1 at the lowest frequencies. At the high frequency end,
the spectrum is largely featureless, extending almost to the inverse
free-fall time at the innermost boundary. Other than a smooth
steepening, there is no special feature at the orbital frequency
of the marginally stable orbit (marked on the figure). The flow is
unsteady everywhere, even well inside the plunging region.

**Figure 6:**
The mass accretion rate through the
inner radial boundary as a function of time. To capture the
time-variability accurately, the accretion rate was measured every 10
timesteps in the simulation, corresponding to an average interval of
0.0037 in time.

**Figure 7:**
Fourier power density of the accretion rate through the inner boundary
(solid curve) and volume-integrated Maxwell stress (dashed curve).
Frequency is in code units; the orbital frequency (i.e., the inverse of
the orbital period at the marginally stable
orbit is marked with an arrow.

Figure 6 provides some idea of how nonsteady the accretion is;
the mass flux is highly spatially inhomogeneous as well. To gain some
sense of the magnitude of the spatial variations, we show the accretion
rate as a function of position on a cylinder at *R*= 5
(Figure 8). As this figure shows, the fluctuations are very
large. If we suppose that the mean accretion rate of 5 is evenly
distributed over the range of altitudes -l < |z| < 1 the mean
mass flux would be ~0.1 As the figure makes plain, there are several
places where the mass flux is five times as great as this; there are
also several places where the mass flux is ~0.3-0.4 *outward*.
Relatively large fluctuations compared to the mean are to be expected
when the transport of angular momentum in the disk is by turbulence.

**Figure 8:**
The accretion rate at cylindrical radius *R*=5; inward mass flow
is given a positive sign. Note that the
* local* mass flux is almost as likely to be outward as inward.

As we have stressed all along, one of the most striking features
of this simulation is the presence of large amplitude nonaxisymmetric
modulations in the density, mass flux, and magnetic field strength.
We have already discussed how
* (,t)*
spacetime plots enable us to
measure their angular velocity. It is possible to learn more
by constructing analogous plots in
(*R*,*t*). For example, in Figure 9
we show the history of the azimuthally- and vertically-averaged
mass flux (eq. )
in this simulation. The diagonal stripes show the radial motion
through the disk of fluctuations in the mass
flux, their slopes indicating their radial group speeds. Typically
these disturbances begin near
*R* ~ 5
and travel both outward and
inward. In the outward direction, they move at a very nearly constant
speed
~0.7 while in the inward direction, their speed is at
least three times faster. Because the outward speed very nearly
matches the magnetosonic speed, we infer that these disturbances
propagate as magnetosonic waves in the outward direction. On the other
hand, inward-going mass flux perturbations are simply advected along
with the dynamical inflow.

We conjecture that their origin lies in the way in which material enters the region of unstable orbits. When a fluid element begins to accelerate inward, it creates a rarefaction wave that moves outward at the magnetosonic speed, while traveling azimuthally at the (faster) orbital speed. On the other hand, the radial inflow speed near and inside the marginally stable orbit is comparable to the magnetosonic speed, so fluctuations travel inward rather more quickly. These fluctuations are never smoothed out into a completely steady flow because individual fluid elements move into the plunging region as a result of the continually changing torques exerted by the MHD turbulence (see the next subsection). Thus, at this level of detail, there is no way to achieve a stationary state.

**Figure 9:**
Vertically- and azimuthally-averaged mass flux as a function of
radius *R* and time. The diagonal stripes denote radially-propagating
waves. Reds and yellows correspond to accretion, blues to outflow.

Although strong nonaxisymmetric pressure waves were generated by MHD turbulence in the shearing box simulations, they had no significant impact on the evolution of those simulations. In this global simulation, however, as Figure 9 shows, they play a significant role in the creation of large fluctuations in the accretion rate.

The total stress must satisfy the equation of angular momentum
conservation. In a steady-state disk, this means that locally the
torque must be equal to the mass accretion rate (which is constant)
times the difference between the local specific angular momentum
and the specific angular momentum carried into the black hole
.
Thus, the stress is equal to

where is the total mass accretion rate, the rotational frequency (for the pseudo-Newtonian potential, ), and is the specific angular momentum ( ). In the time-steady, axisymmetric ``Novikov-Thorne" disk, it is assumed that because the stress is assumed to be zero for all .

For the purpose of studying angular momentum transport, it is useful
to look at the azimuthally- and vertically-integrated magnetic stress,
,
as this is the quantity most directly
related to transfer of the *z*-component of angular momentum. In
particular,
we compare this quantity to the total stress per unit area
that would be expected in the steady state ``Novikov-Thorne" disk. As
discussed above, in the inner part of the disk the mean accretion rate
is about 5 in code units. In Figure 10
(for
and
)
is
contrasted with
averaged
from *t*=750 to *t*=1500. As that figure
shows, in the main body of the disk (approximately
),
the average magnetic stress is almost exactly equal to the expected
stress.
At large radius the two curves diverge because
equation (16) assumes a
time-steady disk,
whereas ours
has only a finite mass, so the behavior at its radial exterior is
quite
different. At the inner edge of the disk, in contrast to the
prediction of the conventional model, the azimuthally- and
vertically-averaged
magnetic stress does not
diminish, but, if anything, increases slightly inward (the final
upward
spike is an artifact of the boundary condition).
Note, however, that the local value of the ratio
in the plunging region varies widely:
from as much as ~10
several scale-heights out of the equatorial
plane to as little as ~10^{-4}
at some locations in the plane.

**Figure 10:**
The azimuthally-averaged and vertically-integrated magnetic
stress averaged from *t*=750 until the end of the simulation
(solid curve) and nominal mean stress per unit area (dashed curve).

In the context of time-steady and axisymmetric disks, the stress is
often discussed in terms of the dimensionless parameter introduced by
Shakura & Sunyaev (1973),

(17) |

It is possible to similarly calibrate the stress in this simulation if
is generalized to be the *local* ratio of stress to
pressure
,
which we will approximate by supposing that
.
Unlike ,
which, in most studies, is taken to be a constant independent of time and
location,
in this simulation is a dynamical and highly variable quantity
that is determined by local and transient magnetic dynamics. For this
reason, the absolute value is significant; although
is positive most places, there are counter-examples.

To give a sense of its scale and range of variation, we first examine
,
the azimuthally-averaged version of this
ratio, as found in the late-time snapshot (Figure 11). In the
main portion of the disk,
ranges between ~10^{-2} and ~10^{-1}.
However, it rises dramatically to ~10
at several scale-heights above and below the mid-plane. The
same tendency has been seen (with much better vertical resolution)
in stratified shearing box simulations (Miller & Stone 2000).
also increases sharply at small radius.
In the mid-plane,
at *R*=3,
but rises to
at the inner radial edge of the simulation.
Moreover, the vertical
rise in
also becomes sharper at small
radius, so that at *R*=3,
at |z| as small as 0.5.

**Figure 11:**
Absolute value of
,
in a logarithmic scale in the snapshot at *t*=1500.
This quantity increases with decreasing radius, and
with increasing distance above or below the disk midplane

In addition to these systematic trends,
is also subject to
sizable fluctuations. Measured in the same way as for the magnetic
energy density and the gas density (i.e., in terms of fluctuations in
azimuth at fixed *R* and *z*), the *rms* fractional
fluctuation amplitude in the bulk of the disk ranges from
to ,
but at large altitude above the disk midplane, the fluctuations can
be much greater. The character of these fluctuations is illustrated
by a different projection of ,
which we call :

(18) |

the density-weighted vertical average of . In this projection (Figure 12), we clearly see the nature of the azimuthal fluctuations: just as for the density fluctuations (shown on smaller scale in Figure 2), the magnetic fluctuations are sheared filaments of limited angular coherence length.

**Figure 12:**
Absolute value of the vertically-averaged ,
plotted
on a logarithmic scale for the snapshot at *t*=1500. Note the
large fluctuations, even when values are compared at fixed radius.
is large at both large and small *R*, but for different
reasons (see text).

These results highlight the limitations of the primary assumption made in many previous studies of disk structure, that is a universal constant, independent of time and location. By contrast, we find that the ratio of stress to pressure varies strongly from place to place, both in the sense of having large fluctuations and in the sense of possessing strong systematic gradients. In many places the gas pressure is relatively smooth while the magnetic stress varies considerably. In other places the variation in is simply due to changes in the background pressure. For example, the magnetic stress does not drop with height as rapidly as the pressure does.

For purposes of comparison with standard disk models, it is possible to
compute (albeit with some trepidation) the value of
found
in this simulation. This quantity, defined in terms of the Maxwell
stress, is shown as a function of radius in Figure 13 for the disk at
the end of the simulation. It is
through most of the disk,
but, starting just outside the marginally stable orbit, it rises inward.
It reaches a peak
near the inner radial boundary of the
simulation. The rise is primarily due to the rapid drop in gas
pressure as the accretion flow accelerates inward. Also shown is the
ratio of the Maxwell stress to the magnetic pressure (dashed line).
This
is 0.2-0.3 throughout most of the disk, and
increases inside *r*_{ms}.

**Figure 13:**
The Shakura-Sunyaev alpha parameter as a function of radius (solid
line) constructed from vertical and azimuthally averaging the Maxwell
stress and the gas pressure. Also shown is the averaged magnetic alpha
as a function fo radius (dashed line), constructed using the magnetic
pressure in place of the gas pressure. The magnetic alpha is
relatively constant in the disk, whereas the Shakura-Sunyaev alpha
varies with the gas pressure.

Whether the Shakura-Sunyaev alpha
is meaningful depends on how it is used. As the
previous discussion emphasized, the local ratio of magnetic stress to
pressure has strong vertical gradients that are completely obscured
when the vertically-integrated version is used. Similarly, there are
large fluctuations in azimuth and as a function of time. For these reasons,
it is useless as a local dynamical variable. However,
the approximate constancy of alpha
as a function of radius *in the disk proper* means that it can
be used as a
rough estimator of
the surface density in the body of the disk, provided its use is not
extended to local properties nor is it thought to be accurate at better
than ``factor of a few" accuracy. This conclusion is entirely
consistent with many standard uses of
and within the spirit
with which it was originally introduced by Shakura and Sunyaev (1973)
as a scaling parameter, not a fundamental variable.

Another significant diagnostic of the flow is the mean specific angular
momentum distribution. There are a variety of ways to characterize the
specific angular momentum; one is by vertical and
azimuthal average,
,
as defined in equation
(9). This function was discussed in Hawley
(2000) for the model
GT4 (see Fig. 14 in Hawley 2000); the present model is very
similar.
At late time the average specific angular momentum is slightly
sub-Keplerian throughout most of the disk. The most striking feature
is that
continues to decline inside *r*_{ms}
rather than leveling off. If all stress ceased at the marginally
stable orbit, material would travel inward from there retaining
precisely its energy and angular momentum. Conversely, the degree to
which fluid changes its specific energy and angular momentum is a
measure of the stresses that exist in the plunging region. The mean
specific angular momentum diminishes from 2.60 at
*r*_{ms} to 2.48 at
the innermost radius. Thus, on average in this snapshot, the accreted
matter has lost 5% of the angular momentum it had when it
was at the marginally stable orbit. It is likely (for reasons
explained in §5) that this number is an underestimate, and in
reality
the specific angular momentum at the innermost radius is rather
smaller.

The total rate at which angular momentum is carried across the inner radial boundary is slightly smaller than these numbers would suggest because there is also a small outward electromagnetic angular momentum flux. This electromagnetic flux is due to a generic property of magnetic fields in shearing plasma: when the orbital frequency decreases outward, the field is swept back so as to transport angular momentum outward, i.e., . However, near the inner radial boundary this is a small effect, the electromagnetic angular momentum flux there is only about 2% of the angular momentum flux carried by the matter.

Further complexity is revealed by a closer examination of the
(*R*,*z*)
distribution of the late time azimuthally-averaged angular momentum
inside the plunging region. In fact, the specific angular momentum
in the cells lying along the equatorial plane
is nearly constant inside *r*_{ms}. However, the
value of
drops with *z* above and below the equator. This is consistent
with the distribution of magnetic field which, of course, accounts for
the nonzero stress. The low value of field near the equator is in part
a consequence of the assumed symmetry of the initial distribution of
field. The initial field loops are symmetric with respect to the
equator (*B*_{r} and
are zero there), and throughout the run the
strongest fields tend to be located above and below the equatorial
plane.

Just as for every other quantity, there is also considerable spiral structure. Figure 14 shows the density-weighted vertically-averaged specific angular momentum. It makes plain how nonaxisymmetric the distribution of torque is in the plunging region.

**Figure 14:**
The density-weighted vertically-averaged specific angular momentum
in the inner portion of the accretion flow.

For a number of reasons, the accreted energy is harder to determine
than the accreted angular momentum (see further discussion in
§5.2).
To introduce this topic, we begin by displaying (Figure 15) the
specific energy *e* in rest-mass units,

in the equatorial plane at late time. The specific energy for a
circular orbit precisely at *r*_{ms} is -0.0625;
here,
is a bit greater than that ()
because of
the thermal, magnetic, and turbulent kinetic energy contributions to
the total energy. In the equatorial plane,
is almost constant from *R*=3 to just outside *R* = 2, but
actually *rises* as *R* decreases from there to *R*=1.6,
reaching a peak that is very nearly unbound: -0.01. The
azimuthally-averaged energy in the equatorial plane then falls sharply
to -0.07. Even this wild swing underestimates the fluctuations: the
smallest specific energy found in the equatorial plane is -0.19. Thus,
we are confident in asserting that there is a large amount of energy
transfer within the plunging region. However, given these large
fluctuations, we do not feel that we can compute a well-defined mean
energy for the accreted matter. It is likely, moreover, that at least
some of the very large fluctuations in the innermost radial zones are
not physical, but due to the inner boundary condition.

**Figure: 15**
The specific energy of matter (in rest-mass units) in the
equatorial plane within *R*=4 (the high-resolution segment of the
simulation). Note the large amplitude spiral fluctuations. The sharp
features immediately surrounding the inner radial boundary are
artifacts.

In fact, the situation is somewhat worse than this. The reason is that dissipation is crucial for computing the accreted energy, and our simulation contains no real dissipation physics. Absent numerical losses (see §5.1, 5.2), energy taken from plunging matter may go several places. One possibility is to build magnetic field energy in the plunging region itself. If the field energy is advected into the black hole without dissipation and radiative losses, the transfer of energy from fluid to magnetic field makes no difference to the accreted energy. It is also possible for magnetic stresses to transfer energy all the way back to the disk proper (and indeed we see significant magnetic stress throughout the plunging region). Again, without dissipation, the energy transferred in that way stays in the fluid, whether in the form of kinetic or magnetic energy. Only in the event that the energy is transferred all the way out to non-accreting matter would it be possible for that energy to avoid being carried back in. Thus, to the degree that this simulation avoids numerical losses, it is very difficult for it to show any significant transfer of energy out of the plunging region. By contrast, in a real disk, dissipation changes energy into heat, which may then be radiated if the thermal timescale is shorter than the inflow timescale; in our simulation, neither physical dissipation nor radiation is present.