3. Results

3.1 Quasi-stationary state

The initial torus was in hydrodynamical equilibrium, but it evolves rapidly as the simulation proceeds. The overall evolution is very similar to that of the GT4 of Hawley 2000 (see Fig. 13 of that paper). The radial field is sheared by the differential rotation, creating toroidal field. Turbulence develops within the disk with the onset of the magnetorotational instability. The resulting Maxwell stresses drive angular momentum transfer; the disk begins to accrete into the central hole and expands radially outward. Relatively rapidly the inner portion of the disk achieves an approximate steady-state. By $t \simeq 700$, the mass accretion rate reaches roughly the long-time mean; after $t \simeq 1000$, the shape of the disk's average radial density distribution $\langle \rho\rangle (R)$ no longer changes appreciably. However, it must be borne in mind that there are always significant fluctuations around all these mean states. Thus, we observe the disk in a quasi-stationary state for several dozen orbital periods at the marginally stable orbit. In much of what follows, we will present a detailed analysis of the state of the disk during this time frame, focusing particularly on the end of the simulation, using that ``snapshot" as a fairly typical sample of conditions in the quasi-stationary disk.

3.2 Gas density

A contour plot of the azimuthally-averaged gas density at the initial and final times is presented in Figure 1. As might be expected, the density diminishes inside the marginally stable orbit and outside the initial outer radial boundary of the disk, while also decreasing with increasing altitude away from the equatorial plane. Although the vertical height of the highest contour levels at the final time are roughly independent of radius within the disk proper, the azimuthally-averaged exponential scale-height of the pressure, H, increases very nearly linearly with radius from $R \simeq 2$ out to $R \simeq 20$; at the marginally stable orbit (R = 3), $H \simeq 0.4$, while at R = 10, $H \simeq 1.5$. This particular disk shape is almost certainly a consequence of the equation of state, in which all fluid elements retain the same specific entropy except for shock heating. In a real disk, one in which turbulent dissipation introduces heat and radiation vents it, the shape might be rather different.

Figure 1 Figure 1: The azimuthally-averaged gas density distribution at the beginning and the end (t=1500) of the simulation. There are 14 evenly spaced contours. The maximum density at t=0 is 34.2; at t=1500, it is 14.3.

Density gradients in radius and altitude are determined primarily by the overall dynamical effects of approximate hydrostatic equilibrium. No such requirements exist for azimuthal fluctuations, whose amplitudes characterize the turbulence. We quantify the azimuthal fluctuation level by defining

\begin{displaymath} {\delta \rho \over \rho}\left(R,z\right) = {1 \over \langle... ...ft[\rho - \langle \rho \rangle_\phi \right]^2 \right\}^{1/2} , \end{displaymath} (15)

where ${\langle X \rangle}_\phi$ refers to the azimuthal average at fixed R and z of the quantity X. Measured in this way, $\delta\rho/\rho \simeq 0.2$ - 0.6 over most of the problem volume, with smaller values typically nearer the equatorial plane. In a few places, $\delta \rho/\rho$ rises to be greater than unity. Thus, quite sizable density fluctuations are generic.

Figure 2 Figure 2: Density in the equatorial plane for R < 10 at t=1500.

These fluctuations are primarily sheared strips whose radial extent is about 0.1R, but extend ~1 radian in azimuthal angle (Figure 2). This pattern (tightly wrapped spiral features) is characteristic of fluctuations in almost all quantities. Part of their nature is revealed by constructing spacetime diagrams in which the surface density at a specific radius is plotted as a function of t and $\phi $. From this visualization we are able to measure the angular speed of these features; it is generally quite close to the local orbital frequency $\Omega$. We can also measure their persistence at any particular radius; it is usually less than an orbital period. In this sense, they can be regarded as fluctuations that merely ``ride" with the local orbital motion, but see further discussion in §3.5.

3.3 Magnetic field

