Hawley and Krolik: Discussion

Hawley & Krolik: Simulations of the Plunging Region

Title Page | 4. Initially Toroidal Field | 6. Conclusion

5. Discussion

5.1.1 Stress distribution

We suggest that the contrasting behavior seen in the toroidal and poloidal field simulations illustrates the effects of different inflow and magnetic dissipation times. The ratio between these two quantities distinguishes the disk proper from the plunging region. In the disk proper, the inflow time is long compared to the dissipation time so that turbulent dissipation controls the magnetic stress; on the other hand, in the plunging region the inflow time is short compared to the dissipation time so that the field and stress evolve primarily through flux-freezing.

The shearing-box approximation corresponds to the disk proper because the inflow time is effectively infinite. In those simulations (Hawley, Gammie & Balbus 1995, 1996), the magnetic field tends to saturate at a roughly constant value with respect to the local pressure. The stress is then a fixed fraction of the field energy density. We find just this sort of behavior in the portion of our simulations well outside r_ms, with -0.3 in the initially-poloidal simulation and somewhat smaller in the initially-toroidal case.

However, in our simulations increases steadily toward unity from to the inner boundary, indicating a qualitative change in the nature of the accretion flow at this radius. At small radii, the inflow time diminishes while the dissipation time (we expect) does not change greatly. Where the inflow time is shorter than the dissipation time, the field evolves via shear amplification and flux-freezing and is decoupled from the local pressure. The dominance of coherent shear over turbulence accounts for the increase in $\alpha _{mag}$ .

The contrast between the initially-poloidal and initially-toroidal cases arises from the fact that the inflow time at any given radius of the disk proper is longer in the toroidal simulation than in the poloidal due to the weaker magnetic stress. As a result, the zone in which the field energy density roughly tracks the pressure extends further inward when the initial field is purely toroidal. The magnetic stress in the toroidal run begins to evolve via flux-freezing at a smaller radius than in the poloidal run and (crudely) follows the diminishing pressure between R=6 and R=3, whereas in the poloidal run the magnetic stress continues increasing inward even when the pressure is decreasing.

In evaluating these contrasts, however, it is essential to remember that the dissipation rate in these simulations is controlled by numerical, not physical effects (see further discussion in §5.2); consequently, our results may not be quantitatively reliable guides to the behavior of real disks.

5.1.2 Which field topology is ``natural"

From these simulations it appears that the level of stress at and inside the marginally stable orbit depends significantly upon whether the initial field topology is toroidal or poloidal. Dependence of the level of MHD turbulence on initial field topology has been observed in previous simulations, both global and local. The linear analysis of the MRI shows that the fastest growing wavelengths of vertical fields are axisymmetric and have larger vertical wavelengths, while the toroidal field instability is nonaxisymmetric, favors small wavelengths, and grows transiently as time-dependent radial wavenumbers swing from leading to trailing. Given these differences it is important to ask which field topology better represents the state of an actual accretion flow.

A result that is common to all of the simulations done to date, both global and local, is that the toroidal field energy dominates the final turbulent state, whatever the initial field configuration. This does not mean, however, that the initial purely toroidal field evolution is the best model; it is in many ways a singular state, and there are sharp differences between it and models that have even a small amount of local net poloidal field. For example, Hawley et al. (1995) carried out several local shearing box runs in which a very much weaker vertical field was added to boxes with strong toroidal fields. The result was greatly enhanced turbulence. As noted by Hawley (2000) even if global simulations begin with no overall net poloidal field, but include some poloidal field in the form of field loops, local regions of the disk still behave like local simulations with poloidal fields, that is they generate strong turbulence and $\alpha \sim 0.1$ .

Hence, if an accretion flow begins from a source with even a small amount of poloidal field, and that field remains in or is amplified by the resulting inflow, one would expect local regions of the disk to behave more like the poloidal field simulation than like the pure toroidal field case. Given the strong influence of even a weak poloidal field, the inner disk dynamics could well be variable over long timescales due to variations in the net poloidal field carried in by accretion from the outer disk. Nonetheless, in a time-average sense we expect real disks to resemble the initially-poloidal simulation more than the initially-toroidal simulation.

