text:
# Global MHD Simulation of the Inner Accretion Disk in a Pseudo-Newtonian Potential

## 1 Introduction

For many years, the physical nature of angular momentum transport in accretion disks has been one of the great outstanding questions of high-energy astrophysics. As a result of nearly a decade of effort, it is now becoming increasingly clear that the mechanism is turbulent magnetic stress in which the energy for the magnetic field is largely drawn from the free energy of orbital shear (as reviewed by Balbus & Hawley 1998). Theoretical grounds for this conclusion were first found in the rediscovery of the magneto-rotational instability (MRI) by Balbus & Hawley (1991). Increasingly detailed simulations of local “shearing boxes” have bolstered our confidence that this is indeed the process found in nature (e.g., Hawley, Gammie & Balbus 1995, 1996; Brandenburg et al. 1995; Stone et al. 1996; Miller & Stone 2000).

What has been missing until recently is studies of the radial structure of disks accreting under the influence of magnetohydrodynamic (MHD) turbulent stresses. Hitherto, only a few such three-dimensional global simulations have been presented (Armitage 1998; Matsumoto 1999; Hawley 2000; Machida, Hayashi, & Matsumoto 2000). In global simulations it is difficult to obtain adequate resolution over the wide range of length scales involved, namely the MHD turbulent length scales within the disk, the disk scale height $`H`$, and the radial distance from the central object $`R`$. These length scales are generally quite disparate, but near the inner boundary of an accretion disk they should become comparable. In this paper we present a global simulation of an accretion disk in which we focus on this inner region where adequate numerical resolution can be obtained.

With the results of this new simulation, we can begin to answer a number of questions of central interest to the mechanics of accretion. For example, the most common view of the inner portions of black hole accretion disks is based on the work of Novikov & Thorne (1973) (e.g. Abramowicz & Zurek 1981; Muchotrzeb & Paczyński 1982; Matsumoto et al. 1984). In this picture, the disk is (nearly) time-steady and axisymmetric. Page & Thorne (1974) argued on heuristic grounds that the $`R`$-$`\varphi `$ component of the stress (the one responsible for angular momentum transport) should go to zero at the marginally stable orbit; later work by Abramowicz and Kato (1989) showed that if the stress is a (constant) small fraction of the local pressure, it automatically approaches zero near that location.

However, because the MHD turbulence is the result of a general, local, rapidly-growing instability, assumptions such as stationarity and axisymmetry are likely inappropriate. One might also question whether there is any reason for the stress to diminish at the disk’s inner edge, given that the growth rate of the MRI does not appreciably diminish near the marginally stable orbit (cf. Krolik 1999; Gammie 1999; see also the footnote in Page & Thorne 1974 in which they speculate that their heuristic arguments might fail if strong magnetic fields were present). Further, because magnetic stresses do not automatically maintain a fixed ratio to the pressure, conclusions based upon $`\alpha `$ parameterizations (e.g., Abramowicz & Kato 1989) may not directly apply. This issue is important because continued stress at the marginally stable orbit could alter both the energy and angular momentum with which matter arrives at the black hole; that is, this boundary condition determines the efficiency of accretion and the rate at which the black hole is spun up.

Because the greatest luminosity should be released near the inner edge of the accretion disk, the detailed behavior and structure of the region near the marginally stable orbit is crucial to understanding the observations of black hole systems. For example, the luminosity of essentially all accreting black holes exhibits sizable fluctuations with a very broad-band distribution of fluctuation power with timescale (e.g., as discussed by Sunyaev & Revnivtsev 2000); can we identify the specific dynamics that drive these variations? Occasionally, the lightcurve exhibits quasi-periodic oscillations (QPOs); can we either identify the mechanism responsible for these, or test some of those mechanisms that have been proposed?

In the following section we will present a technical description of the three-dimensional global simulation we report. We will then set out our results in §3, with their qualitative implications outlined in §4. In §5 we will discuss the degree to which limitations of our simulation hamper direct application of our results, and suggest in which directions the greatest improvement might soon be made. We summarize our conclusions in §6.

## 2 Numerical Method

### 2.1 Equations

An investigation of the inner edge of a disk orbiting a black hole would be best carried out with a fully general relativistic simulation. However, we have not yet written a three dimensional general relativistic MHD code. In the present study we instead make use of the existing three dimensional Newtonian MHD code (Hawley 2000). The code evolves the equations of ideal MHD, i.e.,

$$\frac{\rho }{t}+(\rho 𝐯)=0$$
(1)
$$\rho \frac{𝐯}{t}+(\rho 𝐯)𝐯=\left(P+𝒬+\frac{B^2}{8\pi }\right)\rho \mathrm{\Phi }+\left(\frac{𝐁}{4\pi }\right)𝐁$$
(2)
$$\frac{\rho ϵ}{t}+(\rho ϵ𝐯)=(P+𝒬)𝐯$$
(3)
$$\frac{𝐁}{t}=\times \left(𝐯\times 𝐁\right)$$
(4)
where $`\rho `$ is the mass density, $`ϵ`$ is the specific internal energy, $`𝐯`$ is the fluid velocity, $`P`$ is the pressure, $`\mathrm{\Phi }`$ is the gravitational potential, $`𝐁`$ is the magnetic field vector, and $`𝒬`$ is an explicit artificial viscosity of the form described by Stone & Norman (1992a). The global disk code is written in cylindrical coordinates, $`(R,\varphi ,z)`$. To model important effects associated with the relativistic Schwarzschild metric, we employ the pseudo-Newtonian potential of Paczyński & Wiita (1980). The pseudo-Newtonian potential has the form

