text:
# The Dynamical Structure of Nonradiative Black Hole Accretion Flows

## 1 Introduction

Many black hole accretion sources are remarkably underluminous, and pose a stern challenge for accretion theory. Generally, these sources emit a significant fraction of their radiated energy in the form of X-rays (Elvis et al. 1994), and standard, optically thick Keplerian disk theory cannot account for this type of spectrum. In a classical Keplerian accretion disk (Shakura & Sunyaev 1973), dissipated mechanical energy is lost as local blackbody radiation, and there is a simple relationship between the accretion rate $`\dot{M}`$ and the luminosity $`L`$. There is little freedom to adapt the model to more complex spectra.

Another class of models is characterized by inefficient cooling and high temperatures (e.g., Shapiro, Lightman & Eardley 1976; Ichimaru 1977). One particular model that of late has enjoyed considerable attention is known as an Advection Dominated Accretion Flow, or ADAF for short. The epithet ADAF originated with the influential paper of Narayan & Yi (1994; hereafter NY94), which presented a one-dimensional, time-stationary, self-similar solution (fixed power-law dependencies in $`r`$ for dynamical variables). All velocities —radial, nonradial and thermal—are assumed to have the same Keplerian $`r`$ dependence, $`r^{1/2}`$. Angular momentum transport is effected by a Navier-Stokes viscosity along the lines of the Shakura & Sunyaev (1973) and Lynden-Bell & Pringle (1974) $`\alpha `$ formalism. The gas is hot, with the sound speed comparable to the orbital speed. The internal energy is advected along with the mass into the central black hole; very little is radiated. An ADAF is more like a Bondi flow than a Keplerian disk. The quasi-spherical nature of the flow is often cited as a defining characteristic feature. ADAFs are sufficiently distinct from a Keplerian disk that a transition radius $`r_t`$ is invoked to define the location separating the (outer) Keplerian flow from the (inner) ADAF.

There are two dynamical features of the NY94 self-similar ADAF solution that have attracted special attention. The first is that the specific energy flux of the gas (the Bernoulli parameter) is everywhere positive, the second is that entropy increases inwards, rendering the accretion convectively unstable by the Schwarzschild criterion. We review these in turn.

In a nonradiative accretion flow a positive Bernoulli parameter is a consequence of the fact that the stress transports energy, as well as angular momentum, outwards. In a standard thin disk solution, this transported energy is radiated on its outward journey, so that the local radiated flux (outside of the inner boundary layer) is always greater by a factor of 3 than the energy released locally. By contrast, in a nonradiative accretion flow the transported energy is retained, locally raising the specific energy. If the flow originates at large radius as a marginally bound fluid, the Bernoulli parameter will accordingly be everywhere positive. NY94 noted that a positive Bernoulli parameter might make outflows possible, but since an outflow is very distinct from their ADAF solution, the significance of this observation is unclear. Blandford & Begelman (1999) argued that this property is likely to change fundamentally the dynamical flow properties of the nonradiative solution, and proposed that not only were outflows possible, but that the bulk of the gas and energy would be carried off by a wind. They refer to this concept as an Adiabatic Inflow-Outflow Solution, or ADIOS.

More recently, attention has focused upon the increasing-inward entropy profile that is a straightforward consequence of viscous dissipation, and is a property of any nonadiabatically heated, nonradiative inflow. Narayan & Yi (1995) argued that radial convection would inevitably arise, and serve as a source of outward angular momentum transport—either as the self-consistent source of $`\alpha `$, or at a minimum, as a supplement to it. Narayan, Igumenshchev, & Abramowicz (2000) later explored the consequences of convection moving angular momentum either outward or inward. In one scheme, inward angular momentum transport by convection is envisioned to be capable of balancing the outward $`\alpha `$ transport, creating zero (net) mass and angular momentum flux. Quataert & Gruzinov (2000) argued that convection would produce a state of near marginal convective stability rather like stellar convective zones, described by the classical Høiland criteria. Quataert & Gruzinov referred to this new variation as a Convective Dominated Accretion Flow (CDAF), though formally there is no mass accretion at all!

We regard the defining property of a specific flow to be the fate of the accreting mass and energy. In the Shakura-Sunyaev thin disk model, liberated mechanical energy is both transported outwards and radiated locally. In an ADAF, this liberated energy is advected with the flow through the event horizon of the black hole. In a CDAF, the accretion flow is effectively stifled with no significant mass inflow or outflow. What little energy is extracted from the inner region of the flow is simply transported outward by convection to large radius. In an ADIOS the net accretion is very small because systematic outflows remove the bulk of the mass and energy. The existence of what seem to be under-luminous X ray sources lends credence to the existence of nonradiative flows. The question is, are they best described as ADAFs, CDAFs, ADIOS, or something else?

While flow energetic properties are likely to be complex and (for the moment) highly model-dependent, gross dynamical features—inflow versus outflow, rotational versus pressure support—should be more robust. In this work, we report the results of three-dimensional magnetohydrodynamic (MHD) simulations that follow the accretion of rotating, nonradiating, magnetized gas onto a black hole. We find that the bulk of the mass in nonradiative flows is not slowly rotating and quasispherical like an ADAF, but is much closer to Keplerian: rotationally supported and disk-like. Most of the liberated energy is carried outward by streaming backflow and turbulent transport. The nature of the flow is entirely determined by MHD turbulence, triggered by the magnetorotational instability (MRI; Balbus & Hawley 1991).

An outline of the paper is as follows. In §2 we discuss the numerical model for the simulations. In §3 we present two such simulations, one with resistive dissipation and one without. In §4 we summarize the main features of our model, and discuss the implications of the simulation for the accreting black hole in Sgr A\*. Our conclusions are given in §5.

## 2 Numerical MHD Simulations: Background

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

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

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

### 2.1 Equations

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

$$\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+𝒬)𝐯+\eta J^2,$$
(3)
$$\frac{𝐁}{t}=\times \left(𝐯\times 𝐁\eta _i𝐉\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, $`𝐉`$ is the current, $`𝒬`$ is an explicit artificial viscosity of the form described by Stone & Norman (1992a), and $`\eta `$ is an anomalous resistivity of the form used by SP, namely

$$\eta _i=C_{res}(\mathrm{\Delta }x_i)^2|J_i|/\sqrt{\rho }.$$
(5)
The constant $`C_{res}`$ needs to be large enough to spread a current sheet out over a few zones, but not so large as to turn the overall flow into a resistive one. We use a resistivity constant of $`C_{res}=0.1`$, as did SP.

The form of $`\mathrm{\Phi }`$ is the pseudo-Newtonian gravitational potential introduced by Paczyński & Wiita (1980),

$$\mathrm{\Phi }=GM/(rr_g),$$
(6)
where $`M`$ is the mass of the central black hole, and $`r_g`$ is the gravitational radius, equivalent to the Schwarzschild radius in general relativity. With this potential, the angular momentum of a circular orbit is

$$l_{Kep}=(GMr)^{1/2}r/(rr_g),$$
(7)
and the binding energy is

$$e=\left(\frac{GM}{2rc^2}\right)\left(\frac{[r2r_g]r}{[rr_g]^2}\right)c^2.$$
(8)
The pseudo-Newtonian potential mimics the dynamically important marginally stable circular orbit of the full Schwarzschild metric (defined by $`dl_{Kep}/dr=0`$) at $`r=r_{ms}=3r_g`$.

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

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

### 2.2 Physical Units and Code Units

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

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

$$r_g=3\times 10^5\left(M/M_{}\right)\mathrm{cm}$$
(9)
where $`M/M_{}`$ is the black hole mass. We also set $`GM=1`$, and the speed of light is $`c=(2GM/r_g)^{1/2}`$, or in code units $`\sqrt{2}`$. The orbital time at $`100r_g`$ is

$$2\pi \times 10^3\sqrt{2}r_g/c.$$

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

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

$$L_{Edd}=1.3\times 10^{38}M/M_{}\mathrm{erg}\mathrm{s}^1,$$
(10)
and the Eddington mass accretion rate of

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

## 3 Three Dimensional NRAF Simulations

A previous paper (Hawley, Balbus, & Stone 2001) presents a prototype MHD NRAF simulation. Here we investigate that simulation (designated F1) in detail, and extend it to later times. We also compare F1 with a simulation that includes resistive dissipation (simulation F1r). The motivation for this is that in CDAF models dissipative heating in the central regions is argued to be responsible for global flow structure (Narayan et al. 2000).

We begin with some technical details. These simulations use $`128\times 32\times 128`$ grid zones in cylindrical coordinates $`(R,\varphi ,z)`$. The radial grid extends from $`R=1.5r_g`$ to $`R=220r_g`$. There are 36 equally-spaced zones inside $`R=15r_g`$, and 92 zones increasing logarithmically in size outside this point. There are 50 $`z`$ zones equally spaced between $`10`$ and $`10`$ with the remainder of the zones logarithmically stretched to the $`z`$ boundaries at $`\pm 60`$. The azimuthal domain is uniformly gridded over $`\pi /2`$ in angle. The radial and vertical boundary conditions are simple zero-gradient outflow conditions; no flow into the computational domain is permitted. The $`\varphi `$ boundary is periodic. The magnetic field boundary condition requires the transverse components of the field to be zero outside the computational domain, while the perpendicular component satisfies the divergence-free constraint.

The initial condition is a constant specific angular momentum ($`l`$) torus with a pressure maximum at $`R=100r_g`$, an inner edge at $`R=75r_g`$, and an outer edge at $`R=153r_g`$. The initial magnetic field consists of poloidal loops lying along isodensity contours. It has a volume-averaged $`\beta `$ value $`(=P_{gas}/P_{mag}`$) of 200. The average specific energy of the magnetofluid in the torus is $`2.3\times 10^3c^2`$, i.e., it is bound. Simulation F1 is run for 6.32 orbits at the initial pressure maximum, or 1,807 orbits at $`r_{ms}=3r_g`$. Simulation F1r begins at a time equal to 3.5 pressure maximum orbits in F1, after a global flow has been established. It continues to a time of 6.38 orbits.

