Hawley & Krolik: Simulations of the Plunging Region

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
,
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
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 .
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 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 -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 .
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 rms (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
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
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-
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 (,
); 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: 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:
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 ``" 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
in non-symmetric 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
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
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 flux-freezing and the
strong correlation created by shear is only weakly diminished by
turbulence.

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

(7)

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