The magnetic field distribution is rather more complicated than the density distribution. A poloidal projection (i.e., averaged azimuthally) shows the magnetic pressure to have its largest value outside the dense regions of the torus, with large variations inside the torus (Figure 3). In the shearing-box simulations of Miller & Stone (2000), the field strength likewise peaked away from the midplane, but we find somewhat greater contrast. Field-strength fluctuations are also larger than those for density: the rms fractional fluctuations in the azimuthal direction are about 1 to 1.5 in the disk body, falling to ~ 0.25 several gas scale-heights immediately above and below the disk. Like the density, the characteristic shape of fluctuations is sheared arcs that run for ~1 radian before losing coherence. However, the peaks in magnetic field strength tend to lie adjacent to, rather than coincident with, peaks in gas density.

Figure 3 Figure 3: The azimuthally-averaged magnetic pressure plotted on a logarithmic scale from 10-5 to 10-2. The mean value of magnetic pressure within the region where the density is greater than 0.01 is 7.1 × 10-4. Using the mean gas pressure this corresponds to an average plasma beta value of 13.

The structure of the field in the inner part of the disk is illustrated by Figure 4, the azimuthally-averaged poloidal magnetic field vectors. There is a complicated field pattern within the disk, which is surrounded by a stronger and more organized field of roughly dipolar character. Figure 4 illustrates the poloidal field, but in general, the toroidal field is several times larger; ratios anywhere from ~1 to ~10 are common.

As one expects from magnetic fields in shearing plasma, there is a mean correlation between BR and $B_\phi$ in the sense that the ratio of the Maxwell stress to the magnetic pressure, $\langle -2 B_R B_\phi /B^2 \rangle > 0$ and has (approximately) constant magnitude $\simeq 0.2$ - 0.3 in the body of the disk (see Figure 12). Inside the marginally stable orbit, this ratio increases slowly in magnitude as radial accretion ``combs'' the field out into radial lengths that are strongly sheared, increasing the correlation between BR and $B_\phi$. In fact, the Maxwell stress is so correlated with the average magnetic pressure $\langle B^2/8\pi\rangle$, that a figure showing the stress is almost indistinguishable from one showing the magnetic pressure.

Figure 4 Figure 4: Azimuthally-averaged poloidal magnetic field vectors in the inner region of the simulation overlaid on density contours. The length of each arrow is proportional to the magnitude of the poloidal field.

3.4 Accretion rate

Perhaps the single most basic quantity of interest is the accretion rate, defined in equation (10). The disk is never perfectly time-steady; vR varies substantially from place to place and from time to time. In fact, disk quantities tend to be quite nonaxisymmetric, and much of the final accretion inside rms takes place in a spiral flow. However, time-averaging can remove much of this variability, and Figure 5 shows the time-averaged accretion rate as a function of R along with the instantaneous values of M-dot at t= 1000 and 1500. The time-averaged accretion rate is nearly constant as a function of radius inside R=10, and equal to 5 in code units. We shall adopt this as the fiducial steady state accretion rate for the inner disk. Between R=9 and 14 there is inflow at a decreasing rate; outside R=14 the net flux is outward. In this simulation, the initial disk is entirely contained within the grid, and the outer part of the disk must move outward in order to ``soak up" the angular momentum received from the inner disk material. Over the course of the simulation the total disk mass decreases by 14%; about 90% of this has gone into the black hole.

Figure 5: The accretion rate (in code units) as a function of radius at times t=1000 and t=1500 (solid lines). The dashed curve is the time-average of the accretion rate between those two times. Note that negative accretion rate means outward flow. Figure 5

Figure 6 shows a detailed time-history of the accretion rate through the inner radial boundary. As can be seen from this figure, the accretion rate shows significant time variations after its initial growth phase. The nature of the time-variability can be studied using the Fourier power spectrum of the accretion rate during the latter portion of the simulation (beyond t=600, i.e., after theinitial growth phase is over). This is shown in Figure 7. It is most simply described as a ``red" power-law; roughly speaking, $P(f) \propto f^{-1.7}$. However, there is also some curvature, in the sense that it steepens with increasing frequency. Because the accretion rate does appear to have a well-defined mean value, $d\ln P/d\ln f$ must be > -1 at the lowest frequencies. At the high frequency end, the spectrum is largely featureless, extending almost to the inverse free-fall time at the innermost boundary. Other than a smooth steepening, there is no special feature at the orbital frequency of the marginally stable orbit (marked on the figure). The flow is unsteady everywhere, even well inside the plunging region.

