Hawley & Krolik: Simulations of the Plunging Region

Title Title Page   |   Method 2. Numerical Method   |   Toroidal Field 4. Initially Toroidal Field

3. Initially Poloidal Magnetic 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 $A_\phi (R,z) = \rho (R,z) - \rho_{min}$, for all $\rho$ greater than a minimum value, $\rho_{min} = 0.1$. The poloidal field is then obtained from ${\bf B} = \nabla \times {\bf A}$. 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., $\beta=P_{\rm gas}/P_{\rm mag} = 100$.

Only three differences distinguish this simulation from the previous ones: the earlier simulations treated a full $2\pi$ 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 $\phi $ resolution; and the new simulation was run somewhat longer, to t=1838 rather than t=1500. This longer time corresponds to 84 orbits at rms, 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 well-resolved 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 $t \simeq 600$. 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 $H/R \simeq 0.21$ at R=10 and $H/R \simeq 0.15$ at R=3.

We begin quantitative consideration of this simulation with its mass accretion history. In the GT4 simulation the mean time-averaged 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, lower-resolution work, the typical peak/trough ratio was $\simeq 1.1$-1.2; in the new one this ratio approaches 2.

As before, we also examine the Fourier power spectrum of the time-varying 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 $\simeq -1.7$. In this new simulation, the power spectrum is quite similar, but exhibits slightly less curvature. From the lowest frequency ( $f \simeq 10^{-3}$ inverse time units) up to roughly the orbital frequency at rms (0.046 inverse time units) $\vert\hat{\dot M}\vert^2 \propto f^{-1}$; from that frequency up to $f \simeq 1$, 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 1a
Figure 1b

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 volume-integrated 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 power-law, but in this case $\vert\widehat{M_{r\phi}}\vert^2 \propto f^{-2}$ at high frequencies, while bending to a steeper slope at low frequencies. In shearing-box simulations, the power spectrum of the fluctuating Maxwell stress is also a power-law 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- $\phi $ 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 azimuthally-averaged energy density in magnetic field increased by 50% in the accreting portion of the disk ($R \leq 10$, $\vert z\vert \le 4$); 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 azimuthally-averaged, vertically-integrated 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 rms), 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 rms. In addition, the pressure gradient is slightly shallower when better resolution is employed.

Figure 2

Figure 2: The azimuthally-averaged magnetic pressure, , measured at the end-time 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 initially-toroidal simulation.

Figure 3

Figure 3: Vertically-integrated and azimuthally-averaged pressure and magnetic stress at time t=1490 in two initially-poloidal 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 dot-dashed curves show the magnetic stress and pressure, respectively, in the new higher-resolution simulation. Analogous data for the initially-toroidal simulation are shown in fig. 12.

A popular way of parameterizing the importance of the magnetic stress is through the Shakura-Sunyaev ``$\alpha$" parameter. In its original definition (Shakura & Sunyaev 1973), the disk is supposed to be azimuthally symmetric and time-steady, and is given by the ratio of the vertically-integrated R- component of the stress tensor to the vertically-integrated 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 $\alpha$ in non-symmetric MHD turbulence. One may compute it as a local quantity at each point in R, $\phi $, 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 $\alpha _{SS}$ as the azimuthal average of the ratio between the vertically-integrated magnetic stress and the vertically-integrated pressure, we find (see fig.  4) that is typically - 0.2 in the accretion portion of the disk that is well outside rms (i.e., ). Inside R=5, it rises sharply, reaching $\simeq 0.5$ 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 $\alpha _{SS}$ 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 $\alpha$ 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 $\alpha$ that yields the most constant value is stress divided by magnetic pressure, $\alpha_{mag} \equiv -2B_R B_\phi/B^2$. 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 $\alpha$ 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 $\langle B_R B_\phi\rangle \neq 0$. The value within the disk is typical for stress arising from MHD turbulence. However, $\alpha _{mag}$ increases through the plunging region as turbulence gives way to more coherent fluid flow. Here the field evolves primarily through flux-freezing and the strong correlation created by shear is only weakly diminished by turbulence.

Figure 4

Figure 4: Azimuthally- and time-averaged ratio of the vertically-integrated stress to: the vertically-integrated 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 time-average runs from t=1500 to the end of the simulation.

These results highlight the limitations of the $\alpha$ picture: the stress doesn't actually correlate particularly well with the pressure. Another drawback to parametrizing the stress as an averaged $\alpha$ 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 azimuthally-averaged stress can vary by more than an order of magnitude within plus or minus two gas scale-heights from the midplane; similarly, the vertically-integrated 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

Figure 5: Logarithm (base 10) of the azimuthally-averaged magnetic stress at a late-time 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 azimuthally-averaged stress is negative, i.e., it is attempting to transport angular momentum inward. However, the magnitude of the largest azimuthally-averaged 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 steady-state disk, the local vertically-integrated R-$\phi $ 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(rg), i.e.,

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

where $\dot M$ is the total mass accretion rate, $\Omega(R)$ 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 azimuthally-averaged magnetic stress very nearly matched the stress predicted by (7) for the time-averaged mass accretion rate in the body of the disk, but, in sharp contrast to the prediction of the conventional zero-stress model, maintained a high level from R = 5 to the inner edge of the simulation at R = 1.5.

In the higher-resolution 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 time-averaged, azimuthally-averaged, and vertically-integrated stress is scaled by $\langle \dot M \rangle\Omega/2\pi$ 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, $\Omega$ is the usual orbital frequency for a circular orbit outside rms, but inside rms it is the actual orbital frequency as determined by the angular momentum and radial position of the material. In the figure, $\Omega$ inside rms 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 time-steady 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 time-steady approximation. Angular momentum is being transferred into the outer disk, which is moving outward in response. By definition, in a time-steady 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 $\propto R^{1.3}$. This continuing importance of magnetic stress in the inner portions of the disk is a counterpart to the increase in $\alpha _{SS}$ that occurs at the same place.

Figure 6

Figure 6: Azimuthal- and time-averaged scaled magnetic stress as a function of radius (solid curve) contrasted with the prediction of the zero-stress 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 10-20%. 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 scale-height roughly twice the gas density scale-height, 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 10-20%. In nonstratified cylindrical disk simulations (Armitage et al. 2001; Hawley 2001) the decline of l inside of rms 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 rms and the present simulation supports that conclusion. Even with nonstratified disks, however, strong fields near rms can produce substantial torques and red uctions in l (Reynolds, Armitage & Chiang 2001).

Figure 7

Figure 7: Mass flux weighted azimuthally- and vertically-averaged 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

Figure 8: Net energy in rest-mass 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 rms is 0.0625.

Figure 8 illustrates another aspect of the character of the flow inside rms: 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

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 Title Page   |   Method 2. Numerical Method   |   Toroidal Field 4. Initially Toroidal Field