Accretion powers a wide range of energetic objects, systems ranging in size from low-mass X-ray binary systems to the disks surrounding supermassive black holes in active galaxies and quasars. Theoretical models of these systems must be, necessarily, highly simplified. One such simplification that is often employed is to assume a time stationary disk. But the rapid time-variability seen in space- and ground-based observations emphasizes that accretion disks are highly dynamic systems. These observations highlight the need to solve the fully time-dependent equations describing accretion disks, and that, in turn, requires numerical techniques.

Recently, time-dependent numerical simulations have been applied to the problem of accretion disk transport. The discovery of the magnetorotational instability (MRI) in accretion disks (Balbus & Hawley 1991) has elucidated the physical basis for this transport: magnetohydrodynamic (MHD) turbulence. At the same time, the development of practical three-dimensional MHD codes brought the study of disk MHD turbulence, transport, and evolution into the computational domain.

To date, most of the numerical simulations of the MRI have been carried out in a local disk approximation known as the shearing box (Hawley, Gammie & Balbus 1995, hereafter HGB). In the shearing box model, one considers a frame corotating with the angular velocity at a fiducial radius Ro. By restricting the computational domain to small excursions from this fiducial radius, one can reduce the geometry (but not the dynamics) to a simple Cartesian system while retaining tidal and Coriolis forces. Local shearing box simulations (e.g., HGB; Hawley, Gammie, & Balbus 1996, hereafter HGB2; Brandenburg et al. 1995; Matsumoto & Tajima 1995; Stone et al. 1996) have demonstrated that MHD turbulence is the nonlinear outcome of the MRI, and that outward angular momentum transport is its natural consequence. The total stress TR phiis dominated by the Maxwell, or magnetic, component rather than the kinematic, or Reynolds stress. In disks the amplitude of this stress is often parameterized with the Shakura-Sunyaev (1973) alpha formulation, TR phi=alpha P, where P is the pressure. In the local simulations typical stress values range from alpha ~ 0.1 to 0.01. These simulations are more fully reviewed in Balbus & Hawley (1998).

An important step in increasing the realism of numerical disk studies is to move from this local approximation to fully global disk simulations. The difficulties of three dimensional global simulations and, in particular, the demands such simulations place on computer hardware, make this step challenging. In general, the time- and length-scale differences between the extended radial scales of an accretion disk and the scale of the MHD turbulence are too great to be fully resolved. Nevertheless, global simulations are now possible for somewhat restricted ranges in radius, height, and angle.

Several global MHD disk simulations have already been done. One of the first, by Matusmoto & Shibata (1997), followed the evolution of a thick torus embedded in a weak vertical magnetic field. What transpires is highly dynamic: the outer layers of the torus slough inward in a process referred to as avalanche accretion. This avalanche accretion is the global consequence of the radial streaming motion (the ``channel solution'') which is the nonlinear manifestation of the vertical field instability in the local limit (Hawley & Balbus 1992; Goodman & Xu 1994). The vertical field winds up as it is twisted by the torus rotation, and material is accelerated upward and outward along the field lines.

In subsequent work, Matsumoto (1999) evolved constant angular momentum tori containing an initially toroidal field. The nonaxisymmetric MRI rapidly leads to turbulence. Magnetic energy is amplified and approaches a value of beta = P/(B2/8 pi) ~ 10. Again the onset and appearance of the nonlinear stage of the instability are understandable in light of what has been seen in local simulations. In this simulation, however, the MRI has a global consequence: the resulting MHD stress transports angular momentum outward, rapidly changing the angular momentum distribution from constant with radius to nearly Keplerian.

Recently, Armitage (1998) performed a global calculation of a cylindrical (no vertical structure) Keplerian disk using the ZEUS code (Stone & Norman 1992a,b). The initial magnetic field was Bz = sin(R)/R, with the sine function chosen to vary once over the central part of the disk. This initial state is unstable, and turbulence rapidly develops along with significant angular momentum transport; alpha values are on order 0.1. As in local simulations, the Maxwell stress is several times the Reynolds stress. The magnetic field in the turbulent disk is dominated by the toroidal component, and the turbulence itself is characterized by structures with small azimuthal wavenumbers m. Hawley & Balbus (1999) similarly computed cylindrical models of Keplerian disks and constant angular momentum tori. They found in all cases that the MRI rapidly developed leading directly to turbulence and significant angular momentum transport.

Global simulations such as these have become possible only recently with advances in available computer speed and memory. As computer capacities continue to increase, ever more ambitious simulations will become possible. The aim of this paper is to begin a systematic effort to develop such three-dimensional global MHD disk simulations. To a large extent, this paper will focus on the technical aspects of such simulations, and the gross dynamics of the MRI in global disks. What type of simulations are currently practical? What sort of diagnostics are useful? What are the advantages and limitations of simplifications such as axisymmetry or the cylindrical approximation? As these questions are addressed, however, several significant conclusions about accretion torus dynamics will be obtained.

The plan of this paper is as follows: Section 2 discusses technical aspects of the global simulations, including the equations, the numerical algorithm, test problems, initial conditions, and diagnostic procedures. In Section 3 the results from a range of simulations are presented. Two important simplifications are investigated: two-dimensional axisymmetry, and a cylindrical gravitational potential. These approximations reduce the computational complexity of global simulations, and a comparison with the full three-dimensional treatment allows their limitations to be assessed. Simulations of fully global three dimensional tori follow. Finally Section 4 summarizes the conclusions from this work.

Next: Global Simulations

Back to title page