Hawley & Krolik: Simulations of the Plunging Region

Title Title Page   |   Introduction 1. Introduction   |   Poloidal Field 3. Initially Poloidal Magnetic Field

2. Numerical Method

As in past work (H00; HK01) we evolve the equations of Newtonian MHD in cylindrical coordinates $(R,\phi,z)$, namely
(1)


\begin{displaymath} \rho {\partial{\bf v} \over \partial t} + (\rho {\bf v}\cdot... ...la \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\epsilon {\bf v}) = -(P+{\mathcal Q}) \nabla \cdot {\bf v} \end{displaymath} (3)


\begin{displaymath} {\partial{\bf B}\over \partial t} = \nabla\times\left( {\bf v} \times {\bf B} \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, and is an explicit artificial viscosity of the form described by Stone & Norman (1992a). To model a black hole gravitational field we use the pseudo-Newtonian potential of Paczynski & Wiita (1980) which is

\begin{displaymath} \Phi = - {G M \over r-r_g}, \end{displaymath} (5)

where r is spherical radius, and $r_g \equiv 2GM/c^2$ is the ``gravitational radius,'' akin to the black hole horizon. For this potential, the Keplerian specific angular momentum (i.e., that corresponding to a circular orbit) is

\begin{displaymath} l_{kep} = (GMr)^{1/2} {r \over r-r_g} , \end{displaymath} (6)

and the angular frequency $\Omega = l/R^2$. The innermost marginally stable circular orbit is located at rms=3 rg. We use an adiabatic equation of state, , where P is the pressure, $\rho$ is the mass density, $\epsilon$ is the specific internal energy, K is a constant, and $\Gamma = 5/3$. Radiation transport and losses are omitted. Since there is no explicit resistivity or physical viscosity, the gas can heat only through adiabatic compression or by artificial viscosity which acts in shocks.

The code employs time-explicit Eulerian finite differencing. The numerical algorithm is that of the ZEUS code for hydrodynamics (Stone & Norman 1992a) and MHD (Stone & Norman 1992b; Hawley & Stone 1995). We set GM=1 and rg = 1 (so that $c = \sqrt{2}$), thus establishing the units of time and velocity. The circular orbital period at a radius r is $P_{orb} = 2\pi \Omega^{-1} = 2\pi r^{3/2} (r-r_g)/r$.

In this paper we increase the overall resolution within the disk itself and in the inflow region above what was used in H00 and HK01. The computational grid is laid out in cylindrical coordinates, with 256 x 64 x 192 zones in $R \times \phi \times z$. This represents only a 50% increase in the total number of zones used compared to HK01. In the present simulations, however, the zones are concentrated to increase the effective resolution in the most important regions of the flow. We locate more of the zones near the marginally stable orbit and around the equator, and double the angular resolution while decreasing the angular extent. The radial inner boundary is moved in to $R_{\rm min} =1.25$ and there are 110 equally spaced zones out to R=4. Compared to our earlier simulation, this scheme decreases the zone size $\Delta R$ in the inner region by a factor of 3.3. Beyond R=4, $\Delta R$ gradually increases; the remaining 146 zones extend out to R=36. The z coordinate is centered on the equatorial plane, and runs from -11 to +11. From z=-1 to 1 there are 76 equally spaced zones; again comparing to the earlier simulation, the $\Delta z$ around the equator is smaller by a factor of 2.4. Beyond $z=\pm 1$, $\Delta z$ gradually increases out to the top and bottom boundaries.

The angle $\phi $ spans the range from 0 to $\pi/2$ in 64 equally spaced zones; $\Delta \phi$ is half the size used in H00 and HK01. Although the resolution is improved over H00 and HK01, the domain is only one quarter as large. However, experiments (Hawley 2001) with ``cylindrical'' disks (no vertical gravity) found that reducing the angular domain from $2\pi$ to $\pi/2$ does not alter the qualitative features of the evolution, although it lowered the energy and stress levels by about 10%. Since the practical advantage of limiting the angular domain is great, we use it here and assume that the quantitative effects will be small.

The boundary conditions on the grid are simple zero-gradient outflow conditions; no flow into the computational domain is permitted. The magnetic field boundary condition is set by requiring the transverse components of the field to be zero outside the computational domain, while the perpendicular component satisfies the divergence-free constraint. The $\phi $ direction is, of course, periodic.

The initial condition for the simulations is the same torus used in HK01 and for model GT4 of H00. This is a moderately thick torus ( $H/R \simeq 0.12$ at the pressure maximum) with an angular velocity distribution $\Omega \propto R^{-1.68}$, slightly steeper than Keplerian. The angular momentum within the torus is equal to the Keplerian value at the torus pressure maximum at R=10. As before, the pressure and density at R=10 are Pmax=0.036 and $\rho_{max}=34$, while Porb(R=10)=179. For reference, the orbital period at the marginally stable orbit is Porb =21.8.

We consider two different initial field configurations: poloidal loops, as in HK01 and H00, and a purely toroidal field. These models are discussed in turn in §3 and §4.

Title Title Page   |   Introduction 1. Introduction   |   Poloidal Field 3. Initially Poloidal Magnetic Field