$$\mathrm{\Phi }=\frac{GM}{rr_g}$$
(5)
where $`r`$ is spherical radius, and $`r_g2GM/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

$$\mathrm{}_{kep}=(GMr)^{1/2}\frac{r}{rr_g},$$
(6)
and the angular frequency $`\mathrm{\Omega }=\mathrm{}/R^2`$. The innermost marginally stable circular orbit is located at $`r_{ms}=3r_g`$.

We use an adiabatic equation of state, $`P=\rho ϵ(\mathrm{\Gamma }1)=K\rho ^\mathrm{\Gamma }`$, with $`\mathrm{\Gamma }=5/3`$, and ignore radiation transport and losses. Since there is no explicit resistivity or physical viscosity, the gas can heat only through adiabatic compression or by the action of the artificial viscosity which acts in shocks. Equations (1)-(4) are solved using time-explicit Eulerian finite differencing. The numerical algorithm is that employed by the ZEUS code for hydrodynamics (Stone & Norman 1992a) and MHD (Stone & Norman 1992b; Hawley & Stone 1995). We adopt units where $`GM=1`$ and $`r_g=1`$ (so that $`c=\sqrt{2}`$).

In this paper we repeat—with significantly better resolution—simulation GT4 from Hawley (2000), whose initial state was a thick torus with an angular velocity distribution $`\mathrm{\Omega }R^{1.68}`$, slightly steeper than Keplerian. The angular momentum within the torus is equal to the Keplerian value at $`R=10`$, which determines the location of the pressure maximum in the torus. In somewhat arbitrary units, the pressure at this point is equal to 0.036 and the density $`\rho _{max}=34`$. The initial magnetic field is constructed by setting the toroidal component of the vector potential $`A_\varphi (R,z)=\rho (R,z)\rho _{min}`$, for all $`\rho `$ greater than a minimum value, $`\rho _{min}=0.1`$. Poloidal field components are then constructed from $`𝐁=\times 𝐀`$. This procedure produces large-scale poloidal field lines that follow the equal-density contours inside the torus. The magnetic field is initially zero everywhere $`\rho \rho _{min}`$. As discussed in Hawley (2000), this configuration is motivated by the desire to have well-resolved unstable modes of the poloidal field instability, and to have the initial field fully contained within the torus. The magnetic field level is renormalized so that the total magnetic energy in the torus is 0.01 the total integrated gas pressure. This corresponds to an overall value of $`\beta =P_{\mathrm{gas}}/P_{\mathrm{mag}}=100`$.

The aim of this initial condition is not so much to follow the evolution of this particular torus, as it is to construct a disk that is nearly Keplerian with an inner boundary that is initially close to, but still outside, the marginally stable orbit. To that end, the torus is perturbed with low-amplitude random adiabatic pressure fluctuations. These seed the magneto-rotational instability whose rapid evolution leads to the desired state. Further disk evolution results in inflow through $`r_{ms}`$; it is this resulting accretion flow that will be studied in detail.

The computational grid has $`128\times 128\times 128`$ zones. This is the same number of zones as used in GT4 of Hawley (2000), but here, to increase the resolution within the disk itself and in the inflow region, we locate more of the zones near the marginally stable orbit and around the equator. The center of the coordinate system is excised, i.e., the radial coordinate begins at a nonzero $`R_{\mathrm{min}}=1.5`$; this choice avoids both the coordinate singularity associated with the axis and the singularity at $`r_g`$. The outer radius is set at $`R=31.5`$. The $`z`$ coordinate is centered on the equatorial plane, and runs from -10 to $`+10`$. The angle $`\varphi `$ is periodic and covers the full $`2\pi `$. In the radial direction, there are 30 equally spaced zones in $`R`$ from 1.5 to 4; the radial zones are then graded logarithmically from 4 to the outer boundary at 31.5. In the vertical direction, half of the $`z`$ zones are equally spaced around the equator from $`2`$ to 2, then graded from those locations out to the vertical boundaries at $`z=\pm 10`$. The grid cells are equally spaced in $`\varphi `$. The simulation GT4 from Hawley (2000) used equally spaced zones in $`R`$ and $`z`$, and the outer radial boundary was set at 21.5. We have also run this initial torus at lower resolution for comparison, using $`64^3`$ grid zones.

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. Although this produces a nonzero stress at the boundary, it works well in preventing an artificial buildup of field at the boundary. Most of the time, for $`|z|1`$ along the inner radial boundary and $`R2`$ in the equatorial plane, the poloidal velocity is greater than the magnetosonic speed, so that the influence of the boundary condition on the flow upstream is limited.

Timescales of interest in the simulation are set by the circular orbital period $`P_{orb}=2\pi \mathrm{\Omega }^1=2\pi r^{3/2}(rr_g)/r`$ at significant points within the disk. In units where $`GM=r_g=1`$, $`P_{orb}`$ at the pressure maximum ($`R=10`$) is $`179`$, and at the marginally stable orbit, where the orbital frequency $`\mathrm{\Omega }=1/(2\sqrt{3})`$, the orbital period is 21.8. The local accretion timescale is set by the effective stress, namely $`t_{acc}=R/v_R\mathrm{}/W_{R\varphi }(P_{orb})`$, where $`W_{R\varphi }`$ is the $`R`$$`\varphi `$ component of the stress, $`W_{R\varphi }=\delta v_R\delta v_\varphi v_{AR}v_{A\varphi }`$, and $`v_A`$ is an Alfvén velocity. The goal is to evolve long enough to observe a number of accretion times near the marginally stable orbit. Since the code is time-explicit and therefore Courant-limited, this is difficult to achieve over large radial extents. Here we ran the simulation out to time $`t=1500`$, which is 69 orbits at $`r_{ms}`$.

### 2.2 Diagnostics

The quantity of data associated with a three dimensional simulation is daunting, and this makes analyzing the outcome a challenge. Diagnostics can be computed during the simulation and in post-processing. The best procedures for this are still being developed. Here we briefly describe some of the diagnostic procedures used in the present simulation.

The time-evolution of the disk can be visualized through animation sequences. Animation frames of the $`\varphi =0`$ slice and the equatorial ($`z=0`$) slice are dumped every 4 units in time. Full 3D density information is saved every 10 units of time, and these are used to make an animation of a three dimensional volumetric rendering.

The time-evolution can also be studied in a manageable way through space-time $`(R,t)`$ studies of azimuthally- and vertically-integrated quantities. Examples include the averaged mass density

$$\rho =\frac{\rho R𝑑\varphi 𝑑z}{R𝑑\varphi 𝑑z},$$
(7)
the average angular momentum,

$$\rho \mathrm{}=\frac{\rho \mathrm{}R𝑑\varphi 𝑑z}{R𝑑\varphi 𝑑z},$$
(8)
and average specific angular momentum,

$$\mathrm{}=\frac{\rho \mathrm{}}{\rho }.$$
(9)
The net radial mass flux is

$$\dot{M}=\rho v_R=\rho v_RRd\varphi dz,$$
(10)
here defined to be positive for net accretion (inflow). Similarly, the angular momentum flux,

$$\rho v_R\mathrm{}=\rho v_R\mathrm{}R𝑑\varphi 𝑑z,$$
(11)
can be computed.

Radial angular momentum transport in the disk is due to the $`R`$$`\varphi `$ component of the stress,

$$T^{R\varphi }=\rho \delta v_R\delta v_\varphi \frac{B_RB_\varphi }{4\pi },$$
(12)
where $`\delta v_R`$ and $`\delta v_\varphi `$ are the turbulent portions of the velocity, that is, the departures from the mean flow. The kinematic portion is the Reynolds stress and the magnetic portion is the Maxwell stress. In the simulation we have no way of knowing what the mean flow is. For example, the typical radial velocity is much larger than the time-averaged net radial drift velocity, and there are significant differences between the instantaneous angular velocity and the value associated with a circular orbit (which one might assume should be the mean orbital flow). As a result, there is no unique prescription to compute the Reynolds component of the total stress. Hawley (2000) used a definition for the perturbed orbital velocity $`\delta v_\varphi `$ in terms of the difference between the total instantaneous angular momentum flux, and the mass flux times the average specific angular momentum. The difference then represents the excess or deficit angular momentum transport due to orbital velocity fluctuations compared to the mean. Specifically,

$$\rho \delta v_R\delta v_\varphi =\rho v_Rv_\varphi \rho v_R\mathrm{}/R.$$
(13)
In contrast, we can easily compute the stress due solely to magnetic forces, the second term on the right-hand-side of equation (12); we denote this quantity by $`M^{R\varphi }`$. The vertically- and azimuthally-averaged Maxwell stress is

$$M^{R\varphi }=\frac{(B_RB_\varphi /4\pi )R𝑑\varphi 𝑑z}{R𝑑\varphi 𝑑z}.$$
(14)
We do not distinguish the total magnetic field from the fluctuating field because, as we shall show later, where the field is dynamically important, the magnitude of the fluctuating portion is large compared to any time-averaged mean field.

In addition to these vertically- and azimuthally-averaged quantities, we examine full data dumps of the simulation at specific moments in time. From these we can examine not only the full three dimensional structure, but also azimuthally-averaged and vertically-averaged slices in $`(R,z)`$ and $`(R,\varphi )`$.

## 3 Results

### 3.1 Quasi-stationary state

The initial torus was in hydrodynamical equilibrium, but it evolves rapidly as the simulation proceeds. The overall evolution is very similar to that of GT4 in Hawley (2000; see Fig. 13). The radial field is sheared by the differential rotation, creating toroidal field. Turbulence develops within the disk with the onset of the magnetorotational instability. The resulting Maxwell stresses drive angular momentum transfer; the disk begins to accrete into the central hole and expands radially outward. Relatively rapidly the inner portion of the disk achieves an approximate steady-state. By $`t700`$, the mass accretion rate reaches roughly the long-time mean; after $`t1000`$, the shape of the disk’s average radial density distribution $`\rho (R)`$ no longer changes appreciably. Thus, we observe the disk in a quasi-stationary state for several dozen orbital periods at the marginally stable orbit. In much of what follows, we will present a detailed analysis of the state of the disk during this time frame, focusing particularly on the end of the simulation, using that “snapshot” as a fairly typical sample of conditions in the quasi-stationary disk. However, it must be borne in mind that there are always significant fluctuations around all these mean states.

### 3.2 Gas density

A contour plot of the azimuthally-averaged gas density at the initial and final times is presented in Figure 1. As might be expected, the density diminishes inside the marginally stable orbit and outside the initial outer radial boundary of the disk, while also decreasing with increasing altitude away from the equatorial plane. Although the vertical heights of the highest contour levels at the final time are roughly independent of radius within the disk proper, the azimuthally-averaged exponential scale-height of the pressure, $`H`$, increases very nearly linearly with radius from $`R2`$ out to $`R20`$; at the marginally stable orbit ($`R=3`$), $`H0.4`$, while at $`R=10`$, $`H1.5`$. This particular disk shape is almost certainly a consequence of the equation of state, in which all fluid elements retain the same specific entropy except for shock heating. As a result of this special choice of equation of state, the temperature rises inward as the gas is compressed; by coincidence, this occurs at such a rate as to keep the ratio of sound speed to orbital speed, and therefore $`H/R`$, almost constant at $`0.15`$. In a real disk, one in which turbulent dissipation introduces heat and radiation vents it, the shape might be rather different.

Density gradients in radius and altitude are determined primarily by the overall dynamical effects of approximate hydrostatic equilibrium. No such requirements exist for azimuthal fluctuations, whose amplitudes characterize the turbulence. We quantify the azimuthal fluctuation level by defining

$$\frac{\delta \rho }{\rho }(R,z)=\frac{1}{\rho _\varphi }\left\{\frac{1}{2\pi }𝑑\varphi \left[\rho \rho _\varphi \right]^2\right\}^{1/2},$$
(15)
where $`X_\varphi `$ refers to the azimuthal average at fixed $`R`$ and $`z`$ of the quantity $`X`$. Measured in this way, $`\delta \rho /\rho 0.2`$ – 0.6 over most of the problem volume, with smaller values typically nearer the equatorial plane. In a few places, $`\delta \rho /\rho `$ rises to be greater than unity. Thus, quite sizable density fluctuations are generic.

These fluctuations are primarily sheared strips whose radial extent is $`<0.1R`$, but extend $`1`$ radian in azimuthal angle (Figure 2). This pattern (tightly wrapped spiral features) is characteristic of fluctuations in almost all quantities. Part of their nature is revealed by constructing spacetime diagrams in which the surface density at a specific radius is plotted as a function of $`t`$ and $`\varphi `$. From this visualization we are able to measure the angular speed of these features; it is generally quite close to the local orbital frequency $`\mathrm{\Omega }`$. We can also measure their persistence at any particular radius; it is usually less than an orbital period. In this sense, they can be regarded as fluctuations that merely “ride” with the local orbital motion, but see further discussion in §3.5.

### 3.3 Magnetic field

The magnetic field distribution is rather more complicated than the density distribution. A poloidal projection (i.e., averaged azimuthally) shows the magnetic pressure to have its largest value immediately outside the dense regions of the torus, with large variations inside the torus (Figure 3). In the shearing-box simulations of Miller & Stone (2000), the field strength likewise peaked away from the midplane, but we find somewhat greater contrast. Field-strength fluctuations are also larger than those for density: the rms fractional fluctuations in $`|\stackrel{}{B}|`$ in the azimuthal direction are $`1`$ – 1.5 in the disk body, falling to $`0.25`$ several gas scale-heights immediately above and below the disk. Another measure of the magnitude of the fluctuations is $`|B_i_\varphi |/|B_i|_\varphi `$, the ratio of the mean value of an individual field component to the mean of its absolute value. This quantity is typically $`0.2`$ for all components wherever the gas density is significant, indicating that the field is extremely variable as a function of azimuth (and therefore as a function of time at fixed azimuth, as orbital rotation carries the field around). Only where the gas density is very low does the field become relatively smooth and steady (a result also qualitatively similar to the findings of Miller & Stone 2000). Like the density, the characteristic shape of fluctuations is sheared arcs that run for $`1`$ radian before losing coherence. However, the peaks in magnetic field strength tend to lie adjacent to, rather than coincident with, peaks in gas density.

To illustrate some of the field structure, we show in Figure 4 poloidal slices taken at $`\varphi =\pi /2`$. Panel a) is a close-up of the inner, high-resolution region; panel b) shows (almost) the entire disk. This view, of course, suppresses the toroidal component (which is typically larger than the poloidal component by a factor $`1`$ – 10); in addition, in order to make the figures visually comprehensible, we have omitted small-scale field variations (smaller than 0.125 length units in panel a), smaller than 0.5 length units in panel b). Even with this smoothing, it is apparent that the field is very tangled. The only field property exhibiting any stability is the roughly dipolar pattern seen in the low density regions. This is a remnant of the initial condition. Because it only persists in these very low density regions, and they do not contribute much to the stress (see §3.6), we do not regard this remnant artifact as dynamically significant.

As one expects from magnetic fields in shearing plasma, $`B_R`$ and $`B_\varphi `$ are correlated in the sense that the mean ratio of the Maxwell stress to the magnetic pressure, $`2B_RB_\varphi /B^2>0`$ and has (approximately) constant magnitude $`0.2`$ – 0.3 in the body of the disk. Inside the marginally stable orbit, this ratio increases slowly in magnitude as radial accretion “combs” the field out into radial loops that are strongly sheared, increasing the correlation between $`B_R`$ and $`B_\varphi `$. In fact, the Maxwell stress is so correlated with the average magnetic pressure $`B^2/8\pi `$, that a figure showing the stress is almost indistinguishable from one showing the magnetic pressure.

### 3.4 Accretion rate

Perhaps the single most basic quantity of interest is the accretion rate, $`\dot{M}`$, defined in equation (10). The disk is never perfectly time-steady; $`\rho v_R`$ varies substantially from place to place and from time to time. In fact, disk quantities tend to be quite nonaxisymmetric, and much of the final accretion inside $`r_{ms}`$ takes place in a spiral flow. However, time-averaging can remove much of this variability, and Figure 5 shows the time-averaged accretion rate as a function of $`R`$ along with the instantaneous values of $`\dot{M}(R)`$ at $`t=1000`$ and 1500. The time-averaged $`\dot{M}(R)`$ is nearly constant as a function of radius inside $`R=10`$, with $`\dot{M}=5`$. We shall adopt this as the fiducial steady state accretion rate for the inner disk. Between $`R=9`$ and 14 there is inflow at a decreasing rate; outside $`R=14`$ the net flux is outward. In this simulation, the initial disk is entirely contained within the grid, and the outer part of the disk must move outward in order to “soak up” the angular momentum received from the inner disk material. Over the course of the simulation the total disk mass decreases by 14%; about 90% of this has gone into the black hole.

Figure 6 shows a detailed time-history of the accretion rate through the inner radial boundary. As can be seen from this figure, the accretion rate shows significant time variations after its initial growth phase. The nature of the time-variability can be studied using the Fourier power spectrum of the accretion rate during the latter portion of the simulation (beyond $`t=600`$, i.e., after the initial growth phase is over). This is shown in Figure 7. It is most simply described as a “red” power-law; roughly speaking, $`P(f)f^{1.7}`$. However, there is also some curvature, in the sense that it steepens with increasing frequency. Because the accretion rate does appear to have a well-defined mean value, $`d\mathrm{ln}P/d\mathrm{ln}f`$ must be $`>1`$ at the lowest frequencies. At the high frequency end, the spectrum is largely featureless, extending almost to the inverse free-fall time at the innermost boundary. Other than a smooth steepening, there is no special feature at the orbital frequency of the marginally stable orbit (marked on the figure). The flow is unsteady everywhere, even well inside the plunging region.

Figure 6 provides some idea of how nonsteady the accretion is; the mass flux is highly spatially inhomogeneous as well. To gain some sense of the magnitude of the spatial variations, we show the accretion rate as a function of position on a cylinder at $`R=5`$ (Figure 8). As this figure shows, the fluctuations are very large. If we suppose that the mean accretion rate of 5 is evenly distributed over the range of altitudes $`1z1`$, the mean mass flux would be $`0.1`$. As the figure makes plain, there are several places where the mass flux is five times as great as this; there are also several places where the mass flux is $`0.3`$ – 0.4 outward. Relatively large fluctuations compared to the mean are to be expected when the transport of angular momentum in the disk is by turbulence.

### 3.5 Propagation of fluctuations

As we have stressed all along, one of the most striking features of this simulation is the presence of large amplitude nonaxisymmetric modulations in the density, mass flux, and magnetic field strength. We have already discussed how $`(\varphi ,t)`$ spacetime plots enable us to measure their angular velocity. It is possible to learn more by constructing analogous plots in $`(R,t)`$. For example, in Figure 9 we show the history of the azimuthally- and vertically-averaged mass flux (eq. 10) in this simulation. The diagonal stripes show the radial motion through the disk of fluctuations in the mass flux, their slopes indicating their radial group speeds. Typically these disturbances begin near $`R5`$ and travel both outward and inward. In the outward direction, they move at a very nearly constant speed $`0.07`$, while in the inward direction, their speed is at least three times faster. Because the outward speed very nearly matches the magnetosonic speed, we infer that these disturbances propagate as magnetosonic waves in the outward direction. On the other hand, inward-going mass flux perturbations are simply advected along with the dynamical inflow.

We conjecture that their origin lies in the way in which material enters the region of unstable orbits. When a fluid element begins to accelerate inward, it creates a rarefaction wave that moves outward at the magnetosonic speed, while traveling azimuthally at the (faster) orbital speed. On the other hand, the radial inflow speed near and inside the marginally stable orbit is comparable to the magnetosonic speed, so fluctuations travel inward rather more quickly. These fluctuations are never smoothed out into a completely steady flow because individual fluid elements move into the plunging region as a result of the continually changing torques exerted by the MHD turbulence (see the next subsection). Thus, at this level of detail, there is no way to achieve a stationary state.

Although strong nonaxisymmetric pressure waves were generated by MHD turbulence in the shearing box simulations, they had no significant impact on the evolution of those simulations. In this global simulation, however, as Figure 9 shows, they play a significant role in the creation of large fluctuations in the accretion rate.

### 3.6 $`R`$$`\varphi `$ Stress

The total stress must satisfy the equation of angular momentum conservation. In a steady-state disk, this means that locally the torque must be equal to the mass accretion rate (which is constant) times the difference between the local specific angular momentum $`\mathrm{}`$ and the specific angular momentum carried into the black hole $`\mathrm{}(r_g)`$. Thus, the stress is equal to

$$\overline{S}=\frac{\dot{M}\mathrm{\Omega }(R)}{2\pi }\left[1\frac{\mathrm{}(r_g)}{\mathrm{}(R)}\right],$$
(16)
where $`\dot{M}`$ is the total mass accretion rate, $`\mathrm{\Omega }(R)`$ the rotational frequency (for the pseudo-Newtonian potential, $`\mathrm{\Omega }(R)=1/[R^{1/2}(r1)]`$), and $`\mathrm{}`$ is the specific angular momentum ($`\mathrm{}R^2\mathrm{\Omega }`$). In the time-steady, axisymmetric “Novikov-Thorne” disk, it is assumed that $`\mathrm{}(r_g)=\mathrm{}(r_{ms})`$ because the stress is assumed to be zero for all $`Rr_{ms}`$.

For the purpose of studying angular momentum transport, it is useful to look at the azimuthally- and vertically-integrated magnetic stress, $`𝑑zM^{R\varphi }_\varphi `$, as this is the quantity most directly related to transfer of the $`z`$-component of angular momentum. In particular, we compare this quantity to the total stress per unit area that would be expected in the steady state “Novikov-Thorne” disk. As discussed above, in the inner part of the disk the mean accretion rate is about 5 in code units. In Figure 10 $`\overline{S}`$ (for $`\dot{M}=5`$ and $`\mathrm{}(r_g)=\mathrm{}(r_{ms})`$) is contrasted with $`𝑑zM^{R\varphi }_\varphi `$ averaged from $`t=750`$ to $`t=1500`$. As that figure shows, in the main body of the disk (approximately $`4R15`$), the average magnetic stress is almost exactly equal to the expected stress. At large radius the two curves diverge because equation (16) assumes a time-steady disk, whereas ours has only a finite mass, so the behavior at its radial exterior is quite different. At the inner edge of the disk, in contrast to the prediction of the conventional model, the azimuthally- and vertically-averaged magnetic stress does not diminish, but, if anything, increases slightly inward (the final upward spike is an artifact of the boundary condition). Note, however, that the local value of the ratio $`|M^{R\varphi }/(\rho v_Rv_\varphi )|`$ in the plunging region varies widely: from as much as $`10`$ several scale-heights out of the equatorial plane to as little as $`10^4`$ at some locations in the plane.

In the context of time-steady and axisymmetric disks, the stress is often discussed in terms of the dimensionless parameter introduced by Shakura & Sunyaev (1973),

$$\alpha _{SS}\frac{𝑑zT^{R\varphi }}{𝑑zp}_\varphi .$$
(17)
In many papers following this seminal work, the “$`\alpha `$ model” is generalized so that the local ratio of stress to pressure is measured by the quantity $`\alpha |T^{R\varphi }|/p`$. Often, it is also assumed that the stress is due to some variety of turbulent viscosity. However, here we follow the spirit of the original Shakura & Sunyaev work by supposing that $`\alpha `$ is merely a convenient dimensionless way of measuring the stress in units of the pressure. We do not particularly care whether the stress is due to chaotic or regular mechanisms. In this spirit, we define $`\alpha |M^{R\varphi }|/p`$. Unlike $`\alpha _{SS}`$, which, in most studies, is taken to be a constant independent of time and location, $`\alpha `$ in this simulation is a dynamical and highly variable quantity that is determined by local and transient magnetic dynamics. For this reason, the absolute value is significant; although $`M^{R\varphi }`$ is positive most places, there are counter-examples.

To give a sense of its scale and range of variation, we first examine $`\alpha _\varphi `$, the azimuthally-averaged version of this ratio, as found in the late-time snapshot (Figure 11). In the main portion of the disk, $`\alpha _\varphi `$ ranges between $`10^2`$ and $`10^1`$. However, it rises dramatically (to $`10`$) at several scale-heights above and below the mid-plane. The same tendency has been seen (with much better vertical resolution) in stratified shearing box simulations (Miller & Stone 2000). $`\alpha _\varphi `$ also increases sharply at small radius. In the mid-plane, $`\alpha _\varphi 0.2`$ at $`R=3`$, but rises to $`3`$ at the inner radial edge of the simulation. Moreover, the vertical rise in $`\alpha _\varphi `$ also becomes sharper at small radius, so that at $`R=3`$, $`\alpha _\varphi 10`$ at $`|z|`$ as small as 0.5.

In addition to these systematic trends, $`\alpha `$ is also subject to sizable fluctuations. Measured in the same way as for the magnetic energy density and the gas density (i.e., in terms of fluctuations in azimuth at fixed $`R`$ and $`z`$), the rms fractional fluctuation amplitude in the bulk of the disk ranges from $`1.5`$ to $`3`$, but at large altitude above the disk midplane, the fluctuations can be much greater. The character of these fluctuations is illustrated by a different projection of $`\alpha `$, which we call $`\overline{\alpha }`$:

$$\overline{\alpha }\frac{𝑑z\rho \alpha }{𝑑z\rho },$$
(18)
the density-weighted vertical average of $`\alpha `$. In this projection (Figure 12), we clearly see the nature of the azimuthal fluctuations: just as for the density fluctuations (shown on smaller scale in Figure 2), the magnetic fluctuations are sheared filaments of limited angular coherence length.

These results highlight the limitations of the primary assumption made in many previous studies of disk structure, that $`\alpha _{SS}`$ is a universal constant, independent of time and location. By contrast, we find that the ratio of stress to pressure varies strongly from place to place, both in the sense of having large fluctuations and in the sense of possessing strong systematic gradients. In many places the gas pressure is relatively smooth while the magnetic stress varies considerably. In other places the variation in $`\alpha `$ is simply due to changes in the background pressure. For example, the magnetic stress does not drop with height as rapidly as the pressure does.

For purposes of comparison with standard disk models, it is possible to compute (albeit with some trepidation) the value of $`\alpha _{SS}`$ found in this simulation. This quantity, defined in terms of the Maxwell stress, is shown as a function of radius in Figure 13 for the disk at the end of the simulation. It is $`0.1`$ through most of the disk, but, starting just outside the marginally stable orbit, it rises inward. It reaches a peak $`1`$ near the inner radial boundary of the simulation. The rise is primarily due to the rapid drop in gas pressure as the accretion flow accelerates inward. Also shown is the ratio of the Maxwell stress to the magnetic pressure (dashed line). This $`\alpha _{mag}`$ is 0.2–0.3 throughout most of the disk, and increases inside $`r_{ms}`$.

Whether $`\alpha _{SS}`$ is meaningful depends on how it is used. As the previous discussion emphasized, the local ratio of magnetic stress to pressure has strong vertical gradients that are completely obscured when the vertically-integrated version is used. Similarly, there are large fluctuations in azimuth and as a function of time. For these reasons, $`\alpha _{SS}`$ is useless as a local dynamical variable. However, the approximate constancy of $`\alpha _{SS}`$ as a function of radius in the disk proper means that it can be used as a rough estimator of the surface density in the body of the disk, provided its use is not extended to local properties nor is it thought to be accurate at better than “factor of a few” accuracy. This conclusion is entirely consistent with many standard uses of $`\alpha `$, and within the spirit with which it was originally introduced by Shakura and Sunyaev (1973) as a scaling parameter, not a fundamental variable.

### 3.7 Accreted angular momentum and energy

Another significant diagnostic of the flow is the mean specific angular momentum distribution. There are a variety of ways to characterize the specific angular momentum; one is by vertical and azimuthal average, $`\mathrm{}`$, as defined in equation (9). This function was discussed in Hawley (2000) for the model GT4 (see Fig. 14 in Hawley 2000); the present model is very similar. At late time the average specific angular momentum is slightly sub-Keplerian throughout most of the disk. The most striking feature is that $`\mathrm{}`$ continues to decline inside $`r_{ms}`$ rather than leveling off. If all stress ceased at the marginally stable orbit, material would travel inward from there retaining precisely its energy and angular momentum. Conversely, the degree to which fluid changes its specific energy and angular momentum is a measure of the stresses that exist in the plunging region. The mean specific angular momentum diminishes from 2.60 at $`r_{ms}`$ to 2.48 at the innermost radius. Thus, on average in this snapshot, the accreted matter has lost $``$5% of the angular momentum it had when it was at the marginally stable orbit. It is likely (for reasons explained in §5) that this number is an underestimate, and in reality the specific angular momentum at the innermost radius is rather smaller.

The total rate at which angular momentum is carried across the inner radial boundary is slightly smaller than these numbers would suggest because there is also a small outward electromagnetic angular momentum flux. This electromagnetic flux is due to a generic property of magnetic fields in shearing plasma: when the orbital frequency decreases outward, the field is swept back so as to transport angular momentum outward, i.e., $`B_RB_\varphi <0`$. However, near the inner radial boundary this is a small effect, the electromagnetic angular momentum flux there is only about 2% of the angular momentum flux carried by the matter.

Further complexity is revealed by a closer examination of the $`(R,z)`$ distribution of the late time azimuthally-averaged angular momentum inside the plunging region. In fact, the specific angular momentum in the cells lying along the equatorial plane is nearly constant inside $`r_{ms}`$. However, the value of $`\mathrm{}`$ drops with $`z`$ above and below the equator. This is consistent with the distribution of magnetic field which, of course, accounts for the nonzero stress. The low value of field near the equator is in part a consequence of the assumed symmetry of the initial distribution of field. The initial field loops are symmetric with respect to the equator ($`B_R`$ and $`B_\varphi `$ are zero there), and throughout the run the strongest fields tend to be located above and below the equatorial plane.

Just as for every other quantity, there is also considerable spiral structure. Figure 14 shows the density-weighted vertically-averaged specific angular momentum. It makes plain how nonaxisymmetric the distribution of torque is in the plunging region.

For a number of reasons, the accreted energy is harder to determine than the accreted angular momentum (see further discussion in §5.2). To introduce this topic, we begin by displaying (Figure 15) the specific energy $`e`$ in rest-mass units,

$$ec^2=\frac{1}{2}v^2+\frac{(\mathrm{\Gamma }1)P}{\rho }+\frac{B^2}{8\pi \rho }\frac{1}{(rr_g)},$$
(19)
in the equatorial plane at late time. The specific energy for a circular orbit precisely at $`r_{ms}`$ is -0.0625; here, $`e_\varphi `$ is a bit greater than that ($`0.055`$) because of the thermal, magnetic, and turbulent kinetic energy contributions to the total energy. In the equatorial plane, $`e_\varphi `$ is almost constant from $`R=3`$ to just outside $`R=2`$, but actually rises as $`R`$ decreases from there to $`R=1.6`$, reaching a peak that is very nearly unbound: -0.01. The azimuthally-averaged energy in the equatorial plane then falls sharply to -0.07. Even this wild swing underestimates the fluctuations: the smallest specific energy found in the equatorial plane is -0.19. Thus, we are confident in asserting that there is a large amount of energy transfer within the plunging region. However, given these large fluctuations, we do not feel that we can compute a well-defined mean energy for the accreted matter. It is likely, moreover, that at least some of the very large fluctuations in the innermost radial zones are not physical, but due to the inner boundary condition.

In fact, the situation is somewhat worse than this. The reason is that dissipation is crucial for computing the accreted energy, and our simulation contains no real dissipation physics. Absent numerical losses (see §5.1, 5.2), energy taken from plunging matter may go several places. One possibility is to build magnetic field energy in the plunging region itself. If the field energy is advected into the black hole without dissipation and radiative losses, the transfer of energy from fluid to magnetic field makes no difference to the accreted energy. It is also possible for magnetic stresses to transfer energy all the way back to the disk proper (and indeed we see significant magnetic stress throughout the plunging region). Again, without dissipation, the energy transferred in that way stays in the fluid, whether in the form of kinetic or magnetic energy. Only in the event that the energy is transferred all the way out to non-accreting matter would it be possible for that energy to avoid being carried back in. Thus, to the degree that this simulation avoids numerical losses, it is very difficult for it to show any significant transfer of energy out of the plunging region. By contrast, in a real disk, dissipation changes energy into heat, which may then be radiated if the thermal timescale is shorter than the inflow timescale; in our simulation, neither physical dissipation nor radiation is present.

## 4 Discussion

### 4.1 Black hole spin-up

We are now ready to begin answering some of the questions posed in the Introduction. Most significantly, this simulation has shown that the stress due to MHD turbulence does not automatically cease at the marginally stable orbit. In fact, it continues to be present throughout the plunging region. As an immediate consequence, matter continues to lose angular momentum as it travels from the marginally stable orbit toward the event horizon. Because it then enters the black hole with smaller angular momentum than if no stresses had been present, the rate at which the black hole is spun up is slower than it would have been otherwise (Agol & Krolik 2000; Gammie 1999). It further follows that estimates of the black hole spin distribution are shifted toward slower rotation (e.g., Moderski & Sikora 1996). Although in this simulation the magnitude of this effect is small ($`5\%`$), it is possible that more realistic conditions will increase its size (see §5.1).

### 4.2 Additional dissipation in the inner disk

We have already remarked on the fact that our simulation does not directly compute the dissipation rate. In this section we attempt to estimate what dissipation might be associated with dynamics of the sort seen here. In conventional disk models, the local dissipation rate is directly tied to the local stress:

$$Q=\left|\frac{d\mathrm{\Omega }}{d\mathrm{ln}r}\right|𝑑zT^{R\varphi },$$
(20)
where $`Q`$ is the dissipation rate per unit area. This relationship is predicated upon the notion that the same stress producing local transport also produces local dissipation. This remains true for the mean flow dynamics of MHD turbulence (Balbus & Papaloizou 1999); in fact, numerical dissipation mimicking the physical dissipation at the short wavelength end of the turbulent cascade is the way in which matter in this simulation loses energy as it moves slowly inward in the body of the disk. However, this identification is not necessarily correct as a prediction for the instantaneous rate of dissipation in specific fluid elements. It is especially questionable in the plunging region. The stress is not necessarily turbulent, and there may not be enough time for the energy deposited in the fluid by the stress to be transferred to small enough lengthscales for it to be dissipated. Nonetheless, in this subsection we adopt the speculative assumption that the vertically-integrated dissipation rate is proportional to the vertically-integrated stress in order to speculate about the consequences for disk radiation that may result.

One possible consequence is an altered radial distribution of the disk surface brightness. As shown in Figure 10, the azimuthally-averaged stress in the simulation is nearly constant inside $`R5`$. If the relation of equation (20) holds, the total dissipation rate outside $`R=3`$ in this simulation would be $`50\%`$ greater than in a disk with the same accretion rate, but zero stress at the marginally stable orbit.

If this were to happen, there would be a number of effects of considerable phenomenological interest (Agol & Krolik 2000). Because this additional heat is released in a relatively small area of the disk, it leads to relatively high temperatures, and therefore contributes to the highest-frequency portion of the spectrum. Another consequence stems from the fact that the dissipation is located in the innermost part of the disk, where relativistic effects are strongest. The light that is radiated due to this extra heating should be comparatively strongly beamed into the equatorial plane by a combination of special relativistic Doppler beaming driven by the orbital motion and general relativistic trajectory-bending. Trajectory-bending also forces a large fraction of the photons to return to the disk at larger radius, to be reprocessed there into lower-frequency light.

Whether or not the enhanced stresses in the plunging region lead to additional radiative losses depends on the dissipation and photon diffusion rates. In the disk proper, dissipation takes place in the turbulent cascade, in small lengthscale fluctuations that are short-lived, exhibit little spatial coherence, but are smoothly distributed in a statistical sense. In the time-varying accretion of the plunging region this turbulent cascade may not have time to go to completion. On the other hand, there might be dissipative events, mediated by reconnection for example, that involve structures with relatively large spatial coherence lengths and which are not necessarily smoothly distributed throughout the region. As we will discuss in §5.1, our limited spatial resolution makes it difficult for us to estimate fairly the rate of such events in the plunging region. Moreover, assessing the radiative properties of such dissipative events is beyond the the capabilities of this simulation. However, if events of this sort do occur, significant luminosity might result because the reduced optical depth of the plunging matter (as compared to the disk proper) also greatly reduces the thermal time. Note that this mechanism results in an enhanced radiative efficiency without the transfer of any energy from the plunging region to the disk outside the marginally stable orbit.

Finally, although the azimuthally-averaged stress matches conventional expectations in much of the disk, there are large azimuthal variations in the form of sheared filaments (Figure 12). Translated into a dissipation rate, these spiral waves of enhanced (or diminished) stress become spirals of enhanced (or diminished) heating. Whether this concentration of heating has an effect on the emitted radiation depends on the characteristic timescale of these fluctuations as seen in a specific fluid element. If a particular fluid element has an elevated heating rate for as long as its thermal timescale, the effective radiating area of the disk is reduced to the area of the regions where the heating is greatest, and the typical effective temperature is raised accordingly. On the other hand, if over a thermal timescale individual fluid elements see both positive and negative fluctuations, the effects will average out. We have already estimated that the characteristic speed of fluctuations through the disk is of order the magnetosonic speed; on the other hand, the thermal timescale is $`\alpha _{SS}^1`$ times the orbital period. Consequently, the ratio of the coherence time of fluctuations to the thermal time should be, very roughly, $`\alpha _{SS}v_{orb}/(v_A+c_s)`$; in this simulation this ratio is $`1`$.

### 4.3 Fluctuations in light output

We have emphasized throughout this paper the highly non-stationary behavior that is intrinsic to disk behavior as found in this simulation. One form these fluctuations can take is the coherent waves that are continually radiated outward through the disk from the region of the marginally stable orbit. Even if, on average, dynamics in the plunging region change the mean accretion efficiency by only a small amount, the large fluctuations we observe in the simulation should cause substantial time-dependent variations in the light output.

Once again, because there is no direct connection between the dynamics traced by this simulation and dissipation, we cannot directly predict the time variation of the disk luminosity. In fact, to do so would also require calculating the photon diffusion rate from inside the (optically thick) disk. However, some indication of what may occur may be gleaned from Figures 6 and 7. Some of the fluctuations in the accretion rate (particularly those at the highest frequencies) are due to dynamics in the plunging region; these may or may not lead to dissipation (see §4.2), and therefore may or may not contribute to fluctuations in the light output. For this reason, and following the ansatz of equation (20), in Figure 7 we also plot the Fourier power spectrum of the volume-integrated magnetic stress. Whether the accretion rate or the stress is more closely related to the dissipation rate, the actual output is the convolution of fluctuations in the heating rate with the probability distribution for the photon escape time. In frequency space, this amounts to a (position-dependent) low-pass filter. When the heating occurs deep inside the disk, the cut-off frequency for this filter is $`\alpha _{SS}\mathrm{\Omega }`$; fluctuations in the light due to heating closer to the surface are subject to a less stringent frequency cut-off (in the optically thin limit, there is no cut-off at all). Consequently, we expect the Fourier power spectrum of the disk luminosity to be rather strongly cut off at frequencies above $`10^2`$, but the power spectrum of fluctuations associated with coronal dynamics may be different.

The observation of QPOs in black hole candidates has motivated a number of studies of possible sources of systematic time variability in accretion disks. In a series of papers (Nowak & Wagoner 1991, 1992; Perez et al. 1997) Wagoner and collaborators have explored the linear theory of fluid oscillations in the relativistic portions of stationary accretion disks. In these disks only hydrodynamic effects are considered, and angular momentum transport is modeled by a pseudo-viscosity parameterized by a constant $`\alpha `$. These disks possess a variety of normal modes, including $`p`$-modes that are trapped between $`r_{ms}`$ and the radius where the epicyclic frequency is maximum, and $`g`$-modes, whose greatest amplitude is found under the peak of the epicyclic frequency curve. In the Paczyński-Wiita potential we use, the position of the maximum epicyclic frequency is $`R3.6`$; with the additional correction terms introduced by Nowak and Wagoner (1991), it moves outward to $`R=4`$. We would not expect to find exactly the same mode frequencies as they because of the difference in the potentials used; however, the mode frequencies should in general be $`\mathrm{\Omega }(r_{ms})`$. A perusal of Figure 7 fails to find any evidence for special frequencies in this range.

There are several reasons why trapped disk oscillations do not appear here. First, the normal modes discussed in the previous paragraphs are predicted on the basis of a purely hydrodynamic analysis; no account is made of magnetic forces. Here, however, magnetic forces play a major role in determining the dynamics of fluid elements. For example, if one repeated the derivation of the “diskoseismic” dispersion relation, but allowed for a weak magnetic field, one would rediscover the magneto-rotational instability, thus vitiating the hydrodynamic “normal modes”. More fundamentally, our simulation indicates that there is no underlying quiet flow to support coherent oscillations. The disk is so turbulent that there is no steady equilibrium against which linear perturbations can develop. Put another way, we find that the typical radial motion timescale in the vicinity of $`r_{ms}`$ (where these normal modes are supposed to be concentrated) is comparable to the predicted mode frequencies, so the fundamental assumption that these modes are perturbations to circular orbits is inappropriate.

## 5 Limitations

The results presented here and in Hawley (2000) obviously represent only the first steps in detailed simulations of global accretion disks. In this section we outline how the approximations we have made may impact our conclusions. Of course the only definitive way to investigate the effects of these approximations would be with additional simulations, and this discussion should be viewed as a guide for future work.

### 5.1 Effects of resolution

Foremost among the limitations is the numerical resolution of the simulation. Although we placed 30 radial zones between $`R=1.5`$ and $`R=4`$, that is probably insufficient. There are two reasons for believing this to be the case. First, we may apply the standard test for convergence of comparing two simulations with the same physics but different resolution, and ask whether they show significant differences. Simulation GT4, reported in Hawley (2000), and a new $`64^3`$ run (hereafter referred to as run LR for Lower Resolution) serve as good comparisons. Both GT4 and LR share identical initial and boundary conditions. GT4 had equally spaced grid zones over a distance of 20 in both $`R`$ and $`z`$ for $`\mathrm{\Delta }R=\mathrm{\Delta }z=0.156`$. Simulation LR used equally spaced radial zones over a range of 30 for $`\mathrm{\Delta }R=0.47`$, and used half of the $`z`$ grid zones centered around the equator within $`z=\pm 2`$, for $`\mathrm{\Delta }z=0.125`$. By contrast, here $`\mathrm{\Delta }R=0.0833`$ and $`\mathrm{\Delta }z=0.0625`$ in the inner region.

All three simulations are qualitatively similar in global properties, as measured by volume integration, but differ in certain other respects that are important to physical interpretation. The time-averaged accretion rate in GT4 is about 40% smaller than in the work reported here (although the time-averaged accretion rate in LR is only about 2% smaller). However, the most striking contrast is in magnetic field strength. Higher resolution results in larger magnetic fields, particularly in the poloidal components. Contrasting maps of the magnetic field intensity at similar late time for GT4 and our new simulation (Figure 16), we find that the field is consistently stronger in the new simulation everywhere inside $`R=5`$, rising to typically a factor of 3 greater in magnitude inside $`R=4`$. Significantly, outside $`R=5`$, the field is generally stronger in the old simulation (sometimes by more than a factor of 10), but that, too, is in keeping with the sense of the resolution contrast: GT4 has better resolution in the outer regions.

The field intensity is controlled by a trade-off between the turbulent MHD dynamo driven by the magneto-rotational instability and numerical dissipation (recall that there is no explicit resistivity in this simulation). Whatever the real physics of resistivity in accretion disks may be, it is unlikely that it is as large as the effective numerical resistivity induced by our modest spatial resolution. The contrast between the field intensity in the two simulations is therefore very likely explained by the smaller numerical resistivity of the better-resolved simulation.

Another piece of evidence supporting the belief that the magnetic field is weaker than it would be with better resolution is presented in Figure 17, where we display the magnetic field structure in the equatorial plane within $`R=4`$, plotted for every fourth grid zone. The planar component of the field is organized into long sheared loops. Because the loops are wrapped around one another, there are many places where the field direction reverses within 4 – 8 grid cells. This means that with a modest amount of flow convergence numerical effects could cause the field to decay on a relatively short timescale.

Finally there is the issue of the computational domain employed, in particular the location of the boundaries. Very little outflow leaves through the $`z`$ boundaries (typically only a few percent of the accretion rate). The outer radial boundary was extended in the present simulation compared to GT4, which had the effect of keeping more of the disk on the grid, but otherwise had no obvious effect on the evolution. We expect that there is also little dependence upon the location of the inner boundary because the flow across that boundary is generally supermagnetosonic. Although we made no tests here of the effect of varying $`r_{in}`$, tests of this sort were made on the similar simulations in Hawley (2000), and it appeared to make little difference.

### 5.2 Energy conservation and physical dissipation rates

The artificial dissipation caused by insufficient resolution also has consequences for energy conservation. To very good accuracy, angular momentum is conserved in this simulation. The primary source of numerical dissipation, the truncation error associated with transport terms in the equation, does not alter the total angular momentum, although it may artificially change its distribution. Total energy, however, is different because the energy equation that is solved is the internal energy equation. Numerical dissipation can destroy kinetic energy in regions of strong velocity gradients and magnetic energy in regions of strong magnetic field gradients, but there is no provision in this simulation for putting that energy back into heat (except for dissipation in shocks through the use of an artificial viscosity). Both sorts of dissipation can also occur in real systems, but at rates that might be very different from what happens in the simulation.

In the body of the disk, the problems that numerical dissipation can cause are limited. The single largest items in the energy budget are potential energy and rotational kinetic energy; neither one is significantly affected by numerical effects. The thermal energy is only a small fraction of those energies; this is a prime motivation for solving an internal rather than total energy equation. Artificial energy losses occur in the dissipation of relatively small amounts of energy bound up in small-scale random fluid motions and in magnetic field. In fact, in this regard, numerical dissipation acts in a way that is consistent with the physical dissipation of energy at the short lengthscale end of the turbulent cascade. In both cases, the losses occur on the shortest lengthscales, and in both cases (assuming the disk radiates efficiently) energy is lost from the disk on timescales short compared to the inflow time.

However, near and inside the marginally stable orbit, detailed questions of energy transfer become significant. The contributions to the energy budget due to both turbulent motions and magnetic field grow. At the same time, their importance also grows. For example, it would be of great interest to determine whether plunging fluid elements can, by exerting magnetic stresses, transfer energy back to the disk proper (cf. §3.7). As we have seen, significant magnetic stresses exist throughout this region, so there is some possibility that this happens. However, the rate at which this may occur is regulated by the magnetic field strength and topology, and these are subject to significant artificial numerical dissipation. In addition, any energy that is delivered is likely to take the form of MHD turbulence; whether this turbulence is dissipated (and the energy radiated) quickly or slowly compared to the inflow is critical for evaluating the radiative efficiency that results. Because the inflow, thermal, and dissipation timescales all change significantly in the vicinity of $`r_{ms}`$, it is unlikely that in this region numerical dissipation serves as a successful stand-in for physical dissipation and radiative losses in the way it does at larger radii.

### 5.3 Equation of state

Dissipation and radiation are, of course, key elements in the gas thermodynamics. In the present simulation, we describe the gas’s internal energy in only an approximate way: we assume it is adiabatic except for heating due to artificial viscosity. Explicit treatment of dissipation and radiation would fill in the missing terms in the gas energy equation.

In the absence of a full treatment of the energy equation, we can at best make educated speculations about the impact of our assumption that the gas behaves adiabatically. Some consequences were already discussed in the previous subsection. There are several more, not as closely related to resolution issues.

As we have remarked earlier (§3.7), one of the most dramatic features of this simulation is the coherent disturbances that propagate as magnetosonic waves from the vicinity of the marginally stable orbit outward through the disk. Their behavior depends strongly, of course, on the sound speed in the disk. However, in this simulation, in which all gas is given the same specific entropy initially and then evolved adiabatically, the sound speed could be quite different from what it would be in a more realistic simulation that includes heating from dissipation of MHD turbulence and cooling via radiation. We expect, therefore, that properties that depend on pressure forces (such as the propagation of these waves, or the vertical thickness of the disk) are not well-described in this simulation.

Another consequence of our assumed equation of state, combined with the particular value of specific entropy we assigned the gas, is that the disk is modestly thick. As a result, it begins from a non-Keplerian, partially pressure-supported state and radial pressure gradients drive some of its initial evolution. Magnetic stresses rapidly change the angular momentum distribution within the torus to nearly Keplerian, but without cooling, the internal pressure continues to push matter away from the pressure maximum. Some accretion, therefore, is driven not by magnetic stresses but simply by the pressure gradient. At late time, however, the inner disk establishes a quasi-steady state and nearly all the accretion can be accounted for by the Maxwell stress. It is likely that the observed dynamics are reasonably representative of generic inner disk accretion.

Finally, although the present disk is moderately thin, with $`H/R0.15`$ in the inner region, it is conceivable that new phenomena might appear in a still thinner disk. Although our results provide no particular reason to predict this, it is possible that in a disk with still smaller $`H/R`$ there might be modes with wavelength large compared to $`H`$ yet still small compared to $`R`$. Questions such as these could be investigated by contrasting the present work with simulations of the inner regions of more extended, thinner, initially Keplerian disks.

A possible example of this sort of comparison is provided by the work of Armitage et al. (2000). Like us, they also simulated the inner region of an accretion disk in a pseudo-Newtonian potential using a three dimensional MHD code. They find a smaller reduction of $`\mathrm{}`$ inside $`r_{ms}`$ than we find in the simulation reported here, an effect they conjecture could be due to the smaller temperature and scale height in their disk. However, their scale height is $`0.08`$, not much smaller than our typical value of $`0.15`$. Moreover, there are other differences between the two simulations that may also account for quantitative contrasts: in their stratified simulation, only a limited vertical extent is simulated; their equation of state is isothermal; and they simulate a $`\pi /3`$ wedge, rather than a full $`2\pi `$ in azimuth. Because in our simulation fluctuations with azimuthal wavelength $`1`$ radian or more have considerable amplitude, the limitation to a wedge could have significant consequences. In any event, as they emphasize, the differences between their simulations and ours are simply quantitative, not qualitative. They, too, find a nonzero Maxwell stress inside $`r_{ms}`$, but its amplitude is lower. Further work is required to determine what sets the stress levels in a variety of different situations.

Thus, for all these reasons, it will be important for future simulations to incorporate improved treatment of gas heating and cooling. It should not be difficult to integrate a total energy equation along with the gas momentum (i.e., force) equation. Heating could be treated in a controlled way by inserting an explicit resistivity into the magnetic field evolution equation and an explicit viscosity into the fluid force equation, and placing corresponding terms in the internal energy equation. Radiative losses could be described with either an escape probability formalism, or, in the longer term, solution of the radiation transfer problem.

A simulation with both heating and radiation would simulate the behavior of thin, radiative disks. There has been no previous work of this sort; the significance of such an effort is obvious.

Heating alone would, in effect, lead to a simulation of the popular advection dominated accretion flow (ADAF) model. Some previous simulations of this model have used a phenomenological prescription for viscosity (e.g., Igumenshchev & Abramowicz 1999; Igumenshchev, Abramowicz, & Narayan 2000; Stone, Pringle & Begelman 1999). Others have computed the stress from a global magnetic field (in axisymmetry: Uchida & Shibata 1985 and Stone & Norman 1994; in three dimensions: Matsumoto 1999). Still another approach has been to study the stress created by MHD turbulence including resistive dissipation, but in axisymmetry (Stone & Pringle 2000). Addition of explicit dissipation to the simulation presented here would, in effect, extend the Stone & Pringle work to three dimensions.

### 5.4 Duration

The timestep of the simulation is limited by the Courant condition. The most stringent limitation comes from the high orbital velocities and the radial plunging speeds in the inner part of the grid where the resolution is reasonably fine. The total duration of the simulation was 1500 time units, requiring almost 400,000 timesteps. As discussed in §3.1, this is long enough for a quasi-steady state to be established in the inner part of the disk. Since $`\alpha 0.1`$, the simulation covers several accretion times in the inner part of the disk. Whether or not the results obtained are truly representative, must, however, be tested with simulations of longer duration. An evolution beginning with a nearly Keplerian disk that extends out to large radius could allow the inner disk region to establish itself self-consistently from a long-term accretion flow.

### 5.5 Field topology

What are the limitations due to the choice of initial disk and field topology? The initial poloidal field loops used here are amenable to rapid development of the MRI on large scales. Further, the presence of significant initial radial field quickly generates through simple shear amplification a toroidal field of near-equipartition strength. Thus, the initial field topology guarantees magnetically-driven evolution within a few orbits. In earlier work (Hawley 2000) it was shown that in this regard (generating magnetically-driven inflow after a few orbits), initially poloidal and initially toroidal field configurations are essentially equivalent. However, the symmetry inherent in this choice of initial field topology is partly responsible for the low field intensity near the equatorial plane and the weak dipolar structure characteristic of the low-density regions (§3.3).

Whether the field is contained entirely within the problem area, or, equivalently, whether there is zero net flux through the box, has a greater qualitative impact. Shearing box simulations (Hawley, Gammie & Balbus 1995, 1996) also found strong similarities between simulations with initial toroidal field and simulations with random initial fields. However, they found dramatically different behavior when there was net vertical field. Such field topologies tend to produce larger net angular momentum transport and more violent evolution. These effects have also been observed in global simulations (Matsumoto 1999).

### 5.6 Displacement current

The equations of MHD as we have written them (see especially equation 2) omit the displacement current term. This approximation is valid so long as the Alfvén speed $`v_A<c`$. This condition is satisfied almost everywhere in this simulation, with the only exception a small region near the inner radial boundary ($`1.5R2`$ and $`1|z|2`$). In the equatorial plane, where most of the mass flows, the largest value of $`v_A0.3`$, and that occurs only at the inner radial boundary; everywhere outside $`R=2`$ in the midplane, $`v_A0.15`$. Future simulations in which the displacement current term is included should be somewhat more accurate, but we doubt that this correction will cause any major changes.

### 5.7 Relativity

Obviously, it would be preferable to run the simulation using genuine general relativistic dynamics, rather than the approximation of Newtonian physics in a Paczyński-Wiita potential. We might hope that this approximation is qualitatively reasonable when attempting to mimic dynamics around a non-rotating black hole; after all, the most important aspect of this case is the existence of dynamically unstable orbits inside a critical radius and enforced inward radial motion within a still smaller distance from the center. Further, Gammie & Popham (1998) compute the growth rate of the MRI in full relativity and have shown that it remains undiminished $`r_{ms}`$. The Paczyński-Wiita potential, however, has no representation of the frame-dragging created by a rotating black hole, so it cannot be a good approximation to accretion dynamics around a Kerr black hole. Gammie (1999) and Gammie & Popham (1998) examine some of the interesting possibilities created by the Kerr geometry, but examination of those effects must await a code with true relativistic dynamics.

## 6 Conclusions

The simulation reported here represents another step towards more realistic simulations of accretion disk dynamics. By improving the grid resolution in the inner region, we have shown that magnetic effects become increasingly important near and inside the radius of the marginally stable orbit. Contrary to the most widely-held assumption, the $`R`$$`\varphi `$ component of the stress tensor does not go to zero at $`r_{ms}`$; averaged over azimuth, height, and time, it is almost constant as a function of $`R`$ for $`R2r_{ms}`$. As a result, matter in the region of unstable orbits does not follow simple energy- and angular momentum-conserving free-fall trajectories, although we are not yet ready to make quantitative predictions about how much the energy and angular momentum change.

The ratio $`T^{R\varphi }/P\alpha `$ is commonly used as a measure of the magnitude of the stress, and almost as commonly is thought to be a constant, independent of location and time. Our simulation showed clearly that while it can provide a good qualitative measure of the stress when averaged over several orbits in time, at any particular instant $`\alpha `$ has large systematic gradients (as a function of both $`R`$ and $`z`$) and substantial temporal fluctuations. On average the value of $`\alpha `$ in the disk is $`0.1`$.

This simulation has shown that substantial time- and space variations are the rule. The inner regions of accretion disks around black holes are highly turbulent, non-steady systems. Large-amplitude fluctuations sheared into spiral fragments come and go, many of them initiated in the region immediately outside the marginally stable orbit. As a result, there are continuing fluctuations in the accretion rate; we expect that some of these fluctuations will be mirrored in the light curves of accreting black holes.

Future simulations, guided by these results, will probe the beyond the limitations of the present simulation. Within the context of the same pseudo-Newtonian potential and the present adiabatic equation of state, we can explore the effects of different initial disk configurations and magnetic field strengths and topologies, at higher resolution and for longer time. Beyond this, it will be important to incorporate improved physics into the models. These include general relativistic dynamics and improved treatments of the equation of state and dissipative processes. This future work should bring us still closer to a quantitative understanding of the dynamics of accretion onto black holes.

###### Acknowledgements.

JHK would like to thank Eric Agol for numerous helpful conversations, and Nancy Levenson for invaluable instruction in the wiles of IDL. He would also like to thank Giora Shaviv for comments on the difference between heating rate fluctuations and luminosity fluctuations. JFH thanks Charles Gammie for useful comments on the manuscript. This work was supported by NSF grant AST-0070979, and NASA grants NAG5-9266 and NAG5-7500 to JFH, and NSF Grant AST-9616922 and NASA grant NRA-99-01-ATP-031 to JHK. Simulations were carried out on the Cray T3E system of the San Diego Supercomputer Center of the National Partnership for Advanced Computational Infrastructure, funded by the NSF.

Figure Captions


format: text

text:
# Measurements of the Mass, Total Width and Two-photon Partial Width of the 𝜂_𝑐 Meson

## Abstract

Using 13.4 $`\mathrm{fb}^1`$ of data collected with the CLEO detector at the Cornell Electron Storage Ring, we have observed 300 events for the two-photon production of ground-state pseudo-scalar charmonium in the decay $`\eta _cK_S^0K^{}\pi ^\pm `$. We have measured the $`\eta _c`$ mass to be (2980.4 $`\pm `$ 2.3 (stat) $`\pm `$ 0.6 (sys)) MeV and its full width as (27.0 $`\pm `$ 5.8 (stat) $`\pm `$ 1.4 (sys)) MeV. We have determined the two-photon partial width of the $`\eta _c`$ meson to be (7.6 $`\pm `$ 0.8 (stat) $`\pm `$ 0.4 (sys) $`\pm `$ 2.3 (br)) $`\mathrm{keV}`$, with the last uncertainty associated with the decay branching fraction.

In this Letter, we report a study of two-photon production of the ground-state pseudo-scalar charmonium, i.e., $`\gamma \gamma \eta _c`$. The two space-like photons are radiated by $`e^+`$ and $`e^{}`$ beams, each at an energy of approximately 5.3 GeV. The charmonium spectrum is an ideal testing ground for quantum chromodynamics (QCD) calculations, and producing C-even charmonium states through $`\gamma \gamma `$ fusion provides a clean environment for this purpose.

The two-photon partial width of the $`\eta _c`$ meson can be expressed in next-to-leading order (NLO) perturbative QCD (PQCD), in terms of the $`e^+e^{}`$ partial width of the $`J/\psi `$ meson, as

$$\frac{\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}}{\mathrm{\Gamma }_{ee}^\psi }=\frac{4}{3}(1+1.96\alpha _s/\pi )\times \frac{|\mathrm{\Psi }_{\eta _c}(0)|^2}{|\mathrm{\Psi }_\psi (0)|^2}.$$
(1)
Using the world average value of $`\mathrm{\Gamma }_{ee}^\psi `$, a value of the strong coupling constant $`\alpha _s`$ evaluated at the charm mass scale of (0.28 $`\pm `$ 0.02), and the assumption that the two 1S wave functions, $`\mathrm{\Psi }`$, are the same at the origin, this relationship predicts $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}=(8.2\pm 0.6)`$ keV.

