format: text

text:
# Global General Relativistic Magnetohydrodynamic Simulations of Accretion Tori

## 1 Introduction

Numerical experiments are playing an increasingly central role in investigating accretion into compact objects. A particularly important application is accretion into black holes, either solar-mass black holes in X-ray binaries, or supermassive black holes in active galactic nuclei, quasars, and the like. The investigation of such systems requires the development of algorithms and simulation codes that solve the equations of magnetohydrodynamics (MHD) in full general relativity (GR). To date, however, the only global three-dimensional MHD black hole accretion simulations have been done using a pseudo-Newtonian potential. The pseudo-Newtonian potential is $`\mathrm{\Phi }=GM/(rr_g)`$, where the gravitational radius, $`r_g`$, is equated with the Schwarzschild radius (Paczyński & Wiita 1980). This potential reproduces several important aspects of the Schwarzschild potential, in particular the qualitative behavior associated with the innermost stable orbit. But the system is still fundamentally Newtonian, and cannot account for the effects of a rotating black hole. Fully relativistic simulations are necessary, and in this paper we present the first such three-dimensional MHD black hole accretion models.

Since the pseudo-Newtonian model is substantially simpler than full relativity, it has been the natural starting point for accretion simulations. We begin by briefly reviewing some of those results. Hawley (2000; hereafter H00) computed fully global, three-dimensional prototype models of magnetized accretion tori using Newtonian MHD. H00 compared a variety of approximations: the “cylindrical disk” limit (in which the vertical component of gravity is ignored), two-dimensional versus three-dimensional simulations, and Newtonian versus pseudo-Newtonian potentials. This work demonstrated in a global context that the magnetorotational instability (MRI; Balbus & Hawley 1991) drives the evolution of accretion disks by generating turbulence and transporting angular momentum through significant Maxwell stresses. Tori with near-constant angular momentum distributions, which served as initial conditions, are particularly unstable, and quickly evolve to near-Keplerian disks. It was further found that the Maxwell stress did not go to zero at the radius of the marginally stable orbit, $`r_{ms}`$, the point traditionally regarded as the effective inner boundary of the disk. Instead, the specific angular momentum of the gas continued to decline inside $`r_{ms}`$.

Questions of large-scale disk structure and evolution were addressed with global MHD simulations by Hawley, Balbus & Stone (2001) and Hawley & Balbus (2002). This work examined the accretion of gas from a torus at $`100r_g`$, and the subsequent formation of a disk interior to that point. In this model the accretion flow was assumed to be nonradiating. As a consequence, the resulting disk was vertically thick, with an angular momentum distribution that had a Keplerian slope. The flow also featured a substantial outward flowing magnetized corona. Further, the accretion rate into the near hole region was larger than the accretion rate through the equatorial opening in the centrifugal barrier. Gas “backed up” at this point, creating an hot torus in the inner region near the hole, and an outflow along the centrifugal wall surrounding the evacuated funnel. The primary elements of this model—thick, nearly Keplerian disk, magnetized coronal backflow, hot inner torus—appear to be generic features of nonradiative accretion flows.

The properties of the inner region of an accretion disk are worth particularly close scrutiny since the most extreme relativistic effects and the greatest energy release happen near $`r_{ms}`$. The properties of the inner edge were the focus of the work of Hawley & Krolik (2001, 2002) which used relatively high resolution simulations to study the inflow past $`r_{ms}`$ and into the black hole. Cylindrical disk studies investigating this region were also done by Armitage, Reynolds, & Chiang (2001). Krolik & Hawley (2002) examined the implications of the simulations for the inner edge of a black hole accretion disk. The essential points are as follows: an accretion flow must transition from turbulence within the disk to a plunging infall into the hole, and simulations find that the location of this transition is highly time-variable. Generally the transition begins outside of $`r_{ms}`$, but the full transonic plunging inflow is not fully established until some point inside of $`r_{ms}`$. The inner disk edge is not abrupt: average flow variables, including the stress, are relatively smooth, albeit time-varying, functions of radius. Magnetic stress continues to be important even within the plunging region. In fact, the magnetic stress becomes proportionally more important as flux-freezing increases the field strength relative to the gas pressure, and increases the correlation between the radial and azimuthal field components.

As remarked, all the above work was done within the context of a pseudo-Newtonian potential. Just how well does that potential model the most important behavior in accretion flows? When is that approximation sufficient and when does it fail? What unique properties in the accretion flow are created by the spin of the hole and frame dragging, and do these produce unique observational signatures? Answering questions such as these requires a fully relativistic treatment.

The development of new codes for general relativistic MHD is ongoing. Recent work includes that of Koide, Shibata & Kudoh (1999), Gammie, McKinney, & Tóth (2003), and our own effort (De Villiers & Hawley 2003; hereafter DH03). In DH03 we describe the development and testing of the general relativistic MHD code used in this paper. The hydrodynamic portion of this code has already been employed (De Villiers & Hawley 2002; hereafter DH02) to study the evolution of accretion tori that are subject to the Papaloizou-Pringle instability (Papaloizou & Pringle 1984) in the full Kerr metric. Here we turn our attention to the first in a series of global simulations of accretion tori in the Kerr metric using the equations of MHD. Our aim with these models is to lay the groundwork for future, more detailed, black hole accretion simulations.

In this paper we examine a relativistic version of the constant angular momentum torus model GT1 studied in H00. In §2 we write the equations of general relativistic MHD in the form that the code uses. We describe the constant angular momentum tori that serve as initial conditions for the models in this paper. We list a number of the physics diagnostics used in the simulations. In §3 we present the results of the three-dimensional global MHD simulations, specifically a Schwarzschild hole, a prograde, and a retrograde Kerr hole with near maximal spin. We contrast three-dimensional simulations with two-dimensional axisymmetric models. We also compare a simulation that includes the full range in $`\varphi `$ ($`0\varphi 2\pi `$) with the results from the more economical quarter-plane simulation ($`0\varphi \pi /2`$), to examine the emerging impression that quarter-plane simulations give qualitatively the same results, and are adequate for many studies (Hawley 2001; Nelson & Papaloizou 2003). We contrast the results from the Schwarzschild simulation with an equivalent pseudo-Newtonian model. In §4 we review and summarize our findings.

## 2 Equations of General Relativistic Magnetohydrodynamics

We are studying the evolution of a magnetized fluid in the background spacetime of a Kerr (rotating) black hole. We adopt the familiar Boyer-Lindquist coordinates, $`(t,r,\theta ,\varphi )`$, for which the line element has the form,

$$ds^2=g_{tt}dt^2+2g_{t\varphi }dtd\varphi +g_{rr}dr^2+g_{\theta \theta }d\theta ^2+g_{\varphi \varphi }d\varphi ^2.$$
(1)
We use the metric signature $`(,+,+,+)`$, along with geometrodynamic units where $`G=c=1`$. Time and distance are in units of the black hole mass, $`M`$. The determinant of the 4-metric is $`g`$, and $`\sqrt{g}=\alpha \sqrt{\gamma }`$, where $`\alpha `$ is the lapse function, $`\alpha =1/\sqrt{g^{tt}}`$, and $`\gamma `$ is the determinant of the spatial $`3`$-metric.

The equations of general relativistic ideal MHD are the equation of baryon conservation, the conservation of the fluid energy-momentum tensor, and the induction equation. While the equations themselves are determined, there are a large number of ways that they can be expressed and solved numerically. In our methods paper (DH03) we derive useful forms for these equations, along with a set of primary and secondary variables especially suited to numerical evolution. The first of these is the equation of mass conservation,

$$_tD+\frac{1}{\sqrt{\gamma }}_j(D\sqrt{\gamma }V^j)=0,$$
(2)
where $`D=\rho W`$ is the auxiliary density variable, $`\rho `$ the density, and $`V^j=U^j/U^t`$ is the transport velocity, $`U^\mu `$ is the four-velocity, $`W`$ is the relativistic gamma-factor, and $`U^t=W/\alpha `$. Spatial indices are indicated by roman characters $`i,j=1,2,3`$.

The internal energy equation,