Our simulations, with or without resistive heating, show three principal flow components: (1) a hot, but rotationally-supported, disk extending down to the marginally stable orbit; (2) an extended, low density coronal backflow enveloping the disk; (3) a distinctive, jet-like flow near the hole, that emerges unambiguously in momentum plots. These features are illustrated schematically in Figure 1. We believe that these structures are real and robust, and that they are fundamental generic properties of NRAFs. Physical arguments underlying the origin of NRAF structure are presented in §4.1.

### 3.1 Simulation F1: Ideal MHD

The evolution begins as a thick torus, with $`H/R=0.2`$. The complete evolution is represented in Figure 2 as a series of contour plots in azimuthally-averaged density. During the first orbit, the field in the torus is amplified both by the MRI and by shear, and the flow evolves rapidly. The MRI acts most effectively near the equatorial plane, where the field is predominantly vertical. Here the long-wavelength, nearly axisymmetric modes of the MRI grow rapidly. Significant accretion begins shortly after one orbit, when the magnetic energy in the torus has increased to $`\beta 2`$–10.

Between orbit 1 and 2, low-$`m`$ (azimuthal wavenumber) spiral arms of gas accrete from the inner edge of the torus, forming a vertically thin and very nonuniform “daughter disk.” Strong magnetic fields surround the gas. Due to the initial field topology, a current sheet forms near the equator, which proves unstable to vertical oscillations. At later times, as more gas accretes from the initial torus, the disk fills out and thickens. Inside of $`R=100r_g`$ the NRAF forms a modestly thick, nearly Keplerian disk. Low density material is stripped off this evolving disk to create a backflow. Inside of $`R=100r_g`$, the net mass flux is inward; for $`R>100r_g`$ there is net outflow.

By orbit 2, the mass accretion rate into the central black hole has approached its long-term average value, which is only few percent of the rate at which the torus feeds mass into the inner region. There are two requirements for the gas to accrete into the central black hole, and neither is easy to meet: (1) the angular momentum of the gas must be reduced to a value close to that of the marginally stable circular orbit; (2) the gas also must be relatively cool, or pressure and centrifugal forces will drive the gas back out. In this simulation this is much in evidence: accreting gas must pass through a narrow gap in the equatorial “waist” of the centrifugal barrier. When the gas is hot and geometrically thick, this is most difficult. Instead of accreting, much of the gas splashes off the centrifugal barrier near the hole, either contributing to the formation of a hot, intermittent torus inside of $`R=10r_g`$, or creating a coronal backflow that adds to a growing low density envelope around the Keplerian core disk.

As time advances, the equatorial inflow follows a nearly Keplerian angular momentum distribution, $`l_{Kep}R^{1/2}`$ (see fig. 3). The value of $`l`$ actually remains everywhere slightly sub-Keplerian, except at the innermost radii where it is slightly super-Keplerian (in the hot inner torus). At the end of the evolution, $`l`$ is about $`\mathrm{}<`$ 5% below $`l_{Kep}`$ between $`R30r_g`$ and $`100r_g`$, but for $`R10r_g`$ $`l`$ drops to as much as 15% below the Keplerian value.

Between orbits 4 and 5, the accretion flow matures and becomes reasonably steady, forming a moderately thick disk surrounded by a low density, highly magnetized atmosphere. Between orbit 5 and the end of the run (6.3 orbits), the properties in most of the disk do not change by large amounts. The exception to this general trend is the inner torus ($`R10r_g`$), which sporadically depletes and reforms.

Near the equator, thermal pressure dominates magnetic, i.e., $`\beta >1`$, whereas the surrounding region is strongly magnetized, $`\beta <1`$ (fig. 4). The thermal pressure exponential scale height is $`H7`$–10$`r_g`$, decreasing rapidly inside $`R=10r_g`$. The total pressure (thermal plus magnetic) is much smoother, and has an exponential scale height $`H20r_g`$ at $`R=100r_g`$; $`H`$ decreases slowly inward. The gas density, vertically averaged over one gas pressure scale height, is nearly constant from $`R=50r_g`$ down to $`R=10r_g`$, where it increases rapidly to a peak at $`R=5.5r_g`$. Inside of $`R=10r_g`$ the disk resembles a thickened torus.

The disk and coronal regions can also be distinguished in Figure 3b which shows contours of specific angular momentum. Within the disk the gas is approximately rotating on cylinders, characteristic of a nearly barotropic equation of state. In the corona, surfaces of constant specific angular momentum are swept back, consistent with conserved $`l`$ along outflowing streamlines.

Figure 5 shows contours of entropy, $`S\mathrm{ln}(P\rho ^{5/3})`$, and the entropy along the equatorial plane at the end of the F1 simulation. Note that because of nonadiabatic heating by the artificial viscosity, the entropy increases inward even in the absence of resistive heating.

Magnetic fields in the disk produce a significant Maxwell stress which results in turbulent angular momentum transport. Between orbit 4 and the end of the simulation, the ratio of the vertical- and $`\varphi `$-averaged Maxwell stress to the gas pressure (a form of the $`\alpha `$ stress parameter) ranges from $`0.05`$ to as much as 0.2 inside of $`R=100r_g`$. This is significantly higher than the values found in small local disk simulations with zero net field (e.g., Hawley, Gammie, & Balbus 1996). It is this vigorous outward angular momentum transport that is responsible for building up a Keplerian profile: an initially shallower angular momentum distribution, as is typical in a thick torus equilibrium, exports more and more of its angular momentum to the exterior, until much of the radial pressure support is eliminated. The result is a Keplerian profile.

The mass inflow ($`v_R<0`$) and outflow ($`v_R>0`$) rates through each cylindrical radius are computed as a function of time,

$$\dot{M}=R\rho v_R𝑑z𝑑\varphi .$$
(12)
The instantaneous local inflow and outflow rates are primarily due to the mass flux produced by short term $`\rho `$ and $`v_R`$ fluctuations, and they are nearly equal. This is not a matter of systematic inward mass flow being cancelled by a systematic outward mass flow. Rather, this is intrinsic to the nature of turbulent flow. The long term net accretion rate, $`\dot{M}`$, is always smaller than its rms value. An average over time for the last orbit of the simulation shows that the inflow exceeds the outflow for $`R<90r_g`$ (fig. 6). The ratio of inflow to outflow is about 1.5 between $`R=30r_g`$ and $`90r_g`$, dropping toward unity as one moves inward until $`r_{ms}`$ is reached, at which point the outflow rate plummets to zero. The time-averaged net accretion rate $`\dot{M}`$ is roughly proportional to $`R`$ between $`R=10r_g`$ and $`60r_g`$, i.e., it is not constant. The net $`\dot{M}`$ can vary with $`R`$ if either the flow is not steady, e.g., the mass in the disk is building up with time, or if there is flow out of the grid at the $`z`$ boundaries. If we consider the computational volume bounded by $`R=90r_g`$ on the outside over the last orbit of time, we find that the total mass in this volume declines and that over half of the net loss is through outflow along the $`z`$ boundaries. We conclude that much of the decline of $`\dot{M}`$ with radius can be accounted for by this outflow.

The total mass on the grid has decreased by 8.7% by the end of the simulation. Roughly 51% of this leaves through the outer radial boundary, 44% through the upper and lower $`z`$ boundaries, and the remaining 5% is accreted into the black hole. Integrating the mass flux at $`R=75r_g`$ over time shows that a net 14% of the initial torus mass has moved inward. Only about 3% of the material accreting inward from the initial torus ends up in the central black hole.