The total width of the $`\eta _c`$ meson can be assumed to be dominated by its two-gluon component, i.e., $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}\mathrm{\Gamma }_{gg}^{\eta _c}`$. The ratio of the rates for $`\eta _cgg`$ and $`\eta _c\gamma \gamma `$ is an especially clean prediction of PQCD because the dependencies of these rates on the wave functions and non-perturbative factors are identical in the numerator and denominator. The ratio depends only on the coupling constants and has been calculated in NLO,

$$\frac{\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}}{\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}}\frac{9\alpha _s^2}{8\alpha ^2}\times \frac{(1+4.8\alpha _s/\pi )}{(13.4\alpha _s/\pi )}.$$
(2)
Using the value of $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}`$ estimated in NLO gives $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}`$ as (28 $`\pm `$ 6) MeV; using instead the world average value of $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}`$, one obtains an estimate of $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}`$ as (26 $`\pm `$ 6) MeV. A calculation with fully relativistic decay amplitudes and a sophisticated QCD potential model predicts $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}23`$ MeV. The current world average of $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}=13.2_{3.2}^{+3.8}`$ MeV disagrees with these theoretical expectations. A precise measurement of the full width and two-photon partial width of the $`\eta _c`$ is important for the verification of these PQCD calculations and approximations.

In the two-photon process $`e^+e^{}`$ $``$ $`e^+e^{}\gamma \gamma `$ $``$ $`e^+e^{}\eta _c`$, the photon propagators dictate that the cross section naturally peaks at low momentum transfer, so the photons are almost real (“on shell”). The incident leptons are scattered at very low angles and continue traveling down the beam pipe undetected. Such “untagged” events typically have low net transverse momentum and low visible energy. The production of the $`\eta _c`$ meson in this untagged two-photon process was searched for in the $`K_S^0K^{}\pi ^\pm `$ decay mode.

The data used in this study correspond to an integrated luminosity of 13.4 $`\mathrm{fb}^1`$ and were collected with two configurations (CLEO II and CLEO II.V) of the CLEO detector at the Cornell Electron Storage Ring (CESR). Approximately one third of the data were taken with the CLEO II configuration. The detector components most useful for this study were the concentric tracking devices for charged particles, operating in a 1.5T superconducting solenoid. For CLEO II, this tracking system consisted of a 6-layer straw tube chamber, a 10-layer precision drift chamber, and a 51-layer main drift chamber. The main drift chamber also provided measurements of the specific ionization loss, $`dE/dx`$, used for particle identification. For CLEO II.V, the straw tube chamber was replaced by a 3-layer, double-sided silicon vertex detector and the gas in the main drift chamber was changed from a 50:50 mixture of argon-ethane to a 60:40 helium-propane mixture. These changes gave rise to significant improvements in the momentum and $`dE/dx`$ resolutions for charged tracks. Photons were detected using the high-resolution electromagnetic calorimeter consisting of 7800 CsI crystals. The Monte Carlo simulation of the CLEO detector response was based upon GEANT. Simulated events were processed in the same manner as the data to determine the $`\eta _cK_S^0K^{}\pi ^\pm `$ detection efficiency and the $`K_S^0K^{}\pi ^\pm `$ mass resolution at the $`\eta _c`$ mass.

The $`K_S^0`$ vertex was reconstructed from its decay to $`\pi ^+\pi ^{}`$ and was required to be displaced from the $`e^+e^{}`$ interaction point; the amount of this displacement varied with detector configuration but was $``$1.5 mm. The $`K_S^0`$ candidate was also required to be within 4 standard deviations ($`\sigma `$) of the known $`K_S^0`$ mass; here $`\sigma `$ was determined on an event-by-event basis from the momenta measurements. Furthermore, the $`K_S^0`$ momentum vector was required to point back to the interaction point. Of the two remaining charged tracks, the $`K^{}`$ and $`\pi ^\pm `$ candidates, the one with lower momentum was typically uniquely identified using the particle’s specific ionization ($`dE/dx`$). This fixed the identity of the only remaining unidentified track, because the presence of the $`K_S^0`$ dictated that exactly one of these two be a kaon to conserve strangeness in the event.

The background from processes other than two-photon production was suppressed by requiring that the $`\eta _c`$ candidate have net transverse momentum less than 0.6 GeV/$`c`$ and that visible energy in the event be less than 6 GeV. Also, because the final state had no expected energy deposits in the calorimeter from neutral particles, the total calorimeter energy in the event not matched to charged tracks was required to be less than 0.6 GeV.

For the mass measurement only, we restricted ourselves to those events in which the $`K^{}`$ and $`\pi ^\pm `$ daughters of the $`\eta _c`$ candidate traversed all layers of the tracking volume. The $`K_S^0`$ daughter pions were not required to satisfy the same criterion, in that the kinematic fitting of the $`K_S^0`$ decay corrected for any possible momentum mis-measurement of the daughter tracks. We did not make any such requirements while measuring $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}`$ and $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}`$, because these quantities are relatively insensitive to precise measurements of the track momenta; the distribution of candidate invariant masses for the determination of these two quantities is shown in Fig. 1.

We fitted the background with a power law function ($`AW_{\gamma \gamma }^n`$, with $`W_{\gamma \gamma }`$ the $`K_S^0K^{}\pi ^\pm `$ invariant mass and $`A`$ a multiplicative constant) and the signal with a spin-0 relativistic Breit-Wigner function (describing the natural line shape) convolved with a double Gaussian function (to take into account the detector resolution). The parameters for this double Gaussian were obtained from a Monte Carlo sample that had the $`\eta _c`$ intrinsic width set to zero. We performed a simultaneous, binned, maximum-likelihood fit to the invariant mass distributions for the CLEO II and CLEO II.V datasets, constraining the physical variables $`M_{\eta _c}`$, $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}`$, and $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}`$ in the two datasets to be the same. The constraint on $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}`$ was accomplished by requiring the ratio of the fitted yields to be the same as the ratio of the integrated luminosities of the two data sets times the efficiencies as determined from our simulations. The invariant mass resolution was approximately 9 MeV in CLEO II and 7 MeV in CLEO II.V. The bin width for fitting was chosen as approximately the average of these two resolutions.

