**
Global Simulations
**
*
*

*
Equations and Algorithms
*

The dynamical evolution of an accretion disk is governed by the equations of ideal MHD, i.e.,

where *rho* is the mass density, *epsilon* is the
specific internal energy, **v** is the fluid velocity, *P* is
the pressure, *Phi* is the gravitational potential, and **B**
is the magnetic field vector. In the present study we adopt an
adiabatic equation of state, *P = rho epsilon (Gamma -1) = K
rho ^{Gamma}*, and ignore radiation transport and losses.
Since there is no explicit resistivity or physical viscosity, the gas
can heat only through adiabatic compression, or in shocks through
explicit artificial viscosity,

In some simulations below, the code employs the usual Newtonian
gravitational potential, *Phi = -GM/r*, where *r* is the spherical
radius. To deal with the inner boundary in a more realistic manner,
however, the code also makes use of the pseudo-Newtonian potential
(Paczynsky & Wiita 1980) to model the relativistic effects
associated with a Schwarzschild metric (e.g., minimum stable orbit).
This inner boundary permits supersonic (and super-Alfvenic) accretion
off of the inner radial grid, reducing the likelihood of unphysical
influences from the boundary conditions. The pseudo-Newtonian potential
has the form

(5)

where *r _{g}* is the ``gravitational radius''
(akin to the black hole
horizon), here set equal to 1. For this potential,
the Keplerian specific angular momentum is

(6)

and * Omega R ^{2}= l*. Here we set

Equations (1)-(4) are solved using time-explicit Eulerian finite
differencing. The global disk code is written in in cylindrical
coordinates, *(R, phi,z)*. The center of the coordinate system is
excised, i.e., the radial coordinate begins at a nonzero minimum
*R _{min}*; this avoids the coordinate singularities
associated with the axis. In some simulations, the gravitational
potential is assumed to be cylindrical; this removes the vertical
component of gravity, and reduces the number of

FIGURE 1:A schematic view of the cylindrical coordinate computational domain. The radialRdimension begins at a minimum radiusR, the vertical_{min}zdimension is centered on the equator, and the angular directionphiruns from 0 to2 pi/mwheremis an integer. The simulations in this paper usem=1andm=4.

Most of the three dimensional global simulations are run on parallel computer systems with a version of the global code that uses the Message Passing Interface (MPI) library for interprocess communication. The full computational domain is divided into overlapping subgrids, one for each processor. Boundary zone values are passed between processors using explicit MPI routines. In the description of what follows, however, this subdivision will be ignored; the physical domain size and the number of grid zones discussed will be the totals.

The MHD algorithm now has a long history of use. The magnetic field is
evolved with the constrained transport (CT) approach of Evans and
Hawley (1988) that preserves the constraint *nabla dot B =0*. The
CT scheme was designed to work in any general curvilinear coordinate
system, and here the addition of appropriate geometric terms easily
adapts the procedure to cylindrical coordinates. The algorithm uses
information propagated along Alfven characteristics to solve a
restricted set of characteristic equations for time-advanced fields and
electromotive forces within the CT framework. This approach is known
as the Method of Characteristics Constrained Transport (MOCCT)
algorithm, and its first implementation is described in detail in Stone
& Norman (1992b); the current version is described in Hawley & Stone
(1995).

Considerable effort has gone into determining the most satisfactory
scheme for the individual numerical terms in the present application.
Local conservation of angular momentum is particularly important, and a
number tests using hydrodynamic Keplerian disks on a variety of grids
suggest that some numerical approaches are better than others.
Specifically, it is best to evolve *rho l*, where *l* is the
specific angular momentum, as a fundamental variable (as opposed to
angular velocity, or specific angular momentum alone) and to use
consistent advection (Norman, Wilson & Barton 1980) where the
advection of all variables is tied to the advection of the density. To
minimize errors associated with the coordinate singularity located at
*R=0* we make use of a regularized form of the operator
*R ^{-1} dR * (Evans 1986).

A variety of test simulations were used to verify code performance. The hydrodynamic implementation was tested with Keplerian disks and a series of simulations of constant angular momentum tori similar to those studied by Hawley (1991). Constant and near-constant angular momentum tori are subject to particularly vigorous growth of the Papaloizou-Pringle instability (Papaloizou & Pringle 1984), especially when such tori are ``slender,'' meaning that their cross section is small compared to the radius of the pressure maximum. Slender tori are unstable to the principal mode of the Papaloizou-Pringle instability, with growth rates near the maximum value. These tests gave results that were consistent with past simulations.

For an MHD test, the grid was moved out to large radius and an isolated vertical magnetic field flux tube was embedded in a Keplerian flow. This test is a three-dimensional version of the axisymmetric simulations described in Hawley & Balbus (1991). Again the results were fully consistent with the earlier work.

*
*

*
Accretion Tori
*

Most of the simulations discussed below consider the problem of the
nonlinear evolution of the magnetorotational instability in accretion
tori. Tori have significant internal pressure gradients that balance
the gravitational and centrifugal accelerations in hydrostatic
equilibrium. This pressure gives the torus a vertical thickness
*H* that can be comparable to its radius *R*. Significant
departures from Keplerian angular momentum distributions are possible
as the torus becomes ever thicker; the limit is the constant angular
momentum torus. A particularly useful feature of these tori for
numerical simulations is that they are well-defined equilibria that can
be completely contained on the finite domain of a computational
simulation. This allows their evolution to proceed (at least
initially) independent of the grid boundary conditions.