5.2 Influence of resolution

The contrast between our new results and those obtained from the previous simulations of HK01 and H00 provides a measure of the effects of numerical resolution (Table 1). The grid used in HK01 is, in its inner parts, finer by a factor 1.9 in $\Delta R$ and 2.5 in $\Delta z$ relative to the grid of simulation GT4 in H00; our new grid is finer by a factor 3.3 in $\Delta R$ , a factor 2.4 in $\Delta z$ and 2.0 in $\Delta \phi$ relative to the simulation of HK01. We have also done a second initially-toroidal simulation with lower resolution than the one described in §4. It used 128 x 32 x 128 grid zones. The radial grid concentrated 55 zones linearly spaced inside R=4; outside this point $\Delta R$ gradually increased to the outer boundary at R=32. The vertical grid put 40% of the zones within $\pm 2$ of the equator, with a graded mesh beyond that point out to $z= \pm 8.8$ . Thus in the lower resolution toroidal simulation $\Delta R$ is twice as large, and $\Delta z$ three times as large in the inner disk region. This model was run for 400,000 timesteps to time t=3873.

It should be noted that, since this is a study of the disk structure near r_ms, we have focused the resolution improvements within the inner regions of the flow and around the equator. To the extent that the overall evolution depends on the outer parts of the disk, the effective resolution difference between simulations is far less than that of the inner region. The resolution in the outer disk of HK01 is actually poorer than that in the same portion of GT4.

5.2.1 Accretion rate

With two successive improvements in resolution by roughly factors of 3, the mean accretion rate in the initially-poloidal case increased 60% from GT4 to HK01 and a further 30% in the new simulation. The increase in accretion rate between HK01 and the new simulation might have been slightly greater if the new simulation had spanned a full $2\pi$ in azimuthal angle, rather than being restricted to only a single quadrant. Moreover, the fluctuations about the mean also increased substantially with improving resolution, with the peak to trough ratio increasing from $\simeq 1.1$ to $\simeq 2$ .

A similar rise in accretion rate with increasing resolution was seen in the pair of initially-toroidal simulations. At the lower resolution, the time-averaged accretion rate was $\simeq 0.7$ , fluctuating between a high of 1.7 and a low of 0.3. By contrast, with better resolution the accretion rate (after the end of transient effects) had a mean of 1.7 and fluctuated between 0.5 and 2.6.

These numerical experiments show that increasing resolution produces increased stress and accretion rates. For the poloidal field simulation the observed time-steady mean accretion rate is probably within tens of percent of the ``converged'' value, although almost certainly still on the low side. However, we have not yet achieved sufficient resolution to gauge the convergence of the magnitude of the fluctuations around the mean. For the simulation with the initial toroidal field, the increase in $\dot M$ with resolution was sufficiently great that it is not yet possible to estimate what the converged value might be.

5.2.2 Magnetic field

Between GT4 and HK01 there was a dramatic increase in magnetic field energy density in the inner disk: roughly a factor of 10. However, the improvement in resolution between HK01 and the new simulation, although comparable in magnitude to that between GT4 and HK01, changed the magnetic field much less: an increase of about 50% in mean energy density.

Analogous effects can be seen in the toroidal case. Growth to turbulent saturation takes almost twice as long, and the averaged poloidal field energies (a measure of the strength of the turbulence) in the saturated state are reduced by about a factor of 2 in the poorer resolution simulation.

Stronger magnetic field is mirrored in a rising ratio between magnetic stress and gas pressure. In the stable part of the accretion region (i.e., $5 \le R \le 10$ ), $\alpha_{SS} \simeq 0.05$ -0.1 in HK01 and 0.1-0.2 in the new poloidal simulation. The initially-toroidal simulations show an increase from $\simeq 0.01$ -0.028 in the lower-resolution simulation to $\simeq 0.02$ -0.03 in the higher-resolution run.