Figure 6: The mass accretion rate through the inner radial boundary as a function of time. To capture the time-variability accurately, the accretion rate was measured every 10 timesteps in the simulation, corresponding to an average interval of 0.0037 in time. Figure 6

Figure 7: Fourier power density of the accretion rate through the inner boundary (solid curve) and volume-integrated Maxwell stress (dashed curve). Frequency is in code units; the orbital frequency (i.e., the inverse of the orbital period at the marginally stable orbit is marked with an arrow. Figure 7

Figure 6 provides some idea of how nonsteady the accretion is; the mass flux is highly spatially inhomogeneous as well. To gain some sense of the magnitude of the spatial variations, we show the accretion rate as a function of position on a cylinder at R= 5 (Figure 8). As this figure shows, the fluctuations are very large. If we suppose that the mean accretion rate of 5 is evenly distributed over the range of altitudes -l < |z| < 1 the mean mass flux would be ~0.1 As the figure makes plain, there are several places where the mass flux is five times as great as this; there are also several places where the mass flux is ~0.3-0.4 outward. Relatively large fluctuations compared to the mean are to be expected when the transport of angular momentum in the disk is by turbulence.

Figure 8 Figure 8: The accretion rate at cylindrical radius R=5; inward mass flow is given a positive sign. Note that the local mass flux is almost as likely to be outward as inward.

3.5 Propagation of fluctuations

As we have stressed all along, one of the most striking features of this simulation is the presence of large amplitude nonaxisymmetric modulations in the density, mass flux, and magnetic field strength. We have already discussed how (,t) spacetime plots enable us to measure their angular velocity. It is possible to learn more by constructing analogous plots in (R,t). For example, in Figure 9 we show the history of the azimuthally- and vertically-averaged mass flux (eq. ) in this simulation. The diagonal stripes show the radial motion through the disk of fluctuations in the mass flux, their slopes indicating their radial group speeds. Typically these disturbances begin near R  ~   5 and travel both outward and inward. In the outward direction, they move at a very nearly constant speed ~0.7 while in the inward direction, their speed is at least three times faster. Because the outward speed very nearly matches the magnetosonic speed, we infer that these disturbances propagate as magnetosonic waves in the outward direction. On the other hand, inward-going mass flux perturbations are simply advected along with the dynamical inflow.

We conjecture that their origin lies in the way in which material enters the region of unstable orbits. When a fluid element begins to accelerate inward, it creates a rarefaction wave that moves outward at the magnetosonic speed, while traveling azimuthally at the (faster) orbital speed. On the other hand, the radial inflow speed near and inside the marginally stable orbit is comparable to the magnetosonic speed, so fluctuations travel inward rather more quickly. These fluctuations are never smoothed out into a completely steady flow because individual fluid elements move into the plunging region as a result of the continually changing torques exerted by the MHD turbulence (see the next subsection). Thus, at this level of detail, there is no way to achieve a stationary state.

Figure 9 Figure 9: Vertically- and azimuthally-averaged mass flux as a function of radius R and time. The diagonal stripes denote radially-propagating waves. Reds and yellows correspond to accretion, blues to outflow.

Although strong nonaxisymmetric pressure waves were generated by MHD turbulence in the shearing box simulations, they had no significant impact on the evolution of those simulations. In this global simulation, however, as Figure 9 shows, they play a significant role in the creation of large fluctuations in the accretion rate.

3.6 R- Stress

The total stress must satisfy the equation of angular momentum conservation. In a steady-state disk, this means that locally the torque must be equal to the mass accretion rate (which is constant) times the difference between the local specific angular momentum $\ell$ and the specific angular momentum carried into the black hole $\ell(r_{g})$. Thus, the stress is equal to

\begin{displaymath} \bar S = {\dot M \Omega(R) \over 2\pi} \left[1 - {\ell(r_g) \over \ell(R)}\right], \end{displaymath} (16)

where $\dot M$ is the total mass accretion rate, $\Omega(R)$ the rotational frequency (for the pseudo-Newtonian potential, $\Omega(R) = 1/[R^{1/2}(r-1)]$), and $\ell$ is the specific angular momentum ( $\ell \equiv R^2 \Omega$). In the time-steady, axisymmetric ``Novikov-Thorne" disk, it is assumed that $\ell(r_g) = \ell(r_{ms})$ because the stress is assumed to be zero for all $R \leq r_{ms}$.

For the purpose of studying angular momentum transport, it is useful to look at the azimuthally- and vertically-integrated magnetic stress, $\langle \int dz M^{R\phi}\rangle_\phi$, as this is the quantity most directly related to transfer of the z-component of angular momentum. In particular, we compare this quantity to the total stress per unit area that would be expected in the steady state ``Novikov-Thorne" disk. As discussed above, in the inner part of the disk the mean accretion rate is about 5 in code units. In Figure 10 $\bar S$ (for $\dot M = 5$ and $\ell(r_g) = \ell(r_{ms})$) is contrasted with $\langle \int \, dz \, M^{R\phi}\rangle_\phi$ averaged from t=750 to t=1500. As that figure shows, in the main body of the disk (approximately $4 \leq R \leq 15$), the average magnetic stress is almost exactly equal to the expected stress. At large radius the two curves diverge because equation (16) assumes a time-steady disk, whereas ours has only a finite mass, so the behavior at its radial exterior is quite different. At the inner edge of the disk, in contrast to the prediction of the conventional model, the azimuthally- and vertically-averaged magnetic stress does not diminish, but, if anything, increases slightly inward (the final upward spike is an artifact of the boundary condition). Note, however, that the local value of the ratio $\vert M^{R\phi}/(\rho v_R v_\phi)\vert$ in the plunging region varies widely: from as much as ~10 several scale-heights out of the equatorial plane to as little as ~10-4 at some locations in the plane.

Figure 10
Figure 10: The azimuthally-averaged and vertically-integrated magnetic stress averaged from t=750 until the end of the simulation (solid curve) and nominal mean stress per unit area (dashed curve).

In the context of time-steady and axisymmetric disks, the stress is often discussed in terms of the dimensionless parameter introduced by Shakura & Sunyaev (1973),

\begin{displaymath} \alpha_{SS} \equiv \langle {\int \, dz \, T^{R\phi} \over \int \, dz \, p} \rangle_\phi . \end{displaymath} (17)

It is possible to similarly calibrate the stress in this simulation if $\alpha _{SS}$ is generalized to be the local ratio of stress to pressure $\alpha \equiv \vert T^{R\phi}/p\vert$, which we will approximate by supposing that $T^{R\phi} \simeq M^{R\phi}$. Unlike $\alpha _{SS}$, which, in most studies, is taken to be a constant independent of time and location, $\alpha$ in this simulation is a dynamical and highly variable quantity that is determined by local and transient magnetic dynamics. For this reason, the absolute value is significant; although $M^{R\phi}$ is positive most places, there are counter-examples.

To give a sense of its scale and range of variation, we first examine $\langle \alpha \rangle_\phi$, the azimuthally-averaged version of this ratio, as found in the late-time snapshot (Figure 11). In the main portion of the disk, $\langle \alpha \rangle_\phi$ ranges between ~10-2 and ~10-1. However, it rises dramatically to ~10 at several scale-heights above and below the mid-plane. The same tendency has been seen (with much better vertical resolution) in stratified shearing box simulations (Miller & Stone 2000). $\langle \alpha \rangle_\phi$ also increases sharply at small radius. In the mid-plane, $\langle \alpha \rangle_\phi \simeq 0.2$ at R=3, but rises to $\sim 3$ at the inner radial edge of the simulation. Moreover, the vertical rise in $\langle \alpha \rangle_\phi$ also becomes sharper at small radius, so that at R=3, $\langle \alpha \rangle_\phi \sim 10$ at |z| as small as 0.5.

Figure 11: Absolute value of $\langle \alpha \rangle_\phi$, in a logarithmic scale in the snapshot at t=1500. This quantity increases with decreasing radius, and with increasing distance above or below the disk midplane Figure 11

In addition to these systematic trends, $\alpha$ is also subject to sizable fluctuations. Measured in the same way as for the magnetic energy density and the gas density (i.e., in terms of fluctuations in azimuth at fixed R and z), the rms fractional fluctuation amplitude in the bulk of the disk ranges from $\simeq 1.5$ to $\simeq 3$, but at large altitude above the disk midplane, the fluctuations can be much greater. The character of these fluctuations is illustrated by a different projection of $\alpha$, which we call $\bar{\alpha}$:

alpha (18)

the density-weighted vertical average of $\alpha$. In this projection (Figure 12), we clearly see the nature of the azimuthal fluctuations: just as for the density fluctuations (shown on smaller scale in Figure 2), the magnetic fluctuations are sheared filaments of limited angular coherence length.

Figure 12 Figure 12: Absolute value of the vertically-averaged $\bar{\alpha}$, plotted on a logarithmic scale for the snapshot at t=1500. Note the large fluctuations, even when values are compared at fixed radius. $\bar{\alpha}$ is large at both large and small R, but for different reasons (see text).

These results highlight the limitations of the primary assumption made in many previous studies of disk structure, that $\alpha_{SS}$ is a universal constant, independent of time and location. By contrast, we find that the ratio of stress to pressure varies strongly from place to place, both in the sense of having large fluctuations and in the sense of possessing strong systematic gradients. In many places the gas pressure is relatively smooth while the magnetic stress varies considerably. In other places the variation in $\alpha$ is simply due to changes in the background pressure. For example, the magnetic stress does not drop with height as rapidly as the pressure does.

For purposes of comparison with standard disk models, it is possible to compute (albeit with some trepidation) the value of $\alpha_{SS}$ found in this simulation. This quantity, defined in terms of the Maxwell stress, is shown as a function of radius in Figure 13 for the disk at the end of the simulation. It is $\sim 0.1$ through most of the disk, but, starting just outside the marginally stable orbit, it rises inward. It reaches a peak $\sim 1$ near the inner radial boundary of the simulation. The rise is primarily due to the rapid drop in gas pressure as the accretion flow accelerates inward. Also shown is the ratio of the Maxwell stress to the magnetic pressure (dashed line). This $\alpha _{mag}$ is 0.2-0.3 throughout most of the disk, and increases inside rms.

Figure 13: The Shakura-Sunyaev alpha parameter as a function of radius (solid line) constructed from vertical and azimuthally averaging the Maxwell stress and the gas pressure. Also shown is the averaged magnetic alpha as a function fo radius (dashed line), constructed using the magnetic pressure in place of the gas pressure. The magnetic alpha is relatively constant in the disk, whereas the Shakura-Sunyaev alpha varies with the gas pressure. Figure 13

Whether the Shakura-Sunyaev alpha is meaningful depends on how it is used. As the previous discussion emphasized, the local ratio of magnetic stress to pressure has strong vertical gradients that are completely obscured when the vertically-integrated version is used. Similarly, there are large fluctuations in azimuth and as a function of time. For these reasons, it is useless as a local dynamical variable. However, the approximate constancy of alpha as a function of radius in the disk proper means that it can be used as a rough estimator of the surface density in the body of the disk, provided its use is not extended to local properties nor is it thought to be accurate at better than ``factor of a few" accuracy. This conclusion is entirely consistent with many standard uses of and within the spirit with which it was originally introduced by Shakura and Sunyaev (1973) as a scaling parameter, not a fundamental variable.

3.7 Accreted angular momentum and energy

Another significant diagnostic of the flow is the mean specific angular momentum distribution. There are a variety of ways to characterize the specific angular momentum; one is by vertical and azimuthal average, $\langle \ell \rangle$, as defined in equation (9). This function was discussed in Hawley (2000) for the model GT4 (see Fig. 14 in Hawley 2000); the present model is very similar. At late time the average specific angular momentum is slightly sub-Keplerian throughout most of the disk. The most striking feature is that $\langle \ell \rangle$ continues to decline inside rms rather than leveling off. If all stress ceased at the marginally stable orbit, material would travel inward from there retaining precisely its energy and angular momentum. Conversely, the degree to which fluid changes its specific energy and angular momentum is a measure of the stresses that exist in the plunging region. The mean specific angular momentum diminishes from 2.60 at rms to 2.48 at the innermost radius. Thus, on average in this snapshot, the accreted matter has lost $\sim$5% of the angular momentum it had when it was at the marginally stable orbit. It is likely (for reasons explained in §5) that this number is an underestimate, and in reality the specific angular momentum at the innermost radius is rather smaller.

The total rate at which angular momentum is carried across the inner radial boundary is slightly smaller than these numbers would suggest because there is also a small outward electromagnetic angular momentum flux. This electromagnetic flux is due to a generic property of magnetic fields in shearing plasma: when the orbital frequency decreases outward, the field is swept back so as to transport angular momentum outward, i.e., $\langle B_R B_\phi \rangle < 0$. However, near the inner radial boundary this is a small effect, the electromagnetic angular momentum flux there is only about 2% of the angular momentum flux carried by the matter.

Further complexity is revealed by a closer examination of the (R,z) distribution of the late time azimuthally-averaged angular momentum inside the plunging region. In fact, the specific angular momentum in the cells lying along the equatorial plane is nearly constant inside rms. However, the value of $\ell$ drops with z above and below the equator. This is consistent with the distribution of magnetic field which, of course, accounts for the nonzero stress. The low value of field near the equator is in part a consequence of the assumed symmetry of the initial distribution of field. The initial field loops are symmetric with respect to the equator (Br and $B_\phi$ are zero there), and throughout the run the strongest fields tend to be located above and below the equatorial plane.

Just as for every other quantity, there is also considerable spiral structure. Figure 14 shows the density-weighted vertically-averaged specific angular momentum. It makes plain how nonaxisymmetric the distribution of torque is in the plunging region.

Figure 14 Figure 14: The density-weighted vertically-averaged specific angular momentum in the inner portion of the accretion flow.

For a number of reasons, the accreted energy is harder to determine than the accreted angular momentum (see further discussion in §5.2). To introduce this topic, we begin by displaying (Figure 15) the specific energy e in rest-mass units,

\begin{displaymath} ec^2 = {1\over 2} v^2 + {(\Gamma -1)P\over \rho } + {B^2\over8\pi \rho } - {1\over (r-r_g)}, \end{displaymath} (19)

in the equatorial plane at late time. The specific energy for a circular orbit precisely at rms is -0.0625; here, $\langle e \rangle_\phi$ is a bit greater than that ($\simeq -0.055$) because of the thermal, magnetic, and turbulent kinetic energy contributions to the total energy. In the equatorial plane, $\langle e \rangle_\phi$ is almost constant from R=3 to just outside R = 2, but actually rises as R decreases from there to R=1.6, reaching a peak that is very nearly unbound: -0.01. The azimuthally-averaged energy in the equatorial plane then falls sharply to -0.07. Even this wild swing underestimates the fluctuations: the smallest specific energy found in the equatorial plane is -0.19. Thus, we are confident in asserting that there is a large amount of energy transfer within the plunging region. However, given these large fluctuations, we do not feel that we can compute a well-defined mean energy for the accreted matter. It is likely, moreover, that at least some of the very large fluctuations in the innermost radial zones are not physical, but due to the inner boundary condition.

Figure 15 Figure: 15 The specific energy of matter (in rest-mass units) in the equatorial plane within R=4 (the high-resolution segment of the simulation). Note the large amplitude spiral fluctuations. The sharp features immediately surrounding the inner radial boundary are artifacts.

In fact, the situation is somewhat worse than this. The reason is that dissipation is crucial for computing the accreted energy, and our simulation contains no real dissipation physics. Absent numerical losses (see §5.1, 5.2), energy taken from plunging matter may go several places. One possibility is to build magnetic field energy in the plunging region itself. If the field energy is advected into the black hole without dissipation and radiative losses, the transfer of energy from fluid to magnetic field makes no difference to the accreted energy. It is also possible for magnetic stresses to transfer energy all the way back to the disk proper (and indeed we see significant magnetic stress throughout the plunging region). Again, without dissipation, the energy transferred in that way stays in the fluid, whether in the form of kinetic or magnetic energy. Only in the event that the energy is transferred all the way out to non-accreting matter would it be possible for that energy to avoid being carried back in. Thus, to the degree that this simulation avoids numerical losses, it is very difficult for it to show any significant transfer of energy out of the plunging region. By contrast, in a real disk, dissipation changes energy into heat, which may then be radiated if the thermal timescale is shorter than the inflow timescale; in our simulation, neither physical dissipation nor radiation is present.

Title Title Page   |   Method 2. Numerical Method   |   Discussion 4. Discussion