Hawley & Balbus: Three Dimensional NRAF Simulations

Numerical MHD Numerical MHD   |   The NRAF Model 4. The Generic NRAF Model

3. Three Dimensional NRAF Simulations

A previous paper (Hawley, Balbus, & Stone 2001) presents a prototype MHD NRAF simulation. Here we investigate that simulation (designated F1) in detail, and extend it to later times. We also compare F1 with a simulation that includes resistive dissipation (simulation F1r). The motivation for this is that in CDAF models dissipative heating in the central regions is argued to be responsible for global flow structure (Narayan et al. 2000).

We begin with some technical details. These simulations use $128\times32\times 128$ grid zones in cylindrical coordinates $(R,\phi,z)$. The radial grid extends from R=1.5rg to R=220rg. There are 36 equally-spaced zones inside R=15rg, and 92 zones increasing logarithmically in size outside this point. There are 50 z zones equally spaced between -10 and 10 with the remainder of the zones logarithmically stretched to the z boundaries at $\pm 60$. The azimuthal domain is uniformly gridded over $\pi/2$ in angle. The radial and vertical boundary conditions are simple zero-gradient outflow conditions; no flow into the computational domain is permitted. The $\phi$ boundary is periodic. The magnetic field boundary condition requires the transverse components of the field to be zero outside the computational domain, while the perpendicular component satisfies the divergence-free constraint.

The initial condition is a constant specific angular momentum (l) torus with a pressure maximum at R=100rg, an inner edge at R=75rg, and an outer edge at R=153rg. The initial magnetic field consists of poloidal loops lying along isodensity contours. It has a volume-averaged $\beta$ value (=Pgas/Pmag) of 200. The average specific energy of the magnetofluid in the torus is -2.3 x 10-3c2, i.e., it is bound. Simulation F1 is run for 6.32 orbits at the initial pressure maximum, or 1,807 orbits at rms=3rg. Simulation F1r begins at a time equal to 3.5 pressure maximum orbits in F1, after a global flow has been established. It continues to a time of 6.38 orbits.

Our simulations, with or without resistive heating, show three principal flow components: (1) a hot, but rotationally-supported, disk extending down to the marginally stable orbit; (2) an extended, low density coronal backflow enveloping the disk; (3) a distinctive, jet-like flow near the hole, that emerges unambiguously in momentum plots. These features are illustrated schematically in Figure 1. We believe that these structures are real and robust, and that they are fundamental generic properties of NRAFs. Physical arguments underlying the origin of NRAF structure are presented in §4.1.

Figure 1

Figure 1: A schematic diagram of the NRAF showing its major dynamical structures. An (MRI) turbulent, nearly Keplerian hot disk is surrounded by a dynamic low density envelope. Near the marginally stable orbit, the accretion flow thickens into a small inner torus as the pressure rises, but is still primarily supported by rotation. A centrifugally-evacuated funnel lies along the axis and is surrounded by an unbound outflowing jet confined by magnetic pressure in the corona.

3.1 Simulation F1: Ideal MHD

The evolution begins as a thick torus, with H/R = 0.2. The complete evolution is represented in Figure 2 as a series of contour plots in azimuthally-averaged density. During the first orbit, the field in the torus is amplified both by the MRI and by shear, and the flow evolves rapidly. The MRI acts most effectively near the equatorial plane, where the field is predominantly vertical. Here the long-wavelength, nearly axisymmetric modes of the MRI grow rapidly. Significant accretion begins shortly after one orbit, when the magnetic energy in the torus has increased to $\beta \approx 2$-10.

Figure 2

Figure 2: Density contours at initial pressure-maximum orbits 0, 1, 2, 4, 5, and 6. There are 30 contours, linearly spaced between 0 and 1.2. (The maximum density in the initial torus is defined to be 1.) During orbits 1 and 2 the linear modes of the MRI can be clearly seen. Beginning with orbit 4, a thinner disk structure is established; this disk stretches out radially with advancing time. By orbit 6 the hot inner torus can be seen inside of R=10rg.