Much of the flow through the outer radial boundary is a consequence of the increasing specific angular momentum in the outer part of the torus. The flow through the $`z`$ boundaries, however, is driven from the disk by thermal and magnetic pressure. Perhaps the most interesting part of this outflow is an unbound, high temperature hollow conical outflow, confined to the axis region by surrounding magnetic pressure (see fig. 4 and fig. 7). By orbit 5, 60% of the outflow through the $`z`$ boundaries is in this conical wind. The remainder of the $`z`$ outflow is bound, although its specific energy is greater than that of the initial torus. Whether or not this more widespread outflow can be an unbound wind may therefore be a question of initial conditions: if the flow had begun as a marginally bound gas, the outflow would be unbound.

The initial average specific energy for the gas is $`2.3\times 10^3c^2`$. After 5 orbits this has decreased to $`2.5\times 10^3c^2`$ as the gas remaining on the grid becomes more bound. The total thermal, magnetic, poloidal kinetic, and orbital energies of this gas all increase with time. The magnetic energy has increased the most, receiving 42% of the total net energy increase; 95% of this is in toroidal field. The thermal, orbital and kinetic energies comprise 25%, 21%, and 12% of the increase respectively.

### 3.2 Simulation F1r: Resistive Heating

In this calculation, simulation F1 is rerun beginning at orbit 3.5 with the addition of an explicit resistivity,

$$\eta _i=0.1(\mathrm{\Delta }x_i)^2|J_i|/\sqrt{\rho }.$$
(13)
We refer to the resistive simulation as F1r. The resistivity enhances smallscale reconnection and returns the magnetic energy losses as bulk heating. A study of the effect of resistive heating is motivated in part by claims that such heating will lead to convective instabilities. Resistive heating and turbulent dissipation are physical processes hydrodynamical codes try to mimic using a Navier-Stokes viscosity. In hydrodynamic simulations, this $`\alpha `$ viscosity is argued to be the origin of convective heating and additional angular momentum transport (Narayan et al. 2000).

Full MHD simulations reveal almost no difference between the resistive F1r model and the nonresistive F1 model described above. There is no evidence of convective transport due to resistive heating. This result was already obtained by the 2D MHD simulations of SP, who compared runs with and without resistivity. Although the simulation outcomes will always depend somewhat on the specific implementation of the resistivity, grid resolution and the like, there are good reasons to believe that thermal convection per se is generally of little dynamical importance in nonradiating accretion flows (Balbus & Hawley 2002). They are dominated by MHD turbulence.

There are some minor differences between F1 and F1r that show up in the details. Immediately after switching on the resistivity, for example, the magnetic field energy in F1r declines, and there is a corresponding increase in the thermal energy. The changes are small. Integrated over the entire computational domain, the thermal energy increases by a maximum of 5% over run F1 at 4.5 orbits. As the simulation proceeds, however, the total thermal energy in F1r actually lags behind that of F1. By orbit 6, F1 has 17% more thermal energy than F1r. Model F1 builds up a hotter inner torus, in part because of a larger accretion rate into the inner region. Heating via the artificial viscosity in shocks is more important than heating by the explicit resistivity. The shocks that occur in the simulation are not large-scale, coherent, or time-independent. The nonadiabatic heating appears to arise mainly in small-scale, highly time-dependent dissipation of the turbulence.

Although resistive dissipation increases the overall thermal energy, it must do so at the expense of the magnetic field, which is driving the accretion flow in the first place. The integrated magnetic energy of F1 always exceeds that of F1r. It is 5% greater through orbit 4.5, climbs to 20% at orbit 5, and drops back to 5% by the simulation’s end. During the last two orbits, the total magnetic energy ranges from about 50% to 70% of the thermal energy. Finite resistivity decreases the turbulence levels and therefore the accretion rate as well. The result is that at late times, the internal energy in the resistive run can be smaller in isolated patches, while remaining largely unchanged in others. Figure 8 shows the central temperature as a function of radius at the end of F1 and F1r. The temperature is roughly proportional to $`R^{3/2}`$ within the accreting disk for both models.

## 4 The Generic NRAF Model

### 4.1 The Three Component Model

In this section we assemble into a simple model those features suggested by our simulation that we expect will prove to be enduring and representative of radiatively inefficient accretion flows more generally. This model has the virtue of limited flexibility, is falsifiable, and is grounded in uncontroversial dynamics. Whether it enjoys subsequent support or must ultimately be abandoned, it will have served a useful role. It is schematically illustrated in Figure 1.

The structure emerges as follows. Differentially rotating magnetized gas is subject to the MRI. The MRI produces a turbulent $`R\varphi `$ Maxwell stress, and the angular momentum transport results in accretion. (A somewhat smaller $`R\varphi `$ Reynolds stress component is also present.) The sizable outward transport of angular momentum rapidly changes the flow profile to near Keplerian, at which point the supporting pressure gradients are small. The global energy minimum state separates the matter in the disk’s inner region from the angular momentum in its outer regions, as in the classical disk evolution envisioned by Lynden-Bell & Pringle (1974). The bulk of the gas is decidedly disk-like: whether it radiates or not, the main body of the accretion flow is nearly Keplerian. Absent radiation, the disk is simply hot, vertically thickened, and radial pressure gradients are minimized.

The disk temperature increases rapidly inward, with $`TR^{3/2}`$. The vertical structure is dynamic, with gas lofting away from the equator, as pressure and centrifugal accelerations drive the gas out. This produces the hot, dynamic coronal envelope surrounding the disk. Disk-generated magnetic field also rises into the corona, where the resulting magnetic pressure significantly exceeds the gas pressure.

The disk remains vertically thick as it accretes, making an encounter with the centrifugal barrier inevitable, despite the loss of angular momentum. This is the centrifugal funnel wall. It is present because gravity weakens with increasing vertical distance from the central hole while the centrifugal force remains unchanged. Just outside the marginally stable orbit, a small hot torus of gas accumulates. The specific angular momentum in the torus is slightly greater that of the marginally stable orbit. Hot gas, pressed up against the funnel wall, accelerates up along this centrifugal barrier, and is held against it by the magnetic pressure of the surrounding corona. This is the magnetically-confined jet, but note: the magnetic confinement is from the outside medium! The jet apparently is stable.

The size of the inner torus varies depending on the jet flux, the accretion into the hole, and the rate at which matter is supplied from the Keplerian disk. It is highly variable. The final accretion into the black hole takes place only through the opening in the funnel wall at the equator, like threading a needle. A high resolution torus simulation (Hawley & Krolik 2001) found that the magnetic stress can remain large down to and beyond the marginally stable orbit. This effect might increase $`\dot{M}`$ into the hole, but the present simulation is not sufficiently well-resolved to address this.

To summarize: the combination of gravity, radiative inefficiency, angular momentum, and MHD turbulence found in black hole accretion leads to a three component flow structure—a hot Keplerian disk, an extended corona, and a jet-like central outflow.

### 4.2 Application to Sgr A\*

An important application of ADAF-type systems has been to the source Sgr A\* at the Galactic center. The properties of the black hole system at the Galactic center are reviewed by Melia & Falcke (2001). Compelling dynamical evidence suggests the presence of a massive black hole of $`2.6\times 10^6M_{}`$. Observations in X-ray and radio bands reveal a luminosity substantially below Eddington, making this system a prime candidate as an archetype low-radiative-efficiency accretion flow. Recent Chandra observations find a luminosity in 2-10 keV X-rays of $`2\times 10^{33}`$ erg s<sup>-1</sup>, and also an X-ray flare rapidly rising to a level about 45 times as large, lasting for only $`10^4`$s (Baganoff et al. 2001), indicating that the flare must originate near the black hole.

Many aspects of spectral models for Sgr A\* follow from the simple scaling laws of black hole accretion and will be present in any model, regardless of the detailed dynamics. For example, Quataert & Narayan (1999) demonstrate the impressive range of spectra that may be generated with the adjustment of a few free parameters: the ratio of electron to ion temperature, the magnetic pressure, the run of density with radius, and the accretion rate. Our knowledge of the flow is not yet sufficient to tightly constrain these parameters. As the underlying dynamical models become more sophisticated, the spectral models should become more constrained as well.

Although our simulation lacks a formal treatment of the energetics necessary for a detailed application to Sgr A\*, it is possible to look at some radiative properties of the computed flow and compare them with other, more detailed spectral predictions. The aim is to illustrate how the dynamical structures revealed in the simulation can be compared with current spectral models for Sgr A\*.