It appears very likely that improvements in resolution reduce numerical dissipation of magnetic field. Because of the central role of magnetic torques in driving accretion, the stronger magnetic field then increases the accretion rate. On the basis of the decline in the rate of increase of the magnetic field energy density with resolution in going from HK01 to our new poloidal simulation, it appears that, with regard to the mean field intensity, we are approaching convergence for the poloidal case, but it is likely that further refinement would yield at least modest quantitative changes.

We have less confidence that the initially-toroidal case is sufficiently well-resolved. As noted above, the toroidal field MRI favors high wavenumbers, precisely those that are the most challenging to resolve. The resolution comparison carried out here suggests that even with the higher resolution grid the simulation is not yet converged, and the values reported here must be regarded as lower bounds.

5.2.3 Stress distribution

As we discussed at length in §5.1, the ratio between effective magnetic dissipation time and inflow time plays a critical role in governing the magnitude of stresses in the inner disk. Where the dissipation time is short compared to the inflow time, the magnetic field strength is coupled to the local pressure; where the comparison is reversed, the field evolves by flux-freezing. At even our finest resolution, the dissipation rate is exaggerated by numerical effects. Consequently, in these simulations the magnetic field intensity is tied to the pressure further in than would be expected in a real disk where flux-freezing would set in at a larger radius. Thus, the stress observed in these simulations at small radii is most likely still undervalued.

5.2.4 Energy conservation

Energy conservation and energy flow can be difficult for a finite difference simulation to compute accurately when there are significant departures from equipartition in energy components. Here the bulk of the energy is gravitational and orbital; thermal energy is only 2% of the total at the beginning of the simulation, and 4% at the end. Thermal, magnetic, and poloidal kinetic energy combined are only 14% of the circular orbit binding energy at r_ms at the end of the simulation. Small errors in the gravitational and orbital energy can, in principle, then lead to large errors in the other elements of the energy budget. However, to answer the questions of interest, an accurate accounting of magnetic, thermal, and random kinetic energies is important.

In real disks, the energy of matter reaching the event horizon is determined by a competition between several rates: of energy transfer via work done on other matter and magnetic fields; of turbulent dissipation; of magnetic reconnection; of photon creation; and of photon diffusion. The rate of energy transfer via work is determined by the structure of the large-scale magnetic field and fluid motions. Turbulent dissipation is due to a two-part process: nonlinear wave interactions transfer energy from long wavelengths to short, where they can in turn be dissipated into heat by resistivity and viscosity. Magnetic reconnection is controlled by the rate at which fluid motions force together fluid elements with oppositely directed magnetic field. When the gas is heated, a variety of processes (bremsstrahlung, Compton scattering, etc.) transfer its heat into photon energy. Finally, energy leaves the disk as photons diffuse away.

Generally speaking, in these simulations we account reasonably accurately for large-scale energy transfer, overall disk evolution, and the long-wavelength portion of the turbulent cascade. However, numerical effects substitute for small-scale dissipation and enhance reconnection. Because the nonlinear wave interaction rate on long lengthscales is very rapid compared to the numerical dissipation rate, the effective turbulent dissipation rate in the simulation is controlled by numerical effects. Radiation we ignore entirely, with regard to both photon creation and photon diffusion. The absence of a real treatment of dissipation and photon losses hinders our ability to estimate the energy released by accretion. For example, in HK01 we argued that the observed change in energy between R=3 and R_min was likely to be an underestimate because, in the absence of radiative losses, some of the work done by plunging matter on gas farther out would eventually be brought back in.

