An investigation of the inner edge of a disk orbiting a black hole
would be best carried out with a fully general relativistic
simulation. However, we have not yet written a three dimensional
general relativistic MHD code. In the present study we instead make
use of the existing three dimensional Newtonian MHD code (Hawley
2000).
The code evolves the equations of ideal MHD, i.e.,

where
is the mass density,
is the specific internal
energy, **v**
is the fluid velocity, *P* is the pressure,
is the gravitational potential, and **B** is the magnetic field
vector. The global disk code is written in cylindrical coordinates
(*R*,
,*z*)
To model important effects associated with the relativistic
Schwarzschild metric, we employ the pseudo-Newtonian potential of
Paczynski & Wiita (1980). The pseudo-Newtonian potential has the
form

where *r* is spherical radius,
and
is the ``gravitational radius,''
akin to the black hole horizon. For this potential,
the Keplerian specific angular momentum (i.e., that
corresponding to a circular orbit) is

and
.
The innermost marginally stable circular
orbit is located at *r*_{ms}=3*r*_{g}.

We use an adiabatic equation of state,
We use an adiabatic equation of state,
,
with ,
and ignore radiation
transport and losses. Since there is no explicit resistivity or
physical viscosity, the gas can heat only through adiabatic
compression, or by the action of an explicit artificial viscosity (of
the form described by Stone & Norman 1992a) which acts in shocks.
Equations (1)-(4) are solved using time-explicit
Eulerian finite differencing. The numerical algorithm is that employed
by the ZEUS code for hydrodynamics (Stone & Norman 1992a) and MHD
(Stone & Norman 1992b; Hawley & Stone 1995). We adopt units where
*GM*=1 and *r*_{g} = 1 (so that
).

In this paper we repeat--with significantly better
resolution--simulation GT4 from Hawley (2000), whose initial state was
a thick torus with an angular velocity distribution
,
slightly steeper than Keplerian. The angular momentum
within the torus is equal to the Keplerian value at *R*=10, which
determines the location of the pressure maximum in the torus. In
somewhat arbitrary units, the pressure at this point is equal to 0.036
and the density
.
The initial magnetic field is
constructed by setting the toroidal component of the vector potential
equal to the density inside the disk,
,
for all
greater than a minimum value (
).
The total
integrated magnetic energy is then normalized to a value of
,
using the total integrated gas pressure within the torus.
The torus is perturbed with low-amplitude random adiabatic pressure
fluctuations. The aim of this initial condition is not so much to
follow the evolution of this particular torus, as it is to construct a
disk that is nearly Keplerian with an inner boundary that is initially
close to, but still outside, the marginally stable orbit. Evolution will
rapidly lead to inflow through *r*_{ms};
it is this resulting accretion
flow that will be studied in detail.

The computational grid has
128 x 128 x 128 zones. This is the
same number of zones as used in GT4 of Hawley (2000) except that here,
to increase the resolution within the disk itself and in the inflow
region, we locate more of the zones near the marginally stable orbit
and around the equator. The center of the coordinate system is
excised, i.e., the radial coordinate begins at a nonzero minimum
;
this avoids both the coordinate singularity
associated with the axis, and the
singularity at *r*_{g}. The outer
radius is set at *R*=31.5. The *z* coordinate is centered on the
equatorial plane, and runs from -10 to +10. The angle
is periodic and covers the full
.
In the radial direction, there
are 30 equally spaced zones in *R* from 1.5 to 4; the radial zones are
then graded logarithmically from 4 to the outer boundary at 31.5. In
the vertical direction, half of the *z* zones are equally spaced around
the equator from -2 to 2, then graded from those locations out to the
vertical boundaries at
.
The grid cells are equally spaced in
.
The simulation GT4 from Hawley (2000) used equally spaced
zones in *R* and *z*, and the outer radial boundary was set at 21.5.
We have also run this initial torus at lower resolution for comparison,
using 64^{3} grid zones.

The boundary conditions on the grid are simple zero-gradient outflow conditions; no flow into the computational domain is permitted. The magnetic field boundary condition is set by requiring the transverse components of the field to be zero outside the computational domain, while the perpendicular component satisfies the divergence-free constraint. Although this produces a nonzero stress at the boundary, it works well in preventing an artificial buildup of field at the boundary. Most of the time, for along the inner radial boundary and in the equatorial plane, the poloidal velocity is greater than the magnetosonic speed, so that the influence of the boundary condition on the flow upstream is limited.

Timescales of interest in the simulation are set by the circular
orbital period
at
significant points within the disk. In units where
*GM*=*r*_{g}=1,
*P*_{orb} at the pressure maximum (*R*=10) is
179, and at the
marginally stable orbit, where the orbital frequency
,
the orbital period is 21.8. The local accretion
timescale is set by the effective stress, namely
,
where
is
the *R*-
component of the stress,
,
and *v*_{A} is an
Alfvén velocity. The goal is to evolve long enough to observe a
number of accretion times near the marginally stable orbit. Since the
code is time-explicit and therefore Courant-limited, this is difficult
to achieve over large radial extents. Here we ran the simulation out
to time *t*=1500, which is 69 orbits at
*r*_{ms}.

The quantity of data associated with a three dimensional simulation is daunting, and this makes analyzing the outcome a challenge. Diagnostics can be computed during the simulation and in post-processing. The best procedures for this are still being developed. Here we briefly describe some of the diagnostic procedures used in the present simulation.

The time-evolution of the disk can be visualized through animation
sequences. Animation frames of the
slice and the equatorial
(*z*=0) slice are dumped every 4 units in time. Full 3D density
information is saved every 10 units of time, and these are used to
make
an animation of a three dimensional volumetric rendering.

The time-evolution can also be studied in a manageable way through
space-time (*R*,*t*) studies of azimuthally- and
vertically-integrated
quantities.
Examples include the averaged mass density

the average angular momentum,

and average specific angular momentum,

The net radial mass flux is

here defined to be positive for net accretion (inflow).
Similarly, the angular momentum flux,

can be computed.

Radial angular momentum transport in the disk is due to
the *R*-
component of the stress,

where
and
are the *turbulent* portions
of the velocity, that is, the departures from the mean flow. The
kinematic portion is the Reynolds stress and the magnetic portion is
the Maxwell stress. In the simulation we have no instantaneous way of
knowing what the mean flow is. For example, the typical radial
velocity is much larger than the time-averaged net radial drift
velocity, and there are significant differences between the
instantaneous angular velocity and the value associated with a
circular
orbit (which one might assume should be the mean orbital flow). As a
result, there is no unique prescription to compute the Reynolds
component of the total stress. Hawley (2000) used a definition for
the
perturbed orbital velocity
in terms of the difference
between the total instantaneous angular momentum flux, and the mass
flux times the average specific angular momentum. The difference then
represents the excess or deficit angular momentum transport due to
orbital velocity fluctuations compared to the mean. Specifically,

In contrast, we can easily compute the stress due solely to magnetic
forces, the second term on the right-hand-side of equation (12);
we denote this quantity by .
The
vertically- and azimuthally-averaged Maxwell stress is

In addition to these vertically- and azimuthally-averaged quantities,
we examine full data dumps of the simulation at specific moments in
time. From these we can examine not only the full three dimensional
structure, but also azimuthally-averaged and vertically-averaged
slices in (*R*,*z*) and (*R*,
).