We must first translate between computational and physical units, as discussed in §2.3. The black hole mass is $`2.6\times 10^6M_{}`$, which gives a Schwarzschild radius of $`7.8\times 10^{11}`$ cm. To keep things as general as possible, let the code value $`n=1`$ be equal to a physical value $`n_o\mathrm{cm}^3`$. With this parameterization, the initial torus mass is $`6.8\times 10^{17}n_o`$ gm. The average accretion rate from the inner edge of the initial torus is roughly 2% of the torus mass per orbit (see fig. 6); in physical units this becomes $`6\times 10^{10}n_o`$ gm s<sup>-1</sup>, or $`1.6\times 10^{13}n_o\dot{M}_{Edd}`$. The average accretion rate into the central hole is about a factor of 10 smaller than this.

The accretion rate for Sgr A\* is uncertain. Coker & Melia (1997) estimate a rate of $`10^{22}`$ gm s$`{}_{}{}^{1}=0.03\dot{M}_{Edd}`$ from Bondi-Hoyle accretion of winds from nearby stars. Reconciling this accretion rate estimate with the low X-ray luminosity is a problem, however. Quataert, Narayan, & Reid (1999) argue that the low luminosity requires that the accretion rate at large radius be substantially sub-Eddington. The best fit spectral model of Melia, Liu, & Coker (2001) has an accretion rate into the central hole of $`10^{16}`$ gm s<sup>-1</sup>.

Figure 8 shows the run of temperatures along the equator in the model at the end time. This is a single fluid calculation with a simple equation of state, i.e., we assume that the electron and ion temperatures are the same. We estimate the total bremsstrahlung emissivity using an approximate form of the relativistic ion-electron bremsstrahlung formula (Svensson 1982)

$$\mathrm{\Lambda }=1.2\times 10^{22}n^2\mathrm{\Theta }\mathrm{ln}\left(\mathrm{\Theta }+1.5\right)\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^3$$
(14)
where $`\mathrm{\Theta }=kT/m_ec^2`$, and for simplicity we use $`Z=1`$ and $`n_e=n_i`$. The emission is calculated for each grid zone the total is obtained by summing up over the entire computational volume. The total bremsstrahlung luminosity on the grid at the end of the simulation is $`3\times 10^{21}n_o^2`$ erg s<sup>-1</sup>. The inner torus dominates the total emission. The total thermal energy is $`10^{37}n_o`$ ergs, so the cooling time is $`3\times 10^{16}n_o^1`$ s. The simulation time is $`10^6`$ s, so the cooling time is substantially longer than the dynamical time for $`n_o<10^{10}\mathrm{cm}^_3`$, or $`\dot{M}10^3\dot{M}_{Edd}`$. We note that the formal electron-ion equilibration time from Coulomb collisions (Spitzer 1962) is long compared to the flow time for $`n_o10^8\mathrm{cm}^3`$.

We can develop a qualitative sense of the bremsstrahlung spectrum by computing the nonrelativistic value $`n^2T^{1/2}e^{h\nu /kT}`$ for each grid zone. At the end of the simulation, the bremsstrahlung emission from the hot Keplerian disk inside of $`R=100r_g`$ peaks at a few times $`10^{19}`$Hz. The inner torus contributes the bulk of the highest frequency emission, peaking at $`10^{21}`$Hz. The coronal gas outside of one scale height from the disk emits over a broad range of frequency, but because the density is lower, the total emissivity is about a factor of 10 below that of the Keplerian disk at lower frequencies, and 2.5 orders of magnitude below at the highest energies. Our model is constrained by the observed low X-ray flux in the same way as all accretion models: the net accretion rate must be low. We are aided in achieving this by the coronal backflow, which permits the escape of matter with little additional bremsstrahlung emission, while reducing the hot gas density near the hole.

In theoretical spectral models, the radio and submillimeter emission arises in the innermost regions of the flow from synchrotron emission and Compton scattering. Synchrotron emission is governed by the electron number and energy densities, and the strength of the magnetic field. In some analyses, the field strength $`\beta `$ and electron temperature are treated as free parameters. In the model of Melia et al. (2001) $`\beta `$ comes directly from the connection between the field strength and the magnetic stress that drives the accretion. This is closer in spirit to a direct numerical simulation in which $`\beta `$ emerges self-consistently.

We shall limit our analysis of the sub-mm excess to a calculation of the spatial distribution of the peak synchrotron frequency (Rybicki & Lightman 1979)

$$\nu _c=10^{20}(n_eT^5/\beta )^{1/2}\mathrm{Hz},$$
(15)
postponing a detailed spectral analysis to a future paper. Figure 9 shows that the region of high temperature flow in the inner torus can account for the observed sub-mm emission. The highest peak frequencies are $`2.5\times 10^{11}(n_o/10^8)^{1/2}`$ Hz, and they are found at the inner edge of the disk, where gas is compressed against the centrifugal barrier. The region immediately surrounding the inner disk torus is also an emitting region. The value of $`\beta `$ is $`10`$ at the inner edge of the torus and $`1`$ in the surrounding area. A simple estimate for the synchrotron cooling time (Spitzer 1962) is $`8\times 10^5\beta (10^8/n_o)`$ s, which is longer than the flow time.

It is interesting to compare this result with the spectral model of Melia et al. (2001). Their best-fit model is an accretion flow with $`\dot{M}=10^{16}`$ gm s<sup>-1</sup>, from which they compute the emission emerging from a region inside of $`5r_g`$. The gas has a temperature $`10^{11}`$K, number densities $`10^7`$ cm<sup>-3</sup> and $`\beta 30`$. In our simulation the temperature in the inner torus is $`4\times 10^{10}`$$`10^{11}`$K, $`\beta 10`$, and the maximum number density along the equator lies between $`nn_o`$ and $`10n_o\mathrm{cm}^3`$. The inner torus is time varying over timescales of tens of hours, with fluctuations in temperature of about 50%, and vertically-averaged number density $`n`$ by about a factor of 5 (fig. 10). The dynamics are consistent with significant emission variability. This is encouraging, but given the very simple treatment of energy in our simulation, at present it is only suggestive.

It is interesting to note that key features in our simulation have been suggested independently by others. For example, the jet outflow from the inner torus may also be a source of emission, as in the jet model of Sgr A\* by Falcke & Markoff (2000). The observational evidence for coronal outflows in low luminosity black hole sources has recently emphasized out by Merloni & Fabian (2001). Combining the more sophisticated spectral treatments of these models with the dynamics observed here is an obvious next step.

## 5 Conclusions

In this paper, we have examined two global MHD simulations of a magnetized nonradiative optically thin accretion flow onto a black hole. Thermal conduction is not included, and the gas obeys a $`\gamma =5/3`$ adiabatic equation of state. The initial state is a constant angular momentum torus located at 100 Schwarzschild radii from the hole. Such flow is unstable according to the generalized adiabatic criteria of Balbus (1995), becomes fully turbulent, and in the process transports angular momentum and energy outward and greatly enhanced rates. Several distinctive properties characterize the resulting flow.

Three principal flow structures are apparent: a hot Keplerian disk that extends down to the last marginally stable orbit, a coronal envelope, and a jet along the centrifugal barrier, or funnel wall.

Although the flow is nonradiative and hot, it remains centrifugally supported. The density-weighted specific angular momentum distribution is nearly Keplerian, $`lr^{1/2}`$. There are only small sub-Keplerian departures from this. The near-Keplerian distribution arises naturally from the vigorous MRI-induced turbulence, which transports angular momentum outwards until radial pressure gradients are minimal. No circumstances have been observed in any simulations done to date where a significant non-Keplerian angular momentum distribution (e.g. constant=$`l`$) could be sustained over a large radial extent in the face of the MRI transport.

The only location where radial pressure support becomes somewhat important is in the innermost region of the flow where a hot thickened toroidal structure forms. This is constantly created and destroyed. Such a hot inner torus may be the origin of the submillimeter excess associated with Sgr A.

The entropy and temperature increase inwards. In the simulation the entropy is increased by shock heating (artificial viscosity), and resistivity. Despite this entropy gradient, convection appears to play no role in the dynamics of the flow, whose turbulent properties are determined entirely by the MRI. Marginal stability to the classical Høiland criteria is not achieved, nor could it be achieved as a matter of principle (Hawley et al. 2001). Thus, although little of the mass that accretes at large radius makes it into the hole (in contrast to an ADAF), this is not a CDAF. Mass and energy are carried off by an outflow, more in keeping with the outline of the ADIOS model.

The flow is highly dynamic, inherently multidimensional, and time-variable on all scales. The flow reacts locally to the MRI, which, because it is a fast growing instability, changes on timescales shorter than the local orbital time. Note that even if there exists a time-averaged description of the global flow that is relatively time-steady, this need not correspond to an analytic global equilibrium.

Because these features follow directly from the self-consistent MHD dynamics, we believe that they will be generic to nonradiative accretion flows. Longer evolutions would go farther in establishing that the simulation results have lost memory of the initial conditions. Also, models with a greater radial domain would allow the inflow to proceed over several decades in distance and to liberate greater amounts of gravitational energy. Higher resolution is needed in the inner region to capture the details of the disk boundary and the flow past the marginally stable orbit.

