Title Page

2. Numerical
Method

4.
Initially Toroidal Field
The initial magnetic field topology for the the first of the two simulations reported here is identical to that of the GT4 simulation in H00 and the simulation presented in HK01: poloidal field loops lying along equal density surfaces in the torus. The initial condition for the magnetic field is set by the toroidal component of the vector potential , for all greater than a minimum value, . The poloidal field is then obtained from . The resulting field is fully contained within the torus. The strength of the magnetic field is set so that the ratio of the total integrated gas pressure to magnetic pressure is 100, i.e., .
Only three differences distinguish this simulation from the previous ones: the earlier simulations treated a full in azimuth, whereas this one considers only a single quadrant; the gridding scheme used in the new one offers an improvement of roughly a factor of 3 in R and z resolution, and a factor of 2 in resolution; and the new simulation was run somewhat longer, to t=1838 rather than t=1500. This longer time corresponds to 84 orbits at r_{ms}, and is reached in 560,000 timesteps. Two purposes are served by studying this new simulation: we examine the degree to which the higher resolution of the new simulation allows us to approach numerical convergence, and it provides a standard of comparison for the simulation reported in the next section whose initial magnetic field was purely toroidal.
Just as in earlier simulations, the poloidal field loops of the initial state have unstable MRI wavelengths that are wellresolved on the grid. As a result, the MRI grows rapidly and turbulence develops within the disk. The disk expands as the MHD turbulence redistributes angular momentum. When the inner edge reaches the marginally stable orbit, matter begins to accrete through the inner boundary of the simulation, into the ``black hole". The initial growth phase of total field energy ends at . For the remainder of the simulation, conditions in the disk exhibit a rough stationarity, but with very large fluctuations in almost every quantity. The disk is modestly thick with at R=10 and at R=3.
We begin quantitative consideration of this simulation with its mass accretion history. In the GT4 simulation the mean timeaveraged accretion rate was 4 in code units. In the HK01 simulation, the mean accretion rate was very nearly 5; in the new simulation, the mean accretion rate rises to (fig. 1a). In all cases the accretion rate is characterized by large fluctuations with time. However, the amplitude of the fluctuations is larger in the new simulation. In the previous, lowerresolution work, the typical peak/trough ratio was 1.2; in the new one this ratio approaches 2.
As before, we also examine the Fourier power spectrum of the timevarying accretion rate. We previously found that the power spectrum was a smooth function of frequency that gradually steepened towards higher frequencies; its mean logarithmic slope over two decades in frequency was . In this new simulation, the power spectrum is quite similar, but exhibits slightly less curvature. From the lowest frequency ( inverse time units) up to roughly the orbital frequency at r_{ms} (0.046 inverse time units) ; from that frequency up to , the logarithmic slope is (fig. 1b). A small local peak appears at around a period of 41, which corresponds to the circular orbit frequency at R=4.2. The significance of this peak is unclear; if it is a real feature, it might be interpreted as suggesting that the final ``plunge" into the black hole is controlled by events a short way outside the marginally stable orbit. This would not be surprising as the size of the fluctuations in the MHD turbulence in this region is of the same order as the local scale height and the turbulence timescale is comparable to the infall time. However, it is also possible that the peak is a statistical fluctuation.
Figure 1: Upper panel: Mass accretion rate at the inner edge as a function of time in the initially poloidal simulation. Lower panel: Solid line is Fourier power density per logarithmic frequency interval of the accretion rate into the black hole (i.e. . Dashed line is Fourier power in the same units for the volumeintegrated Maxwell stress. To avoid transients associated with the initial start up and linear growth phase, the spectrum is computed for time units.
The shape of the Maxwell stress power spectrum is different. It, too, may be described as a broken powerlaw, but in this case at high frequencies, while bending to a steeper slope at low frequencies. In shearingbox simulations, the power spectrum of the fluctuating Maxwell stress is also a powerlaw of index 2 over several decades in frequency around the local orbital frequency; it therefore appears that the dominant stress fluctuations in these global simulations are controlled primarily by local effects. Note that the observable fluctuating quantity, the luminosity, may be more closely related to the stress than to the accretion rate because it is tied to the dissipation rate. However, it is further modified by the distribution of photon escape times.
The origin of the larger accretion rate seen in this simulation is a stronger magnetic field (fig. 2) which produces a greater R stress. As we also found in HK01, the field strength is typically somewhat greater near the disk surface than in the equatorial plane. Compared to the HK01 simulation, the azimuthallyaveraged energy density in magnetic field increased by 50% in the accreting portion of the disk (, ); at larger radius the net flow is outward, and there is little mass or field at higher altitude. The stress is also larger than in HK01, both in absolute terms and as a fraction of the disk pressure. Figure 3 compares the azimuthallyaveraged, verticallyintegrated magnetic stress and gas pressure at the endtime of the HK01 simulation to the same quantities in this new simulation. In both cases, although the pressure has a clear maximum as a function of radius (as happens in almost every disk model due to the sharp inward drop in density as the radial velocity grows near r_{ms}), the magnetic stress increases monotonically inward. However, the magnitude of the stress is everywhere greater in the higher resolution simulation; it is 2.8 times larger at r_{ms}. In addition, the pressure gradient is slightly shallower when better resolution is employed.
Figure 2: The azimuthallyaveraged magnetic pressure, , measured at the endtime in the initially poloidal simulation, plotted on a logarithmic scale from 10^{3.5} to 10^{1.5}. Fig. 11 shows the same quantity in the initiallytoroidal simulation.
Figure 3: Verticallyintegrated and azimuthallyaveraged pressure and magnetic stress at time t=1490 in two initiallypoloidal simulations. The solid curve shows the R magnetic stress in the simulation of HK01; the dotted curve shows the pressure in the same simulation. The dashed and dotdashed curves show the magnetic stress and pressure, respectively, in the new higherresolution simulation. Analogous data for the initiallytoroidal simulation are shown in fig. 12.
A popular way of parameterizing the importance of the magnetic stress is through the ShakuraSunyaev ``" parameter. In its original definition (Shakura & Sunyaev 1973), the disk is supposed to be azimuthally symmetric and timesteady, and is given by the ratio of the verticallyintegrated R component of the stress tensor to the verticallyintegrated disk pressure. Although it remains useful to compare observed stress levels to the pressure, there is no unique way (or most physically significant way) to define in nonsymmetric MHD turbulence. One may compute it as a local quantity at each point in R, , and z or one may average it in any of a variety of ways. Different definitions and averaging procedures yield quantitatively distinct results. However, certain trends do show consistent behavior.
For example, as already suggested in Figure 3, there is a general rise toward smaller radii in the importance of magnetic stresses relative to pressure stresses. Defining as the azimuthal average of the ratio between the verticallyintegrated magnetic stress and the verticallyintegrated pressure, we find (see fig. 4) that is typically  0.2 in the accretion portion of the disk that is well outside r_{ms} (i.e., ). Inside R=5, it rises sharply, reaching near R=3, and approaching at the innermost edge of the simulation. In the earlier simulation, 0.08 between R=5 and R=10, rising inside R=3, but always at a lower level than the new simulation.
The ``hump" in around R=20 is also seen in the HK01 simulation, but with smaller amplitude. The hump arises because the magnetic pressure generally has a larger scale height than does the gas pressure. As a result, in the outer portion of the disk the stress also falls less rapidly with R than the gas pressure, creating a peak in outside the gas pressure maximum.
If one chooses to parameterize the stress by the total pressure, gas plus magnetic, the rise is less dramatic. The definition of that yields the most constant value is stress divided by magnetic pressure, . This has a value between 0.2 and 0.3 throughout the main part of the disk, and rises slowly toward 1 inside R=6. This definition of measures the degree of correlation between the toroidal and radial components of the field. Because the shear has a consistent sense, the turbulence is anisotropic and . The value within the disk is typical for stress arising from MHD turbulence. However, increases through the plunging region as turbulence gives way to more coherent fluid flow. Here the field evolves primarily through fluxfreezing and the strong correlation created by shear is only weakly diminished by turbulence.
Figure 4: Azimuthally and timeaveraged ratio of the verticallyintegrated stress to: the verticallyintegrated gas pressure, i.e., (top solid curve); the gas plus magnetic pressure (bottom solid curve); and the magnetic pressure, i.e., (dashed line) in the initially poloidal simulation. The timeaverage runs from t=1500 to the end of the simulation.
These results highlight the limitations of the picture: the stress doesn't actually correlate particularly well with the pressure. Another drawback to parametrizing the stress as an averaged is that doing so obscures its highly variable nature. There is tremendous irregularity in the distribution of stress with altitude and azimuth (fig. 5). At a fixed radius within the accreting portion of the disk, the azimuthallyaveraged stress can vary by more than an order of magnitude within plus or minus two gas scaleheights from the midplane; similarly, the verticallyintegrated stress may vary by comparable amounts as a function of azimuth. Often, but not always, the stress is greater near the disk surface than in its midplane.
Figure 5: Logarithm (base 10) of the azimuthallyaveraged magnetic stress at a latetime in the initially poloidal simulation. The largest, most consistent stresses are associated with the plunging region near and within the marginally stable orbit. In the black regions of the plot, the azimuthallyaveraged stress is negative, i.e., it is attempting to transport angular momentum inward. However, the magnitude of the largest azimuthallyaveraged negative stress is at most only in these units.
Another way to look at the radial stress distribution is
with the ``scaled stress.'' In a steadystate disk, the local
verticallyintegrated R
stress must be equal to the mass
accretion rate times the difference between the local specific
angular momentum l(R)
and the specific angular momentum carried into the black hole
l(r_{g}), i.e.,
where is the total mass accretion rate, the rotational frequency, and l is the specific angular momentum. the quantity in square brackets is sometimes called the ``reduction factor". In HK01, we found that the time and azimuthallyaveraged magnetic stress very nearly matched the stress predicted by (7) for the timeaveraged mass accretion rate in the body of the disk, but, in sharp contrast to the prediction of the conventional zerostress model, maintained a high level from R = 5 to the inner edge of the simulation at R = 1.5.
In the higherresolution simulation the stress behaves in a similar, but not quite identical, fashion. In particular, the finer resolution leads to increased stress near the marginally stable orbit. In figure 6, the timeaveraged, azimuthallyaveraged, and verticallyintegrated stress is scaled by in order to highlight its effective ``reduction factor", i.e. the amount by which the stress is diminished due to the outward angular momentum flux. Note that for this purpose, is the usual orbital frequency for a circular orbit outside r_{ms}, but inside r_{ms} it is the actual orbital frequency as determined by the angular momentum and radial position of the material. In the figure, inside r_{ms} is approximated by assuming that the angular momentum is constant in the plunging region; this approximation (as shown in the following paragraphs) depresses the result by about 5% in the mean.
The scaled stress represents the contrast between the local specific angular momentum and the accreted angular momentum per accreted mass. In a timesteady disk this quantity would approach unity at large radius, but in the simulation it is considerably smaller than 1 there because the finite outer edge in our mass distribution leads to a strong inconsistency there with the timesteady approximation. Angular momentum is being transferred into the outer disk, which is moving outward in response. By definition, in a timesteady disk the scaled stress must be zero at the innermost boundary; the conventional model predicts that it goes to zero at R=3 and stays zero at all smaller radii. In the simulation, the scaled stress tracks the conventional model fairly well between R=12 and R=5, but at smalle r radii, instead of going sharply to zero at R=3, it falls more gradually: between R=1.25 and R=5 it is . This continuing importance of magnetic stress in the inner portions of the disk is a counterpart to the increase in that occurs at the same place.
Figure 6: Azimuthal and timeaveraged scaled magnetic stress as a function of radius (solid curve) contrasted with the prediction of the zerostress boundary condition model (dashed curve). As in Figure 4, the time average starts at t=1500.
The result of this continuing stress is a transfer of both angular momentum and energy from the plunging region inside R=3 to the disk proper at greater radius. As shown in Figure 7, the mean specific angular momentum falls by about 5% inside the marginally stable orbit, while the mean binding energy per unit mass rises by about 1020%. These significant, but modest, changes in the mean hide the much larger amounts of angular momentum and energy transfer that can occur in individual fluid elements. Some individual fluid elements arrive at R=1.3 with binding energy almost twice the mean binding energy at R=3, while others pass the ``event horizon" with slightly positive net energy (fig. 8). Because the magnetic stress has a vertical scaleheight roughly twice the gas density scaleheight, the torque per unit mass felt by matter near the disk surface is greater than that expressed in the midplane; in consequence, the specific angular momentum is systematically smaller off the equatorial plane, often by 1020%. In nonstratified cylindrical disk simulations (Armitage et al. 2001; Hawley 2001) the decline of l inside of r_{ms } can be significantly reduced. A systematic study of the influences of gas pressure (i.e., H/R) and computational domain size carried out by Hawley (2001) suggests that stratification plays a dominant role in determining dl/dR inside of r_{ms} and the present simulation supports that conclusion. Even with nonstratified disks, however, strong fields near r_{ms} can produce substantial torques and red uctions in l (Reynolds, Armitage & Chiang 2001).
Figure 7: Mass flux weighted azimuthally and verticallyaveraged specific binding energy (solid curve) and angular momentum (dashed curve) at the end values at R=3 to emphasize the continuing change at smaller radii. The spikes in the binding energy just outside the inner radius of the simulation are artifacts.
Figure 8: Net energy in restmass units in a slice through the equatorial plane in the plunging region at late time in the initial poloidal field simulation. Note that in these units the binding energy of a circular orbit at r_{ms} is 0.0625.
Figure 8 illustrates another aspect of the character of the flow inside r_{ms}: as gas plunges toward the ``event horizon", transfers of energy and angular momentum between adjacent fluid elements become stronger and stronger. Orbital shear stretches initially coherent regions into spirals stretching roughly a radian in azimuth. Local contrasts become especially strong deep in the plunging region because gas there can no longer exert forces and torques on gas farther out.
Viewed in a poloidal slice (fig. 9), the simulation reveals a different effect: in the plunging region, gas is concentrated where the magnetic shear or current density is particularly great. Interestingly, this concentration is much less noticeable considered as a function of azimuth. Generally speaking the gas and magnetic pressures are anticorrelated in the disk (although within the disk the gas pressure is generally larger). The total pressure is much smoother than either the gas or magnetic pressure separately.
Figure 9: Magnitude of the electric current (color contours) compared with gas density (line contours) in a poloidal slice at , both in a logarithmic scale. The two quantities are very well correlated.
Title Page

2. Numerical
Method

4.
Initially Toroidal Field