$$_tE+\frac{1}{\sqrt{\gamma }}_j(E\sqrt{\gamma }V^j)+P_tW+\frac{P}{\sqrt{\gamma }}_j(W\sqrt{\gamma }V^j)=0,$$
(3)
evolves $`E=ϵD`$, the auxiliary energy variable, where $`ϵ`$ is the specific internal energy. The isotropic pressure $`P`$ is determined by the equation of state of an ideal gas, $`P=\rho ϵ(\mathrm{\Gamma }1)`$, where $`\mathrm{\Gamma }`$ is the adiabatic exponent. For these simulations we take $`\mathrm{\Gamma }=5/3`$.

The equation of momentum conservation,

$$_t\left(S_j\alpha b_jb^t\right)+\frac{1}{\sqrt{\gamma }}_i\sqrt{\gamma }\left(S_jV^i\alpha b_jb^i\right)+\frac{1}{2}\left(\frac{S_ϵS_\mu }{S^t}\alpha b_\mu b_ϵ\right)_jg^{\mu ϵ}+\alpha _j\left(P+\frac{b^2}{2}\right)=0.$$
(4)
is written in terms of the auxillary four-momentum, $`S_\mu =(\rho h+b^2)WU_\mu `$, where $`h=1+ϵ+P/\rho `$ is the relativistic enthalpy, $`b^\mu `$ is the magnetic field four-vector in the rest-frame of the fluid, and $`b^2=g^{\mu \nu }b_\mu b_\nu `$. The momentum is subject to the normalization condition $`g^{\mu \nu }S_\mu S_\nu =(\rho h+b^2)^2W^2`$, which is algebraically equivalent to the more familiar velocity normalization $`U^\mu U_\mu =1`$.

Finally, the induction equation is solved in the form

$$F_{\alpha \beta ,\gamma }+F_{\beta \gamma ,\alpha }+F_{\gamma \alpha ,\beta }=0,$$
(5)
where $`F_{\alpha \beta }`$ is the electromagnetic field strength tensor, written in terms of the constrained transport (CT) magnetic field variables of Evans & Hawley (1988),

$$^r=F_{\varphi \theta },^\theta =F_{r\varphi },^\varphi =F_{\theta r}.$$
(6)
We assume infinite conductivity (the flux-freezing condition), $`F^{\mu \nu }U_\nu =0`$. The induction equation then reads

$`_j\left(^j\right)`$ $`=0`$ $`(\nu =0),`$ (7)
$`_t\left(^i\right)_j\left(V^i^j^iV^j\right)`$ $`=0`$ $`(\nu =i),`$ (8)
where $`V^\mu =U^\mu /U^t`$ is the transport velocity with $`U^t=W/\alpha `$. The CT variables $`^i`$ are related to the magnetic four-vector $`b^\mu `$ by

$$b^i=\frac{^i}{\sqrt{4\pi }\sqrt{\gamma }W}+V^ib^t,$$
(9)
and

$`b^t`$ $`=`$ $`{\displaystyle \frac{W}{\sqrt{4\pi }\alpha ^2\sqrt{\gamma }}}\left[g_{rr}V^r^r+g_{\theta \theta }V^\theta ^\theta +\left(g_{\varphi \varphi }V^\varphi +g_{t\varphi }\right)^\varphi \right].`$ (10)
Note that the factor $`\sqrt{4\pi }`$ is incorporated into the definition of $`b^\mu `$.

For reference, the energy-momentum tensor, written using the primitive variables, has the familiar form

$$T^{\mu \nu }=\left(\rho h+b^2\right)U^\mu U^\nu +\left(P+\frac{b^2}{2}\right)g^{\mu \nu }b^\mu b^\nu ,$$
(11)
which readily identifies the magnetic pressure, $`P_{mag}=b^2/2`$.

The GRMHD code evolves time-explicit, operator-split, finite difference forms of equations (2)—(4) and evolves the magnetic field induction equation (8) using the CT algorithm and a modified method of characteristic technique as described in DH03. In this paper we evolve a fiducial constant-angular momentum torus containing poloidal field loops orbiting around a prograde rotating Kerr hole, a retrograde rotating Kerr hole, and a nonrotating Schwarzschild hole. The fiducial models use a $`128\times 128\times 32`$ grid in $`(r,\theta ,\varphi )`$. The azimuthal grid spans the quarter plane, i.e., $`0\varphi \pi /2`$, with periodic boundary conditions in $`\varphi `$. One test simulation was run with a full $`2\pi `$ azimuthal grid, using $`128\times 128\times 128`$ zones. The outer radial boundary is set to $`r_{out}=120M`$ in all cases; the inner radial boundary is located just outside the horizon, with the specific value depending on the location of the horizon as determined by the Kerr spin parameter $`a`$. The radial grid is spaced logarithmically to maximize the resolution near the inner boundary. The radial boundary conditions are zero gradient boundary conditions, where the contents of the active zones adjacent to the boundary are copied into neighboring ghost zones. The $`\theta `$-grid runs over $`0.02\pi \theta 0.98\pi `$, with reflecting boundary conditions enforced at the polar axis. This avoids the substantial reduction in timestep that occurs if the $`\theta `$ grid runs up to the coordinate singularity at the axis. Test simulations have found no significant differences resulting from the reduced $`\theta `$ grid. Because the accreting gas has angular momentum, it is excluded from the region near the rotation axis.

A density and energy floor is employed to ensure that the dynamic range between the vacuum and the accreting gas does not become too great. The floor is set to be $`10^7`$ the initial density maximum in the torus. The timestep is set to one half the minimum time required for light to cross one grid zone. Because the grid is fixed, the timestep too is fixed throughout the simulation.

### 2.1 Torus Initial State

For the first simulations we choose a particular initial state that has been studied extensively in previous pseudo-Newtonian simulations: the constant-angular momentum torus containing weak poloidal magnetic field loops. We adopt the definition of specific angular momentum $`l=U_\varphi /U_t`$. The initial axisymmetric torus is generated from the constant-$`l`$ axisymmetric GR hydrodynamic equilibrium equations for an adiabatic gas, as described in Hawley, Smarr, & Wilson (1984; see also DH02). A particular torus is determined by the choice of $`l`$, which defines a set of effective potential surfaces, and the specific binding energy of the effective potential corresponding to the torus surface, designated $`\left(U_t\right)_s`$. Alternately, one can specify the radius of the pressure maximum (where $`l`$ is equal to the Keplerian value), and the radius of the inner edge of the torus (which sets $`\left(U_t\right)_s`$).

The initial magnetic field is obtained from the definition of $`F_{\mu \nu }`$ in terms of the $`4`$-vector potential, $`A_\mu `$,

$$F_{\mu \nu }=_\mu A_\nu _\nu A_\mu .$$
(12)
Our initial field consists of axisymmetric poloidal field loops, laid down along isodensity surfaces within the torus by defining $`A_\mu =(A_t,0,0,A_\varphi )`$, where