Our data sets in principle readily lend themselves to observational modeling. A very preliminary application to Sgr A\* relies on low accretion rates into the hole and emission from a hot, time-varying inner region. A detailed follow-up investigation on the radiative properties of our simulated flows is currently being pursued.

We acknowledge support under NSF grant AST-0070979, and NASA grant NAG5-9266. The simulations were carried out on the IBM Bluehorizon system of the San Diego Supercomputer Center. We thank Julian Krolik and James Stone for useful discussions related to this work.


format: text

text:
# Charmed and Bottom Baryons from Lattice NRQCD

## I INTRODUCTION

A comprehensive knowledge about the mass spectrum and spin splittings of heavy baryons is important for our understanding of quantum chromodynamics. However, except for singly heavy charmed baryons and only one singly heavy bottom baryon ($`\mathrm{\Lambda }_b`$), most of the heavy baryon masses have not yet been measured experimentally pdg . On the theoretical side there are many results on heavy baryon masses from different models including, for example, a number of quark model variationshv\_ex1 ; hv\_ex2 ; hv\_ex3 .

Using lattice QCD, substantial work has been done in the heavy meson sector. However, so far only very few results have been reported for heavy baryonshv1 ; hv2 ; hv3 ; hv4 , and there is only one work hv2  where heavy (bottom) quarks are treated nonrelativistically. A further study of charmed and bottom heavy baryons on the lattice using nonrelativistic QCD (NRQCD) therefore seems worthwhile.

Extraction of the experimentally observed mass splittings between vector and pseudoscalar mesons remains a challenging problem in lattice QCD; quenched calculations have so far been unable to extract the observed mass splittingsquen\_res , and unquenched studies have not resolved the issueunquen\_res . Therefore, it is natural to ask whether lattice results for baryon mass splittings also exhibit similar suppression compared to experiment.

Empirically, spin splittings in baryons are smaller than those in the meson sector. Moreover, in a lattice simulation the correlators for baryons, particularly for spin 3/2 states, are noisier than those for mesons, and thus, by using lattice QCD it is comparatively difficult to extract a reliable mass spectrum for heavy baryons. In this work we report on the charmed and bottom baryon mass spectrum and mass splittings by using a nonrelativistic heavy quark action and an improved light quark action on anisotropic lattices.

In section II we summarize different charmed and bottom heavy baryons with their relevant quantum numbers and discuss our choice for interpolating fields. Section III presents numerical simulation details. For heavy quarks, we use the nonrelativistic action from Ref. spwave , while a tadpole improved gauge action and an improved Dirac-Wilson action of the D234 type alford  are used for light quarks. Since these actions were previously detailed elsewhere hv3 ; hv4 , we will describe them only in an appendix. The calculations are done on two different anisotropic lattices with the same gauge configurations as were used in Ref. hv4  at $`\beta =2.1`$ and $`\beta =2.3`$.

In section IV we present our results. Masses are calculated using two methods; the first uses the standard NRQCD relation between mass and energy while the second employs mass splittings to calculate masses. As mass splittings can be estimated more accurately than masses, errors in the second method are smaller than those obtained from the first one. The overall systematic uncertainty is estimated by including scale uncertainty, uncertainty due to the choice of a time window for fitting correlation functions, error due to extrapolation to the physical light quark masses, uncertainty in fixing charm and bottom masses and uncertainty from our determination of the lattice anisotropy.

Spin splittings are discussed in section V. From our results, along with other published results, we conclude that the suppression of mass splittings is not present in the baryon sector in the same way as it is in the meson sector. Over the whole mass range where data are available, quenched lattice QCD simulations yield mass differences between spin 3/2 and spin 1/2 baryons which are comparable to or larger than experimental values.

## II CHARMED and BOTTOM BARYONS

Singly and doubly charmed and bottom baryons are summarized in Tables I and II respectively. Table II also includes doubly heavy states containing two different heavy quarks (charmed and bottom quarks together). Quark content, as well as the spin-parity $`J^P`$, the isospin $`I`$, and $`s_l`$ which identifies the total spin of the light quarks (also spin-flavor symmetry: $`s_l=0`$ is symmetric while $`s_l=1`$ is antisymmetric) are shown. Notice that masses for many singly heavy states are not measured yet and there are no data at all on masses for doubly heavy states.

To project out heavy baryon states we use the same interpolating operators as were used in Ref. hv4 . For $`\mathrm{\Sigma }`$-like baryons we choose

$$\mathrm{\Sigma }:ϵ^{abc}[q_a^TC\gamma _5Q_b]q_c$$
(1)
where $`q`$ is a light quark field and $`Q`$ is a heavy quark field. Here $`a,b,c`$ are color indices whereas Dirac indices have been suppressed. For $`\mathrm{\Sigma }_Q`$, $`q`$ is $`u`$ or $`d`$ and for $`\mathrm{\Omega }_Q`$, $`q`$ is $`s`$. For doubly heavy $`\mathrm{\Sigma }`$-like baryons with equal heavy masses, we interchange the role of light and heavy fields, i.e., to get $`\mathrm{\Xi }_{QQ}`$, we change $`qQ`$ and $`Qu`$ or $`d`$. Similarly, for $`\mathrm{\Omega }_{QQ}`$, the change is $`qQ,Qs`$.

The $`\mathrm{\Xi }_Q^{}`$ is $`\mathrm{\Sigma }`$-like but it contains two different light flavors so it is considered separately as