For an adiabatic equation of state, a gravitational potential
*Phi*, and an assumed rotation law * Omega ~ R ^{-q}*,
the density in the torus is determined by

(7)

where *l _{k}* is the Keplerian angular momentum at the pressure
maximum and

Tori are known to be subject to the purely hydrodynamic global instability discovered by Papaloizou & Pringle (1984). The Papaloizou-Pringle instability has been extensively investigated both with perturbation analyses, and with several numerical simulations. The mechanism of the instability was elucidated by Goldreich, Goodman & Narayan (1986) and involves the establishment of a global wave that gains energy through the exchange of angular momentum between the inner and outer region of the disk. Amplification requires wave reflection at at least one disk boundary. The nonlinear evolution of the Papaloizou-Pringle instability has been followed for both slender (Hawley 1987), and relatively wide (Blaes & Hawley 1988; Hawley 1991) tori. In the slender torus case, the instability saturates as a highly nonaxisymmetric orbiting fluid ellipse (the ``planet'' solution; Hawley 1987) which itself proves to be an equilibrium solution (Goodman, Narayan, Goldreich 1987). Wider tori saturate in spiral pressure waves (Blaes & Hawley 1988; Hawley 1991).

Torus linear stability and, as we shall see, nonlinear dynamics are
quite different when a magnetic field is present. The largest growth
rate of the MRI greatly exceeds those typical of the hydrodynamic
Papaloizou-Pringle instability. Linear stability analyses bear this
out. Curry & Pudritz (1996) and Ogilvie & Pringle (1996) have done
global linear MHD stability analyses of cylindrical systems (*1/R*
gravitational potential) with a variety of initial angular momentum
distributions. They find rapidly growing unstable modes for a wide
range of equilibria and field strengths.

Although global linear analyses are most appropriate for strong magnetic fields whose most unstable wavelengths are comparable to the size of the torus, weak fields have unstable wavelengths that are much smaller than the torus dimensions. The essential physics of the MRI in these circumstances is local, and the nonlinear evolution should be similar to that seen in the local shearing box simulations. This is, of course, a testable proposition. In any case, global simulations will provide information regarding the ultimate consequences of the MRI for the evolution, and, indeed, existence of thick tori. Obviously, the initial equilibrium torus cannot persist in the face of vigorous instability and the resulting angular momentum transport. Its overall global evolution, however, may present novel features, and the end state of its evolution may represent a more realistic accretion disk structure.

*
*

*
Diagnostics
*

A global simulation will involve a grid with a million to a billion zones, run for possibly hundreds of thousand of timesteps. A major task is to develop useful diagnostics that concisely and adequately characterize the essential behavior of an evolving disk in such a simulation.

The simplest type of diagnostic is the volume integral of a quantity over the entire grid. For example, the total mass is

(8)

Time histories are computed of this quantity and others such as total kinetic and magnetic energies (by component), and the radial mass flux

(9)

through both the inner and outer boundaries.

Standard steady-state thin disk models are often expressed in terms of
vertically averaged quantities. Here such values can be
computed by integrating
density-weighted quantities over *z* and *phi*. The
vertically- and azimuthally-averaged mass density is

(10)

The averaged mass flux, * M *, is defined by
(9),
and one can construct
similar averages for other quantities. For non-Keplerian disks the
radial dependence of the average specific angular momentum
will be of interest, and this can be computed from

(11)

The stress * T _{R phi}* consists of a magnetic component
(the Maxwell stress)

(12)

and a kinematic component (the Reynolds stress)

(13)

The value of the perturbed angular velocity * delta v _{phi}* is
defined for a Keplerian disk by

(14)

Together the Maxwell and Reynolds stress make up the total stress,
which, in the Shakura-Sunyaev (1973) parameterization, is set equal to
*alpha P*. The stresses observed in the simulations can be similarly
scaled using the vertically averaged pressure to derive an *alpha* as
a function of *R*. In the simulations, *alpha* varies with both time
and space. Its utility is mainly to provide a familiar point of
reference for characterizing the observed stress. The local shearing
box simulations (HGB; HGB2) found that the Maxwell stress, which
dominates the total transport, is highly correlated with the total
magnetic pressure. This suggests an alternate stress parameter, the
``magnetic alpha value'', defined as

*alpha _{m } = 2 B_{R}
B_{phi} / B^{2}* .

In a steady state disk, the total stress, Maxwell plus Reynolds, would be
equal to radial angular momentum flux carried inward by the net radial
drift velocity, minus the angular momentum flux at the inner boundary
of the disk where the stress is presumed to vanish. The latter term is
generally ignored for radii large compared to the inner disk boundary.
In the simulations none of these assumptions hold *a priori*: the
typical torus radius is comparable to the inner disk boundary, the
stress does not vanish at the inner disk boundary, and the disk will
not be in steady state.
In fact, the simulations will prove to be highly dynamic
with rapid variations in both time and space. This is best
appreciated from time-dependent animations which also have been used
as diagnostics for these global simulations. While reporting time- and
space-averaged quantities tends to obscure this feature,
such averages are more practical for a written summary.