As noted above, two separate sets of fits were performed, one for the measurements of $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}`$ and $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}`$ and another for the determination of the $`\eta _c`$ mass. The full width and yield were obtained from the distributions shown in Fig. 1, with the total observed yield being $`N_{obs}=300\pm 32`$. The fit to the width gives $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}=(27.0\pm 5.8)`$ MeV, with the uncertainty being only statistical. The two-photon partial width was determined by first correcting for the detector efficiency, $`ϵ`$, and then dividing by the number of events expected, $`N_1`$, for a two-photon partial width of 1 keV:

$$\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}=N_{obs}/(ϵN_1).$$
(3)
The quantity $`N_1`$ was determined using

$$N_1=_{\eta _c}_{K_S}\sigma _{e^+e^{}e^+e^{}\eta _c};$$
(4)
$``$ is the integrated luminosity. The cross section for $`e^+e^{}e^+e^{}\eta _c`$ was obtained from Monte Carlo simulation, using the formalism of Budnev et al. and setting $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}=1`$ keV; this choice of this value has no effect on the extracted result. Also, $`_{\eta _c}(\eta _cK_s^0K^{}\pi ^\pm )`$ and $`_{K_S}(K_S^0\pi ^+\pi ^{})`$, with the world average values used. This procedure gives $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}=(7.6\pm 0.8)`$ keV, with this uncertainty coming from statistics only.

The mass was obtained from fits to the more restrictive set of events, as described above, with a total signal size of 195 $`\pm `$ 24, yielding $`M_{\eta _c}=(2980.4\pm 2.3)`$ MeV, the uncertainty being statistical only.

Possible sources of systematic uncertainty for the measured mass, full width, and two-photon partial width were studied. The results are summarized in Table I, in which the individual uncertainties are added in quadrature to obtain the total systematic uncertainty.

The mass calibration of our detector was checked by measuring the masses of the well known $`K_S^0`$, $`\varphi `$, $`D`$, and $`J/\psi `$ mesons using decay modes involving only charged tracks. The measurements were found to be in good agreement with their respective world averages when we limited ourselves to events in which the charged tracks traversed all layers of the tracking volume. The fitted mass including events outside this more restrictive detector volume is consistent with the value we quote, but with substantially larger systematic uncertainties. Our particle identification procedure did not introduce any significant systematic bias to the mass measurement. The measured mass was also insensitive to variations in the signal shape used to fit the $`\eta _c`$ resonance.

The systematic uncertainties in the measurement of the full width were dominated by effects due to the mass resolution, the particle identification procedure, and the signal shape used in the fit. Our ability to predict the actual mass resolution was tested by studying the reconstructed $`D`$ mesons in Monte Carlo simulation and data; the agreement was found to be better than 0.1 MeV.

The particle identification procedure was unable to distinguish between a charged kaon and a charged pion if the track momentum was above 0.8 GeV/$`c`$, for which the expected ionization losses are nearly equal for the two species. This led to a broadening of the reconstructed resonance and was taken into account by the wider Gaussian of the double Gaussian resolution function. This effect was limited to less than 5$`\%`$ of the events. We estimated the possible uncertainty due to this mis-assignment of particle species by completely removing the fraction of events having two possible $`\eta _c`$ candidates and assigning the corresponding change in the measured width as the systematic error.

The accuracy of the fitting method in extracting the Breit-Wigner width of the resonance was checked by extracting the $`\eta _c`$ widths from sets of simulation events generated with different intrinsic widths. We varied the parameters of the signal shape used to fit the $`\eta _c`$ resonance within their uncertainties, derived from a comparison of the fit to $`D`$ meson decays in Monte Carlo simulation and data, to estimate the effect on the measured width. The measured width from the more restrictive set of events used for mass measurement was within 0.2 MeV of the corresponding measurement using the full sample, and no significant correlation was found between the measured mass and full width.

There were several sources of uncertainty for the estimation of the efficiency, which in turn affected the measurement of $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}`$. These were dominated by the uncertainties in the tracking and trigger efficiencies. The effect of a possible presence of resonant substructure ($`K^{}\overline{K}`$) in $`\eta _c`$ decay was studied and was found to give an insignificant variation in the detection efficiency. We estimated the uncertainty in the measured partial width from this effect by considering the possibility that all the $`\eta _c`$ mesons decay through $`K^{}\overline{K}`$. Our initial investigation showed roughly a third of the $`\eta _cK_s^0K^{}\pi ^\pm `$ events proceed via $`K^{}`$(1430) with no evidence for any $`K^{}`$(892); a detailed analysis of this substructure is beyond the scope of this Letter.