However, at some level numerical effects mimic radiation by causing energy to be lost when converging fluid elements carry with them oppositely directed velocity or magnetic field. We can estimate the relative importance of numerical energy losses by a simple argument. If energy were exactly conserved locally, the only change in the integrated energy over the problem area would be due to material that flows off the grid. Here that effectively means material accreted because, for example, over the course of the poloidal simulation 19% of the initial mass passes through the inner boundary while only 1.5% leaves through the outer boundary. The binding energy of the accreted matter is, on average, about 0.118 in our units, so that mass-loss alone should increase the energy by 0.0224 if the initial mass is defined as unity. Starting from a total energy of -0.042, that means the final energy should be -0.0196. However, we actually find a final energy of -0.0272. The difference (a loss of 0.0076) can be attributed to numerical losses. Because the numerical dissipation time is shorter than the inflow time over much of the disk, this total energy loss is very close to the total radiation loss predicted by the conventional model, but there is no reason to believe that the specific location and rate of numerical losses in the simulation accurately predict what would happen in a real disk.

5.3 Stress at the marginally stable orbit

In these simulations we observe stress at the marginally stable orbit, yet it has been a long-standing belief that the stress must be close to zero at that point. As discussed in §1, at least two such arguments have been put forward in support of this notion. One is that the flow in the plunging region would necessarily be low density and would lack the inertia necessary to exert a torque on the disk (Novikov & Thorne 1973; Page & Thorne 1974). However, Page & Thorne point out that this argument doesn't apply if the torques are carried by magnetic fields.

The second argument relies upon the stress model, i.e., that the stress is . This point is worth revisiting in slightly more detail. Equation (7) relates the mass accretion rate to the stress in a steady state accretion flow. Using the $\alpha$ model for the stress, this equation can be rewritten

$\begin{displaymath} \alpha_{SS} (c_s/v_R) (c_s/R\Omega) = \Delta l/l(R) \end{displaymath}$

(8)

where $\Delta l$ is the difference between the local angular momentum l(R) and the angular momentum carried into the hole, l(r_g). Once the inflow passes r_ms, the inflow speed increases and the pressure drops precipitously, and with it the stress according to the $\alpha$ prescription. The $\alpha$ model is implicitly hydrodynamic, so the point at which plunging matter loses contact with the disk is the sonic point. There, by definition, $v_R \sim c_s$ , and $c_s/R\Omega = H/R < 1$ . Since $\alpha _{SS}$ is generally thought to be considerably less than unity, it would seem to follow that l changes very little in the plunging region. This is the basis for the ``proof'' that the stress goes to zero at r_ms.

In the simulations, however, we find a scaled stress, i.e. $\Delta l/l$ , at r_ms that is $\approx 0.05$ - 0.1. Why does the proof fail? Because the reasoning used relies on a strict application of the $\alpha$ relation, that is, that the pressure determines the stress. But the magnetic stress in the plunging region doesn't depend upon the gas pressure; $\alpha$ climbs rapidly inside r_ms as the gas pressure drops and the field evolves mainly by flux freezing. In fact, we find in the initially poloidal simulation that $\alpha _{SS}$ increases by a factor of 100 as matter falls through this zone.

Since the torque is magnetic, can we save relation (8) by replacing gas pressure with magnetic pressure, and setting the torque equal to $\alpha_{mag} \rho v_A^2$ ? Even with this readjustment, $\alpha$ still doesn't determine the stress in a manner that can be used in a proof; $\alpha _{mag}$ is also not a constant in the plunging region. It is roughly constant with value 0.2-0.3 in the disk proper where magnetic turbulence is active, but rises toward 1 in the plunging region. The stress in the plunging region becomes more efficient relative to the magnetic pressure as the radial field is stretched out in the inflow, and its strength relative to the toroidal field grows. Moreover, $v_A/(R\Omega)$ increases slowly inward in the flux-freezing regime, whereas in the conventional $\alpha$ model $c_s/(R\Omega)$ falls rapidly inward.

These problems illustrate the difficulties that can arise when the $\alpha$ formalism is taken too literally. The $\alpha$ parametrization cannot be used as the basis of a proof. At best, it might be possible to use the reformulated version of this argument using $\alpha _{mag}$ to estimate a time-averaged $\Delta l/l$ . In that case, the observed change in l of 5-10% inside the plunging region is consistent with the value of $v_A/(R\Omega)$ that we find ( for the poloidal field case, for the toroidal field case).

Title Page | 4. Initially Toroidal Field | 6. Conclusion