Between orbit 1 and 2, low-m (azimuthal wavenumber) spiral arms of gas accrete from the inner edge of the torus, forming a vertically thin and very nonuniform ``daughter disk.'' Strong magnetic fields surround the gas. Due to the initial field topology, a current sheet forms near the equator, which proves unstable to vertical oscillations. At later times, as more gas accretes from the initial torus, the disk fills out and thickens. Inside of R=100rg the NRAF forms a modestly thick, nearly Keplerian disk. Low density material is stripped off this evolving disk to create a backflow. Inside of R=100rg, the net mass flux is inward; for R>100rg there is net outflow.

By orbit 2, the mass accretion rate into the central black hole has approached its long-term average value, which is only few percent of the rate at which the torus feeds mass into the inner region. There are two requirements for the gas to accrete into the central black hole, and neither is easy to meet: (1) the angular momentum of the gas must be reduced to a value close to that of the marginally stable circular orbit; (2) the gas also must be relatively cool, or pressure and centrifugal forces will drive the gas back out. In this simulation this is much in evidence: accreting gas must pass through a narrow gap in the equatorial ``waist'' of the centrifugal barrier. When the gas is hot and geometrically thick, this is most difficult. Instead of accreting, much of the gas splashes off the centrifugal barrier near the hole, either contributing to the formation of a hot, intermittent torus inside of R=10rg, or creating a coronal backflow that adds to a growing low density envelope around the Keplerian core disk.

Figure 3

Figure 3: Evolution of the specific angular momentum. (a) The mass-weighted, vertical- and azimuthally-averaged specific angular momentum l is plotted as a function of radius at each orbital time at the pressure maximum. For reference, the dashed line depicts the value lKep corresponding to a circular orbit. At t=0 l is a constant (=10.1). After only one orbit of time the flow has established a nearly Keplerian distribution throughout the radial domain, which it maintains throughout the simulation. (b) Contours of azimuthally-averaged specific angular momentum at the end of the simulation.

As time advances, the equatorial inflow follows a nearly Keplerian angular momentum distribution, $l_{Kep} \propto R^{1/2}$ (see fig. 3). The value of l actually remains everywhere slightly sub-Keplerian, except at the innermost radii where it is slightly super-Keplerian (in the hot inner torus). At the end of the evolution, l is about $\mathrel{\rlap{\lower 3pt\hbox{$\mathchar''218$}} \raise 2.0pt\hbox{$<$}}$ 5% below lKep between $R\sim 30r_g$ and 100rg, but for $R\sim 10r_g$ l drops to as much as 15% below the Keplerian value.

Between orbits 4 and 5, the accretion flow matures and becomes reasonably steady, forming a moderately thick disk surrounded by a low density, highly magnetized atmosphere. Between orbit 5 and the end of the run (6.3 orbits), the properties in most of the disk do not change by large amounts. The exception to this general trend is the inner torus ($R \le 10r_g$), which sporadically depletes and reforms.

Near the equator, thermal pressure dominates magnetic, i.e., $\beta > 1$, whereas the surrounding region is strongly magnetized, $\beta < 1$ (fig. 4). The thermal pressure exponential scale height is $H\approx 7$-10rg, decreasing rapidly inside R=10rg. The total pressure (thermal plus magnetic) is much smoother, and has an exponential scale height $H\approx 20r_g$ at R=100rg; H decreases slowly inward. The gas density, vertically averaged over one gas pressure scale height, is nearly constant from R=50rg down to R=10rg, where it increases rapidly to a peak at R=5.5rg. Inside of R=10rg the disk resembles a thickened torus.

Figure 4

Figure 4: Contours of the logarithm of total pressure, magnetic plus gas, at the end of the F1 simulation. The shaded regions overlaid on the contours show where gas pressure exceeds magnetic, i.e., beta > 1. The bulk of the coronal envelope is magnetically dominated. Gas pressure dominates in the disk, the hot inner torus, and in the jet along the funnel wall.

