5. Limitations

The results presented here and in Hawley (2000) obviously represent only the first steps in detailed simulations of global accretion disks. In this section we outline how the approximations we have made may impact our conclusions. Of course the only definitive way to investigate the effects of these approximations would be with additional simulations, and this discussion should be viewed as a guide for future work.

5.1 Effects of resolution

Foremost among the limitations is the numerical resolution of the simulation. Although we placed 30 radial zones between R=1.5 and R=4, that is probably insufficient. There are two reasons for believing this to be the case. First, we may apply the standard test for convergence of comparing two simulations with the same physics but different resolution, and ask whether they show significant differences. Simulation GT4, reported in Hawley (2000), and a new 643 run (hereafter referred to as run LR for Lower Resolution) serve as good comparisons. Both GT4 and LR share identical initial and boundary conditions. GT4 had equally spaced grid zones over a distance of 20 in both R and z for $\Delta R = \Delta z = 0.156$. Simulation LR used equally spaced radial zones over a range of 30 for $\Delta R = 0.47$, and used half of the z grid zones centered around the equator within $z = \pm 2$, for $\Delta z= 0.125$. By contrast, here $\Delta R = 0.0833$ and $\Delta z = 0.0625$ in the inner region.

All three simulations are qualitatively similar in global properties, as measured by volume integration, but differ in certain other respects that are important to physical interpretation. The time-averaged accretion rate in GT4 is about 40% smaller than in the work reported here (although the time-averaged accretion rate in LR is only about 2% smaller). However, the most striking contrast is in magnetic field strength. Higher resolution results in larger magnetic fields, particularly in the poloidal components. Contrasting maps of the magnetic field intensity at similar late time for GT4 and our new simulation (Figure 16), we find that the field is consistently stronger in the new simulation everywhere inside R=5, rising to typically a factor of 3 greater in magnitude inside R=4. Significantly, outside R=5, the field is generally stronger in the old simulation (sometimes by more than a factor of 10), but that, too, is in keeping with the sense of the resolution contrast: GT4 has better resolution in the outer regions.

Figure 16 Figure 16: Logarithm of the ratio of the density-weighted mean magnetic field intensity in the new simulation to the same quantity in GT4 at analogous late times. Only the region inside R=10 is shown. Because these are different simulations, one cannot expect features to match up precisely. However, there is a clear general tendency for the field in the new simulation to be stronger than in the old simulation inside R=5 and weaker outside.

The field intensity is controlled by a trade-off between the turbulent MHD dynamo driven by the magneto-rotational instability and numerical dissipation (recall that there is no explicit resistivity in this simulation). Whatever the real physics of resistivity in accretion disks may be, it is unlikely that it is as large as the effective numerical resistivity induced by our modest spatial resolution. The contrast between the field intensity in the two simulations is therefore very likely explained by the smaller numerical resistivity of the better-resolved simulation.

Another piece of evidence supporting the belief that the magnetic field is weaker than it would be with better resolution is presented in Figure 17, where we display the magnetic field structure in the equatorial plane within R=4, plotted for every fourth grid zone. The planar component of the field is organized into long sheared loops. Because the loops are wrapped around one another, there are many places where the field direction reverses within 4 - 8 grid cells. This means that with a modest amount of flow convergence numerical effects could cause the field to decay on a relatively short timescale.

Figure 17 Figure 17: The planar component of the magnetic field in the equatorial slice of the high-resolution region of the simulation. The length of each arrow is proportional to the intensity of the magnetic field. For clarity, the field is shown in every fourth cell.

Finally there is the issue of the computational domain employed, in particular the location of the boundaries. Very little outflow leaves through the z boundaries (typically only a few percent of the accretion rate). The outer radial boundary was extended in the present simulation compared to GT4, which had the effect of keeping more of the disk on the grid, but otherwise had no obvious effect on the evolution. We expect that there is also little dependence upon the location of the inner boundary because the flow across that boundary is generally supermagnetosonic. Although we made no tests of the effect of varying rin, tests of this sort were made on the similar simulations in Hawley (2000), and it appeared to make little difference.

