2. Numerical Method

2.1 Equations

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.,

eqn1 (1)


eqn 2 (2)


eqn3 (3)


eqn4 (4)

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. The global disk code is written in cylindrical coordinates (R, phi ,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

\begin{displaymath} \Phi = - {G M \over r-r_g} \end{displaymath} (5)

where r is spherical radius, and $r_g \equiv 2GM/c^2$ 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

\begin{displaymath} \ell_{kep} = (GMr)^{1/2} {r \over r-r_g} , \end{displaymath} (6)

and $\Omega R^2 = \ell$. The innermost marginally stable circular orbit is located at rms=3rg.

We use an adiabatic equation of state, We use an adiabatic equation of state, $P=\rho\epsilon(\Gamma -1) = K\rho^\Gamma$, with $\Gamma = 5/3$, 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 rg = 1 (so that $c = \sqrt{2}$).

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 $\Omega \propto R^{-1.68}$, 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 $\rho_{max}=34$. The initial magnetic field is constructed by setting the toroidal component of the vector potential equal to the density inside the disk, $A_\phi = \rho (R,z)$, for all $\rho$ greater than a minimum value ( $\rho_{min} = 0.1$). The total integrated magnetic energy is then normalized to a value of $\beta=100$, 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 rms; 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 $R_{\rm min} =1.5$; this avoids both the coordinate singularity associated with the axis, and the singularity at rg. 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 $\phi $ is periodic and covers the full $2\pi$. 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 $z = \pm 10$. The grid cells are equally spaced in $\phi $. 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 643 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 $\vert z\vert \leq 1$ along the inner radial boundary and $R \leq 2$ 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 $P_{orb} = 2\pi \Omega^{-1} = 2\pi r^{3/2} (r-r_g)/r$ at significant points within the disk. In units where GM=rg=1, Porb at the pressure maximum (R=10) is 179, and at the marginally stable orbit, where the orbital frequency $\Omega = 1/(2\sqrt 3)$, the orbital period is 21.8. The local accretion timescale is set by the effective stress, namely $t_{acc} = R/\langle v_R \rangle \sim \ell/W_{R\phi} ( \gg P_{orb})$, where $W_{R\phi}$ is the R-$\phi $ component of the stress, $W_{R\phi} = \langle \delta v_R \, \delta v_\phi - {v_{AR} v_{A\phi} \rangle}$, and vA 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 rms.

2.2 Diagnostics

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 $\phi = 0$ 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

\begin{displaymath} \langle \rho\rangle = {\int \rho R d\phi dz \over \int R d\phi dz}, \end{displaymath} (7)

the average angular momentum,

\begin{displaymath} \langle \rho \ell \rangle = {\int \rho \ell R d\phi dz \over \int R d\phi dz}, \end{displaymath} (8)

and average specific angular momentum,

\begin{displaymath} \langle \ell \rangle= {\langle \rho \ell \rangle \over \langle\rho\rangle }. \end{displaymath} (9)

The net radial mass flux is

\begin{displaymath} \langle \dot M\rangle = \langle -\rho v_R \rangle = \int -\rho v_R R d\phi dz, \end{displaymath} (10)

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

\begin{displaymath} \langle -\rho v_R \ell \rangle = \int -\rho v_R \ell R d\phi dz, \end{displaymath} (11)

can be computed.

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

\begin{displaymath} T^{R\phi} = \rho \delta v_R \, \delta v_\phi - {B_R B_\phi \over 4\pi}, \end{displaymath} (12)

where $\delta v_R$ and $\delta v_\phi$ 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 $\delta v_\phi$ 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,

\begin{displaymath} \langle \rho \delta v_R \delta v_\phi \rangle = \langle \rho... ... \rangle - \langle \rho v_R \rangle \langle \ell \rangle / R . \end{displaymath} (13)

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 $M^{R\phi}$. The vertically- and azimuthally-averaged Maxwell stress is

\begin{displaymath} \langle M^{R\phi}\rangle = {\int (-B_R B_\phi / 4\pi) R d\phi dz \over \int Rd\phi dz} . \end{displaymath} (14)

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, $\phi$).

Title Title Page   |   Introduction 1. Introduction   |   Results 3. Results