The disk and coronal regions can also be distinguished in Figure 3b which shows contours of specific angular momentum. Within the disk the gas is approximately rotating on cylinders, characteristic of a nearly barotropic equation of state. In the corona, surfaces of constant specific angular momentum are swept back, consistent with conserved l along outflowing streamlines.

Figure 5 shows contours of entropy, $S \propto \ln (P\rho^{-5/3})$, and the entropy along the equatorial plane at the end of the F1 simulation. Note that because of nonadiabatic heating by the artificial viscosity, the entropy increases inward even in the absence of resistive heating.

Figure 5

Figure 5: Contours of azimuthally-averaged entropy (top), and along the equator (bottom) at the end of the F1 simulation. The values are normalized to the initial entropy of the torus, S(t)-Si.

Magnetic fields in the disk produce a significant Maxwell stress which results in turbulent angular momentum transport. Between orbit 4 and the end of the simulation, the ratio of the vertical- and $\phi$-averaged Maxwell stress to the gas pressure (a form of the $\alpha$ stress parameter) ranges from $\sim 0.05$ to as much as 0.2 inside of R=100rg. This is significantly higher than the values found in small local disk simulations with zero net field (e.g., Hawley, Gammie, & Balbus 1996). It is this vigorous outward angular momentum transport that is responsible for building up a Keplerian profile: an initially shallower angular momentum distribution, as is typical in a thick torus equilibrium, exports more and more of its angular momentum to the exterior, until much of the radial pressure support is eliminated. The result is a Keplerian profile.

The mass inflow (vR < 0) and outflow (vR > 0) rates through each cylindrical radius are computed as a function of time,

\begin{displaymath} \dot M = \int \int R \rho v_R dz d\phi . \end{displaymath} (12)

The instantaneous local inflow and outflow rates are primarily due to the mass flux produced by short term $\rho$ and vR fluctuations, and they are nearly equal. This is not a matter of systematic inward mass flow being cancelled by a systematic outward mass flow. Rather, this is intrinsic to the nature of turbulent flow. The long term net accretion rate, $\dot M$, is always smaller than its rms value. An average over time for the last orbit of the simulation shows that the inflow exceeds the outflow for R<90rg (fig. 6). The ratio of inflow to outflow is about 1.5 between R=30rg and 90rg, dropping toward unity as one moves inward until rms is reached, at which point the outflow rate plummets to zero. The time-averaged net accretion rate $\dot M$ is roughly proportional to R between R=10rg and 60rg, i.e., it is not constant. The net $\dot M$ can vary with R if either the flow is not steady, e.g., the mass in the disk is building up with time, or if there is flow out of the grid at the z boundaries. If we consider the computational volume bounded by R=90rg on the outside over the last orbit of time, we find that the total mass in this volume declines and that over half of the net loss is through outflow along the z boundaries. We conclude that much of the decline of $\dot M$ with radius can be accounted for by this outflow.

Figure 6

Figure 6: Accretion rate averaged over the last orbit of time. The solid line is the mass inflow rate, the dotted line the mass outflow rate, and the dashed line the difference, i.e., the net accretion rate. Between R=10rg and R=60rg the net accretion rate is nearly proportional to R.

The total mass on the grid has decreased by 8.7% by the end of the simulation. Roughly 51% of this leaves through the outer radial boundary, 44% through the upper and lower z boundaries, and the remaining 5% is accreted into the black hole. Integrating the mass flux at R=75rg over time shows that a net 14% of the initial torus mass has moved inward. Only about 3% of the material accreting inward from the initial torus ends up in the central black hole.

Much of the flow through the outer radial boundary is a consequence of the increasing specific angular momentum in the outer part of the torus. The flow through the z boundaries, however, is driven from the disk by thermal and magnetic pressure. Perhaps the most interesting part of this outflow is an unbound, high temperature hollow conical outflow, confined to the axis region by surrounding magnetic pressure (see fig. 4 and fig. 7). By orbit 5, 60% of the outflow through the z boundaries is in this conical wind. The remainder of the z outflow is bound, although its specific energy is greater than that of the initial torus. Whether or not this more widespread outflow can be an unbound wind may therefore be a question of initial conditions: if the flow had begun as a marginally bound gas, the outflow would be unbound.

