Hawley & Balbus: Nonradiative Black Hole Accretion

Introduction | 3. Three Dimensional NRAF Simulations

2. Numerical MHD Simulations: Background

Rotating, non-radiative accretion flows (NRAFs) can be numerically simulated. A number of axisymmetric hydrodynamical simulations have been performed (e.g., Igumenshchev & Abramowicz 1999, 2000; Stone, Pringle, & Begelman 1999), that drive accretion by employing an explicit kinematic viscosity, $\nu$ . The results depend strongly upon both the specific recipe adopted for $\nu$ , and its magnitude. For example, Igumenshchev & Abramowicz (1999, 2000) found that large viscosity flows accrete directly into the central hole. High viscosity accretion flows are therefore associated with the quasi-spherical ADAF flows. The low viscosity simulations (Igumenshchev & Abramowicz 1999, 2000; Stone et al. 1999), on the other hand, have a much more difficult time reaching the central hole. These flows have been interpreted as CDAFs.

MHD fluids, however, are fundamentally different from unmagnetized fluids, and essential features cannot be modeled with a Navier-Stokes viscosity formalism. Accretion simulations must be MHD.

The first NRAF MHD simulations were done by Stone & Pringle (2001; hereafter SP). These simulations were axisymmetric, and this is a limitation. First, the anti-dynamo theorem (e.g., Moffatt 1978) prevents the indefinite maintenance of a poloidal magnetic field in the face of dissipation. Indeed, toward the end of the SP simulations, the turbulence begins to die down, persisting only in flow close to the black hole. Second, axisymmetric simulations tend to over-emphasize the ``channel" mode (Hawley & Balbus 1992), which produces coherent streaming in the disk plane rather than the more generic MHD turbulence. Finally, the toroidal field MRI cannot be simulated in axisymmetry. Consequently, a fully self-consistent accretion simulation requires three dimensional MHD.

2.1 Equations

The simulations described in this paper evolve the three-dimensional equations of MHD: the continuity equation, the equation of motion, an internal energy equation, and the induction equation:

$\begin{displaymath} {\partial\rho\over \partial t} + \nabla\cdot (\rho {\bf v}) = 0, \end{displaymath}$

(1)