The selection requirements on total visible energy and unmatched energy clusters were shown by simulation to be essentially 100% efficient for our signal process and free of systematic bias. Possible bias from the transverse momentum requirement was investigated by changing the nature of the form factor in the simulation and shown to also be negligible.

In our analysis, we investigated the possible effects of interference between the $`K_S^0K^{}\pi ^\pm `$ resonant and non-resonant final states. From the distribution of net transverse momentum for events in the sidebands of the signal lineshape, we estimated that one third of the background events were not of the type $`\gamma \gamma K_S^0K^{}\pi ^\pm `$; these included events with at least one missing particle ($`\pi _0`$, $`\gamma `$) as well as events of the type $`e^+e^{}`$ hadrons, $`e^+e^{}\tau ^+\tau ^{}`$, and $`\gamma \gamma \tau ^+\tau ^{}`$. Such events cannot have interference with our signal $`K_S^0K^{}\pi ^\pm `$ events. Study of the helicity angle distributions in the $`\eta _c`$ rest frame indicated that the sideband events predominantly have $`J=2`$ (or higher) while our signal has $`J=0`$. Due to the preferential production of states with natural parity in two-photon un-tagged processes, a majority of the remaining background events were expected to have the natural spin-parity ($`0^+`$) compared to the unnatural spin-parity ($`0^{}`$) of the signal events. The acceptance of our detector is symmetric in polar angle and uniform in azimuth, making the interference between these states of opposite parity vanish. We did not include any possible effects due to interference on the measured mass, full width and two-photon partial width. Further, the signal shape showed no distortions and the goodness of fit to the hypothesis that ignored interference was very good, as shown in Fig. 1.