5.2 Energy conservation and physical dissipation rates

The artificial dissipation caused by insufficient resolution also has consequences for energy conservation. To very good accuracy, angular momentum is conserved in this simulation. The primary source of numerical dissipation, the truncation error associated with transport terms in the equation, does not alter the total angular momentum, although it may artificially change its distribution. Total energy, however, is different because the energy equation that is solved is the internal energy equation. Numerical dissipation can destroy kinetic energy in regions of strong velocity gradients and magnetic energy in regions of strong magnetic field gradients, but there is no provision in this simulation for putting that energy back into heat (except for dissipation in shocks through the use of an artificial viscosity). Both sorts of dissipation can also occur in real systems, but at rates that might be very different from what happens in the simulation.

In the body of the disk, the problems that numerical dissipation can cause are limited. The single largest items in the energy budget are potential energy and rotational kinetic energy; neither one is significantly affected by numerical effects. The thermal energy is only a small fraction of those energies, and this is a prime motivation for solving an internal rather than total energy equation, although our assumption of an adiabatic equation of state is somewhat artificial. Artificial energy losses occur in the dissipation of relatively small amounts of energy bound up in small-scale random fluid motions and in magnetic field. In fact, in this regard, numerical dissipation acts in a way that is consistent with dissipating energy at the short lengthscale end of the turbulent cascade.

However, near and inside the marginally stable orbit, detailed questions of energy transfer become significant. The contributions to the energy budget due to both turbulent motions and magnetic field grow. At the same time, their importance also grows. For example, it would be of great interest to determine whether plunging fluid elements can, by exerting magnetic stresses, transfer energy back to the disk proper. As we have seen, significant magnetic stresses exist throughout this region, so there is some possibility that this happens. However, the rate at which this may occur is regulated by the magnetic field strength and topology, and these are subject to significant artificial numerical dissipation. In addition, any energy that is delivered is likely to take the form of MHD turbulence; whether this turbulence is dissipated (and the energy radiated) quickly or slowly compared to the inflow is critical for evaluating the radiative efficiency that results. Because the inflow, thermal, and dissipation timescales all change significantly in the vicinity of rms, it is unlikely that numerical dissipation serves as a successful stand-in for physical dissipation and radiative losses in the way it does at larger radii. Introducing physical dissipation and radiation rates, and controlling extraneous dissipation, are therefore important improvements to be added to future simulations.

5.3 Duration

The timestep of the simulation is limited by the Courant condition. The most stringent limitation comes from the high orbital velocities and the radial plunging speeds in the inner part of the grid where the resolution is reasonably fine. The total duration of the simulation was 1500 time units, requiring almost 400,000 timesteps. As discussed in §3.1, this is long enough for a quasi-steady state to be established in the inner part of the disk. Since $\alpha \sim 0.1$, the simulation covers several accretion times in the inner part of the disk. Whether or not the results obtained are truly representative, must, however, be tested with simulations of longer duration. An evolution beginning with a nearly Keplerian disk that extends out to large radius could allow the inner disk region to establish itself self-consistently from a long term accretion flow.

5.4 Field Topology

What are the limitations due to the choice of initial disk and field topology? The initial poloidal field loops used here are amenable to rapid development of the MRI on large scales. Further, the presence of significant initial radial field quickly generates through simple shear amplification a toroidal field of near-equipartition strength. Thus, the initial field topology guarantees magnetically-driven evolution within a few orbits. In earlier work (Hawley 2000) it was shown that in this regard (generating magnetically-driven inflow after a few orbits), initially poloidal and initially toroidal field configurations are essentially equivalent. However, the symmetry inherent in this choice of initial field topology is partly responsible for the low field intensity near the equatorial plane.