The initial average specific energy for the gas is -2.3 x 10-3c2. After 5 orbits this has decreased to -2.5 x 10-3c2 as the gas remaining on the grid becomes more bound. The total thermal, magnetic, poloidal kinetic, and orbital energies of this gas all increase with time. The magnetic energy has increased the most, receiving 42% of the total net energy increase; 95% of this is in toroidal field. The thermal, orbital and kinetic energies comprise 25%, 21%, and 12% of the increase respectively.

Figure 7

Figure 7: Azimuthally-averaged momentum vectors in the inner region, R < 20rg, at the end time of simulation F1. The overlaid contour is density, showing the shape of the inner torus. The outflow along the funnel wall is clearly evident. This outflow is confined by magnetic pressure in the corona.

3.2 Simulation F1r: Resistive Heating

In this calculation, simulation F1 is rerun beginning at orbit 3.5 with the addition of an explicit resistivity,

\begin{displaymath} \eta_i = 0.1 (\Delta x_i)^2 \vert J_i\vert/\sqrt{\rho}. \end{displaymath} (13)

We refer to the resistive simulation as F1r. The resistivity enhances smallscale reconnection and returns the magnetic energy losses as bulk heating. A study of the effect of resistive heating is motivated in part by claims that such heating will lead to convective instabilities. Resistive heating and turbulent dissipation are physical processes hydrodynamical codes try to mimic using a Navier-Stokes viscosity. In hydrodynamic simulations, this $\alpha$ viscosity is argued to be the origin of convective heating and additional angular momentum transport (Narayan et al. 2000).

Full MHD simulations reveal almost no difference between the resistive F1r model and the nonresistive F1 model described above. There is no evidence of convective transport due to resistive heating. This result was already obtained by the 2D MHD simulations of SP, who compared runs with and without resistivity. Although the simulation outcomes will always depend somewhat on the specific implementation of the resistivity, grid resolution and the like, there are good reasons to believe that thermal convection per se is generally of little dynamical importance in nonradiating accretion flows (Balbus & Hawley 2002). They are dominated by MHD turbulence.

There are some minor differences between F1 and F1r that show up in the details. Immediately after switching on the resistivity, for example, the magnetic field energy in F1r declines, and there is a corresponding increase in the thermal energy. The changes are small. Integrated over the entire computational domain, the thermal energy increases by a maximum of 5% over run F1 at 4.5 orbits. As the simulation proceeds, however, the total thermal energy in F1r actually lags behind that of F1. By orbit 6, F1 has 17% more thermal energy than F1r. Model F1 builds up a hotter inner torus, in part because of a larger accretion rate into the inner region. Heating via the artificial viscosity in shocks is more important than heating by the explicit resistivity. The shocks that occur in the simulation are not large-scale, coherent, or time-independent. The nonadiabatic heating appears to arise mainly in small-scale, highly time-dependent dissipation of the turbulence.

Although resistive dissipation increases the overall thermal energy, it must do so at the expense of the magnetic field, which is driving the accretion flow in the first place. The integrated magnetic energy of F1 always exceeds that of F1r. It is 5% greater through orbit 4.5, climbs to 20% at orbit 5, and drops back to 5% by the simulation's end. During the last two orbits, the total magnetic energy ranges from about 50% to 70% of the thermal energy. Finite resistivity decreases the turbulence levels and therefore the accretion rate as well. The result is that at late times, the internal energy in the resistive run can be smaller in isolated patches, while remaining largely unchanged in others. Figure 8 shows the central temperature as a function of radius at the end of F1 and F1r. The temperature is roughly proportional to R-3/2 within the accreting disk for both models.

Figure 8

Figure 8: Temperature at the equator as a function of radius at the endpoints of simulation F1 (solid line) and F1r (dot-dash line). For reference the binding energy is also plotted (dashed line). The maximum binding energy of 0.0625 mp c2 occurs at R=rms =3rg.

Numerical MHD Numerical MHD   |   The NRAF Model 4. The Generic NRAF Model