In summary, we have measured the mass, full width, and two-photon partial width of the $`\eta _c`$ produced in two-photon collisions. The mass measurement of (2980.4 $`\pm `$ 2.3 (stat) $`\pm `$ 0.6 (sys)) $`\mathrm{MeV}`$ compares well with the world average of (2979.8 $`\pm `$ 2.1) $`\mathrm{MeV}`$. The measured total width of (27.0 $`\pm `$ 5.8 (stat) $`\pm `$ 1.4 (sys)) $`\mathrm{MeV}`$ disagrees with the world average of ($`13.2_{3.2}^{+3.8}`$) MeV, which consists of measurements with large relative uncertainties (40$``$100$`\%`$). Our measured width is consistent with theoretical expectations. The measured two-photon partial width of the $`\eta _c`$ meson of (7.6 $`\pm `$ 0.8 (stat) $`\pm `$ 0.4 (sys)) keV agrees well with the world average of $`(7.5_{1.4}^{+1.6})`$ keV and theoretically expected values, and is a significant improvement in terms of experimental precision. We use the world average of the $`\eta _cK\overline{K}\pi `$ branching fraction of $`(5.5\pm 1.7)\%`$. The uncertainty in $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}`$ due to the uncertainty in this branching fraction is $`\pm `$2.3 keV and is stated separately from the other contributions. From the ratio of our measured full width and two-photon partial width, we have extracted $`\alpha _s`$ at the charm mass scale to be $`0.285\pm 0.025`$, for which we have added our sources of uncertainty in quadrature. We have used the NLO calculation in Eq. 2 to estimate $`\alpha _s`$, thus making the result dependent on renormalization scheme and scale; we have not included such theoretical uncertainties in our quoted value.

Our measurements of $`\mathrm{\Gamma }_{\mathrm{tot}}^{\eta _c}`$ and $`\mathrm{\Gamma }_{\gamma \gamma }^{\eta _c}`$ show that PQCD calculations are able to reliably predict the ratios of the decay rates of a heavy quarkonium system, where non-perturbative effects cancel.

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. This work was supported by the National Science Foundation, the U.S. Department of Energy, the Research Corporation, the Natural Sciences and Engineering Research Council of Canada, the A.P. Sloan Foundation, the Swiss National Science Foundation, the Texas Advanced Research Program, and the Alexander von Humboldt Stiftung.