Whether the field is contained entirely within the problem area, or, equivalently, whether there is zero net flux through the box, has a greater qualitative impact. Shearing box simulations (Hawley, Gammie & Balbus 1995; 1996) also found strong similarities between simulations with initial toroidal field and simulations with random initial fields. However, they found dramatically different behavior when there was net vertical field. Such field topologies tend to produce larger net angular momentum transport and more violent evolution. These effects have been observed in both shearing box (Hawley, Gammie & Balbus 1995) and global simulations (Matsumoto 1999).

5.5 Displacement current

The equations of MHD as we have written them (see especially equation 2) omit the displacement current term. This approximation is valid so long as the Alfvén speed vA < c. This condition is satisfied almost everywhere in this simulation, with the only exception a small region near the inner radial boundary ( $1.5 \leq R \leq 2$ and $1 \leq \vert z\vert \leq 2$). In the equatorial plane, where most of the mass flows, the largest value of $v_A \simeq 0.3$, and that occurs only at the inner radial boundary; everywhere outside R = 2 in the midplane, $v_A \leq 0.15$. Future simulations in which the displacement current term is included should be somewhat more accurate, but we doubt that this correction will cause any major changes.

5.6 Relativity

Obviously, it would be preferable to run the simulation using genuine general relativistic dynamics, rather than the approximation of Newtonian physics in a Paczynski-Wiita potential. We might hope that this approximation is qualitatively reasonable when attempting to mimic dynamics around a non-rotating black hole; after all, the most important aspect of this case is the existence of dynamically unstable orbits inside a critical radius and enforced inward radial motion within a still smaller distance from the center. However, it has no representation of the frame-dragging created by a rotating black hole, so it cannot be a good approximation to accretion dynamics around a Kerr black hole. Gammie (1999) and Gammie & Popham (1998) examine some of the interesting possibilities created by the Kerr geometry, but examination of those effects must await a code with true relativistic dynamics.

5.7 Equation of state

As we have remarked earlier (§3.7), one of the most dramatic features of this simulation is the coherent disturbances that propagate as magnetosonic waves from the vicinity of the marginally stable orbit outward through the disk. Their properties depend strongly, of course, on the sound speed in the disk. However, in this simulation, in which all gas is given the same specific entropy initially and then evolved adiabatically, the sound speed could be quite different from what it would be in a more realistic simulation. In a real disk, there is substantial heating from dissipation of MHD turbulence, and also substantial cooling via radiation. Neither is considered here. We expect, therefore, that properties that depend on pressure forces (such as the propagation of these waves, or the vertical thickness of the disk) are not well-described in this simulation.

This lack can, in principle, be remedied in future simulations. It would not be difficult to integrate a total energy equation along with the gas momentum (i.e., force) equation. Heating could be treated in a controlled way by inserting an explicit resistivity into the magnetic field evolution equation and an explicit viscosity into the fluid force equation, and placing corresponding terms in the internal energy equation. Radiative losses could be described with either an escape probability formalism, or, in the longer term, solution of the radiation transfer problem.

Another consequence of our assumed equation of state, combined with the particular value of specific entropy we assigned the gas, is that the disk is modestly thick. As a result, it begins from a non-Keplerian, partially pressure-supported state and radial pressure gradients drive some of its initial evolution. Magnetic stresses rapidly change the angular momentum distribution within the torus to nearly Keplerian, but without cooling the internal pressure continues to push matter away from the pressure maximum. Some accretion, therefore, is driven not by magnetic stresses but simply by the pressure gradient. At late time, however, the inner disk establishes a quasi-steady state and nearly all the accretion can be accounted for by the Maxwell stress. It is likely that the observed dynamics are reasonably representative of generic inner disk accretion.

Finally, although the present disk is moderately thin, with $H/R
\simeq 0.15$ in the inner region, it is conceivable that new phenomena might appear in a still thinner disk. Although our results provide no particular reason to predict this, it is possible that in a disk with still smaller H/R there might be modes with wavelength large compared to H yet still small compared to R. Questions suchas these could be investigated by contrasting the present work with simulations of the inner regions of more extended, thinner, initially Keplerian disks.

Title Title Page   |   Discussion 4. Discussion   |   Conclusion 6. Conclusions