advertisement

Draft version February 1, 2008 Preprint typeset using LATEX style emulateapj v. 11/12/01 NUMERICAL MODELS OF VISCOUS ACCRETION FLOWS NEAR BLACK HOLES Jonathan C. McKinney1 and Charles F. Gammie1,2 arXiv:astro-ph/0204045v1 3 Apr 2002 jcmcknny@uiuc.edu, gammie@uiuc.edu Draft version February 1, 2008 ABSTRACT We report on a numerical study of viscous fluid accretion onto a black hole. The flow is axisymmetric and uses a pseudo-Newtonian potential to model relativistic effects near the event horizon. The numerical method is a variant of the ZEUS code. As a test of our numerical scheme, we are able to reproduce results from earlier, similar work by Igumenshchev and Abramowicz and Stone et al. We consider models in which mass is injected onto the grid as well as models in which an initial equilibrium torus is accreted. In each model we measure three “eigenvalues” of the flow: the accretion rate of mass, angular momentum, and energy. We find that the eigenvalues are sensitive to rin , the location of the inner radial boundary. Only when the flow is always supersonic on the inner boundary are the eigenvalues insensitive to small changes in rin . We also report on the sensitivity of the results to other numerical parameters. Subject headings: accretion disks, black hole physics, hydrodynamics, turbulence, galaxies: active The α model artfully avoids the question of the origin and nature of turbulence in accretion disks. This allows useful estimates to be made absent the solution to a difficult, perhaps intractable, problem. Recently, however, significant progress has been made in understanding the origin of turbulence in accretion flows. It is now known that, in the magnetohydrodynamic (MHD) approximation, an accreting, differentially rotating plasma is destabilized by a weak magnetic field (Balbus & Hawley 1991; Hawley & Balbus 1991). This magneto-rotational instability (MRI) generates angular momentum transport under a broad range of conditions. Numerical work has shown that in a plasma that is fully ionized, which is likely the case for the inner regions of most black hole accretion flows, the MRI is capable of sustaining turbulence in the nonlinear regime (Hawley & Balbus 1991; Hawley et al. 1995; Hawley 2000; Hawley & Krolik 2001). Studies of unmagnetized disks have greatly reduced the probability that a linear or nonlinear hydrodynamic instability drives disk turbulence. While there are known global hydrodynamic instabilities that could in principle initiate turbulence, these have turned out to saturate at low levels or require conditions that are not relevant to an accretion disk near a black hole. As of this writing, no local, linear or nonlinear hydrodynamic instabilities that transport angular momentum outwards are known to exist in Keplerian disks (Balbus & Hawley 1998). Work on magnetized disks has now turned to global numerical models. These are possible thanks to advances in computer hardware and algorithms. Recent work by Hawley (2000, 2001), Stone & Pringle (2001), and Hawley & Krolik (2001) considers the evolution of inviscid, nonrelativistic MHD accretion flows in two or three dimensions. Some of this work uses a pseudo-Newtonian, or Paczynski & Wiita (1980), potential as a model for the effects of strong-field gravity near the event horizon. Other work on global models has considered the equations of viscous, compressible fluid dynamics as a model for 1. introduction Black hole accretion flows are the most likely central engine for quasars and active galactic nuclei (AGN) (Zel’dovich 1964; Salpeter 1964). As such they are the subject of intense astrophysical interest and speculation. Recent observations from XMM-Newton, Chandra, Hubble, VLBA, and other ground- and space-based observatories have expanded our understanding of the time variability, spectra, and spatial structure of AGN. Radio interferometry, in particular, has been able to probe within a few hundred gravitational radii (GM/c2 ) of the central black hole, e.g. Lo et al. (1998); Junor et al. (1999); Doeleman et al. (2001). Despite these observational advances, only instruments now in the concept phase will have sufficient angular resolution to spatially resolve the inner accretion disk (Rees 2001). And so there remain fundamental questions that we can only answer by folding observations through models of AGN structure. All black hole accretion flow models require that angular momentum be removed from the flow in some way so that material can flow inwards. In one group of models, angular momentum is removed directly from the inflow by, e.g., a magneto-centrifugal wind (Blandford & Payne 1982). Here we will focus on the other group of models in which angular momentum is diffused outward through the accretion flow. It has long been suspected that the diffusion of angular momentum through an accretion flow is driven by turbulence. The α model (Shakura & Sunyaev 1973) introduced a phenomenological shear stress into the equations of motion to model the effects of this turbulence. This shear stress is proportional to αP , where α is a dimensionless constant and P is the (gas or gas + radiation) pressure. This shear stress permits an exchange of angular momentum between neighboring, differentially rotating layers in an accretion disk. In this sense it is analogous to a viscosity (see also Lynden-Bell & Pringle (1974)) and is often referred to as the “anomalous viscosity.” 1 Physics Department and Center for Theoretical Astrophysics, University of Illinois at Urbana-Champaign, Loomis Laboratory, 1110 W.Green St. Urbana,IL 61801 2 Department of Astronomy and National Center for Supercomputing Applications (NCSA) 1 2 McKinney & Gammie the accreting plasma (Igumenshchev & Abramowicz 1999; Stone et al. 1999; Igumenshchev & Abramowicz 2000; Igumenshchev et al. 2000). The viscosity is meant to mock up the effect of small scale turbulence, presumably generated by magnetic fields, on the large scale flow. In light of work on numerical MHD models, this may seem like a step backwards. The MHD models, however, are computationally expensive and introduce new problems with respect to initial and boundary conditions. It therefore seems reasonable to investigate the less expensive α based viscosity models. In this paper we investigate axisymmetric, numerical, viscous inflow models. This work was motivated by the earlier work of Igumenshchev & Abramowicz (1999, 2000) and Stone et al. (1999), hereafter referred to as IA99, IA00, and SPB99, respectively (IA99 and IA00 are collectively referred to as IA, in which case SPB99 is simply referred to as SPB). These authors studied similar viscous inflow models yet found different radial scaling laws for radial velocity, density, and angular momentum. They also found different values for the accretion rate of mass and angular momentum. They used different experimental designs, however. We set out to discover whether the results from these authors differed due to numerical methods or model parameters. Along the way, we took a systematic approach to studying numerical parameters and boundary conditions. One particular point of concern, which will be described in greater detail below, is the inner boundary condition. This lies in the energetically dominant portion of the flow, so errors there can potentially corrupt the entire model. In this paper we show that aspects of results presented by other researchers are sensitive to model and numerical parameters. These results should be useful to others contemplating large-scale numerical models of black hole accretion. The paper is organized as follows. In §2 we discuss our models. In §3 we discuss numerical methods. In §4 we discuss a fiducial solution and results from a survey of other solutions. In §5 we summarize our results. 2. model We are interested in modeling the plasma within a few hundred GM/c2 of a black hole. We will consider only axisymmetric models (the work of Igumenshchev et al. (2000) suggests that 2D and 3D viscous models give similar results). Throughout we use standard spherical polar coordinates r, θ, and φ. We solve numerically the axisymmetric, nonrelativistic equations of compressible hydrodynamics in the presence of an anomalous stress Π, which is meant to model the effects of small-scale turbulence on the mean flow. The governing equations then express the conservation of mass Dρ + ρ(∇ · v) = 0, (1) Dt momentum, Dv = −∇P − ρ∇Ψ − ∇ · Π, (2) ρ Dt and energy, Du = −(P + u)(∇ · v) + Φ. (3) Dt Here, as usual, D/Dt ≡ ∂/∂t + v · ∇ is the Lagrangian time derivative, ρ is the mass density, v is the velocity, u is the internal energy density, P is the pressure, and Ψ is the gravitational potential. The dissipation function Φ is given by the product of the anomalous stress tensor Π with the rate-of-strain tensor e (given explicitly in the Appendix) Φ = Πij eij , (4) (sum over indices) where the anomalous stress tensor is the term-by-term product Πij = −2ρνeij Sij , (5) (no sum over indices) where Sij is a symmetric matrix filled with 0 or 1 that serves as a switch for each component of the anomalous stress. The equation of state is P = (γ − 1)u. (6) For the gravitational potential we use the pseudoNewtonian potential of Paczynski & Wiita (1980): Ψ = −GM/(r − rg ) (here rg ≡ 2GM/c2 ). This potential reproduces features of the orbital structure of a Schwarzschild spacetime, including an innermost stable circular orbit (ISCO) located at r = 6GM/c2 . In a few cases we use the Newtonian potential Ψ = −GM/r for comparison with others’ work (IA and SPB exclusively use a Newtonian potential). We must now make some choices for the anomalous stress tensor. One might argue on very general grounds for a Navier-Stokes prescription, and indeed IA and we use a prescription where all elements of S are 1 (the “IA prescription”). SPB, on the other hand, use the NavierStokes prescription with all components zero except Srφ , Sθφ , Sφr , and Sφθ (the “SPB prescription”). SPB justify this choice by arguing that it more appropriately models MHD turbulence; this was later supported by results presented in Stone & Pringle (2001). We must also choose a viscosity coefficient. We consider three different prescriptions: one similar to those chosen by IA; a second viscosity coefficient similar to that chosen by SPB; and a third, similar form that vanishes rapidly near the poles. Explicitly, ν = α(c2s /ΩK ), (7) ν = α(ρ/ρ0 )Ω0 r02 , (8) ν = α(c2s /ΩK )sin3/2 (θ), (9) are the IA,pSPB, and MG prescriptions, respectively, where cs = γP/ρ is the sound speed, and Ω2K ≡ GM 1 ∂Ψ = r ∂r r(r − rg )2 (10) is the “Keplerian” angular velocity. Here ρ0 and Ω0 are values of ρ and Ω at a fiducial radius r0 . The choice of viscosity coefficient for IA and MG is based, as usual, on dimensional arguments. SPB99’s choice focuses the viscosity where most of the matter is, a numerical convenience. Our MG prescription is a small modification of the IA prescription to concentrate the viscosity toward the equator. These choices are to a large extent arbitrary, although one might attempt to motivate the choice by comparison with MHD simulations, as do Stone & Pringle (2001). Nevertheless, some dynamical Numerical Viscous Accretion Flows properties of MHD turbulence, such as the elastic properties that produce magnetic tension and hence Alfven waves, can never be modeled with a viscosity. The model has boundaries at θ = ±π/2 and at r = rin , rout . At the θ boundaries we use the usual polar axis boundary conditions. At the radial boundaries we use “outflow” boundary conditions; ideally these boundary conditions should be completely transparent to outgoing waves. Fuel for the accretion flow must be provided either in the initial conditions or continuously over the model evolution. Some global numerical accretion flow models have started with an equilibrium torus. Examples include Hawley (1991), Hawley et al. (1995), Hawley (2000), Hawley (2001), Hawley et al. (2001), and Hawley & Krolik (2001). Others have started with an initial configuration of matter that is not in equilibrium. For example, matter may be placed in orbit about the black hole, but with sub-Keplerian angular momentum, so that once the simulation commences it immediately falls toward the hole. Examples of this approach include Hawley et al. (1984), Koide et al. (1999), Koide et al. (2000), and Meier et al. (2001). This approach may enhance transients associated with the choice of initial conditions, although it can also be physically well motivated, as in studies of core-collapse supernovae. One can also inject fluid continuously onto the computational grid over the course of the evolution. Examples of this include IA99 and IA00. One might also use an inflow boundary condition at the outer radial edge of the grid, as in Blondin et al. (2001). The main advantage of injection models is that they allow one to achieve a steady, or statistically steady, state. In this paper we will consider only the equilibrium tori and on-grid injection models. The equilibrium tori (Papaloizou & Pringle 1984) are steady-state solutions to the equations of inviscid hydrodynamics. They assume a polytropic distribution of mass and internal energy, P = Kργ , and a power-law rotation profile, Ω ∝ (r sin θ)−q . The radial and meridional components of the velocity vanish. There are 5 parameters that describe the torus: (1) the location of the torus pressure maximum r0 ; (2) the location of the innermost edge of the torus, rt,in < r0 ; (3) the maximum value of the density, ρ0 = ρ(r0 ); (4) the angular velocity gradient q = d ln Ω/d ln R; and (5) the entropy constant K. On-grid injection adds matter to the model at a constant rate. The matter is injected with a non-zero specific angular momentum and with zero radial or meridional momentum in a steady pattern ρ̇(r, θ) which is typically symmetric about the equator. Parameters for this scheme include: (1) a characteristic radius for injection rinj ; (2) the rate of mass injection Ṁinj ; (3) the specific angular momentum of the injected fluid, vφ = f1 rΩK . We usually set f1 = 0.95, so that the fluid circularizes near rinj . This restricts transients associated with circularization to the outer portions of the computational domain; (4) the internal energy of the injected fluid, u = f2 ρΨ. We always set f2 = 0.2 so that the fluid is marginally bound, i.e. has Bernoulli parameter Be ≡ (1/2)v 2 + c2s /(γ − 1) + Ψ < 0. One must also choose the injection pattern ρ̇(r, θ). IA99 choose a radially narrow region, but do not explicitly give ρ̇(r, θ). We use a Gaussian, but found that none of the re- 3 sults are sensitive to the precise profile. The accretion rate of mass, energy, and angular momentum are insensitive to large changes in the size of the injection region except for the extreme cases of filling the entire θ width or injecting in 2 locations. These extreme cases are sufficiently different to be referred to as a completely different model; they lead to a qualitative change in the flow. For example, a full range θ injection region has matter that will collide with any outflow at the poles. A bipolar injection leads to an equatorial outflow. Our models have radial width σr = 0.05(rin + rout )/2 and σθ = π/8. 3. numerical methods Our numerical method is based on ZEUS-2D (Stone & Norman 1992) with the addition of an explicit scheme for the viscosity. ZEUS is an operator-split, finite-difference algorithm on a staggered mesh that uses an artificial viscosity to capture shocks (in addition to the anomalous viscosity in equations [2]). This algorithm guarantees that momentum and mass are conserved to machine precision. Total energy is conserved only to truncation error, so total energy conservation is useful in assessing the accuracy of the evolution. The inner and outer radial boundary conditions are implemented by copying primitive variable values (ρ, u, and v) from the last zone on the grid into a set of “ghost zones” immediately outside the grid. Inflow from outside the grid is forbidden; we set vr (rin ) = 0 if vr (rin ) > 0 and vr (rout ) = 0 if vr (rout ) < 0. Since we expect inflow on the inner boundary, this switch should seldom be activated. We have found that frequent activation of the switch is usually an indication of a numerical problem. We use a radial grid uniform in log(r − rg ). We require that dr(r)/(rin − rg ) ≤ 1/4 so that the structure of the pseudo-Newtonian potential is well resolved. The grid is uniform in θ. The grid has Nr × Nθ zones. 3.1. Numerical Treatment of Low Density Regions Like many schemes for numerical hydrodynamics, ZEUS can tolerate only a limited dynamic range in density. It is therefore necessary to impose a density minimum ρf l to avoid small or negative densities. Our procedure for imposing the floor is equivalent to adding a small amount of mass to the grid every time the floor is invoked. Mass is added in such a way that momentum is conserved. To monitor the effect of the density floor, we track the rate of change of total mass and total energy (from kinetic energy change) due to this procedure, Ṁf l and Ėf l . We set ρf l = 10−10 Ṁinj c3 /(GM )2 for injection runs and ρf l = 10−5 ρ0 for torus runs. Lower values for ρf l do not lead to a significant change in the solution. Larger values of ρf l give Ṁf l ∼ Ṁ , the accretion rate through the inner boundary. The atmosphere also becomes more massive and begins to affect torus stability– vertical oscillations are excited in the inner disk by a Kelvin-Helmholtz like instability. This should be avoided. We must also surround the torus in a low density atmosphere in the initial conditions. The density of the atmosphere is ρf l and the internal energy density is u = Uo ρΨ, where Uo is a constant fraction of order unity (e.g. IA and we choose Uo = 0.2). The addition of the atmosphere has no effect on the solution since the mass source’s evolution 4 McKinney & Gammie eventually dominates the flow everywhere. SPB99 choose a different method of constructing the initial atmosphere but obtain late-time results that are similar to ours. It is also necessary to impose a floor uf l on the internal energy density. This we take to be the minimum value of u in the initial atmosphere. As for the mass, we track the rate of change of total energy due to the internal energy floor, that along with the kinetic energy is included in Ėf l . 3.2. Diagnostics Global numerical simulations of accretion flows are complicated; it is possible to measure many quantities associated with the flow. Some are astrophysically relevant, and some are not. In our view particular interest attaches to the time-averaged flux of mass, energy, and angular momentum through the inner boundary. Physically, these are directly related to the luminosity of the accretion flow and the rate of change of mass and angular momentum of the central black hole. As described by Narayan & Popham (1993), these are in a sense the nonlinear “eigenvalues” of the model. The mass accretion rate is Z Ṁ = ρv · dS, (11) S where S is the inner radial surface of the computational domain. The total energy accretion rate is Z 1 Ė = (( v 2 + h + Ψ)ρv + Π · v) · dS, (12) 2 S where h = (u + P )/ρ = γu/ρ is the specific enthalpy with our equation of state. The angular momentum accretion rate is Z (13) L̇ = r sin θ(ρvvφ + Π · φ̂) · dS. S It is also sometimes useful to focus on the reduced eigenvalues l = L̇/Ṁ and e = Ė/Ṁ . These value of mass, energy, and angular momentum are recorded at about 2 grid zones away from Rin . This avoids any error that may occur when evaluating directly on the boundary where the inflow boundary condition is applied. We also track volume-integrated quantities, the flux of mass, energy, and angular momentum across all boundaries, and floor added quantities in order to evaluate the consistency of the results. Mass and angular momentum are conserved to machine precision, although “machine precision” implies a surprisingly large random walk in the integrated quantities over the full integration because the calculation requires millions of timesteps. Total energy is conserved to truncation error, not machine precision, and thus is a useful check on the quality of the simulation. Total energy conservation implies Ėerr = Ėvol + Ė + Ėout − Ėf l , (14) where Ėvol is the rate of change of the volume integrated total energy, Ėout is the flux of total energy through the outer radial boundary, and Ėf l is the rate of total energy added due to the kinetic energy change (because of the mass density floor) and internal energy density floor. Ideally, Ėerr = 0. Truncation errors can (and do) lead to cumulative, rather than random, changes in the total energy. A useful gauge of the magnitude of these errors is Ėerr /E. For all runs we performed the error rate is within 10% of 10−5 c3 /GM for a torus run and within 10% of 10−4 c3 /GM for a viscous injection run. 3.3. Code Tests Our version of ZEUS reproduces all hydrodynamic test results from Stone & Norman (1992), including their spherical advection and Sod shock tests. We also find excellent agreement with steady spherical accretion solutions, i.e. Bondi flow /citepbondi52. An inviscid equilibrium torus run also persists for many dynamical times with insignificant deviations from the initial conditions. We have parallelized our code using the MPI message passing library. On the Origin 2000 at NCSA we are able to achieve about 2.5 × 107 zone updates per second using 240 CPUs, or about 35 GFLOPs. This is 159 times faster than the single CPU speed, which represents a parallel efficiency of 66%. 4. results The initial motivation for undertaking this calculation was to understand differences between results reported in IA and SPB. Using our code, which is based on the same algorithm used in SPB’s calculations, we ran a series of tests attempting to reproduce SPB’s results. These test calculations used all of SPB’s model choices, including SPB’s viscosity prescription, a Newtonian potential, and a torus for the mass source. We were able to reproduce most quantitative results reported in SPB99’s torus calculations. This includes their radial scaling laws. For example, in a model that is identical to SPB99’s Run B, we find Ṁ ∝ r, ρ ∝ r0 , c2s ∝ r−1 , vφ ∝ r−1/2 , and |vr | ∝ r−1 . These power law slopes are identical to those reported by SPB99. As another example, we found Ṁ = 1.23 × 10−3 torus masses per torus orbit at the pressure maximum for a model identical to SPB’s Model A (their fiducial model); SPB report Ṁ = 1.0 × 10−3 in the same units. Given the fluctuations in mass accretion rate, our value and SPB’s value are fully consistent. We were even able to reproduce certain numerical artifacts associated with the inner radial boundary, such as a density drop and temperature spike near the inner boundary. Recall that SPB evolve an initial torus and allow it to accrete; IA use a different experimental design in which matter is steadily injected onto the grid. They also use a different viscosity prescription. We ran a second series of test calculations attempting to reproduce IA’s results. These test calculations used all of IA’s model choices, including viscosity prescription, Newtonian potential, etc. We were able to reproduce all of IA99’s calculations except those that include thermal conduction (which we did not attempt to reproduce). In each case we found that the qualitative nature of the flow is similar to that described in IA99. In particular, we agree on which models are stable and unstable and which models exhibit outflows. We also find qualitative agreement with their contour plots of, e.g., density pressure, mass flux, and Mach number. We also find qualitative agreement with their radial run of cs /VK and specific angular momentum. Our results do not agree precisely, but this is likely due to small differences in mass Numerical Viscous Accretion Flows injection scheme (because IA99 do not give their ρ̇(r, θ)). Finally, we can also reproduce the radial scalings given in IA00 for their model A. The most significant difference between the results of IA and SPB was due to the choice of anomalous stress prescription, as might have been anticipated. Qualitatively, the stress components that are included in IA and not SPB tend to smooth the flow and suppress turbulence. Thus SPB99’s simulations result in more vigorous convection than simulations performed by IA. The choice of mass supply (torus vs. injection) also leads to a significant difference between IA and SPB’s results; this is discussed in more detail below. The fact that we can reproduce both SPB’s and IA’s results using a single code is consistent with the hypothesis that differences between their reported results (e.g. the lower degree of convection reported in IA than SPB, and the differences in radial scaling laws) is due to differences in experimental design and viscosity prescription rather than numerical methods. While we cannot completely rule out the possibility that SPB, IA, and we have made identical experimental errors, this seems unlikely. This comparison thus lends credibility to SPB, IA, and our numerical results. 4.1. Fiducial Model Evolution We now turn from reproducing earlier viscosity models to considering new aspects of our own models. First, consider the evolution of a “fiducial” model (Run A in Table A1 and Table A2). The fiducial model has rin = 2.7GM/c2 , rout = 600GM/c2, rinj = 495GM/c2, γ = 3/2, α = 0.1, Nr = 108, Nθ = 50. It uses a pseudo-Newtonian potential, mass is supplied by injection, and the viscosity prescription follows IA. It was run from t = 0 to t = 7.3 × 105 GM/c3 . Run A is similar to IA99’s “Model 5”, except that it uses a pseudo-Newtonian potential. In a statistically steady state the flow is characterized by a quasi-periodic outflow. Hot bubbles form at the interface between bound (Bernoulli parameter Be = (1/2)v 2 + c2s /(γ − 1) + Ψ < 0) and unbound (Be > 0) material. These hot bubbles are buoyant and move away from the black hole. This appears to be a low-frequency, low wavenumber convective mode (IA refer to it as a “unipolar outflow”). Higher wavenumber convective modes are evidently suppressed by the viscosity. Figure A1 shows time-averaged plots of various quantities in the fiducial run. The time average is performed from 200 equally spaced data dumps from t = 4.3 × 105 GM/c3 to t = 7.3 × 105 GM/c3 . We show only the region rin < r < 30GM/c2 , whereas the computational domain is much larger: rin < r < 600GM/c2. The injection point is located far outside the plotted domain at r = 496GM/c2 . Because of the strong time-dependence of the flow in the fiducial run, the flow at any instant may look very different from these time averaged plots. The upper left panel in Figure 1 shows the average density; notice that the density is not symmetric about the equator. This is because the flow involves long-timescale quasi-periodic variations which are not quite averaged out over the course of the run. The upper right corner shows 3 5 the Bernoulli parameter Be. Dotted lines are negative; notice that there is a substantial amount of fluid near the equator that is unbound in the sense that Be > 0. Nevertheless, this material is still flowing inward in a nearly laminar fashion. Near the poles, the time-averaged Be < 0, but this region experiences large fluctuations. Polar outflows are typically associated with positive fluctuations in Be. The lower left panel shows the scaled mass flux r2 sin θ(ρv). Notice that much of the mass flux is along the surfaces of the inflow rather than along the equator. The lower right panel shows the scaled angular momentum flux r3 sin2 θ(ρvvφ + Π · φ̂). As for the mass flux, most of the activity is along the surface of the flow. Figure A2 shows the time series of the reduced eigen˙ values: Ṁ /Ṁinj , e = Ė/(Ṁ c2 ), and l = Lc/(GM Ṁ ). Also shown as dashed lines are the thin disk values for e and l. These assume a thin, cold disk terminating at r = 6GM/c2 . The low value of l is due to two effects. First, the disk is already sub-Keplerian by the time the flow reaches the innermost stable circular orbit. In addition, there are residual viscous torques in the plunging region that lower the specific angular momentum of the accreted material (see Figure A3, below). Notice that Figure A2 shows a smooth evolution that varies on a timescale τ ≈ 4 × 104 at late time. The largest variations in mass accretion rate are related to the appearance of large convective bubbles. Figure A3 shows the θ and time averaged run of several quantities with radius. The averages are taken over 4.3×105GM/c3 < t < 7.3×105GM/c3 and |θ−π/2| < π/6 3 . The upper left plot shows the run of density. Notice that here, as for the other quantities, there is a spike near rinj , an intermediate region, and then an inner, roughly power-law region. The upper right plot shows (cs /c)2 ; the lower left shows |vr |/c. Notice that the radial velocity exceeds the speed of light at the inner boundary. Similarly the azimuthal velocity vφ /c shown in the lower right panel approaches the speed of light. Also shown in that panel is the circular velocity (dashed line). Evidently the flow is slightly sub-Keplerian at most radii. The radial run of flow quantities in the inner regions can be fit by power laws, as done by IA and SPB. Our best fit power laws for the fiducial model (Run A) over 2.7GM/c2 < r < 20GM/c2 are ρ ∝ r−0.6 , cs ∝ r−0.5 , |vr | ∝ r−2 , and vφ ∝ r−0.8 . Between 2.7GM/c2 < r < 6GM/c2 , vφ is best fit by vφ ∝ r−0.9 , which is nearly, but not exactly, consistent with conservation of fluid specific angular momentum (vφ ∝ r−1 ). Angular momentum is not exactly conserved at r < 6GM/c2 because of viscous torques. The careful reader may notice that the power law slopes quoted in the last paragraph are not consistent with a constant mass accretion rate. This is because the power laws are derived from averages over |θ − π/2| < π/6, following IA. If one averages over all θ and time, then the resulting profiles are consistent with constant mass, energy, and angular momentum accretion rates, as they must be for a flow that is steady when averaged over large times. 4.2. Dependence on Inner Boundary Location and Gravitational Potential Averaging over |θ − π| < π/36 produces nearly identical results, but we have chosen to use IA’s range in θ. 6 McKinney & Gammie Having established that the differences between SPB and IA’s models are due to model choices rather than numerics, we were also interested in studying whether any features of global viscous accretion models are strongly dependent on numerical parameters. The first parameter we considered was the location of the inner boundary. In models that use a Newtonian gravitational potential (such as IA and SPB) the location of the inner boundary is not an interesting parameter in the sense that there is no physical lengthscale that one can compare rin to: it is simply a scaling parameter. In models that use a pseudoNewtonian potential, however, there is a feature (a “pit”) in the potential on a lengthscale GM/c2 . Starting with our fiducial model, then, what is the effect of shifting rin ? Our fiducial run has rin = 2.7GM/c2 . This may be compared with Run B, which has rin = 6GM/c2 . Figure A4 compares the accretion rates in the two runs. Evidently there are two changes in the solution. First, the time-averaged accretion rates differ by a large factor. The mean mass accretion rate is factor of 3 lower in Run B than Run A. The reduced eigenvalues e and l also differ by about 50% (see Table 2). Second, the time variation of the accretion rates differs, with Run B showing far more short-timescale variations. The short-timescale variations are due to the interaction of unstable convective modes with the boundary conditions. Inspection of the runs reveals an enhancement of convection and turbulence near rin in Run B. The differences between Run B and Run A are caused by the boundary location. Gradual variation of rin (in models not discussed in detail here) reveals that if the flow on the inner boundary is everywhere and always supersonic, then the solution is similar to Run A. If the flow is subsonic, then the solution exhibits artifacts like those seen in Run B. Evidently forcing the flow to be supersonic on the inner boundary causally disconnects the flow from the boundary. 4 This eliminates nonphysical reflection of linear and nonlinear waves from the boundary and renders the precise implementation of the numerical boundary conditions irrelevant. We do not want the reader to think that this problem arises because we happened to choose the wrong numerical implementation of the boundary conditions. Our implementation is the standard ZEUS outflow boundary condition, and it is widely used in astrophysical problems. While it may be possible to implement more transparent boundary conditions in the context of other numerical schemes, a survey of the numerical literature shows that in multiple dimensions this is an area of active research (Roe 1989; Karni 1991; Dedner et al. 2001; Bruneau & Creusé 2001), and that no general solution to the problem has been found. Furthermore, a simple example shows that no local extrapolation scheme can work for all accretion problems. Consider a numerical model of a steady spherical inflow (Bondi flow) in a gravitational potential Ψ(r). Let us suppose that we are primarily interested in accurately measuring Ṁ . We know from the analytic solution of the problem that Ṁ depends on the shape of the potential everywhere outside the sonic point. If we place the inner boundary rin outside the sonic point and use a local extrapolation scheme, we won’t always get the correct answer because the local extrapolation doesn’t have any information about the shape of the potential between rin and the sonic point. Put differently, one can’t determine a global solution from local extrapolation at the boundary. The key point is that, while aspects of the solution may be accurate, Ṁ (and L̇ and Ė) are sensitive to the boundary conditions. It is worth noting that the outer boundary is always in causal contact with the flow, but does not cause the same type of artifacts as the inner boundary. Experiments show that the flow is qualitatively insensitive to the location and implementation of the outer boundary condition. The time averaged Ṁ , however, is sensitive to both rout and rinj /rout . The time averaged l and e scale out this mass dependence and so are qualitatively and quantitatively insensitive to both rout and rinj /rout . It is also worth noting the effects of changing the gravitational potential. Run C (identical to IA99 Model 5) is identical to Run B except that the potential is now Newtonian. It is qualitatively similar to Run B, but Ṁ is now a factor of 5 lower than Run A. Run C also has the property that l oscillates about 0.0. This is a problem if the focus of the simulation is measuring Ṁ or L̇. To summarize: the location of the inner radial boundary can determine the character of the flow. If the flow is everywhere and always supersonic (or super-fastmagnetosonic in MHD) on the inner boundary then boundary-related corruption of the flow is impossible. Since it is computationally expensive to place the inner boundary very deep in the potential (for our model, the time step dt ∼ (rin − 2GM/c2 )), the optimal location for the inner boundary is just inside the radius where the radial Mach number always exceeds 1. The results of IA and SPB do not focus on the timedependence of the accretion values, so much of their discussion is unaffected by their treatment of the inner radial boundary. As discussed below, there are small changes related to the appearance of outflows. 4.3. Comparison of Torus and Injection Models The torus and injection methods represent sharply different approaches to studying accretion flows. The equilibrium torus presents a physically well-posed problem, but the accretion flow is transient: no steady state can be achieved. The injection method reaches a quasi-steady state, but much of the computational domain is wasted on evolving the injection region, which has no astrophysical analog: it is nonphysical. It is natural to ask whether these two widely used schemes for supplying mass can be made comparable or used to measure any of the same quantities. We selected two runs, E (torus) and F (injection), that had similar mass distributions in an evolved state. The torus run was studied at a time when Ṁ was close to its maximum. Run F’s Ṁ is a factor of 10 larger than Run E’s. This difference might have been anticipated from the sensitivity of the injection run to rout and rinj /rout : runs in which mass is concentrated closer to the outer boundary tend to have lower accretion rates because more of the 4 Although the viscous fluid equations of motion are not hyperbolic, and the flow in the supersonic region is in principle in causal contact with the rest of the flow, the coupling is exponentially weak. Numerical Viscous Accretion Flows mass escapes through the outer boundary. The time averaged mass accretion rate is therefore strongly dependent on the method of mass supply. The energy and angular momentum accretion rates per unit mass are, however, insensitive to the experimental design. We find that e and l differ by less than 3% in Runs E and F (see Table A1 and Table A2). These quantities are apparently set by conditions near the inner boundary (the ISCO for the Pseudo-Newtonian potential), and can be measured in either type of experiment. 4.4. Other Parameters We have varied rinj and rout /rinj , and as reported above, these strongly affect the time-averaged value of Ṁ . The sense of the effect is that a simulation with a larger rout /rinj loses less matter through the outer boundary, and this results in more matter streaming back into the black hole (by up to a factor of 10). The qualitative nature of the flow, however, is roughly independent of rout /rinj in that, e.g., the temporal power spectrum of Ṁ is similar. The qualitative nature of the flow is dependent on rinj . If one fixes rout /rinj and all other parameters, the range in α where unipolar outflows are observed tends to become smaller and disappears altogether for rinj as small as 40GM/c2 . The dependence of accretion models similar to ours on α has already been investigated by IA. They find that the flow changes from turbulent to laminar as α is increased and the higher viscosity damps modes of increasing lengthscale. IA find that α . 0.03 the flow is turbulent, and for α & 0.3 the flow is laminar. For 0.03 < α < 0.3 the flow exhibits a “unipolar” outflow. Our results are in agreement with IA. However, our models with a pseudo-Newtonian potential and super-sonic flow at the inner radial boundary show a slight shift in the values of α that exhibit unipolar outflows. We did find a critical value of α ≈ 0.5 above which no supersonic flow at rin could be achieved due to viscous heating, at least for γ = 3/2 and γ = 5/3 and for rin ≥ 2.1GM/c2 . Smaller values of rin were not computationally practical. This high α is typically associated with a bipolar outflow, as seen by IA. Even in this case, however, a choice of rin = 2.1GM/c2 instead of rin = 6GM/c2 leads to a qualitatively different profile for the flow. The flow with smaller rin = 2.1GM/c2 has a bipolar outflow starting at larger radius (10GM/c2 ) rather than immediately on the boundary as with rin = 6GM/c2 , and the mass accretion rate increases by a factor of 3. Finally, we studied the dependence of the results on numerical resolution. We find that Nr × Nθ = 108 × 50 is sufficient at α = 0.1 to resolve the shortest wavelength convective mode. Also, we chose our value of r0 to agree with SPB99’s torus models. We experimented with varying r0 and find that, all things being equal, smaller r0 gives more laminar flow. 5. summary Work in this field will shortly focus on global MHD models in pseudo-Newtonian potentials and in full general relativity. In our view it is useful to understand the solution space for physically and numerically simpler viscous models before turning to MHD. It is even possible, as Stone & 7 Pringle (2001) have claimed, that viscous hydrodynamics provides a crude approximation to the MHD results. In any event, this investigation provides a preview of some of the experimental issues that will play a role in most future global numerical investigations of accretion flows. This investigation was initially motivated by a desire to understand whether the differences between earlier global viscous hydrodynamics simulations performed by IA and SPB were caused by differences in experimental design or numerical method. IA and SPB reported different degrees of convective turbulence in their models and found different radial scalings for vertically averaged quantities such as temperature and density. Using a single code, we were able to reproduce both sets of results. We conclude that the differences are due to experimental design. We also found, while reproducing IA and SPB’s results, that some aspects of our solutions were sensitive to the numerical treatment of the region close to the inner boundary in models that use a pseudo-Newtonian potential. In particular, l and e, the specific angular momentum and energy of accreted material, are strongly dependent on rin , the location of the inner boundary. When the flow is supersonic at rin the location of the boundary does not affect l and e. But when the flow is subsonic at rin the flow interacts strongly with the numerical boundary condition. This produces spurious outflow events and makes l and e dependent on rin . Evidently for accurate measurement of these quantities it is necessary to isolate the numerical boundary condition behind a sonic transition that is located within the computational domain; one must place the inner boundary condition inside a “sound horizon”. We are not saying that all models that lack a sonic transition in the computational domain are fatally flawed. Whether the treatment of the inner boundary condition is problematic or not depends on what is being measured. For the nonlinear eigenvalues L̇, Ė, and Ṁ that we have focused on here, however, the treatment of the inner boundary condition is crucial. Furthermore, the only guarantee that the inner boundary condition is not governing the solution is to isolate it behind a sonic transition; this is the only completely safe choice. The location of the inner boundary may prove even more crucial in MHD models. Much of the character of the flow is determined in the turbulent, energetically important region of the flow just outside the fast magnetosonic transition, just as the region immediately outside the sonic transition determines the nature of the viscous flows described in this paper. We have performed some preliminary numerical tests and find that, as in the viscous flow, l and e for an MHD flow are sensitive to the treatment of the inner boundary. For various reasons it may prove difficult to achieve a fast magnetosonic transition in the computational domain; accurate treatment of the inner boundary condition may require fully (general) relativistic MHD. We have also compared two commonly used experimental designs for black hole accretion flow studies: models that begin with an equilibrium torus, and models that continuously inject fluid onto the grid. The choice between these models is to some degree a matter of taste. We find the equilibrium torus slightly easier to initialize and analyze. Remarkably, the two different approaches produce indistinguishable measurements for l and e, the specific 8 McKinney & Gammie angular momentum and energy of the accreted material. A parallel, viscous, axisymmetric hydrodynamics code based on that used in this paper can be found at http://kerr.physics.uiuc.edu. This work was supported in part by a GE fellowship and NASA GSRP Fellowship Grant NGT5-50343 to JCM, an NCSA Faculty Fellowship for CFG, the UIUC Research Board, and NSF Grant AST-0093091. Computations were done in part under National Computational Science Alliance grants AST010012N and AST010009N using the Origin 2000 and posic Linux cluster at NCSA. We thank Stu Shapiro and the referee for comments that substantially improved this paper. APPENDIX rate of strain tensor The rate of strain tensor e is a symmetric tensor that in spherical polar coordinates has ∂vr 1 err = − (∇ · v), ∂r 3 1 ∂vθ vr 1 eθθ = + − (∇ · v), r ∂θ r 3 vθ 1 1 ∂vφ vr + cot θ − (∇ · v) + , eφφ = r r 3 r sin θ ∂φ 1 ∂vr 1 ∂ vθ (r ( ) + ), 2 ∂r r r ∂θ 1 ∂vr 1 ∂ vφ ), = (r ( ) + 2 ∂r r r sin θ ∂φ erθ = erφ and eθφ = 1 sin θ ∂ vφ 1 ∂vθ ( ( )+ ). 2 r ∂θ sin θ r sin θ ∂φ (A1) (A2) (A3) (A4) (A5) (A6) REFERENCES Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214 Balbus, S. A. & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 Blandford, R. & McKee, C. 1982, ApJ, 255, 419 Blandford, R. D. & Payne, D. G., MNRAS, 199, 883 Blondin, J., Borkowski, K., & Reynolds, S. 2001, ApJ, 557, 782 Bondi, H, 1952, MNRAS, 112, 195 Bruneau C. H. & Creusé E. 2001, Int. J. Numer. Meth. Fluids, 36, 807 Dedner A. et al. 2001, J. Comp. Phys., 171, 448 Doeleman, S. S. et al. 2001, ApJ, 121, 2610 Gebhardt, K. et al.2001, ApJ, 543, 5 Hawley, J. F., Wilson, J., & Smarr, L. 1984, ApJS, 55, 211 Hawley, J. F. 1991, ApJ, 381, 496 Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 223 Hawley, J. F., Gammie, C. F.,& Balbus, S. A. 1995, ApJ, 440, 742 Hawley, J. F. 2000, ApJ, 528, 462 Hawley, J. F. 2001, ApJ, 554, 534 Hawley, J. F. & Krolik, J. 2001, ApJ, 548, 348. Hawley, J. F., Balbus, S.A., & Stone, J. 2001, ApJ, 554, 49 Igumenshchev, I. V. & Abramowicz, M. A. 1999, MNRAS, 303, 309 Igumenshchev, I. V. & Abramowicz, M. A. 2000, ApJS, 130, 463 Igumenshchev, I. V., Abramowicz, M. A., & Narayan, R. 2000, ApJ, 537, L27 Junor, W. & Biretta, J. & Livio, M. 1999, Nature, 401, 891 Karni, S. 1991, Int. J. Numer. Meth. Fluids, 13, 201 Koide, S., Shibata, K., & Kudoh, T. 1999, ApJ, 522, 727 Koide, S., Meier, D. L., Shibata, K., & Kudoh, T. 2000, ApJ, 536, 668 Lo, K. Y., Shen, Z., Zhao, J., & Ho, P. T. P. 1998, ApJ, 508, 61 Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603 Meier, D., Koide, S., & Uchida, Y. 2001, Science, 291, 84 Narayan, R. & Popham, R. 1993, Nature, 362, 820 Paczynski, B. & Wiita, P. J. 1980, A&A, 88, 23 Papaloizou, J. & Pringle, J. 1984, MNRAS, 208, 721 Rees, M. J. 2001, preprint (astro-ph/0101266) Roe, P. L. 1989, Comp. Fluids, 17, 221 Salpeter, E. E. 1964, ApJ, 140, 796 Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337 Stone, J. M., Pringle, J. E., & Begelman, M. C. 1999, MNRAS, 310, 1002 Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753 Stone, J. M. & Pringle, J. E. 2001, MNRAS, 322, 461 Zel’dovich Ya. B. 1964, Sov. Phys.-Dokl., 9, 195 Numerical Viscous Accretion Flows Fig. A1.— Time-averaged spatial structure of fiducial run (Run A; α = 0.1, rin = 2.7GM/c2 , and rout = 600GM/c2 ). Shown are the density (upper left), Bernoulli parameter (Be = (1/2)v2 + c2s /(γ − 1) + Ψ) (upper right; dotted line is a negative contour), scaled mass flux r 2 sin θ(ρv) (lower left), and scaled angular momentum flux r 3 sin2 θ(ρvvφ + Π · φ̂) (lower right). The flow is not symmetric about the equator because the flow exhibits long timescale antisymmetric variations. Convective bubbles form at the interface between positive and negative Bernoulli parameter (i.e. unbound and bound matter). 9 10 McKinney & Gammie Fig. A2.— The evolution of Ṁ /Ṁinj , e = Ė/(Ṁ c2 ), and l = L̇ c/(GM Ṁ ) in the fiducial run (Run A). The dotted line indicates the thin disk value. The run has clearly entered a quasi-steady state. The evolution is relatively smooth with a small variation on a timescale τ ≈ 4 × 104 . This is the timescale for convective bubble formation (the low point in mass accretion rate is when bubble forms). For this model the bubble forms at alternate poles. A full cycle requires of order one rotation period at the injection radius. Numerical Viscous Accretion Flows 11 Fig. A3.— The radial run of θ and time averaged quantities from the fiducial run (Run A). Shown are the density (upper left), squared sound speed (upper right), radial velocity (lower left), specific angular momentum (lower right; solid line), and circular orbit specific angular momentum (lower right; dashed line). Crudely speaking, the inner flow is consistent with a radial power law. The best fits to a power law are: ρ ∝ r −0.6 , cs ∝ r −0.5 , |vr | ∝ r −2 , and vφ ∝ r −0.8 . The plots are averaged over θ = π/2 ± π/6. 12 McKinney & Gammie Fig. A4.— The effect of moving the inner boundary on the accretion rates of mass, angular momentum, and energy (Run A vs. Run B). The top panel shows Ṁ /Ṁinj , the middle panel Ė/(Ṁ c2 ), and the bottom panel L̇ c/(GM Ṁ ). The solid curve is Run A, which has rin = 2.7GM/c2 . The dashed curve is Run B, which has rin = 6GM/c2 . Evidently Run B has a different variability structure and different time averaged values for the accretion rates. The relatively rapid and high-amplitude variations in Run B are due to nonphysical interactions with the inner radial boundary. Only by ensuring a supersonic flow (as in Run A) can one avoid these nonphysical effects. Numerical Viscous Accretion Flows 13 Table A1 Parameter List Run Nr A B C D E F 108 80 80 64 128 128 Nθ 50 50 50 40 80 80 Visc. IA IA IA MG MG MG Potential PN PN Newt. PN PN PN Rin /rg 1.35 3 3 1.2 1.4 1.4 Rout /rg 300 300 300 76 21 81 Rinj /rg 248 248 248 62 17 21 γ α 3/2 3/2 3/2 3/2 5/3 5/3 0.1 0.1 0.1 1.0 0.01 0.01 tf (c3 /GM) 7.3 × 105 4.4 × 105 7.3 × 105 3.3 × 103 3.0 × 104 6.8 × 104 Note. — IA and MG are viscosity prescription described in equations 7-9. PN is the pseudo-Newtonian potential of Paczynski & Wiita (1980). Run B uses Run C as initial conditions. rg = 2GM/c2 . Rinj for Run F is the position of the torus density peak ρ0 . 14 McKinney & Gammie Table A2 Results List Run A B C D E F Steady State Time (GM/c3 ) 5 4.3 × 10 2.4 × 105 2.4 × 105 ≥ 3.3 × 103 5.5 × 103 2.0 × 104 Max. Mach at Rin Ṁ -1.4 +0.0 +0.0 -0.4 -3.0 -3.4 5.96 × 10 1.95 × 10−2 9.46 × 10−3 ≥ 3.59 × 10−2 3.48 × 10−2 5.03 × 10−1 −2 −(Ė/(Ṁ c2 )) × 10−2 L̇ c/(GM Ṁ ) 2.06 3.72 6.77 4.64 3.01 3.11 1.75 1.29 .0746 −0.167 3.41 3.35 Note. — Runs A-E are injection runs with mass accretion rate units in Ṁinj and Run F is a torus run with mass accretion rate unit in ρ0 (GM)2 /c3 . Run C’s angular momentum fluctuations are 10 times the average value shown.