$$A_\varphi =\{\begin{array}{cc}k(\rho \rho _{cut})\hfill & \text{for }\rho \rho _{cut}\hfill \\ 0\hfill & \text{for }\rho <\rho _{cut}\hfill \end{array},$$
(13)
where $`\rho _{cut}`$ is a cutoff density corresponding to a particular isodensity surface within the torus. Using (12), it follows that $`^r=_\theta A_\varphi `$ and $`^\theta =_rA_\varphi `$. The constant $`k`$ is set by the input parameter $`\beta `$, the ratio of the gas pressure to the magnetic pressure, using the volume-integrated gas pressure divided by the volume-integrated magnetic energy density in the initial torus. In all calculations reported here the average initial field strength is $`\beta =100`$. The constant $`\rho _{cut}`$ is chosen to keep the initial magnetic field away from the outer edge of the disk. Here we use $`\rho _{cut}=0.5\rho _{max}`$, where $`\rho _{max}`$ is the maximum density at the center of the torus. This choice means that the initial field loops will be confined well inside the torus.

### 2.2 Evolution Diagnostics

One of the challenges of global, three-dimensional simulations of complex phenomena is to extract and distill useful diagnostic information that accurately characterizes the most important physical properties of the system. One computes, of course, a complete set of variables, at all grid zones for all timesteps; choosing from this data what to evaluate and to save is the issue. Given the novelty of three-dimensional MHD GR accretion simulations, there is as yet no definitive set of diagnostic quantities. We anticipate that standards will develop and evolve as our understanding of the simulation dynamics of black hole accretion flows improves. In previous pseudo-Newtonian calculations (e.g. H00), the evolution was characterized by several reduced quantities extracted from the simulation data. The set of diagnostics presented here is larger than those originally discussed in H00, and serves a first attempt at creating a new set of reduced quantities in a relativistic framework.

Entire restart files are saved at periodic intervals, and these can be examined in detail after a simulation. The more difficult choice is what data to save at frequent time intervals as a running history of the simulation. We have used full dumps of individual variables, both fundamental and derived; for the present simulations we saved density, and magnetic and thermal pressure 30 times per orbit. We also computed a more extensive set of azimuthally-averaged variables, variables averaged on spherical shells, and volume-integrated quantities. We define the average of a quantity $`𝒳`$ on a shell at radius $`r`$ as

$$𝒳(r)=\frac{1}{𝒜(r)}𝒳\sqrt{g}𝑑\theta 𝑑\varphi $$
(14)
where the area of a shell is $`𝒜(r)`$ and the bounds of integration range over the $`\theta `$ and $`\varphi `$ grids. For these simulations we compute shell-averaged values of density, $`\rho `$, angular momentum, $`\rho l`$, gas pressure, $`P`$, and magnetic pressure $`b^2`$. From these quantities, we derive quantities such as the density-weighted average specific angular momentum, $`l=\rho l/\rho `$.

Fluxes through the shell are computed in the same manner, but are not normalized with the area. We evaluate the rest mass flux $`\rho U^r`$, energy flux

$$T_{}^{r}{}_{t}{}^{}=\rho U^rU_t+b^2U^rU_tb^rb_t,$$
(15)
and the angular momentum flux $`T_{}^{r}{}_{\varphi }{}^{}`$, saved as two separate parts, i.e., the fluid part

$$T_{}^{r}{}_{\varphi (\mathrm{F})}{}^{}=\rho U^rU_\varphi ,$$
(16)
and the magnetic part,

$$T_{}^{r}{}_{\varphi (\mathrm{M})}{}^{}=b^2U^rU_\varphi b^rb_\varphi .$$
(17)
Again, various quantities can be subsequently derived from these fluxes and shell averages.

Volume-integrated quantities are computed using

$$\left[𝒬\right]=𝒬\sqrt{g}𝑑r𝑑\theta 𝑑\varphi .$$
(18)
The volume-integrated quantities are the total rest mass, $`\left[\rho U^t\right]`$, and total energy $`\left[T^{tt}\right]`$. In all the simulations these integrated values are computed and saved $`30`$ times per orbit at the initial pressure maximum.

## 3 Results

### 3.1 Retrograde, Schwarzschild, and Prograde Black Holes

Our fiducial simulations consider and contrast three initial torus models: a torus orbiting a Schwarzschild black hole, and tori orbiting rapidly rotating retrograde and prograde Kerr holes. Although the initial conditions cannot be made identical, we choose parameters for the initial tori that keep the inner edge of the disk and the location of the initial pressure maximum roughly the same as the black hole rotation parameter is varied. Table 1 lists the general properties of the models, where $`l`$ denotes the specific angular momentum, $`r_{in}`$ the inner edge of the disk (in the equatorial plane), $`(U_t)_s`$ is the binding energy of the surface of the torus, $`r_{Pmax}`$ the location of the pressure maximum (also in the equatorial plane), $`T_{orb}`$ the orbital period at the pressure maximum in units of $`M`$. For reference we also list $`r_{ms}`$, the location of the marginally stable orbit, the values of $`(U_\varphi )_{ms}`$ and $`(U_t)_{ms}`$ for a particle in a circular orbit at $`r_{ms}`$ (Bardeen, Press, & Teukolsky 1972), and $`T_{orb(ms)}`$, the orbital period at the marginally stable orbit in units of $`M`$. All radii in the table are in units of $`M`$. The retrograde SFR model is an extreme Kerr hole, $`a=0.998`$, while the prograde torus orbits a rapidly rotating Kerr hole with the rotation parameter $`a=0.9`$ chosen so that the marginally stable orbit is outside the static limit.

Figure 1 shows the initial density profiles for these fiducial tori, plotted on the same spatial scale to illustrate clearly the relative size of the disks. Overlaid on the tori are contours of the initial magnetic pressure. The choice of vector potential density cut-off (eq. 13), $`\rho _{cut}=0.5\rho _{max}`$, yields initial poloidal loops concentrated in the immediate vicinity of the pressure maximum. Although this keeps the field well away from the steep drop-off in density near the inner edge of the disk, the initial confinement of the field to the core of the disk gives rise to a violent transient as the field is amplified and buoyantly escapes from the core.

The three models, SFR, SF0, and SFP, were evolved for a time equivalent to $`10.0`$ orbits at their respective pressure maxima. Figures 2, 3, and 4 show polar and equatorial slices through the dataset for the density $`\rho `$ taken at $`t=2.0`$, $`4.0`$, $`10.0`$ orbits. During the evolution the tori pass through three phases: (1) the initial rapid growth of the toroidal magnetic field by shear; (2) the initial nonlinear saturation of the poloidal field MRI; and (3) evolution by sustained MHD turbulence.

During the first phase, the shearing of the radial field quickly generates a significant toroidal field. The increase in magnetic pressure causes the inner region of the torus to expand, driving both accretion and outflow. The second phase, the violent nonlinear saturation of the MRI, begins at about orbit 2 and features strong radial motions. Panels (a) and (b) of each figure show plots at $`t=2.0`$ orbits, where this violent transient state is at its peak. A strong, thin current sheet is clearly visible on the equator. A bubble of high magnetic pressure straddles the equator, driving lower density gas towards the black hole, while just above and below this bubble, pincers of denser material are also driven towards the black hole. The equatorial slices at $`t=2.0`$ show a series of spirals inside $`r15M`$, that represent a well-organized flow of material towards the hole. In model SFR these spiral waves reverse direction near the hole due to frame dragging. This flow feature was also observed in the hydrodynamic simulations of DH02.

As gas is driven inward by magnetic pressure, it also expands upward, carrying with it significant magnetic field. This leads to the formation of a highly magnetized outflow which moves out above and below the torus. This backflow is a generic feature of geometrically thick, nonradiating accretion flows, and has been observed in previous pseudo-Newtonian simulations (Hawley & Balbus 2002; Hawley & Krolik 2001, 2002).

Through orbit 2 the flow retains a high degree of symmetry across the equator, but when the pincers reach the black hole, the near-horizon inflow and the inner part of the torus become turbulent, as can be seen at $`t=4.0`$ orbits in panels (c) and (d). At this stage, the outer regions of the torus are still relatively well ordered, and two large outbound magnetic bubbles of low density gas can be seen above and below the equatorial plane, near the outer edge of the disk. Shortly after this, the entire disk becomes turbulent, and panels (e) and (f), taken at $`t=10.0`$ orbits, are typical of this mature phase of the evolution. At this late stage, large spiral waves can also be seen moving through the outer regions of the disks.

The magnetic nature of the low density bubbles at $`t=2.0`$ and $`4.0`$ orbits seen in figures 24 is made clear by Figure 5, which shows plots of the density in the SF0 model at $`t=2.0`$ orbits. Figure 5a is at $`t=2.0`$ orbits and 5b corresponds to 4.07 orbits, a time when the bubbles are especially prominent. The figure shows a color plot of density with logarithmic contours overlaid with contours of constant magnetic pressure (dark lines). To highlight the edge of the bubbles, only contours of low magnetic pressure are shown (these range from $`10^4`$ to $`10^3`$ of the maximum magnetic pressure on the grid). The regions inside the density-minimum bubbles correspond to local magnetic pressure maxima.

The general evolutionary histories of the three models are similar: growth of a strong current sheet along the equator in the first orbit, initial saturation of the poloidal field MRI between the second and fourth orbit, and sustained turbulence thereafter. The three models also have important differences. The most apparent difference lies in the structure of the inflow near the horizon. At one extreme, retrograde model SFR has a slender inflow with decreasing pressure and density, whereas prograde model SFP develops a new pressure maximum region inside $`r5M`$. We refer to this structure as the mini-torus for its appearance when viewed in a polar slice. Model SF0 is an intermediate case. These differences can be understood as arising from the relative distances between the inner edge of the initial torus and the marginally stable orbit. This distance is the greatest for SFR (see Table 1) and the least for SFP.

A more detailed history of the evolution of the tori can be gained from Figure 6, a space-time diagram $`\rho (r,t)`$ for models SFR, SF0, and SFP. The first orbit in each model shows the result of increasing magnetic pressure: gas is driven in toward the hole and the density in the heart of the torus decreases. Outgoing waves result from backflow driven from the inner region of the torus by this increased pressure. After orbit 2 there are strong fluctuations within the torus, and a series of time-variable accretion events. In models SF0 and SFP there is a notable shift of the denser edge of the torus toward their appropriate marginally-stable orbits, and hence toward the black hole. In model SFR the denser material remains outside of $`r=10M`$ except for occasional accretion streams. After orbit 6 the evolution is less violent. The formation of the mini-torus in model SFP, which appears after this point in time as a high-density region inside of $`r5M`$, is especially obvious.

The spacetime diagram also reveals many outgoing waves. The first of these, which appear in the linear MRI growth stage (between 0 and 3 orbits) originate as the backflow from the accreting gas. From the diagram we can measure the speed of propagation of this flow at $`v0.07c`$. A prominent tilted strip, visible in each plot between $`t1`$ and $`3`$ orbits, is associated with outward moving material propagating from inner edge of the inflow along the surface of the disk; this feature is especially prominent in animations of the data. A series of weaker outgoing waves are visible at later times. These are often internal spiral pressure waves rather than bulk outflows. These waves are also visible in animations of the data, and originate deep within the turbulent accretion flow. In model SFP, for example, strong waves originate within the mini-torus.

A quantitative view of the density distribution is provided in Figure 7, which shows $`\rho `$ for each disk model at $`t=10`$ orbits. In each case the density begins to decline outside of the location of the marginally stable orbit. In model SFP $`r_{ms}=2.32M`$, and the density increases roughly as $`1/r`$ inside the initial torus’s inner boundary. The region of peak density, lying between $`r=3M`$ and $`5M`$, is the region referred to as the mini-torus.

Figure 8 shows the accretion rate through the innermost radial zones, $`\dot{M}=\rho U^r(r_{min},t)`$, normalized to the initial torus mass for each of the models. The violent transient that results from the onset of the strong nonlinear phase of the MRI is clearly visible as a sharp spike at $`t=2.2`$ orbits in all three curves. This spike corresponds to a strong radial streaming flow of short duration. This type of flow is typical for the early nonlinear stage of the MRI operating on a background poloidal field. Strong fluctuations continue for two orbits, but by the fifth orbit all models settle into a steadier overall accretion rate with fluctuations about the mean on short time scales. This figure shows that, as a function of initial torus mass, model SFR has the greatest accretion rate, and model SFP the least. Accretion is easier to accomplish in model SFR because the initial torus inner boundary is just outside the location of the marginally stable orbit. In SFP, on the other hand, the disk must work its way down almost to $`r=2M`$ before reaching the marginally stable orbit. This is a large enough distance that the flow must establish a new, extended turbulent disk structure. The mini-torus is the innermost region of that new disk. Comparing the three accretion rate curves also reveals that model SFR has less variability on short timescales compared to the other models. Variability in the accretion rate is determined by the rate at which matter is fed in from the turbulent disk. Fluctuations in $`\dot{M}`$, therefore, originate mainly in the turbulent disk, outside of $`r_{ms}`$, and the resulting time-dependence of $`\dot{M}`$ reflects the turbulent frequencies of the inner region of the disk. These frequencies are naturally lower in the SFR model. The accretion rates in these models were sampled only 30 times per orbit, however, which makes it difficult to quantify the variability, or to distinguish the variability seen in SF0 from that of SFP.

Table 2 accounts for the total mass in each model. The values are normalized to the initial mass of the SF0 torus. The total initial mass, $`M_0`$, and the mass at 10 orbits, $`M_{\mathrm{final}}`$, are obtained by volume integration. The mass that enters the black hole, $`\mathrm{\Delta }M_{\mathrm{in}}`$, and the mass that leaves through the outer boundary, $`\mathrm{\Delta }M_{\mathrm{out}}`$, are obtained by integrating the mass flux over time at the respective radial boundaries; the integral is approximate since the flux is sampled only 30 times per orbit. The fraction of initial disk mass that has left the grid through the inner boundary is given in the final column. Model SFR has lost the greatest percentage into the hole, and model SFP the least. The low mass loss through the outer grid boundary for all models means that the radial location of this boundary has had a minimal effect on the overall simulation. There are significant differences among the three models in the amount of gas moving outward, however. Model SFP has the greatest relative outflow; by orbit 10, 22% of the original torus mass has moved beyond $`r=40M`$, a radius outside the outer edge of the initial torus. In model SF0 this amount is 14%, and in SFR, 10%.

The space-time diagram of gas pressure $`P_{gas}`$ closely resembles that of the averaged density (fig. 6), and is therefore not reproduced here. The pressure decreases with time at the location of the original pressure maximum as the material in the torus is redistributed. The magnetic pressure, $`P_{mag}`$, grows rapidly in the early stages of the simulations ($`t<3`$ orbits), first primarily due to shear of the initial radial field, then due to the strong growth of the MRI. Once MRI-driven turbulence is established, the maximum of magnetic pressure shifts sharply inward with the infalling magnetized gas. In models SF0 and SFP the location of the gas pressure maximum also moves toward the black hole as the turbulent disk extends inward to the radius of the marginally stable orbit.

A quantitative view of the pressure evolution in the three models is given in Figure 9, where $`P_{mag}`$ and $`P_{gas}`$ are plotted as a function of time at the radius of the initial pressure maximum, and also near the horizon. In each case the scale is normalized to the initial pressure maximum. The time-evolution inside the torus is very similar in all three runs, particularly during the initial linear growth phase and nonlinear saturation. This indicates that the evolution of the MRI is mostly independent of the spin of the black hole even relatively close ($`r=15M`$) to the hole. Gas pressure dominates except for a brief period around $`t=2.0`$ orbits where the violent magnetic transient is at its peak. During the extended turbulent phase beyond $`t=4`$ orbits, $`\beta 5`$ to 10. Near the black hole magnetic pressure dominates, especially in model SFR. The two are more nearly equal in model SFP.

Where the flow shifts from gas to magnetic pressure dominance is determined by the location of the marginally stable orbit. This is illustrated in Figure 10, which shows the averaged gas and magnetic pressures as a function of radius at orbit 10. Outside of $`r_{ms}`$, the disk evolves due to turbulent stresses driven by the MRI, and gas pressure continues to dominate. Inside of $`r_{ms}`$, the flow transitions to plunging infall, and the field evolves primarily by flux freezing. The magnetic pressure becomes increasingly important, helping to confine and guide the flow within the plunging region. The point where the averaged magnetic pressure exceeds the gas pressure is located at about $`0.7r_{ms}`$ in both SFR and SF0. In SFP the inner grid boundary radius exceeds $`0.7r_{ms}`$, and the transition to magnetic dominance is not seen.

Figure 11 shows the azimuthally averaged gas pressure, $`P_{gas}`$, magnetic pressure, $`P_{mag}`$, and total pressure, $`P_{tot}`$, at $`t=10.0`$ orbits for model SF0 (the other two models have similar profiles). Within the disk, total pressure is relatively smooth, consistent with a thickened equilibrium torus, and dominated by gas pressure. Magnetic pressure is increasingly important at the inner edge of the disk, and gas and magnetic pressure are comparable in the low-density coronal outflow region surrounding the denser disk. Magnetic pressure dominates in the funnel region near the axis, the pressure contours are spherically symmetric, and the magnetic field is mainly radial. The material in the funnel consists of fast outward-moving gas without angular momentum, several orders of magnitude lower in density than the maximum density in the disk (and hence below the cut-off in the contour plots shown above). It is important to note that the funnel region is initialized to a numerical vacuum, i.e. the numerical density floor of the code, and never rises much above this value, so conclusions about the physical properties of the funnel region must be reached cautiously. The point is that material from the disk, which contains significant angular momentum, is kept well away from the funnel by the centrifugal barrier.

Accretion is driven by the outward transport of angular momentum through the disk via the MRI. Each torus began with constant specific angular momentum, $`l`$, but Figure 12 shows that by orbit 10 the distribution $`l(r)`$ is approaching a more Keplerian slope. For reference, the Keplerian distributions of $`l`$ are shown as a dotted lines. In each model the inner portion of the torus becomes nearly Keplerian after only one orbit of time. The outer regions evolve more slowly. The specific angular momentum at late time is everywhere sub-Keplerian, indicating that some radial support is supplied by gas pressure (e.g. fig. 11). The inflow toward the horizon shows a steady decline in $`l`$ down to edge of the grid. This is indicative of continued magnetic stress acting on the fluid, a result first noted in the first pseudo-Newtonian simulations (e.g., Figure 6 of H00). In the fully relativistic models SFR and SF0 there is a noticeable roll-off in specific angular momentum for $`r<3M`$ (see Fig. 13). This roll-off is especially pronounced in the first orbit during the arrival of low angular momentum fluid that is released from the surface of the disk during a start-up transient. The feature remains visible and well resolved in the late stages as well. Model SFP does not show this abrupt drop in $`l`$, and this is likely due to the lack of an extended plunging region between the marginally stable orbit and the grid inner boundary.

A question of particular physical interest is, what is the total specific angular momentum carried by the accretion flow into the black hole? To answer this we divide the total angular momentum flux, $`T_{}^{r}{}_{\varphi }{}^{}`$, by the mass flux, $`\rho U^r`$, and examine the value at the inner radial boundary as a function of time. This is plotted in Figure 14. The values obtained by time-averages over the last three orbits are $`4.06`$ for SFR, $`3.32`$ for SF0, and $`1.94`$ for SFP. These are all slightly below the value of $`U_\varphi `$ corresponding to a circular orbit at $`r_{ms}`$ (shown as dashed lines in the figure), and considerably below the value in the initial torus.

A related question is the value of the specific energy carried into the hole. This is computed by the ratio of the energy flux to the mass flux, $`T_{}^{r}{}_{t}{}^{}/\rho U^r`$ at the inner radial boundary. In SFR this value is $`0.969`$, in SF0 it is $`0.957`$, and in SFP it is $`0.867`$. These are slightly less bound than the binding energy of the marginally stable circular orbit, $`(U_t)_{ms}`$, although, in all cases, more bound than the value of the initial torus ($`0.98`$). Although we do not reproduce the plot here, the specific binding energy of the inflow into the hole in SFR and SF0 does not vary much with time after about 4 orbits, while the value in SFP exhibits significant fluctuations. This is another consequence of the absence of an established plunging inflow inside of $`r_{ms}`$ for SFP.

The relationship between the magnetic field and the fluid’s specific angular momentum is illustrated in Figure 15, which is an azimuthally-averaged plot of the specific angular momentum in the accretion flow near the black hole for model SF0. Selected contours of azimuthally-averaged magnetic pressure are overlaid on the plot. The specific angular momentum shows considerable variation with angle $`\theta `$. Regions of high magnetic pressure correspond to regions of low specific angular momentum. The fluid with the highest specific angular momentum that reaches the black hole does so in a region of low magnetic pressure. The density-weighted, shell-averaged value $`l`$ cannot completely capture this complexity, but it is clear from examining points at constant $`r`$ near the inner edge of the flow, that it does fairly characterize the reduction in $`l`$ that occurs near the horizon.

Much of the behavior of the three simulations depends on the location of the marginally stable orbit relative to the accretion torus. Krolik & Hawley (2002) defined the “turbulence edge” of a disk as the point where the magnetic field switches from being controlled by MHD turbulence to simple flux-freezing. Concomitant with this is a transition from a region where the velocity field is dominated by fluctuations to one where the velocity corresponds to streaming inflow. In the pseudo-Newtonian simulations, Krolik & Hawley (2002) found this edge to be located at about $`1.2r_{ms}`$. Figure 16 is a density-weighted radial inflow velocity average, $`\rho U^r/\rho `$, as a function of $`r/r_{ms}`$ at $`t=10`$ orbits in the three fiducial models. The shell-averaged velocities within the turbulent portions of the disk are relatively small and variable. At a point just outside of $`r_{ms}`$, however, $`U^r`$ begins a smooth and steady increase with decreasing radius. From $`r_{ms}`$ inward the slopes of the three curves show close agreement.

### 3.2 Comparisons: Axisymmetric and full $`2\pi `$ Evolution

#### 3.2.1 Axisymmetry

Reducing the extent of the azimuthal direction is a very effective way of reducing the computational demands of a given problem. Axisymmetric (two-dimensional) simulations are particularly useful in this regard. But to what extent do axisymmetric simulations convey useful information? The limitations of axisymmetric pseudo-Newtonian simulations have been previously discussed by Hawley, et al. (2001). Here we carry out axisymmetric versions of our three fiducial models to investigate these issues in the general relativistic regime. These are designated SFR-2D, SF0-2D, and SFP-2D. The axisymmetric simulations will also provide a more complete picture of the role of azimuthal modes in the growth and maintenance of the MRI in the Kerr metric.

The axisymmetric simulations are run for a time equivalent to $`10`$ orbits at the pressure maximum. As in the fiducial simulations, the magnetic pressure increases rapidly due to the growth of toroidal field by shear. The toroidal field MRI, however, is inherently nonaxisymmetric, so only the poloidal MRI modes operate. As has previously been noted (e.g., Hawley et al. 2001), one of the main consequences of axisymmetry is that the evolution is dominated by a particularly violent form of the poloidal MRI, the so-called “channel solution” (Hawley & Balbus 1992). Figure 17 shows the density for model SF0-2D at $`t=3.0`$ orbits (models SFR-2D and SFP-2D are similar). The prominent radially-extended features are characteristic of the channel solution. The channel solution is itself subject to an instability which destroys its coherence, but that instability is nonaxisymmetric (Goodman & Xu 1994). Another known consequence of axisymmetry is that magnetic turbulence cannot be sustained indefinitely. In axisymmetry the toroidal field cannot be converted into poloidal field to create a dynamo; this is Cowling’s anti-dynamo theorem.

The combination of the channel solution plus the anti-dynamo theorem paradoxically produces both more violent initial accretion and weaker long-term accretion. In the axisymmetric models the overall accretion totals are reduced from their three-dimensional counterparts (Table 2). However, in axisymmetry the initial infall of matter is more pronounced; there are significant bursts of accretion from the disk to the horizon.

Figure 18 shows the accretion rate into the hole for model SFP-2D compared with that in SFP. The most obvious feature is that accretion in SFP-2D is highly time-variable, even episodic. In addition, the initial accretion flow takes longer to reach the hole; strong inflow doesn’t occur until 4.5 orbits, when there are sharp, high-amplitude spikes in $`\dot{M}`$. After this, accretion continues at a steady low rate with periodic brief increases. In this model, matter accumulates in a mini-torus near the black hole, much like it does in three dimensions. At late time the mini-torus continues to accrete into the hole, but inflow from the main disk to the mini-torus is greatly reduced. Figure 19 compares the averaged density $`\rho `$ at orbit 10 in model SFP and SFP-2D. One can see that the mini-torus is more or less distinct from the main torus, indicating that the accretion flow between these two regions has been greatly reduced.

In models SFR-2D and SF0-2D (not shown) there are strong bursts of accretion between orbits 3 and 5, and beyond that time the accretion rate declines steadily, punctuated by an occasional flare. By the end of the three axisymmetric models, the accretion rate has quieted considerably and much of the matter remains out in the region of the initial torus, relatively undisturbed.

#### 3.2.2 Full $`2\pi `$ plane

It has been noted in pseudo-Newtonian studies that quarter-$`\varphi `$-plane simulations are often adequate for capturing the essential aspects of accretion disk evolution (Hawley 2001; Nelson & Papaloizou 2003). Hawley (2001) found good qualitative agreement between the two types of grids, with an overall 10% reduction in field and stress levels in the restricted-$`\varphi `$ simulations.

To examine the differences between a full- and a quarter-plane simulation in the present context, we repeat SFP on a $`128^3`$ grid with the full range in the $`\varphi `$ coordinate ($`0\varphi 2\pi `$). This simulation is labeled SFP-2$`\pi `$. Since we simply increase the number of $`\varphi `$ zones by a factor of 4, the grid zone size, $`\mathrm{\Delta }\varphi `$, remains the same. Of the models run here, the prograde torus provided the greatest radial extent of turbulent disk evolution, and provides the most interesting comparison. Polar and equatorial slices through the dataset for the density $`\rho `$ are similar to Figure 4, so the gross overall dynamics of the disk are largely unaffected by the choice of $`\varphi `$ grid.

There are differences, however. A comparison of the density at orbit 10 in SFP-2$`\pi `$ and SFP (fig. 19) shows that the density in the mini-torus is substantially reduced in SFP-2$`\pi `$. Although the mini-torus is initially comparable in size in the two models, it seems to partially dissipate in the late stages of SFP-2$`\pi `$. In SFP the mini-torus appears to dissipate at $`t5`$ orbits, but then reconstitutes itself. This same behavior was seen in the pseudo-Newtonian simulation of Hawley & Balbus (2002); an inner torus formed, dissipated, and reformed during the length of the simulation. In SFP-2$`\pi `$ the density within the mini-torus declines after its formation until the end of the run. How can we account for this? First, the mini-torus forms, in part, because the accretion rate into the near-hole region doesn’t necessarily match the accretion flow that passes beyond the marginally stable orbit and into the hole itself. The initial nonlinear MRI, for example, can rapidly dump considerable material into the near-hole region. At later times inflow from the extended disk may be reduced, while an increase in turbulence within the mini-torus can increase the rate of its accretion into the hole, thereby yielding the observed descrease in density within the SFP-2$`\pi `$ mini-torus. It would seem, therefore, that greater stress levels in SFP-2$`\pi `$ have a noticeable impact on disk structure.

The long-term behavior of the accretion rate may constitute another significant difference between the models. Figure 18 compares the accretion rates in SFP-2$`\pi `$ and SFP. The fluxes are similar in magnitude, but the full-grid simulation is smoother, and has a slightly greater mass flux after $`t=7`$ orbits. Table 2 shows that the total integrated accretion in SFP-2$`\pi `$ is about 6% greater than that of SFP. This would indicate that the low-order azimuthal modes, whose existence is precluded in the quarter-grid simulations, do play a role in the accretion process and can increase the effective stress and angular momentum transport. The reduction in the level of the fluctuations, particularly when one also considers the axisymmetric result, also shows the role of the nonaxisymmetric modes in decreasing the coherence of the axisymmetric MRI modes (i.e., the channel solution). It is also possible that a reduction in $`\mathrm{\Delta }\varphi `$ would similarly reduce the fluctuations in $`\dot{M}`$, even with a quarter-grid domain. This hypothesis remains to be tested with higher resolution simulations.

### 3.3 Comparison with a Pseudo-Newtonian Simulation

To compare the results from the relativistic code with those that use a pseudo-Newtonian potential, we ran model PN0, an analogue of SF0. These two models have several differences in addition to Newtonian versus relativistic physics. The PN0 code is the same as used in earlier studies (e.g., Hawley & Krolik 2001, 2002) and, although the numerics in the two codes are based on the same principles of time-explicit Eulerian differencing, they are, of course, distinct algorithms. The PN0 code employs a cylindrical grid, $`(R,\varphi ,z)`$, with a unit of radius twice that of the Schwarzschild metric, i.e., $`R=2M`$. The radial grid has 128 zones and is equivalent to that of SF0, while the $`z`$ grid has 128 zones, with 64 of the zones located between $`z=\pm 5`$. In PN0 the initial torus has a pressure maximum at $`R=7.65`$ and an inner edge at $`R=4.75`$. The angular momentum is $`l=3.18`$ (unit value smaller by $`2^{1/2}`$ compared to SF0). The initial magnetic field consists of toroidal loops with $`\beta =100`$, set up in the torus according to eqn. (13). The orbital period at the pressure maximum is $`T_{orb}=115.6`$. The pseudo-Newtonian orbital period differs from that of SF0 by the radial unit difference $`2^{3/2}`$ ($`P_{orb}\mathrm{\Omega }^1r^{3/2}`$), plus a factor due to time-dilation, $`(12/r)`$. In the GR code time is measured in terms of the clock of the observer at infinity, whereas there is no distinction between local and global time in the Newtonian computation. When comparing the results from the two simulations we normalize time to orbits at the initial pressure maximum.

Figure 20 compares the history of the accretion into the hole for runs PN0 and SF0. The plots are very similar in overall shape. In both models the initial phase of magnetic pressure-driven accretion, the onset of the MRI saturation phase, and the transition to more or less steady turbulent accretion occur at the same time, normalized by the orbital period at the pressure maximum. Table 2 gives the mass budget for model PN0; the values are normalized to the initial torus mass. The mass loss through the outer boundary includes both the upper and lower $`z`$ boundaries as well as the outer cylindrical radial boundary. PN0 has lost a comparable amount of mass as the relativistic SF0 model.

Figure 21 shows the averaged specific angular momentum $`l`$ at 10 orbits in models SF0 and PN0. Note that in PN0 the average is taken over cylinders rather than spherical shells. However, since it is a density-weighted average, and most of the mass is near the equator, this difference should not be too great. One sees that there is substantial agreement between the two curves, including a continued decline inside of the marginally stable orbit. PN0 also shows a slight downturn in slope near the horizon, although not as pronounced as in SF0.

One can compute a mass averaged inflow velocity by dividing $`\dot{M}`$ by $`\rho `$. This averaged inflow velocity at 10 orbits in PN0 is larger than in SF0 by about 20% at $`r_{ms}`$, and increases more rapidly inward. The PN0 velocity reaches a maximum of about 2.4$`c`$ at the horizon while SF0 remains subluminal and approaches $`c`$. If accurate velocities are important to the analysis, a relativistic treatment is clearly preferred.

## 4 Discussion

In this paper, we present a preliminary survey of the properties of accreting tori in the Kerr metric. This work is the relativistic analogue of earlier work by Hawley (2000) which used a pseudo-Newtonian potential. We carried out three-dimensional simulations of an accretion torus around prograde and retrograde Kerr holes with large $`a`$ values, as well as around a nonrotating Schwarzschild hole. These simulations are contrasted with axisymmetric simulations and a pseudo-Newtonian simulation with equivalent initial conditions.

In comparing the three fiducial models, we find that the most important effect of black hole rotation is in setting the location of the marginally stable orbit. In retrograde model SFR, which has the outermost $`r_{ms}`$, the initial torus begins with its inner edge almost at the marginally stable orbit. As pressure builds and angular momentum is transported out through the torus, a slender accretion flow is set up which remains relatively smooth in the late stages of the simulation. The ease with which the matter is accreted is indicated by the observation that over the course of the simulation, $`42\%`$ of the initial disk mass enters the black hole in 10 orbits. At the other extreme, in the prograde model, SFP, the gas must slowly accrete from the initial torus down to $`r_{ms}=2.32M`$, a radius just outside the static limit. This accretion is driven by angular momentum transport from the MHD turbulence, and the flow develops into a dense, slightly sub-Keplerian disk. The inner edge of this disk features a local pressure maximum inside $`r5M`$ which we refer to as the “mini-torus.” All of these factors reduce the total accretion through the inner radial boundary into the hole over the course of the 10 orbit evolution. Only 18% of the initial torus mass is lost into the hole. A greater fraction of the mass is lifted to radii outside the outer boundary of the initial torus. The Schwarzschild model SF0 is an intermediate case which develops both a thick, turbulent accretion flow and a plunging inflow. SF0 loses $`31\%`$ of its initial mass to the black hole.

In all three models the accretion rate into the hole is highly time variable. SFR shows less variability on short time scales. The accretion rate into the hole is determined by the rate at which material from the turbulent disk is fed into the plunging flow from the disk’s “turbulence edge,” a point located near, but generally slightly outside, the radius of the marginally stable orbit. The accretion rate variability will be roughly determined, therefore, by the orbital frequency at the turbulence edge. Assuming that fluctuations in $`\dot{M}`$ have observable consequences (e.g., as a possible origin of high-frequency QPOs), then these simulations confirm the idea that the hole’s spin parameter $`a`$ helps determine the variability frequency. The intrinsic variability of the turbulence itself, however, prevents the relationship between $`a`$ and variability frequency from ever being more than approximate.

As in our earlier hydrodynamic simulations, we observe the effect of frame dragging most prominently in the retrograde model SFR, where the spiral accretion stream shows the characteristic “dog leg” pattern near the hole as the stream is dragged into co-rotation with the black hole inside $`r3M`$. The effects of prograde frame dragging are much more difficult to discern directly.

We find that, as in earlier pseudo-Newtonian simulations, an extended magnetized corona forms where the magnetic pressure is comparable to the gas pressure. Much of this corona is flowing radially outward. The largest outflows occur in SFP; again, because it has higher specific angular momentum relative to the marginally stable orbit, the gas is not as easily lost to the hole as in SFR.

In all the models the initial magnetic field corresponded to a volume-averaged strength of $`\beta =100`$. In the subsequent evolution the disk remains dominated by gas pressure, with $`\beta 10`$, while the coronal regions are more near equipartition, with $`\beta 1`$. Magnetic pressure dominates within the centrifugal funnel region near the black hole. However, there is no significant matter density within the funnel; the specific angular momentum of the gas prevents the accreted gas from approaching the axis.

An observation that clearly emerges from the simulations is the dominant role of the magnetic fields within the plunging inflow inside the marginally stable orbit. In models SF0 and SFR there is a noticeable transition from a gas-pressure dominated to a magnetic-pressure dominated flow inside $`r0.7r_{ms}`$. For model SFP, however, this transition does not occur; the marginally stable orbit lies too close to the inner radial boundary. Instead we see near-equality of the two pressures at the inner radial boundary. A higher-resolution simulation of model SFP which includes more zones within the ergosphere is required to resolve the plunging inflow.

The magnetic fields within the plunging region continue to exert a stress on the fluid. Regions of strong field correspond to the lowest specific $`l`$ values in the gas entering the hole. In all three models the specific angular momentum flux entering the hole, $`T_{}^{r}{}_{\varphi }{}^{}/\dot{M}`$, is slightly less than the value associated with the marginally stable orbit: 96% for SF0 and SFR, and 92 % for SFP.

We also compared the fiducial simulations to axisymmetric simulations. The overall initial evolution of the axisymmetric and three-dimensional simulations are similar. The magnetic pressure increases, causing the torus to expand and creating outflows, while the angular momentum within the torus is redistributed to a more Keplerian form. However, the duration of MRI-driven turbulence is short-lived in the axisymmetric case since the azimuthal modes that sustain the MRI are absent under axisymmetry. Further, as the SFP-2D simulation demonstrates, the turbulence can die out in one region of the flow while continuing in others, therefore significantly altering the global evolutionary properties.

As noted earlier by Hawley & Balbus (2002), axisymmetric simulations tend to be dominated by the poloidal field MRI channel solution which is characterized by radially extended filamentary structure and large angular momentum transport. The channel solution is also present in the three-dimensional simulations, but it is short-lived due to nonaxisymmetric modes, and its breakup heralds the onset of MRI-driven turbulence. These axisymmetric MRI modes are highly time-variable, causing accretion to occur in bursts. Comparing the results of the full $`2\pi `$ plane simulation SFP-2$`\pi `$ with the standard quarter-plane simulation shows that with greater axial extent, the accretion flow into the hole becomes smoother in time. This points to the role of nonaxisymmetric modes in moderating the influence of the axisymmetric channel solution.

Since axisymmetric simulations cannot model essential aspects of the accretion process, their utility is somewhat limited. Used with care they can provide an inexpensive way to evolve an idealized initial condition, such as a constant-$`l`$ torus with loops of field, into a more extended near-equilibrium structure that then serves as the starting point for a full three-dimensional simulation. They also provide a valuable tool for surveys of various initial conditions. Since axisymmetric simulations can be much more highly resolved, they also should continue to be useful as numerical experiments to investigate specific physical processes in detail.

In comparing full- and quarter-plane simulations we find that, while there are noticeable differences, the quarter-plane simulations capture much of the qualitative behavior. This is the consistent with the conclusion from previous pseudo-Newtonian studies. Their affordability makes them very useful for running a range of models, such as the present work. The accretion rate in the full-plane simulation is less strongly variable, presumably due to a reduction in the power of the axisymmetric poloidal field MRI. The time-averaged accretion rate, however, is higher, and this influences the size of the mini-torus that forms in the prograde runs. The optimal strategy may be to run models on restricted grids until the startup transients are past and the turbulent disk is well-established. The $`\varphi `$-grid can then be increased to get an improved representation of the approximately steady-state disk.

A direct comparison of relativistic model SF0 and pseudo-Newtonian model PN0 reveals good agreement: the two simulations show similar late-time angular momentum distributions as well as comparable accretion histories through the inner radial boundary. In the regions outside the marginally stable orbit relativistic and pseudo-Newtonian results overlap, which should not be surprising. This is a region of weaker gravitational field and the pseudo-Newtonian potential is a very good approximation to the Schwarzschild metric. The differences between the two simulations arise in the region inside $`r3M`$, where general relativistic effects become prominent. This is also the region where the velocities in pseudo-Newtonian models are likely to be the most inaccurate. In relativity velocities are limited to less than $`c`$, whereas there is no such limit in a pseudo-Newtonian simulation. We further note that there is no unique time-normalization that is appropriate for all locations in the flow when comparing a relativistic simulation with a Newtonian one. If details of the time-variability of the flow in the near hole region are important, relativistic calculations are required. Simulating accretion into Kerr holes also requires a relativistic treatment.

As the time-dependent models improve, it is important to connect the results to observable phenomena. Recently, Armitage & Reynolds (2003) took a first step toward this goal with a study of the effects of time variability on certain emission processes in the inner regions of a simulated black hole accretion disk. Machida & Matsumoto (2003) have evolved a global disk to search for regions within the plunging flow that might account for X-ray flares. These studies must ultimately be done with full relativity. In the present work we see both qualitative and quantitative differences among the models with different $`a`$ values in the plunging region inside the marginally stable orbit, where effects of black hole rotation are increasingly important. More detailed analyses of this plunging region will be the subject of subsequent papers.

This work was supported by NSF grant AST-0070979 and NASA grant NAG5-9266. The simulations were carried out on the Origin 2000 system at NCSA, and the Bluehorizon system of NPACI.


text:
# Effect of Antiferromagnetic Interlayer Coupling on Current-Assisted Magnetization Switching

## Abstract

We compare magnetization switching in Co/Cu/Co nanopillars with uncoupled and dipole-field coupled Co layers. In uncoupled nanopillars, current-driven switching is hysteretic at low magnetic field H and changes to reversible, characterized by telegraph noise, at high H. We show that dipolar coupling both affects the switching current and causes the switching to become reversible at small H. The coupling thus changes the switching to reversible, hysteretic, and then reversible again as H increases. We describe our results in terms of current-assisted thermal activation.

Observations of current-driven magnetization direction switching in Co/Cu/Co cornellorig ; cornellapl ; grollier ; cornelltemp ; cornellquant ; wegrowe ; sun2  and Py/Cu/Py \[Py = Permalloy = Ni84Fe16urazhdin  nanopillars have generated great interest, both for science—studies of magnetic systems driven far out of equilibrium, and technology—in magnetic random access memory (MRAM) this effect might eliminate the need for magnetic field-driven switching. In this Letter we show that dipolar coupling between the magnetic layers can affect the switching. Specifically, we show that the switching current is not determined solely by the switching mechanism cornellquant ; grollier , but varies with dipolar coupling, and that sufficiently strong antiferromagnetic (AF) coupling leads to reversible (non-hysteretic) switching at small magnetic field H. This effect may find application in high-sensitivity field sensors, e.g. read-heads of computer hard drives.

Our devices were nanofabricated using the following steps. First, a Cu(80)/Co(20)/Cu(6-10)/Co(2.5)/Cu(5)/Au(15) multilayer was sputtered onto an oxidized Si wafer in Ar pressure of 2 mTorr. All thicknesses in this Letter are in nanometers. An Al(50) nanopillar with lateral dimensions of about 70 nm by 130 nm was then formed by a combination of e-beam lithography and evaporation. The Al was used as a mask for ion-milling the multilayer into a nanopillar. Dipolar coupling between the two Co layers was minimized by timing the ion-milling to stop in the Cu layer above the bottom Co(20) layer. When desired, AF dipolar coupling was achieved by continuing the ion-milling about half way through the bottom Co(20) layer. Magnetic poles at the edges of the two patterned Co layers then AF-couple to minimize the magnetic energy. Without breaking the vacuum, a SiO(25-40) layer was deposited for electrical insulation between the device electrodes. The Al mask was removed by ion-milling with the ion-beam parallel to the sample surface, followed by wet etching. The ion-milling removed metals back-sputtered onto the Al mask. Finally, after a short ion-milling to clean the surface, a Au(150) top contact was sputtered onto the top Au layer. All measurements were performed at room temperature (295 K). Differential resistances dV/dI were measured using a standard four-probe setup with lock-in detection, adding an ac current of amplitude 20 $`\mu `$A at 8 kHz to the dc current I. Most uncoupled devices tested had room temperature resistances of about 1.5 $`\mathrm{\Omega }`$ and magnetoresistances (MR) of about 5%, similar to the best MR values reported by others cornellorig . The coupled devices had larger resistances (because of the additional interface and Co layer in the nanopillar) but similar MRs. We define positive current to flow from the bottom (extended) to the top Co layer.

Typical data for an uncoupled sample are shown in Fig. 1, and comparative data for a coupled sample are shown in Fig. 2. We begin with the uncoupled sample, the data for which are consistent with prior studies of uncoupled Co/Cu/Co cornellapl  and Py/Cu/Py urazhdin . Fig. 1(a) shows that the MR at I = 0 is symmetrically hysteretic in H and the increase in resistance from the low resistance parallel (P) orientation of magnetic moments to the high resistance anti-parallel (AP) orientation occurs in a single sharp step. The current-driven switching is hysteretic at H = 0 Oe (solid curve in Fig. 1(b)), with switching to the high resistance anti-parallel (AP) state at positive $`I_s^{PAP}`$, and to the low resistance parallel (P) state at negative $`I_s^{APP}`$. The resistances and their changes in Fig. 1(b) are close to those in Fig. 1(a). The dotted curve in Fig. 1(b) shows that at large enough H the hysteretic step in dV/dI turns into a non-hysteretic (reversible) switching peak.

Fig. 1(c) shows the switching diagram obtained from MR data such as those in Fig. 1(a). We attribute the slight H-asymmetry of the diagram to a combination of the current-induced Oersted field and sample shape asymmetry. The P$``$AP switching field $`H_s^{PAP}70`$ Oe is independent of I over the range $`1.7`$ mA$`<I<1.7`$ mA. We attribute the P$``$AP switching in this range to reversal of the extended Co layer, unaffected by the small current density in that layer. We associate the AP$``$P transition with reversal of the thin patterned Co layer at $`H_s^{APP}`$ determined by its shape anisotropy. At $`I>1.7`$ mA, $`H_s^{PAP}`$ strongly varies with I. In this regime, we attribute the P$``$AP transition to reversal of the patterned Co layer, induced by $`I>0`$. As the extended Co layer reverses at $`H70`$ Oe, the patterned Co layer reverses simultaneously to keep the AP configuration favored by $`I>0`$. Such simultaneous reversal sometimes produces a weak feature in dV/dI at small H.

At $`I>4`$ mA, the data in Fig. 1(c) become non-hysteretic. Time resolved measurements show that the non-hysteretic switching peak in dV/dI is characterized by telegraph noise switching between the P and AP states, with random distribution of dwell times in the P state ($`\tau _P`$) and AP state ($`\tau _{AP}`$urazhdin . When I is increased, and H adjusted to keep $`\tau _P=\tau _{AP}`$, the average noise period decreases approximately exponentially (Fig. 1(d)). The switching diagram, Fig. 1(c), is asymmetric with respect to the current direction; switching is hysteretic at $`I<0`$, and non-hysteretic only at large $`I>0`$. This difference occurs because application of H favors the P state, while $`I>0`$ favors the AP state.

Fig. 2 shows data similar to those in Fig. 1, but for a sample with strong AF dipolar coupling between the two Co layers, produced by partial patterning of the extended Co(20) layer. Again, the field-driven MR (Fig. 2(a)) is hysteretic and approximately symmetrically about $`H=0`$. But now the MR contains multiple steps, likely because the highly inhomogeneous dipole field favors nonuniform magnetization states of the nanopillar. We determine the switching points from the jump into or from the lowest resistance state. As H is reduced from a large positive value (solid line in Fig. 2(a)), coupling between the two Co layers causes the thinner Co (2.5) layer to switch to the AP state at $`H=0.16`$ kOe, then at small $`H<0`$ both the thicker and thinner Co layers flip together to stay in the AP state, and finally at $`H0.8`$ kOe the Co(2.5) layer reverses again to return to the P state.

At small H =50 Oe the behavior as a function of I in Fig. 2(b) is quite different from that in Fig. 1(b); instead of hysteretic switching, the data show a non-hysteretic, downward peak. This behavior is the most important feature of the data for coupled samples. In Fig. 2(b), the switching is non-hysteretic for small H, becomes hysteretic for intermediate H=0.3 kOe, and then non-hysteretic again for large enough H. Fig. 2(c) shows the switching diagram. In contrast to Fig. 1(c), all switching points now represent reversal of the the Co(2.5) layer and switching becomes reversible at large enough magnitude of current in both directions, at small H for negative I, and large H for positive I. The symmetry between effects of positive and negative current is extended further by time-resolved measurements of telegraph noise, which is present close to the reversible switching points for both current directions. Fig. 2(d) shows that as the magnitude of I is increased, while adjusting H to keep $`\tau _P=\tau _{AP}`$, the telegraph noise period decreases approximately exponentially at similar rates for both current directions. Not all samples with partially patterned bottom Co layer exhibited non-hysteretic switching at small H. Due to weaker AF-coupling, some samples exhibited only a dip in $`I_s^{PAP}`$ at small H, resulting in switching diagrams intermediate between Fig. 1(c) and Fig. 2(c).

We interpret the data in Figs. 1,2 in terms of thermal activation over an effective switching barrier wegrowe ; cornelltemp ; urazhdin . In Co/Cu/Co nanopillars, both the semi-classical spin-torque slonczewski  and quantum threshold berger ; tsoi  models predict magnetic excitation in the P state at $`I>0`$ and in the AP state at $`I<0`$. Such current-driven excitation can be described in terms of a current-dependent effective switching barrier zhang , or a current-dependent effective magnetic temperature $`T_m`$, with an effective barrier modified by the current only through the variation of magnetization with temperature wegrowe ; urazhdin . Here $`T_m`$ depends on the magnetic configuration, $`T_m^{AP}(I)T_m^P(I)`$. These alternative approaches give mathematically similar results, but differ in details that still need experimental testing. We use the latter model, which we find more transparent.

In magnetically uncoupled samples, the magnetization orientation of the patterned layer is bistable at $`H=0`$, $`I=0`$: at T=295 K, the barrier height for switching between the two magnetization orientations significantly exceeds the thermal energy. In AF coupled samples at small H, the P$``$AP switching barrier is reduced by the dipolar field, leading to a thermally activated P$``$AP transition. The reverse AP$``$P transition cannot be thermally activated because the corresponding switching barrier is significantly higher. Thus, at I=0, H=0, AP is the only stable orientation of the nanopillar, as seen in Fig. 2(b,c). Current-driven magnetic excitation at $`I<0`$ increases $`T_m^{AP}`$. At sufficiently high $`I<0`$ the AP$``$P transition becomes thermally activated, leading to telegraph noise switching between the AP and P states (Fig. 3(a)). At larger H=0.3 kOe in Fig. 2(b), H nearly compensates the dipole field, leading to hysteretic switching similar to that in Fig. 1(a) at H=0 in an uncoupled sample, and illustrated in Fig. 3(b). As H is further increased, both the uncoupled and coupled samples behave similarly; As shown in Fig. 3(c), the AP$``$P transition becomes thermally activated. At large enough $`I>0`$ the P$``$AP transition also becomes activated due to current-induced increase of $`T_m^P`$, leading to telegraph noise both in uncoupled (Fig. 1(d)) and AF-coupled (Fig. 2(d)) nanopillars.

In summary, we have shown that dipolar AF coupling between magnetic layers leads to reversible current-driven magnetization switching at small H. Similarly, reversible switching at small H should be induced by the opposite current direction in ferromagnetically exchanged-coupled samples. We will demonstrate such behavior elsewhere urazhdin2 . Thus, the switching current and hysteresis are not intrinsic characteristics of the current-driven switching mechanisms, as they are strongly affected by the coupling between the magnetic layers.

Fig. 2(c) shows that, in the reversible switching regime at small H, the switching field of AF-coupled magnetically nanopillars can be adjusted by changing the applied current. In this regime, magnetically coupled nanopillars with exchange-biased extended magnetic layer and zeroed-out switching field may find application as high sensitivity field sensors.

We acknowledge helpful discussions with Norman O. Birge, important contributions to the development of the sample preparation technique by K. Eid and J. Caballero, and support from the MSU CFMR, CSM, the MSU Keck Microfabrication facility, the NSF through Grants DMR 02-02476 and 98-09688, and Seagate Technology.