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 rhoGamma, 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, Q. Artificial viscosity appears as an additional pressure term in equation (3), i.e., -(P+Q) nabla dot v (Stone & Norman 1992a). The absence of physical shear viscosity in these equations is worth noting. This is the actual state of the plasma: the molecular viscosity is utterly negligible. Accretion disks evolve due to transport resulting from what is normally referred to as anomalous viscosity. But this is not an additional term in the equations; transport results from the dynamics contained within the above equations.

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 rg 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 R2= l. Here we set rg= 1, and, for both gravitational potentials, GM=1. This determines the units of time in the simulation with Omega =1 at R=1 for a Newtonian gravitational potential. All times given will be reported in these units; the corresponding orbital periods at locations in the tori will also be given where appropriate.

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 Rmin; 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 z grid zones required. When the vertical component of gravity is included, the z coordinate is centered on the equatorial plane. The angle phi is periodic and covers the full 2 pi or some integer fraction of 2 pi. In some simulations, only the two-dimensional axisymmetric system (R,z) will be computed. A schematic of the computational domain is shown in Figure 1.

FIGURE 1: A schematic view of the cylindrical coordinate computational domain. The radial R dimension begins at a minimum radius Rmin, the vertical z dimension is centered on the equator, and the angular direction phi runs from 0 to 2 pi/m where m is an integer. The simulations in this paper use m=1 and m=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 lk is the Keplerian angular momentum at the pressure maximum and C is a constant of integration that establishes the zero pressure surface of the torus. A given torus is specified by the angular velocity gradient q (with q=2 corresponding to constant angular momentum), the radial location of the pressure maximum Rkep, and the radial location of the inner edge of the torus, Rin} (which determines C). Outside of the torus the grid is filled with a cold gas with low density. A floor value of the density and internal energy is applied to ensure that these zones do not become evacuated, or develop extremely high sound speeds. The floor value of density should be set low enough so that its contribution to the total mass or energy is not significant at any place on the grid. Here the floor value is 10-3 and the typical torus initial density maximum is ~ 10.

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 TR 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 vphi is defined for a Keplerian disk by delta vphi = vphi - R OmegaKep. For other angular momentum distributions, including average angular momentum distributions that change with time, we adopt an alternative definition for the perturbed orbital velocity delta vphi in terms of the difference between the total instantaneous angular momentum flux, and the mass flux times the average angular momentum. The difference then represents the excess or deficit angular momentum transport due to orbital velocity fluctuations compared to the mean. Specifically,

(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

alpham = 2 BR Bphi / B2 .

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.