$`\mathrm{\Xi }^{}`$ $`:`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left\{ϵ^{abc}[q_a^TC\gamma _5Q_b]q_c+ϵ^{abc}[q_a^TC\gamma _5Q_b]q_c^{}\right\},`$ (2)
with $`q=u`$ or $`d`$ and $`q^{}=s`$.

The $`\mathrm{\Lambda }`$-like baryons involve three distinct flavors. A simple choice is the $`\mathrm{ℎ𝑒𝑎𝑣𝑦}\mathrm{𝑙𝑎𝑚𝑏𝑑𝑎}`$ :

$$\mathrm{\Lambda }:ϵ^{abc}[q_a^TC\gamma _5q_b^{}]Q_c.$$
(3)
where for $`\mathrm{\Lambda }_Q`$, $`q=u`$, $`q^{}=d`$, and for $`\mathrm{\Xi }_Q`$, $`q=u`$, $`q^{}`$ = s. A more symmetrical choice would be the octet lambda

$`\mathrm{\Lambda }_o`$ $`:`$ $`{\displaystyle \frac{1}{\sqrt{6}}}ϵ^{abc}\{2[q_a^TC\gamma _5q_b^{}]Q_c+[q_a^TC\gamma _5Q_b]q_c^{}`$ (4)
$`[q_a^TC\gamma _5Q_b]q_c\},`$
with the same flavor assignment as for the heavy lambda. One can use either of these $`\mathrm{\Lambda }`$ states as they give consistent resultshv4 . We choose the octet-lambda $`(\mathrm{\Lambda }_o)`$ for this work. For spin 3/2 baryons we choose the following interpolating field:

$$\mathrm{\Sigma }^{}:ϵ^{abc}[q_a^TC\gamma _\mu q_b^{}]Q_c,$$
(5)
where for $`\mathrm{\Sigma }_Q^{}`$, $`q=q^{}`$ is $`u`$ or $`d`$ and for $`\mathrm{\Omega }_Q^{}`$, $`q=q^{}`$ is $`s`$. To get $`\mathrm{\Xi }_Q^{}`$, one needs to consider $`q=u`$ or $`d`$ and $`q^{}=s`$. Similarly, to get the doubly heavy states with equal heavy masses one needs to interchange the role of light and heavy fields. For example, to get $`\mathrm{\Sigma }_{QQ}^{}`$, one needs to change $`q,q^{}Q`$ and $`Qu`$ or $`d`$, whereas $`\mathrm{\Omega }_Q^{}`$ requires $`q,q^{}Q`$ and $`Qs`$.

The operator in Eq. (5) has both spin 1/2 and spin 3/2 states. At zero momentum the corresponding correlation function can be written as benmer :

$`C_{ij}(t)=(\delta _{ij}{\displaystyle \frac{1}{3}}\gamma _i\gamma _j)C_{3/2}(t)+{\displaystyle \frac{1}{3}}\gamma _i\gamma _jC_{1/2}(t),`$ (6)
where $`i,j`$’s are spatial Lorentz indices and $`C_{3/2(1/2)}`$ are the spin projections for spin 3/2 (1/2) states. By choosing different Lorentz components the spin 3/2 part, $`C_{3/2}(t)`$, is extracted and used to calculate the mass of the spin 3/2 baryons.

Operators for baryons with two unlike heavy flavors may be constructed from the above interpolating operators by interchanging the role of heavy and light fields. For example, $`\mathrm{\Xi }_{QQ^{}}^{}`$ and $`\mathrm{\Omega }_{QQ^{}}^{}`$ can be obtained from Eq. (5) by letting $`q`$, $`q^{}Q,Q^{}`$ and $`Qq`$ with $`q=u`$ or $`d`$ and $`q=s`$ respectively. For $`\mathrm{\Xi }_{QQ^{}}^{}`$ and $`\mathrm{\Omega }_{QQ^{}}^{}`$ we use the symmetrical form again, as given by Eq. (3), making the same replacements i.e., $`q`$, $`q^{}Q,Q^{}`$ and $`Qq`$ with $`q=u`$ or $`d`$ and $`q=s`$ respectively. Finally, $`\mathrm{\Xi }_{QQ^{}}`$ and $`\mathrm{\Omega }_{QQ^{}}`$ are the doubly heavy analogs of $`\mathrm{\Lambda }`$ and they can be obtained from Eqs. (3) and (4) as previously.

## III NUMERICAL SIMULATION

### III.1 Actions

The gauge action as well as the heavy quark NRQCD action used for this work are described in detail in Ref. spwave . The gauge action is tadpole improved and the leading classical error is quartic in lattice spacing. The Hamiltonian corresponding to the NRQCD action is complete to $`𝒪(1/M^3)`$ in the classical continuum limit. For light quarks we use a Dirac-Wilson action of the D234 type alford  which has been used previously and detailed in Refs. hv3 ; hv4 . Its leading classical errors are cubic in lattice spacing. All these actions are summarized in the appendix.

### III.2 Simulation Details

This work is done with two sets of quenched gauge configurations (at $`\beta =2.1`$ and $`2.3`$) on anisotropic lattices with a bare aspect ratio $`a_s/a_t=2`$, where the spatial lattice spacing varies from about 0.22 to 0.15 fm.

The renormalized anisotropy is obtained from

$$\xi =\frac{a_s}{a_t}=\frac{a_sV(r_2)a_sV(r_1)}{a_tV(r_2)a_tV(r_1)},$$
(7)
where $`V(r)`$ is the potential between a static quark-antiquark pair with separation $`r`$, and is extracted from an exponential fit to a sequence of Wilson loops. In the numerator of Eq. (7) the sequence of Wilson loops extends in a coarsely-spaced direction, and in the denominator the sequence extends in the finely-spaced direction. The separation $`r`$ may be along a lattice axis or off-axis, and various possibilities were included in the calculation. However, the separation $`r`$ never includes the finely-spaced direction, so that $`V(r)`$ itself is independent of $`a_s/a_t`$. It is convenient to avoid using the largest values of $`r`$, where the exponential fit becomes noisier and the periodicity of the lattice can affect a determination of the anisotropy. Our results are

$`\xi ={\displaystyle \frac{a_s}{a_t}}=\{\begin{array}{c}1.96(2),\mathrm{for}\beta =2.1\hfill \\ 1.99(3),\mathrm{for}\beta =2.3.\hfill \end{array}`$ (10)
We used fixed time boundaries to construct quark propagators, and gauge fields were generated using a pseudo-heat-bath Monte Carlo algorithm with 400 ($`\beta =2.1`$) to 800 ($`\beta =2.3`$) sweeps between saved configurations. For $`\beta =2.1`$, we use 720 configurations and for $`\beta =2.3`$ the number of configurations is 442. Two sets of bare masses are used for each heavy quark while four sets of hopping parameters are used for the light one. Bare masses for heavy quarks are chosen to surround the physical value so that an interpolation can be used. For example, at $`\beta =2.1`$, the charm mass is in between 1.2 and 1.5 and the bottom mass is in between 5.0 and 6.0. The charm mass is fixed by setting the $`\eta _c`$ mass to its experimental value, whereas the $`B^0`$ mass is used to fix the bottom mass. The hopping parameter corresponding to the strange quark is fixed from the $`D_s`$ meson mass. The temporal lattice spacing and correspondingly the scale is fixed by setting the $`\rho `$-meson mass to its experimental value. Summaries of lattice parameters as well as hopping parameters for heavy and light fields are given in Tables III and IV respectively.

Correlation functions are calculated using interpolating operators in local form at both source and sink. In addition to that we use a gauge invariant smearing for quark propagators at the sink using the smearing function from Eq. (13) of Ref. smear . These local and sink-smeared correlators are fitted simultaneously to obtain hadron masses. The required correlations among different quantities are taken into account by covariant matrices obtained from singular value decomposition, and the statistical error is estimated from bootstrapping the fitting procedure. As in Ref. hv4 , local correlators are fitted with two exponential functions ($`A\mathrm{exp}(m_1x)+B\mathrm{exp}(m_2x)`$), while the sink-smeared correlation function is fitted with a single exponential ($`C\mathrm{exp}(m_1x)`$). The mass parameter for the sink-smeared fit is constrained to be the same as the lowest mass of the fit to the local correlator. The time window for the fit is chosen in a way such that the ending time is large and the fit is stable under variation of both starting and ending time by a few time steps.

Light quark extrapolation is done by extrapolating the hadron masses extracted at four light quark masses with the form $`c_0+c_2m_\pi ^2+c_3m_\pi ^3`$, where $`m_\pi `$ is the pion mass. In most of the cases the cubic ($`m_\pi ^3`$) contributions are small and they are included only to get systematic errors.

Figs. 1-3 show some representative examples of our simulation results. We plot the effective mass for different heavy baryons versus time $`t`$, where the effective mass is defined to be $`M(t)=ln(g(t)/g(t+1))`$ with $`g(t)`$ being the zero-momentum time correlation function of baryon fields. Open symbols in these figures are for calculations with a correlation function with local source and sink, filled symbols are for local source and smeared sink. There is good agreement between local and smeared results at large times.

It should be noted that the actual fits to determine the masses are performed directly with the correlation functions, and not on the effective masses plotted in Figs. 1-3, but the plots provide an indication of the quality of our data. Although the sink-smeared results appear to be somewhat noisy they are quite helpful in constraining the two-exponential fit of the local correlation function.

### III.3 Mass Extraction

The kinetic mass of a nonrelativistic state can be extracted from the usual NRQCD relation spwave

$`M_{kin}`$ $`=`$ $`{\displaystyle \frac{2\pi ^2}{N_s^2\xi ^2a_t[a_t(E_pE_0)]}},`$ (11)
which is derived from $`E=\frac{p^2}{2M_{kin}}`$. Here $`N=La_s`$, with $`L`$ being the lattice size and $`a_s`$ the spatial lattice spacing. $`\xi =a_t/a_s`$ is the anisotropy whereas $`E_0`$ and $`E_p`$ are simulation energies corresponding to the ground state and the state with momentum $`p=\frac{2\pi }{La_s}`$, respectively. Mass differences between two states ($`H_1`$ and $`H_2`$) with the same heavy quark can be obtained by taking the difference of their zero momentum simulation energies:

$`M_{H_1}M_{H_2}`$ $`=`$ $`E_{sim}^1(0)E_{sim}^2(0),`$ (12)
which follows from the lattice NRQCD expression for the hadron mass

$`M_H`$ $`=`$ $`E_{sim}(0)+ZM_QE_{shift},`$ (13)
where $`E_{sim}`$ is the simulation energy at zero momentum and the last two terms represent the renormalized heavy quark mass. The bare quark mass $`M_Q`$ has both a multiplicative ($`Z`$) and additive ($`E_{shift}`$) renormalizationdavthack  which should be independent of hadronic state. A more precise result is obtained for heavy hadron masses with the heaviest light quark ($`\kappa =0.229`$) rather than with a lighter light quark ($`\kappa =0.233`$ and higher).

Moreover, mass differences (Eq. 10) can be calculated more precisely than masses (Eq. 9). Therefore, for example, one can calculate a meson mass from the relation

$`M(q_l,Q)=M(q_h,H)\mathrm{\Delta }M=M(q_h,H)\mathrm{\Delta }E,`$ (14)
where

$`\mathrm{\Delta }M=\mathrm{\Delta }E=E(q_h,Q)E(q_l,Q).`$ (15)
Here $`q_h`$ and $`q_l`$ denote the heaviest light quark and a lighter one respectively, and $`M(q_h,H)`$ is extracted by using Eq. (9). Eq. (13) is valid as long as $`Z`$ in Eq. (11) is the same i.e., both states consist of the same heavy quark $`Q`$.

Similarly, masses of singly and doubly heavy baryons can be extracted from meson masses by using

$`M(q_1q_2,Q)`$ $`=`$ $`M(q_h,Q)\mathrm{\Delta }E_{sh},`$ (16)
$`M(q_1,QQ)`$ $`=`$ $`M(QQ)\mathrm{\Delta }E_{dh},`$ (17)
where

$`E_{sh}`$ $`=`$ $`E(q_h,Q)E(q_1q_2,Q),`$ (18)
$`E_{dh}`$ $`=`$ $`E(QQ)E(q_1,QQ).`$ (19)
For example, the $`\mathrm{\Sigma }_{c(b)}`$ mass is extracted by taking its difference (at each $`\kappa `$) with the $`D(B^0)`$ mass ($`m`$) at $`\kappa =0.229`$ and then subtracting that from $`m`$. Masses extracted by using Eq. (9) and Eqs. (12-17) are consistent with each other. However, errors in the second method are smaller than the previous one.

## IV RESULTS

The mass spectrum and spin splittings of heavy quark baryons have been computed on an anisotropic lattice using the NRQCD heavy quark action. Results are summarized in Table V, where the first error is the statistical error obtained from a bootstrap analysis with a bootstrap sample size equal to the configuration sample size. The second error is an overall systematic error due to scale and anisotropy uncertainties, the uncertainty due to choosing a time window, the light quark extrapolation error and the strange quark mass uncertainty. Mesons, singly heavy baryons and doubly heavy baryons are separated into different groups by horizontal lines. In Table VI we have compared our results with those obtained by using a relativistic (D234) heavy quark action hv4  and experimental numbers (where available). One can notice that the NRQCD results and D234 results are consistent with each other. Results are also consistent with a previous NRQCD calculation hv2 . As in Ref. hv4 , it is found that the suppression of spin splittings is not present in the baryon sector, although such a suppression is known to be characteristic of the heavy meson sector. One can also notice that the spin splittings for doubly heavy baryons are as large as their singly heavy counterparts.

## V DISCUSSION and SUMMARY

In order to put the results of the present calculation into perspective it is useful to consider spin splittings over the whole range of available quark masses. We start with mesons where it has been known for a long time that the squared mass difference $`M_V^2M_P^2`$ for vector and pseudoscalar mesons is approximately constant for all mesons of the form $`Q\overline{q}`$, where $`q`$ is up or down and $`Q`$ is any light or heavy flavor. This relation is illustrated in Fig. 4 for the mass pairs $`(\rho ,\pi ),`$ $`(K^{},K),`$ $`(D^{},D)`$ and $`(B^{},B).`$ Also shown are the results of lattice simulations including the present work. The tendency for quenched lattice QCD to underestimate the spin splittings relative to experimental values is clearly visible.

In Ref. hv4  we showed that it is useful to consider the behavior of the spin splittings in the baryon sector as a function of quark mass also in terms of the mesonic average mass $`(M_V+M_P)/2`$. The results for the baryon pairs $`(\mathrm{\Delta },N)`$, $`(\mathrm{\Sigma }^{},\mathrm{\Sigma })`$, $`(\mathrm{\Sigma }_c^{},\mathrm{\Sigma }_c)`$ and $`(\mathrm{\Sigma }_b^{},\mathrm{\Sigma }_b)`$ are shown in Fig. 5. It is a remarkable empirical fact that the baryon spin splitting scales almost exactly like the inverse of the average meson mass. The implication is that the ratio of meson to baryon spin splittings is almost constant. This was discussed in Ref. hv4  and was to some extent anticipated by Lipkin Lipkin  from the point of view of the quark model (see also Lipkin and O’Donnell LipO ).

The results of quenched lattice calculations are also shown in Fig. 5. The suppression of spin splittings relative to experiment, visible for mesons, is not seen for baryons. The results of the present lattice NRQCD calculation support this conclusion in the charm and bottom sectors. It is clear that a definitive measurement of $`\mathrm{\Sigma }_b`$ and $`\mathrm{\Sigma }_b^{}`$ masses would be highly desirable to extend the experimental comparison to larger mass values.

From the point of view of lattice NRQCD our results present an interesting challenge. As is well known, the spin splittings of both charmonium trot  and heavy-light mesons quen\_res ; spwave  are clearly underestimated by lattice NRQCD simulations. Up to now, these simulations have used couplings modified only by mean-field tadpole factors. It is tempting to speculate that there are additional large perturbative corrections to these couplings. In particular it might seem that a correction to the quark coupling with the chromomagnetic field (the $`c_4`$-term in Eq. (A5)) has the potential to cure the spin splitting deficiency for both charmonium and heavy-light mesons. However one has to be cautious in wishing for such a cure as it would upset the already reasonable values for spin splittings in the baryon sector.

To summarize, we have calculated the masses of baryons containing one or two heavy quarks using quenched lattice QCD. NRQCD is used to describe charm and bottom quarks. In the charm sector the results of this work are compatible with those obtained previously where a Dirac-Wilson action of the D234 type was used for the heavy quark. No suppression of the spin splittings observed in lattice NRQCD simulations of heavy-light mesons is seen in the heavy baryon sector.

This and our previous work hv4  leave a number of difficult open questions. One would like to be able to improve the lattice calculations of baryons to reduce the uncertainties to the same level achievable in mesons. Also how (and whether) the addition of dynamical quarks to the simulations will solve the dilemma of spin splittings has yet to be understood. A phenomenological issue is to understand the remarkable constancy in the meson to baryon spin splitting ratio over the whole available quark mass range. On the experimental side it will be a significant challenge to extend baryon mass measurements in the bottom and doubly heavy sectors.

###### Acknowledgements.

This work was supported in part by the Natural Sciences and Engineering Research Council of Canada. NM is also thankful to CCS of University of Kentucky for its support. Some of the computing was done on hardware funded by the Canada Foundation for Innovation, with contributions from Compaq Canada, Avnet Enterprise Solutions, and the Government of Saskatchewan.

## Appendix A Details of Actions

### A.1 NRQCD Action

The heavy quark action is nonrelativistic and is discretized to give the following Green’s function propagation:

$`G_{\tau +1}`$ $`=`$ $`\left(1{\displaystyle \frac{a_tH_B}{2}}\right)\left(1{\displaystyle \frac{a_tH_A}{2n}}\right)^n`$ (20)
$`{\displaystyle \frac{U_4^{}}{u_t}}\left(1{\displaystyle \frac{a_tH_A}{2n}}\right)^n\left(1{\displaystyle \frac{a_tH_B}{2}}\right)G_\tau ,`$
The nonrelativistic Hamiltonian is complete to $`O(1/M^3)`$ in the classical continuum limit :

$`H`$ $`=`$ $`H_0+\delta H,`$ (21)
$`H_0`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }^{(2)}}{2M}},`$ (22)
$`\delta H`$ $`=`$ $`\delta H^{(1)}+\delta H^{(2)}+\delta H^{(3)}+O(1/M^4)`$ (23)
$`\delta H^{(1)}`$ $`=`$ $`{\displaystyle \frac{c_4}{u_s^4}}{\displaystyle \frac{g}{2M}}𝝈\stackrel{~}{𝐁}+c_5{\displaystyle \frac{a_s^2\mathrm{\Delta }^{(4)}}{24M}},`$ (24)
$`\delta H^{(2)}`$ $`=`$ $`{\displaystyle \frac{c_2}{u_s^2u_t^2}}{\displaystyle \frac{ig}{8M^2}}(\stackrel{~}{𝚫}\stackrel{~}{𝐄}\stackrel{~}{𝐄}\stackrel{~}{𝚫})`$ (25)
$`{\displaystyle \frac{c_3}{u_s^2u_t^2}}{\displaystyle \frac{g}{8M^2}}𝝈(\stackrel{~}{𝚫}\times \stackrel{~}{𝐄}\stackrel{~}{𝐄}\times \stackrel{~}{𝚫})`$
$`c_6{\displaystyle \frac{a_s(\mathrm{\Delta }^{(2)})^2}{16n\xi M^2}},`$
$`\delta H^{(3)}`$ $`=`$ $`c_1{\displaystyle \frac{(\mathrm{\Delta }^{(2)})^2}{8M^3}}{\displaystyle \frac{c_7}{u_s^4}}{\displaystyle \frac{g}{8M^3}}\{\stackrel{~}{\mathrm{\Delta }}^{(2)},𝝈\stackrel{~}{𝐁}\}`$ (26)
$`{\displaystyle \frac{c_9ig^2}{8M^3}}𝝈\left({\displaystyle \frac{\stackrel{~}{𝐄}\times \stackrel{~}{𝐄}}{u_s^4u_t^4}}+{\displaystyle \frac{\stackrel{~}{𝐁}\times \stackrel{~}{𝐁}}{u_s^8}}\right)`$
$`{\displaystyle \frac{c_{10}g^2}{8M^3}}\left({\displaystyle \frac{\stackrel{~}{𝐄}^2}{u_s^4u_t^4}}+{\displaystyle \frac{\stackrel{~}{𝐁}^2}{u_s^8}}\right)c_{11}{\displaystyle \frac{a_s^2(\mathrm{\Delta }^{(2)})^3}{192n^2\xi ^2M^3}}.`$
Here a tilde signifies discretization errors have been removed. In particular,

$`\stackrel{~}{E}_i`$ $`=`$ $`\stackrel{~}{F}_{4i},`$ (27)
$`\stackrel{~}{B}_i`$ $`=`$ $`{\displaystyle \frac{1}{2}}ϵ_{ijk}\stackrel{~}{F}_{jk},`$ (28)
$`\stackrel{~}{F}_{\mu \nu }(x)`$ $`=`$ $`{\displaystyle \frac{5}{6}}F_{\mu \nu }(x){\displaystyle \frac{1}{6u_\mu ^2}}U_\mu (x)F_{\mu \nu }(x+\widehat{\mu })U_\mu ^{}(x)`$ (29)
$`{\displaystyle \frac{1}{6u_\mu ^2}}U_\mu ^{}(x\widehat{\mu })F_{\mu \nu }(x\widehat{\mu })U_\mu (x\widehat{\mu })(\mu \nu ).`$
The various spatial lattice derivatives are defined as follows:

$`a_s\mathrm{\Delta }_iG(x)`$ $`=`$ $`{\displaystyle \frac{1}{2u_s}}[U_i(x)G(x+\widehat{ı})`$ (30)
$`U_i^{}(x\widehat{ı})G(x\widehat{ı})],`$
$`a_s\mathrm{\Delta }_i^{(+)}G(x)`$ $`=`$ $`{\displaystyle \frac{U_i(x)}{u_s}}G(x+\widehat{ı})G(x),`$ (31)
$`a_s\mathrm{\Delta }_i^{()}G(x)`$ $`=`$ $`G(x){\displaystyle \frac{U_i^{}(x\widehat{ı})}{u_s}}G(x\widehat{ı}),`$ (32)
$`a_s^2\mathrm{\Delta }_i^{(2)}G(x)`$ $`=`$ $`{\displaystyle \frac{U_i(x)}{u_s}}G(x+\widehat{ı})2G(x)`$ (33)
$`+{\displaystyle \frac{U_i^{}(x\widehat{ı})}{u_s}}G(x\widehat{ı}),`$
$`\stackrel{~}{\mathrm{\Delta }}_i`$ $`=`$ $`\mathrm{\Delta }_i{\displaystyle \frac{a_s^2}{6}}\mathrm{\Delta }_i^{(+)}\mathrm{\Delta }_i\mathrm{\Delta }_i^{()},`$ (34)
$`\mathrm{\Delta }^{(2)}`$ $`=`$ $`{\displaystyle \underset{i}{}}\mathrm{\Delta }_i^{(2)},`$ (35)
$`\stackrel{~}{\mathrm{\Delta }}^{(2)}`$ $`=`$ $`\mathrm{\Delta }^{(2)}{\displaystyle \frac{a_s^2}{12}}\mathrm{\Delta }^{(4)},`$ (36)
$`\mathrm{\Delta }^{(4)}`$ $`=`$ $`{\displaystyle \underset{i}{}}\left(\mathrm{\Delta }_i^{(2)}\right)^2.`$ (37)
### A.2 Gauge Field Action

The leading classical errors of the gauge field action are quartic in lattice spacing. The action is

$`S_G(U)`$ $`=`$ $`{\displaystyle \frac{5\beta }{3}}[{\displaystyle \frac{1}{u_s^4\xi }}{\displaystyle \underset{\mathrm{ps}}{}}(1{\displaystyle \frac{1}{3}}\mathrm{ReTr}U_{\mathrm{ps}})`$ (38)
$`{\displaystyle \frac{1}{20u_s^6\xi }}{\displaystyle \underset{\mathrm{rs}}{}}\left(1{\displaystyle \frac{1}{3}}\mathrm{ReTr}U_{\mathrm{rs}}\right)`$
$`+{\displaystyle \frac{\xi }{u_s^2u_t^2}}{\displaystyle \underset{\mathrm{pt}}{}}\left(1{\displaystyle \frac{1}{3}}\mathrm{ReTr}U_{\mathrm{pt}}\right)`$
$`{\displaystyle \frac{\xi }{20u_s^4u_t^2}}{\displaystyle \underset{\mathrm{rst}}{}}\left(1{\displaystyle \frac{1}{3}}\mathrm{ReTr}U_{\mathrm{rst}}\right)`$
$`{\displaystyle \frac{\xi }{20u_s^2u_t^4}}{\displaystyle \underset{\mathrm{rts}}{}}(1{\displaystyle \frac{1}{3}}\mathrm{ReTr}U_{\mathrm{rts}})],`$
where anisotropic ratio $`\xi a_s/a_t`$ and $`\beta `$ is the lattice gauge field coupling constant.
ps : spatial plaquettes
rs : spatial planar 1$`\times `$2 rectangles,
pt : plaquettes in the temporal-spatial plane,
rst(rts) : rectangles with the long side in a spatial(temporal) direction.

### A.3 Light Quark Action

For light quarks, we used a D234 action alford ; hv4  with parameters set to their tadpole-improved classical values. Its leading classical errors are cubic in lattice spacing and the action can be written as

$`S_F(\overline{q},q;U)`$ $`=`$ $`{\displaystyle \frac{4\kappa }{3}}{\displaystyle \underset{x,i}{}}\left[{\displaystyle \frac{1}{u_s\xi ^2}}D_{1i}(x){\displaystyle \frac{1}{8u_s^2\xi ^2}}D_{2i}(x)\right]`$ (39)
$`+{\displaystyle \frac{4\kappa }{3}}{\displaystyle \underset{x}{}}\left[{\displaystyle \frac{1}{u_t}}D_{1t}(x){\displaystyle \frac{1}{8u_t^2}}D_{2t}(x)\right]`$
$`+{\displaystyle \frac{2\kappa }{3u_s^4\xi ^2}}{\displaystyle \underset{x,i<j}{}}\overline{\psi }(x)\sigma _{ij}F_{ij}(x)\psi (x)`$
$`+{\displaystyle \frac{2\kappa }{3u_s^2u_t^2\xi }}{\displaystyle \underset{x,i}{}}\overline{\psi }(x)\sigma _{0i}F_{0i}(x)\psi (x){\displaystyle \underset{x}{}}\overline{\psi }(x)\psi (x),`$
where

$`D_{1i}(x)`$ $`=`$ $`\overline{\psi }(x)(1\xi \gamma _i)U_i(x)\psi (x+\widehat{i})`$ (40)
$`+\overline{\psi }(x+\widehat{i})(1+\xi \gamma _i)U_i^{}(x)\psi (x),`$
$`D_{1t}(x)`$ $`=`$ $`\overline{\psi }(x)(1\gamma _4)U_4(x)\psi (x+\widehat{t})`$ (41)
$`+\overline{\psi }(x+\widehat{t})(1+\gamma _4)U_4^{}(x)\psi (x),`$
$`D_{2i}(x)`$ $`=`$ $`\overline{\psi }(x)(1\xi \gamma _i)U_i(x)U_i(x+\widehat{i})\psi (x+2\widehat{i})`$ (42)
$`+\overline{\psi }(x+2\widehat{i})(1+\xi \gamma _i)U_i^{}(x+\widehat{i})U_i^{}(x)\psi (x),`$
$`D_{2t}(x)`$ $`=`$ $`\overline{\psi }(x)(1\gamma _4)U_4(x)U_4(x+\widehat{t})\psi (x+2\widehat{t})`$ (43)
$`+\overline{\psi }(x+2\widehat{t})(1+\gamma _4)U_4^{}(x+\widehat{t})U_4^{}(x)\psi (x),`$
$`gF_{\mu \nu }(x)`$ $`=`$ $`{\displaystyle \frac{1}{2i}}\left(\mathrm{\Omega }_{\mu \nu }(x)\mathrm{\Omega }_{\mu \nu }^{}(x)\right)`$ (44)
$`{\displaystyle \frac{1}{3}}\mathrm{Im}\left(\mathrm{Tr}\mathrm{\Omega }_{\mu \nu }(x)\right),`$
$`\mathrm{\Omega }_{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{4}}[U_\mu (x)U_\nu (x+\widehat{\mu })U_\mu ^{}(x+\widehat{\nu })U_\nu ^{}(x)`$ (45)
$`+U_\nu (x)U_\mu ^{}(x\widehat{\mu }+\widehat{\nu })U_\nu ^{}(x\widehat{\mu })U_\mu (x\widehat{\mu })`$
$`+U_\mu ^{}(x\widehat{\mu })U_\nu ^{}(x\widehat{\mu }\widehat{\nu })U_\mu (x\widehat{\mu }\widehat{\nu })U_\nu (x\widehat{\nu })`$
$`+U_\nu ^{}(x\widehat{\nu })U_\mu (x\widehat{\nu })U_\nu (x+\widehat{\mu }\widehat{\nu })U_\mu ^{}(x)].`$