$\begin{displaymath} \rho {\partial{\bf v} \over \partial t} + (\rho {\bf v}\cdot... ...a \Phi + \left( {{\bf B}\over 4\pi}\cdot \nabla\right){\bf B}, \end{displaymath}$

(2)

$\begin{displaymath} {\partial\rho\epsilon\over \partial t} + \nabla\cdot (\rho\e... ... {\bf v}) = -(P+{\mathcal Q}) \nabla \cdot {\bf v} + \eta J^2, \end{displaymath}$

(3)

$\begin{displaymath} {\partial{\bf B}\over \partial t} = \nabla\times\left( {\bf v} \times {\bf B} -\eta_i {\bf J} \right), \end{displaymath}$

(4)

where $\rho$ is the mass density, $\epsilon$ is the specific internal energy, ${\bf v}$ is the fluid velocity, P is the pressure, $\Phi$ is the gravitational potential, ${\bf B}$ is the magnetic field vector, ${\bf J}$ is the current, ${\mathcal Q}$ is an explicit artificial viscosity of the form described by Stone & Norman (1992a), and $\eta$ is an anomalous resistivity of the form used by SP, namely

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

(5)

The constant C_res needs to be large enough to spread a current sheet out over a few zones, but not so large as to turn the overall flow into a resistive one. We use a resistivity constant of C_res = 0.1, as did SP.

The form of $\Phi$ is the pseudo-Newtonian gravitational potential introduced by Paczynski & Wiita (1980),

$\begin{displaymath} \Phi = - GM/(r-r_g), \end{displaymath}$

(6)

where M is the mass of the central black hole, and r_{g< /SUB> is the gravitational radius, equivalent to
the Schwarzschild radius in general relativity. With this potential,
the angular momentum of a circular orbit is}

l_Kep = (GMr)^1/2 r/(r-r_g),

(7)

and the binding energy is

$\begin{displaymath} e= -\left({GM\over2rc^2} \right)\left( {[r-2r_g] r \over [r-r_g]^2}\right) c^2. \end{displaymath}$

(8)

The pseudo-Newtonian potential mimics the dynamically important marginally stable circular orbit of the full Schwarzschild metric (defined by dl_Kep/dr = 0) at r=r_ms=3r_g.

The equation of state is $P=\rho\epsilon (\gamma -1)$ , with $\gamma=5/3$ . Radiation transport and losses are by assumption dynamically unimportant in an NRAF, and are omitted. There is no explicit shear viscosity; angular momentum is transported by Maxwell and Reynolds stresses arising from magnetic and velocity correlations in the MRI-induced turbulence. Similarly, the gas is not heated directly by an $\alpha$ -viscosity; nonadiabatic heating comes from the artificial viscosity ${\mathcal Q}$ and the resistivity. Both allow the entropy of the gas to increase.

We evolve the equations using time-explicit Eulerian finite differencing with the ZEUS algorithms (Stone and Norman 1992a,b; Hawley & Stone 1995).

2.2 Physical Units and Code Units

The results presented below are generally presented in terms of scale-free code units. Time is measured in orbital periods at the location of the pressure maximum of the initial torus, R=100r_g. (Here, R is the cylindrical radius.) This is 286 orbital periods at the marginally stable orbit, r_ms.

It is often convenient to have astrophysical scales associated with these values. Following Paczynski & Wiita (1980), we equate the gravitational radius r_g =1 with the Schwarzschild radius,

$\begin{displaymath} r_g = 3 \times 10^{5} \left( M/M_\odot\right) \ {\rm cm} \end{displaymath}$

(9)

where $M/M_\odot$ is the black hole mass. We also set GM=1, and the speed of light is c=(2GM/r_g)^1/2, or in code units $\sqrt 2$ . The orbital time at 100r_g is

$\begin{displaymath} \simeq 2\pi \times 10^3 \sqrt{2} r_g/c. \end{displaymath}$

This is defined to be $2\pi\times 10^3$ code units of time, so that one code time unit is $\sqrt{2}r_g/c\simeq 14 M/M_\odot$ $\mu$ s. For a $2.6\times 10^6 M_\odot$ hole, one orbit at R=100r_g is 2.3 x 10⁵ s; the entire simulation covers about 18 days in the life of the accretion flow.

In the absence of self-gravity and radiation, there is no density scale. The physical number density can be either specified outright, set by the mass of the initial torus, or assigned a value by specifying the accretion rate. The accretion rate can be scaled by the usual Eddington luminosity

$\begin{displaymath} L_{Edd} = 1.3\times 10^{38} M/M_\odot \ {\rm erg\ s^{-1}}, \end{displaymath}$

(10)

and the Eddington mass accretion rate of

$\begin{displaymath} \dot M_{Edd} = L_{Edd}/c^2 = 1.4\times 10^{17} M/M_\odot \ {\rm gm\ s^{-1}}. \end{displaymath}$

(11)

Since a nonradiating rotationally-supported gas is approximately virial in a black hole potential, the gas temperature will generally be of order the binding energy (8), a significant fr action of the rest mass energy (recall that $m_p c^2/k \simeq 10^{13}$ K). In the pseudo-Newtonian potential, the binding energy of the marginally stable orbit is 0.0625c². In our single fluid calculation there is only one temperature, T, and no distinction is made between ion and electron temperatures. A gas simulation such as this provides no constraints on the electron-ion interaction. However, if the ions provide all the dynamical pressure, it would be straightforward to consider the consequences of a two temperature plasma where the electron temperature T_e is some fraction $\delta$ of the ion temperature T_i.

Introduction | 3. Three Dimensional NRAF Simulations