text:
# Numerical Models of Viscous Accretion Flows Near Black Holes

## 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/c^2`$) 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 $`\alpha `$ 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 $`\alpha P`$, where $`\alpha `$ 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.”

The $`\alpha `$ 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 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 $`\alpha `$ 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/c^2`$ 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`$, $`\theta `$, and $`\varphi `$.

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

$$\frac{D\rho }{Dt}+\rho (𝐯)=0,$$
(1)
momentum,

$$\rho \frac{D𝐯}{Dt}=P\rho \mathrm{\Psi }𝚷,$$
(2)
and energy,

$$\frac{Du}{Dt}=(P+u)(𝐯)+\mathrm{\Phi }.$$
(3)
Here, as usual, $`D/Dt/t+𝐯`$ is the Lagrangian time derivative, $`\rho `$ is the mass density, $`𝐯`$ is the velocity, $`u`$ is the internal energy density, $`P`$ is the pressure, and $`\mathrm{\Psi }`$ is the gravitational potential. The dissipation function $`\mathrm{\Phi }`$ is given by the product of the anomalous stress tensor $`𝚷`$ with the rate-of-strain tensor $`𝐞`$ (given explicitly in the Appendix)

$$\mathrm{\Phi }=\mathrm{\Pi }_{ij}e^{ij},$$
(4)
(sum over indices) where the anomalous stress tensor is the term-by-term product

$$\mathrm{\Pi }_{ij}=2\rho \nu e_{ij}S_{ij},$$
(5)
(no sum over indices) where $`S_{ij}`$ 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=(\gamma 1)u.$$
(6)
For the gravitational potential we use the pseudo-Newtonian potential of Paczynski & Wiita (1980): $`\mathrm{\Psi }=GM/(rr_g)`$ (here $`r_g2GM/c^2`$). This potential reproduces features of the orbital structure of a Schwarzschild spacetime, including an innermost stable circular orbit (ISCO) located at $`r=6GM/c^2`$. In a few cases we use the Newtonian potential $`\mathrm{\Psi }=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 $`𝐒`$ are $`1`$ (the “IA prescription”). SPB, on the other hand, use the Navier-Stokes prescription with all components zero except $`S_{r\varphi }`$, $`S_{\theta \varphi }`$, $`S_{\varphi r}`$, and $`S_{\varphi \theta }`$ (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,

$$\nu =\alpha (c_s^2/\mathrm{\Omega }_K),$$
(7)
$$\nu =\alpha (\rho /\rho _0)\mathrm{\Omega }_0r_0^2,$$
(8)
$$\nu =\alpha (c_s^2/\mathrm{\Omega }_K)\mathrm{sin}^{3/2}(\theta ),$$
(9)
are the IA, SPB, and MG prescriptions, respectively, where $`c_s=\sqrt{\gamma P/\rho }`$ is the sound speed, and

$$\mathrm{\Omega }_K^2\frac{1}{r}\frac{\mathrm{\Psi }}{r}=\frac{GM}{r(rr_g)^2}$$
(10)
is the “Keplerian” angular velocity. Here $`\rho _0`$ and $`\mathrm{\Omega }_0`$ are values of $`\rho `$ and $`\mathrm{\Omega }`$ at a fiducial radius $`r_0`$.

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 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 $`\theta =\pm \pi /2`$ and at $`r=r_{in},r_{out}`$. At the $`\theta `$ 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\rho ^\gamma `$, and a power-law rotation profile, $`\mathrm{\Omega }(r\mathrm{sin}\theta )^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 $`r_0`$ ; (2) the location of the innermost edge of the torus, $`r_{t,in}<r_0`$; (3) the maximum value of the density, $`\rho _0=\rho (r_0)`$; (4) the angular velocity gradient $`q=d\mathrm{ln}\mathrm{\Omega }/d\mathrm{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 $`\dot{\rho }(r,\theta )`$ which is typically symmetric about the equator. Parameters for this scheme include: (1) a characteristic radius for injection $`r_{inj}`$; (2) the rate of mass injection $`\dot{M}_{inj}`$; (3) the specific angular momentum of the injected fluid, $`v_\varphi =f_1r\mathrm{\Omega }_K`$. We usually set $`f_1=0.95`$, so that the fluid circularizes near $`r_{inj}`$. This restricts transients associated with circularization to the outer portions of the computational domain; (4) the internal energy of the injected fluid, $`u=f_2\rho \mathrm{\Psi }`$. We always set $`f_2=0.2`$ so that the fluid is marginally bound, i.e. has Bernoulli parameter $`Be(1/2)v^2+c_s^2/(\gamma 1)+\mathrm{\Psi }<0`$.

One must also choose the injection pattern $`\dot{\rho }(r,\theta )`$. IA99 choose a radially narrow region, but do not explicitly give $`\dot{\rho }(r,\theta )`$. We use a Gaussian, but found that none of the results 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 $`\theta `$ 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 $`\theta `$ 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 $`\sigma _r=0.05(r_{in}+r_{out})/2`$ and $`\sigma _\theta =\pi /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 ). 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 ($`\rho `$, $`u`$, and $`𝐯`$) 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 $`v_r(r_{in})=0`$ if $`v_r(r_{in})>0`$ and $`v_r(r_{out})=0`$ if $`v_r(r_{out})<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 $`\mathrm{log}(rr_g)`$. We require that $`dr(r)/(r_{in}r_g)1/4`$ so that the structure of the pseudo-Newtonian potential is well resolved. The grid is uniform in $`\theta `$. The grid has $`N_r\times N_\theta `$ 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 $`\rho _{fl}`$ 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, $`\dot{M}_{fl}`$ and $`\dot{E}_{fl}`$.

We set $`\rho _{fl}=10^{10}\dot{M}_{inj}c^3/(GM)^2`$ for injection runs and $`\rho _{fl}=10^5\rho _0`$ for torus runs. Lower values for $`\rho _{fl}`$ do not lead to a significant change in the solution. Larger values of $`\rho _{fl}`$ give $`\dot{M}_{fl}\dot{M}`$, 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 $`\rho _{fl}`$ and the internal energy density is $`u=U_o\rho \mathrm{\Psi }`$, where $`U_o`$ is a constant fraction of order unity (e.g. IA and we choose $`U_o=0.2`$). The addition of the atmosphere has no effect on the solution since the mass source’s evolution 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 $`u_{fl}`$ 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 $`\dot{E}_{fl}`$.

### 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

$$\dot{M}=\underset{S}{}\rho 𝐯𝑑𝐒,$$
(11)
where $`S`$ is the inner radial surface of the computational domain. The total energy accretion rate is

$$\dot{E}=\underset{S}{}((\frac{1}{2}v^2+h+\mathrm{\Psi })\rho 𝐯+𝚷𝐯)𝑑𝐒,$$
(12)
where $`h=(u+P)/\rho =\gamma u/\rho `$ is the specific enthalpy with our equation of state. The angular momentum accretion rate is

$$\dot{L}=\underset{S}{}r\mathrm{sin}\theta (\rho 𝐯v_\varphi +𝚷\widehat{\varphi })𝑑𝐒.$$
(13)
It is also sometimes useful to focus on the reduced eigenvalues $`l=\dot{L}/\dot{M}`$ and $`e=\dot{E}/\dot{M}`$. These value of mass, energy, and angular momentum are recorded at about 2 grid zones away from $`R_{in}`$. 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

$$\dot{E}_{err}=\dot{E}_{vol}+\dot{E}+\dot{E}_{out}\dot{E}_{fl},$$
(14)
where $`\dot{E}_{vol}`$ is the rate of change of the volume integrated total energy, $`\dot{E}_{out}`$ is the flux of total energy through the outer radial boundary, and $`\dot{E}_{fl}`$ 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, $`\dot{E}_{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 $`\dot{E}_{err}/E`$. For all runs we performed the error rate is within $`10\%`$ of $`10^5c^3/GM`$ for a torus run and within $`10\%`$ of $`10^4c^3/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\times 10^7`$ 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 $`\dot{M}r`$, $`\rho r^0`$, $`c_s^2r^1`$, $`v_\varphi r^{1/2}`$, and $`|v_r|r^1`$. These power law slopes are identical to those reported by SPB99. As another example, we found $`\dot{M}=1.23\times 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 $`\dot{M}=1.0\times 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 $`c_s/V_K`$ and specific angular momentum. Our results do not agree precisely, but this is likely due to small differences in mass injection scheme (because IA99 do not give their $`\dot{\rho }(r,\theta )`$). 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 1 and Table 2). The fiducial model has $`r_{in}=2.7GM/c^2`$, $`r_{out}=600GM/c^2`$, $`r_{inj}=495GM/c^2`$, $`\gamma =3/2`$, $`\alpha =0.1`$, $`N_r=108`$, $`N_\theta =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\times 10^5GM/c^3`$.

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+c_s^2/(\gamma 1)+\mathrm{\Psi }<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 1 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\times 10^5GM/c^3`$ to $`t=7.3\times 10^5GM/c^3`$. We show only the region $`r_{in}<r<30GM/c^2`$, whereas the computational domain is much larger: $`r_{in}<r<600GM/c^2`$. The injection point is located far outside the plotted domain at $`r=496GM/c^2`$. 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 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 $`r^2\mathrm{sin}\theta (\rho 𝐯)`$. 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 $`r^3\mathrm{sin}^2\theta (\rho 𝐯v_\varphi +𝚷\widehat{\varphi })`$. As for the mass flux, most of the activity is along the surface of the flow.

Figure 2 shows the time series of the reduced eigenvalues: $`\dot{M}/\dot{M}_{inj}`$, $`e=\dot{E}/(\dot{M}c^2)`$, and $`l=\dot{Lc}/(GM\dot{M})`$. Also shown as dashed lines are the thin disk values for $`e`$ and $`l`$. These assume a thin, cold disk terminating at $`r=6GM/c^2`$. 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 3, below). Notice that Figure 2 shows a smooth evolution that varies on a timescale $`\tau 4\times 10^4`$ at late time. The largest variations in mass accretion rate are related to the appearance of large convective bubbles.

Figure 3 shows the $`\theta `$ and time averaged run of several quantities with radius. The averages are taken over $`4.3\times 10^5GM/c^3<t<7.3\times 10^5GM/c^3`$ and $`|\theta \pi /2|<\pi /6`$ <sup>1</sup><sup>1</sup>1Averaging over $`|\theta \pi |<\pi /36`$ produces nearly identical results, but we have chosen to use IA’s range in $`\theta `$.. The upper left plot shows the run of density. Notice that here, as for the other quantities, there is a spike near $`r_{inj}`$, an intermediate region, and then an inner, roughly power-law region. The upper right plot shows $`(c_s/c)^2`$; the lower left shows $`|v_r|/c`$. Notice that the radial velocity exceeds the speed of light at the inner boundary. Similarly the azimuthal velocity $`v_\varphi /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/c^2<r<20GM/c^2`$ are $`\rho r^{0.6}`$, $`c_sr^{0.5}`$, $`|v_r|r^2`$, and $`v_\varphi r^{0.8}`$. Between $`2.7GM/c^2<r<6GM/c^2`$, $`v_\varphi `$ is best fit by $`v_\varphi r^{0.9}`$, which is nearly, but not exactly, consistent with conservation of fluid specific angular momentum ($`v_\varphi r^1`$). Angular momentum is not exactly conserved at $`r<6GM/c^2`$ 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 $`|\theta \pi /2|<\pi /6`$, following IA. If one averages over all $`\theta `$ 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

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 $`r_{in}`$ to: it is simply a scaling parameter. In models that use a pseudo-Newtonian potential, however, there is a feature (a “pit”) in the potential on a lengthscale $`GM/c^2`$. Starting with our fiducial model, then, what is the effect of shifting $`r_{in}`$?

Our fiducial run has $`r_{in}=2.7GM/c^2`$. This may be compared with Run B, which has $`r_{in}=6GM/c^2`$. Figure 4 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 $`r_{in}`$ in Run B.

The differences between Run B and Run A are caused by the boundary location. Gradual variation of $`r_{in}`$ (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. <sup>2</sup><sup>2</sup>2Although 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. 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$`\stackrel{´}{e}`$, 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 $`\mathrm{\Psi }(r)`$. Let us suppose that we are primarily interested in accurately measuring $`\dot{M}`$. We know from the analytic solution of the problem that $`\dot{M}`$ depends on the shape of the potential everywhere outside the sonic point. If we place the inner boundary $`r_{in}`$ 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 $`r_{in}`$ 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, $`\dot{M}`$ (and $`\dot{L}`$ and $`\dot{E}`$) 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 $`\dot{M}`$, however, is sensitive to both $`r_{out}`$ and $`r_{inj}/r_{out}`$. The time averaged $`l`$ and $`e`$ scale out this mass dependence and so are qualitatively and quantitatively insensitive to both $`r_{out}`$ and $`r_{inj}/r_{out}`$.

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 $`\dot{M}`$ 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 $`\dot{M}`$ or $`\dot{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-fast-magnetosonic 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(r_{in}2GM/c^2)`$), 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 time-dependence 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 $`\dot{M}`$ was close to its maximum. Run F’s $`\dot{M}`$ is a factor of $`10`$ larger than Run E’s. This difference might have been anticipated from the sensitivity of the injection run to $`r_{out}`$ and $`r_{inj}/r_{out}`$: runs in which mass is concentrated closer to the outer boundary tend to have lower accretion rates because more of the 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 1 and Table 2). 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 $`r_{inj}`$ and $`r_{out}/r_{inj}`$, and as reported above, these strongly affect the time-averaged value of $`\dot{M}`$. The sense of the effect is that a simulation with a larger $`r_{out}/r_{inj}`$ 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 $`r_{out}/r_{inj}`$ in that, e.g., the temporal power spectrum of $`\dot{M}`$ is similar. The qualitative nature of the flow is dependent on $`r_{inj}`$. If one fixes $`r_{out}/r_{inj}`$ and all other parameters, the range in $`\alpha `$ where unipolar outflows are observed tends to become smaller and disappears altogether for $`r_{inj}`$ as small as $`40GM/c^2`$.

The dependence of accretion models similar to ours on $`\alpha `$ has already been investigated by IA. They find that the flow changes from turbulent to laminar as $`\alpha `$ is increased and the higher viscosity damps modes of increasing lengthscale. IA find that $`\alpha 0.03`$ the flow is turbulent, and for $`\alpha 0.3`$ the flow is laminar. For $`0.03<\alpha <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 $`\alpha `$ that exhibit unipolar outflows.

We did find a critical value of $`\alpha 0.5`$ above which no supersonic flow at $`r_{in}`$ could be achieved due to viscous heating, at least for $`\gamma =3/2`$ and $`\gamma =5/3`$ and for $`r_{in}2.1GM/c^2`$. Smaller values of $`r_{in}`$ were not computationally practical. This high $`\alpha `$ is typically associated with a bipolar outflow, as seen by IA. Even in this case, however, a choice of $`r_{in}=2.1GM/c^2`$ instead of $`r_{in}=6GM/c^2`$ leads to a qualitatively different profile for the flow. The flow with smaller $`r_{in}=2.1GM/c^2`$ has a bipolar outflow starting at larger radius ($`10GM/c^2`$) rather than immediately on the boundary as with $`r_{in}=6GM/c^2`$, 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 $`N_r\times N_\theta =108\times 50`$ is sufficient at $`\alpha =0.1`$ to resolve the shortest wavelength convective mode. Also, we chose our value of $`r_0`$ to agree with SPB99’s torus models. We experimented with varying $`r_0`$ and find that, all things being equal, smaller $`r_0`$ 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 & 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 $`r_{in}`$, the location of the inner boundary. When the flow is supersonic at $`r_{in}`$ the location of the boundary does not affect $`l`$ and $`e`$. But when the flow is subsonic at $`r_{in}`$ the flow interacts strongly with the numerical boundary condition. This produces spurious outflow events and makes $`l`$ and $`e`$ dependent on $`r_{in}`$. 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 $`\dot{L},\dot{E}`$, and $`\dot{M}`$ 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 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 A Rate of Strain Tensor

The rate of strain tensor $`𝐞`$ is a symmetric tensor that in spherical polar coordinates has

$$e^{rr}=\frac{v_r}{r}\frac{1}{3}(v),$$
(A1)
$$e^{\theta \theta }=\frac{1}{r}\frac{v_\theta }{\theta }+\frac{v_r}{r}\frac{1}{3}(v),$$
(A2)
$$e^{\varphi \varphi }=\frac{v_r}{r}+\frac{v_\theta }{r}\mathrm{cot}\theta \frac{1}{3}(v)+\frac{1}{r\mathrm{sin}\theta }\frac{v_\varphi }{\varphi },$$
(A3)
$$e^{r\theta }=\frac{1}{2}(r\frac{}{r}(\frac{v_\theta }{r})+\frac{1}{r}\frac{v_r}{\theta }),$$
(A4)
$$e^{r\varphi }=\frac{1}{2}(r\frac{}{r}(\frac{v_\varphi }{r})+\frac{1}{r\mathrm{sin}\theta }\frac{v_r}{\varphi }),$$
(A5)
and

$$e^{\theta \varphi }=\frac{1}{2}(\frac{\mathrm{sin}\theta }{r}\frac{}{\theta }(\frac{v_\varphi }{\mathrm{sin}\theta })+\frac{1}{r\mathrm{sin}\theta }\frac{v_\theta }{\varphi }).$$
(A6)


format: text

text:
# Coupled Replicator Equations for the Dynamics of Learning in Multiagent Systems

## Abstract

Starting with a group of reinforcement-learning agents we derive coupled replicator equations that describe the dynamics of collective learning in multiagent systems. We show that, although agents model their environment in a self-interested way without sharing knowledge, a game dynamics emerges naturally through environment-mediated interactions. An application to rock-scissors-paper game interactions shows that the collective learning dynamics exhibits a diversity of competitive and cooperative behaviors. These include quasiperiodicity, stable limit cycles, intermittency, and deterministic chaos—behaviors that should be expected in heterogeneous multiagent systems described by the general replicator equations we derive.

preprint: Santa Fe Institute Working Paper 02-04-017

Adaptive behavior in multiagent systems is an important interdisciplinary topic that appears in various guises in many fields, including biology Camazine et al. (2001), computer science Simon (1996), economics Young (1998), and cognitive science Hutchins (1996). One of the key common questions is how and whether a group of intelligent agents truly engages in collective behaviors that are more functional than individuals acting alone.

Suppose that many agents interact with an environment and each independently builds a model from its sensory stimuli. In this simple type of coupled multiagent system, collective learning (if it occurs) is a dynamical behavior driven by agents’ environment-mediated interaction Rossler (1987); Taiji and Ikegami (1999). Here we show that the collective dynamics in multiagent systems, in which agents use reinforcement learning Sutton and Barto (1998), can be modeled using a generalized form of coupled replicator equations.

While replicator dynamics were introduced originally for evolutionary game theory Taylor and Jonker (1978), the relationship between reinforcement learning and replicator equations has been developed only recently Borgers and Sarin (1997). Here, we extend these considerations to multiagent systems, introducing the theory behind a previously reported game-theoretic model Sato et al., (2002). We show that replicator dynamics emerges as a special case of the continuous-time limit for multiagent reinforcement learning systems. The overall approach, though, establishes a general framework for dynamical-systems analyses of adaptive behavior in collectives.

Notably, in learning with perfect memory, our model reduces to the form of a multipopulation replicator equation introduced in Ref. Taylor (1979). For two agents with perfect memory interacting via a zero-sum rock-scissors-paper game the dynamics exhibits Hamiltonian chaos Sato et al., (2002). In contrast, as we show here, with memory decay multiagent systems generally become dissipative and display the full range of nonlinear dynamical behaviors, including limit cycles, intermittency, and deterministic chaos.

Our multiagent model begins with simple reinforcement learning agents. To simplify the development, we assume that there are two such agents $`X`$ and $`Y`$ that at each time step take one of $`N`$ actions: $`i=1,\mathrm{},N`$. Let the probability for $`X`$ to chose action $`i`$ be $`x_i(n)`$ and $`y_i(n)`$ for $`Y`$, where $`n`$ is the number of the learning iterations from the initial state at $`n=0`$. The agents’ choice distributions at time $`n`$ are $`𝐱(n)=(x_1(n),\mathrm{},x_N(n))`$ and $`𝐲(n)=(y_1(n),\mathrm{},y_N(n))`$, with $`\mathrm{\Sigma }_ix_i(n)=\mathrm{\Sigma }_iy_i(n)=1`$.

Let $`R_{ij}^X`$ and $`R_{ij}^Y`$ denote the reward for $`X`$ taking action $`i`$ and $`Y`$ action $`j`$ at step $`n`$, respectively. Given these actions, $`X`$’s and $`Y`$’s memories, $`Q_i^X(n)`$ and $`Q_i^Y(n)`$, of the past benefits from their actions are governed by

$`Q_i^X(n+1)Q_i^X(n)`$ $`=`$ $`R_{ij}^X\alpha _XQ_i^X(n)\mathrm{and}`$ (1)
$`Q_i^Y(n+1)Q_i^Y(n)`$ $`=`$ $`R_{ij}^Y\alpha _YQ_i^Y(n),`$
where $`\alpha _X,\alpha _Y[0,1)`$ control each agent’s memory decay rate and $`Q_i^X(0)=Q_i^Y(0)=0`$. The agents chose their next actions according to the $`Q`$’s, updating their choice distributions as follows:

$$x_i(n)=\frac{e^{\beta _XQ_i^X(n)}}{\mathrm{\Sigma }_je^{\beta _XQ_j^X(n)}}\mathrm{and}y_i(n)=\frac{e^{\beta _YQ_i^Y(n)}}{\mathrm{\Sigma }_je^{\beta _YQ_j^Y(n)}},$$
(2)
where $`\beta _X,\beta _Y[0,\mathrm{}]`$ control the learning sensitivity: how much the current choice distributions are affected by past rewards. Using Eq. (2), the dynamic governing the change in agent state is given by:

$$x_i(n+1)=\frac{x_i(n)e^{\beta _X(Q_i^X(n+1)Q_i^X(n))}}{\mathrm{\Sigma }_kx_j(n)e^{\beta _X(Q_k^X(n+1)Q_k^X(n))}},$$
(3)
and similarly for $`y_i(n+1)`$.

Consider the continuous-time limit corresponding to agents performing a large number of actions (iterates of Eqs. (1)) for each choice-distribution update (iterates of Eq. (3)). In this case, we have two different time scales—that for agent-agent interactions and for learning. We assume that the learning dynamics is very slow compared to interactions and so $`𝐱`$ and $`𝐲`$ are essentially constant during the latter. Then, based on Eq. (3), continuous-time learning for agent $`X`$ is governed by

$$\dot{x}_i=\beta _Xx_i(\dot{Q}_i^X\mathrm{\Sigma }_j\dot{Q}_j^Xx_j),$$
(4)
and for the dynamic governing memory updates we have

$$\dot{Q}_i^X=R_i^X\alpha _XQ_i^X,$$
(5)
where $`R_i^X`$ is the reward for $`X`$ choosing action $`i`$, averaged over $`Y`$’s actions during the time interval between learning updates. Putting together Eqs. (2), (4), and (5) one finds

$$\frac{\dot{x_i}}{x_i}=\beta _X[R_i^X\mathrm{\Sigma }_jx_{ij}R_j^X]+\alpha _XI_i^X,$$
(6)
where $`I_i^X\mathrm{\Sigma }_jx_j\mathrm{log}(x_j/x_i)`$ represents the effect of memory with decay parameter $`\alpha _X`$. (The continuous-time dynamic of $`Y`$ follows in a similar manner.) Eq. (6), extended to account for any number of agents and actions, constitutes our general model for reinforcement-learning multiagent systems.

Simplifying again, assume a fixed relationship between pairs $`(i,j)`$ of $`X`$’s and $`Y`$’s actions and between rewards for both agents: $`R_{ij}^X=a_{ij}`$ and $`R_{ij}^Y=b_{ij}`$. Assume further that $`𝐱`$ and $`𝐲`$ are independently distributed, then the time-average rewards for $`X`$ and $`Y`$ become

$$R_i^X=\mathrm{\Sigma }_ja_{ij}y_j\mathrm{and}R_i^Y=\mathrm{\Sigma }_jb_{ij}x_j,$$
(7)
In this restricted case, the continuous-time dynamic is:

$`{\displaystyle \frac{\dot{x}_i}{x_i}}`$ $`=`$ $`\beta _X[(A𝐲)_i𝐱A𝐲]+\alpha _XI_i^X,`$
$`{\displaystyle \frac{\dot{y}_i}{y_i}}`$ $`=`$ $`\beta _Y[(B𝐱)_i𝐲B𝐱]+\alpha _YI_i^Y,`$ (8)
where $`(A)_{ij}=a_{ij}`$ and $`(B)_{ij}=b_{ij}`$, $`(A𝐱)_i`$ is the $`i`$th element of the vector $`A𝐱`$, and $`\beta _X`$ and $`\beta _Y`$ control the time-scale of each agent’s learning.

We can regard $`A`$ and $`B`$ as $`X`$’s and $`Y`$’s game-theoretic payoff matrices for action $`i`$ against opponent’s action $`j`$ 111Eqs. (7) specify the von Neumann-Morgenstern utility (J. von Neumann and O. Morgenstern, Theory of Games and Economic Behavior, (Princeton University Press, 1944)).. In contrast with game theory, which assumes agents have exact knowledge of the game structure and of other agent’s strategies, reinforcement-learning agents have no knowledge of a “game” in which they are playing, only a myopic model of the environment—other agent(s)—given implicitly via the rewards they receive. Nonetheless, a game dynamics emerges—via $`R^X`$ and $`R^Y`$ in Eq. (6)—as a description of the collective’s global behavior.

Given the basic equations of motion for the reinforcement-learning multiagent system (Eq. (8)), one becomes interested in, on the one hand, the time evolution of each agent’s state vector in the simplices $`𝐱\mathrm{\Delta }_X`$ and $`𝐲\mathrm{\Delta }_Y`$ and, on the other, the dynamics in the higher-dimensional collective simplex $`(𝐱,𝐲)\mathrm{\Delta }_X\times \mathrm{\Delta }_Y`$. Following Ref. Hofbauer (1996), we transform from $`(𝐱,𝐲)\mathrm{\Delta }_X\times \mathrm{\Delta }_Y`$ to $`𝐔=(𝐮,𝐯)𝐑^{2(N1)}`$ with $`𝐮=(u_1,\mathrm{},u_{N1})`$ and $`𝐯=(v_1,\mathrm{},v_{N1})`$, where $`u_i=\mathrm{log}(x_{i+1}/x_1)\mathrm{and}v_i=\mathrm{log}(y_{i+1}/y_1),(i=1,\mathrm{},N1)`$. The result is a new version of our simplified model (Eqs. (8)), useful both for numerical stability during simulation and also for analysis in certain limits:

$`\dot{u}_i`$ $`=`$ $`\beta _X{\displaystyle \frac{\mathrm{\Sigma }_j\stackrel{~}{a}_{ij}e^{v_j}+\stackrel{~}{a}_{i1}}{1+\mathrm{\Sigma }_je^{v_j}}}\alpha _Xu_i\mathrm{and}`$
$`\dot{v}_i`$ $`=`$ $`\beta _Y{\displaystyle \frac{\mathrm{\Sigma }_j\stackrel{~}{b}_{ij}e^{u_j}+\stackrel{~}{b}_{i1}}{1+\mathrm{\Sigma }_je^{u_j}}}\alpha _Yv_i,`$ (9)
where $`\stackrel{~}{a}_{ij}=a_{i+1,j}a_{1,j}`$ and $`\stackrel{~}{b}_{ij}=b_{i+1,j}b_{1,j}`$. Since the dissipation rate $`\gamma `$ in $`𝐔`$ is

$$\gamma =\mathrm{\Sigma }_i\frac{\dot{u}_i}{u_i}+\mathrm{\Sigma }_j\frac{\dot{v}_j}{v_j}=(N1)(\alpha _X+\alpha _Y),$$
(10)
Eqs. (8) are conservative when $`\alpha _X=\alpha _Y=0`$ and the time average of a trajectory is the Nash equilibrium of the game specified by $`A`$ and $`B`$, if a limit set exists in the interior of $`\mathrm{\Delta }_X\times \mathrm{\Delta }_Y`$222Cf. P. Schuster et al, Biol. Cybern. 40, 1 (1981).. Moreover, if the game is zero-sum, the dynamics are Hamiltonian in $`𝐔`$ with

$`H`$ $`=`$ $`(\mathrm{\Sigma }_jx_j^{}u_j+\mathrm{\Sigma }_jy_j^{}v_j)`$
$`+`$ $`\mathrm{log}(1+\mathrm{\Sigma }_je^{u_j})+\mathrm{log}(1+\mathrm{\Sigma }_je^{v_j}),`$
where $`(𝐱^{},𝐲^{})`$ is an interior Nash equilibrium Hofbauer (1996).

To illustrate the dynamical-systems analysis of learning in multiagent systems using the above framework, we now analyze the behavior of the two-person rock-scissors-paper interaction 333Such interactions are observed in natural social and biological communities; cf. B. Kerr et al, Nature 418, 171 (2002).. This familiar game describes a nontransitive three-sided competition: rock beats scissors, scissors beats paper, and paper beats rock. The reward structure (environment) is given by:

$$A=\left[\begin{array}{ccc}ϵ_X& 1& 1\\ 1& ϵ_X& 1\\ 1& 1& ϵ_X\end{array}\right]\mathrm{and}B=\left[\begin{array}{ccc}ϵ_Y& 1& 1\\ 1& ϵ_Y& 1\\ 1& 1& ϵ_Y\end{array}\right],$$
(12)
where $`ϵ_X,ϵ_Y[1.0,1.0]`$ are the rewards for ties. The mixed Nash equilibrium is $`x_i^{}=y_i^{}=1/3,(i=1,2,3)`$—the centers of $`\mathrm{\Delta }_X`$ and $`\mathrm{\Delta }_Y`$. If $`ϵ_X=ϵ_Y`$, the game is zero-sum.

In the special case of perfect memory ($`\alpha _X=\alpha _Y=0`$) and with equal learning sensitivity ($`\beta _X=\beta _Y`$), the linear version (Eqs. (8)) of our model (Eq. (6)) reduces to multipopulation replicator equations Taylor (1979):

$$\frac{\dot{x}_i}{x_i}=\left[(A𝐲)_i𝐱A𝐲\right]\mathrm{and}\frac{\dot{y}_i}{y_i}=\left[(B𝐱)_i𝐲B𝐱\right].$$
(13)
The game-theoretic behavior in this case with rock-scissors-paper interactions (Eqs. (12)) was investigated in Sato et al., (2002). Here, before contrasting our more general setting, we briefly recall the behavior in these special cases, noting several additional results.

Figure 1 shows Poincaré sections of Eqs. (13)’s trajectories on the hyperplane $`(\dot{u}_1=0,\dot{v}_1>0)`$ and representative trajectories in the individual agent simplices $`\mathrm{\Delta }_X`$ and $`\mathrm{\Delta }_Y`$. When $`ϵ_X=ϵ_Y=0.0`$, we expect the system to be integrable and only quasiperiodic tori should exist. Otherwise, $`ϵ_X=ϵ_Y>0.0`$, Hamiltonian chaos can occur with positive-negative pairs of Lyapunov exponents Sato et al., (2002). The dynamics is very rich, there are infinitely many distinct behaviors near the unstable fixed point at the center—the classical Nash equilibrium—and a periodic orbit arbitrarily close to any chaotic one. Moreover, when the game is not zero-sum ($`ϵ_Xϵ_Y`$), transients to heteroclinic cycles are observed Sato et al., (2002): On the one hand, there are intermittent behaviors in which the time spent near pure strategies (the simplicial vertices) increases subexponentially with $`ϵ_X+ϵ_Y<0`$ and, on the other hand, with $`ϵ_X+ϵ_Y>0`$, chaotic transients persist; cf. Chawanya (1995).

Our framework goes beyond these special cases and, generally, beyond the standard multipopulation replicator equations (Eqs. (13)) due to its accounting for the effects of individual and collective learning and since the reward structure and the learning rules need not lead to linear interactions. For example, if the memory decay rates ($`\alpha _X`$ and $`\alpha _Y`$) are positive, the system becomes dissipative and exhibits limit cycles and chaotic attractors; see Fig. 2. Figure 3 (top) shows a diverse range of bifurcations as a function of $`ϵ_Y`$: dynamics on the hyperplane ($`\dot{u}_1=0`$, $`\dot{v}_1>0`$) projected onto $`y_1`$. When the game is nearly zero-sum, agents can reach the stable Nash equilibrium, but chaos can also occur, when $`ϵ_X+ϵ_Y>0`$. Figure 3 (bottom) shows that the largest Lyapunov exponent is positive across a significant fraction of parameter space; indicating that chaos is common. The dual aspects of chaos, irregularity and coherence, imply that agents may behave cooperatively or competitively (or switch between both) in the collective dynamics. Such global behaviors ultimately derive from self-interested, myopic learning.

Within this framework a number of extensions suggest themselves as ways to investigate the emergence of collective behaviors. The most obvious is the generalization to an arbitrary number of agents with an arbitrary number of strategies and the analysis of behaviors in thermodynamic limit; see, e.g., Marsili et al. (2000) as an alternative approach. It is relatively straightforward to develop an extension to the linear-reward version (Eqs. (8)) of our model. For three agents $`X`$, $`Y`$, and $`Z`$, one obtains:

$$\frac{\dot{x}_i}{x_i}=\beta _X[\mathrm{\Sigma }_{j,k}a_{ijk}y_jz_k\mathrm{\Sigma }_{j,k,l}a_{jkl}x_jy_kz_l]\alpha _XI_i^X,$$
(14)
with tensor $`(A)_{ijk}=a_{ijk}`$, and similarly for $`Y`$ and $`Z`$. Not surprisingly, this is also a conservative system when the $`\alpha `$’s vanish. However, extending the general collective learning equations (Eq. (6)) to multiple agents is challenging and so will be reported elsewhere.

To be relevant to applications, one also needs to develop a statistical dynamics generalization van Nimwegen et al. (1999) of the deterministic equations of motion to account for finite and fluctuating numbers of agents and also finite histories used in learning. Finally, another direction, especially useful if one attempts to quantify collective function in large multiagent systems, will be structural and information-theoretic analyses Crutchfield and Young (1989) of local and global learning behaviors and, importantly, their differences. Analyzing the stored information in each agent versus that in the collective, the causal architecture of information flow between an individual agent and the group, and how individual and global memories are processed to sustain collective function are projects now made possible using this framework.

We presented a dynamical-systems model of collective learning in multiagent systems, which starts with reinforcement learning agents and reduces to coupled replicator equations, demonstrated that individual-agent learning induces a global game dynamics, and investigated some of the resulting periodic, intermittent, and chaotic behaviors with simple (linear) rock-scissors-papers game interactions. Our model gives a macroscopic description of a network of learning agents that can be straightforwardly extended to model a large number of heterogeneous agents in fluctuating environments. Since deterministic chaos occurs even in this simple setting, one expects that in high-dimensional and heterogeneous populations typical of multiagent systems intrinsic unpredictability will become a dominant collective behavior. Sustaining useful collective function in multiagent systems becomes an even more compelling question in light of these results.

The authors thank J. D. Farmer, E. Akiyama, P. Patelli, and C. Shalizi. This work was supported at the Santa Fe Institute under the Network Dynamics Program by Intel Corporation and under DARPA agreement F30602-00-2-0583. YS’s participation was supported by the Postdoctoral Researchers Program at RIKEN.


text:
# Prospects for Detection of Intracluster Gas Bulk Velocities Through the Sunyaev-Zel’dovich Effect

## 1  Introduction

Clusters of galaxies are believed to form from the infall/merging of smaller scale systems. Memory of the formation process may be probed through the X-ray analysis of the intracluster medium (ICM), where most of the cluster’s visible mass resides. This is due to the fact that the merging of sub-units is predicted to create gas substructures such as temperature and density inhomogeneities, destruction of cooling flows & metal abundance gradients and generation of gas bulk velocities. The comparison of observed structures to numerical hydrodynamical+N body simulations can then provide clues as to the cluster’s evolutionary stage (e.g. Evrard 1990; Katz & White 1993; Navarro, Frenk & White 1995; Evrard, Metzler, & Navarro 1996; Roettiger, Burns & Loken 1993,1996; Schindler & Muller 1993; Pearce, Thomas & Couchman 1994; Gomez et al. 2001 and references therein).

Until very recently the physical characteristics of the ICM had been derived almost exclusively by the analysis of distributions of electron temperature ($`T_e`$) and surface brightness ($`S_X`$) (or gas density), even though in the last few years various numerical simulations of off-center cluster mergers have been published predicting the presence of long-lasting residual intracluster gas velocities approaching a few thousand km s-1 (e.g. Roettiger, Loken & Burns 1997, Ricker 1998; Roettiger, Stone, & Mushotzky 1998; Takizawa & Mineshige 1998; Burns et al. 1999; Takizawa 1999, 2000; Roettiger & Flores 2000). Furthermore, since we are always looking at 2-D projections, velocity maps may be a crucial addition, necessary to break the degeneracies associated with analyses based solely on $`T_e`$ & $`S_X`$ distributions.

Recently, Dupke & Bregman (2001a,b) detected ICM bulk motions in two galaxy clusters using ASCA data: Perseus (Abell 426) and Centaurus (Abell 3526). The intracluster gas velocity distributions detected in these two clusters are roughly consistent with systematic gas bulk rotation with correspondent circular velocities $`>`$ 1000 km s-1, implying that a significant fraction of the intracluster gas energy can be kinetic. Since X-ray measurements are haunted by the need of precise knowledge of instrumental gain 111conversion between pulse height and incoming photon energy these detections are significant only at the (2-4) $`\sigma `$ level. Realistically, the significance of these measurements are not expected to improve substantially within the next several years even with Chandra and XMM-Newton, since the spectral resolution of the imaging spectrometers on-board these satellites, ACIS and EPIC222We are not considering the gratings since they are less suited for this kind of spatially-resolved spectroscopy, is only slightly better than that of the SISs on-board ASCA. Furthermore, it will take a few years before detailed gain variations across those detectors are known well enough for reliable ICM velocity measurements. To compensate for the lack of knowledge of the gain variations in X-ray spectrometers it is necessary to take several consecutive off-center long-exposure observations of a cluster, so that sky regions with discrepant radial velocities can be observed at the same CCD position. However, this kind of measurement is very time consuming if one does not have prior knowledge of the best orientations for velocity measurements.

Alternatively, ICM velocity measurements can be, in principle, corroborated/detected by the use of the kinetic S-Z effect (Sunyaev & Zel’dovich 1970, 1972, 1980). Intracluster gas bulk velocities as high as those detected in A3526 should generate significantly different levels of Comptonization of the cosmic microwave background radiation (CMBR) towards different direction of the cluster (red-shifted and blue-shifted sides). This effect could be detected with current (or in development) instruments, such as the BOLOCAM333www.astro.caltech.edu/l̃gg/bolocam/bolocam.html, or, with smaller spatial resolution, ACBAR444astrophysics.phys.cmu.edu/research/viper/ and the High Frequency Instrument on-board PLANCK555astro.estec.esa.nl/Planck/.

In this paper we realistically predict differential variation of the CMBR temperature as it passes through the intracluster gas of the central regions of the Centaurus cluster, $`(\frac{\mathrm{\Delta }T}{T})_{CMBR}`$ $`10^4`$. We propose that this effect can be observed soon in cool clusters such as Centaurus at frequencies that minimize the thermal component (214 GHz). We suggest that differential S-Z measurements of clusters, coupled with X-ray spatially resolved spectroscopy can provide crucial information not just to cosmology but also to the study of intracluster gas dynamics.

## 2  The Mechanism

The CMBR intensity and temperature variations with respect to the background towards the direction of a rotating cluster has two components; a thermal and a kinetic component. The latter can be associated with some global peculiar velocity and also with internal bulk motions, which is the one that we are interested in this work. If the resulting radial velocity of the kinetic S-Z component at some projected radius “b” is denoted by Vr the thermal and kinetic variation of intensity ($`I`$) and temperature ($`T`$) of the CMB in the non-relativistic approximation can be can be given by (for reviews see, e.g., Rephaeli 1995; Birkinshaw 1999):

$`\mathrm{\Delta }I_\nu `$ $`2{\displaystyle \frac{(kT)^3}{(hc)^2}}{\displaystyle \frac{x^4e^x}{(e^x1)^2}}(K(T_e,\nu )\beta _r(b))\tau (b)`$ (1)
$$(\frac{\mathrm{\Delta }I}{I})_\nu \frac{xe^x}{e^x1}(K(T_e,\nu )\beta _r(b))\tau (b)$$
()
$$(\frac{\mathrm{\Delta }T}{T})_\nu (K(T_e,\nu )\beta _r(b))\tau (b)$$
()
where

$$K(T_e,\nu )=\frac{kT_e}{m_ec^2}(x\frac{e^x+1}{e^x1}4)$$
()
$$\beta _r(b)=\frac{V_r(b)}{c}andx=\frac{h\nu }{kT}$$
()
where $`T_e`$ & $`T`$ are respectively the ICM and CMB temperatures, and the other parameters have their usual meanings. If the gas number density $`n(r)`$ follows a king-like profile $`n(r)=n_0(1+(\frac{r}{r_c})^2)^{\frac{3}{2}\beta }`$, where $`r_c`$ and $`n_0`$ are respectively the core radius and the central density, the optical depth is given as a function of the projected radius “b” by

$$\tau (b)=\sigma _Tn_0r_cB(\frac{1}{2},\frac{3}{2}\beta \frac{1}{2})(1+(\frac{b}{r_c})^2)^{\frac{3}{2}\beta +\frac{1}{2}}$$
()
where $`B(p,q)=_0^{\mathrm{}}x^{p1}(1+x)^{p+q}𝑑x`$ is the Beta function of p,q.

## 3 The Case for Abell 3526

Centaurus (Abell 3526) is a BM type I, nearby (z$``$ 0.0104), X-ray bright cooling flow (accretion rate $`<`$ 30-50 Myr-1), cold (kTe$``$3-3.5 keV)) cluster. Its X-ray emission is relatively smooth except for the central 1 where Chandra found X-ray substructures (Sanders & Fabian 2001), and is peaked on the cD galaxy NGC 4696. The existence of galaxy velocity bi-modality in Centaurus had been shown by Lucey, Currie & Dickens (1986a,b) and was more recently confirmed by Stein et al. (1997). Two galaxy groups are clearly separated: the main group (Cen 30), which is centered on the cD galaxy (NGC 4696) shows an average radial velocity of 3397$`\pm `$139 km s-1 and a velocity dispersion of 933$`\pm `$118 km s-1. The second group (Cen 45) is associated with the galaxy NGC 4709 at $`>`$15 from NGC 4696. It has an average radial velocity of 4746$`\pm `$43 km s-1 and a velocity dispersion of 131$`\pm 4`$3 km s-1. One arcmin at Centaurus distance corresponds to $``$ 19 h$`{}_{}{}^{1}{}_{50}{}^{}`$ kpc. The general interpretation for this bi-modality is that Cen 45 is being accreted by the Centaurus cluster. This explanation has been further strengthened by: 1)the discovery of higher gas temperatures ($`>`$5 keV) associated with Cen 45 by Churazov et al. (1999) and 2) a marginally significant intracluster gas velocity difference of $``$ 1800 km s-1 between Cen 45 and Cen 30 detected by Dupke & Bregman (2001b), hereafter DB01.

More interestingly, DB01 found a small-scale (within the central 10) significant ($`>`$99.8% confidence) intracluster gas velocity gradient more or less symmetric with respect to the cluster’s center, consistent with bulk gas rotation with a correspondent circular velocity of 1.59$`\pm `$0.32 $`\times `$ 103 km s-1. The nature of this velocity gradient is likely to be related to some previous merger on the main body of the Centaurus cluster (Cen 30). The maximum velocity difference is found roughly along an axis at a position angle of $``$ 40 (Figure 1). Therefore, velocity measurements by means of the kinetic S-Z effect on CMBR temperature and intensity variations could be optimized along that axis (e.g. by choosing symmetric subregions for source/background along that position angle, such as the extremes of the box in Figure 1).

In Figure 2 we show the absolute and relative variations of intensity and CMBR temperature towards the directions of maximum (dotted) and minimum (dashed) velocities as observed by ASCA, i.e., 5 away from the cluster’s core (regions P3 & P7 in Figure 1), as a function of wavelength. Both kinetic and thermal effects are added. The gas density profile used was obtained taking into account that we are looking at spatial regions with the same scale as (or smaller than) the cluster’s core radius. Therefore, it is more appropriate to choose surface brightness profile fittings that realistically take into account the central cluster region, where the surface brightness is enhanced and single $`\beta `$ models diverge significantly from the observed profiles. We believe that the double $`\beta `$ fitting profiles used in the extensive ROSAT analysis of Mohr, Mathiesen & Evrard (1999) (hereafter MME99) are the most reliable for that purpose and we use them in our calculations666The choice of MME99 is also a conservative one, since other surface-brightness-derived density profiles available in the literature, e.g. Jones & Forman 1999, typically produce higher values for the optical depth within the spatial regions considered, thus enhancing the magnitude of $`\frac{\mathrm{\Delta }T}{T}`$.

MME99 surface brightness fitting function $`S_X(b)=_{i=1}^2S_{X0_i}(1+(\frac{b}{r_{c_i}})^2)^{3\beta +\frac{1}{2}}`$, assumed a constant $`\beta `$ but different core radii and normalizations for the dominant components near the center and away from the central regions. To deproject that surface brightness fitting and obtain a “physical” density profile MME99 expressed the total gas density as being proportional to the primary density component multiplied by a “fudge” function ($`f(r)`$) that compensated for the distortions caused by the central component. MME99 constrained $`f(r)`$ using the observed surface brightness for each radial bin. Following their methodology we use for our $`\tau (b)`$ the following integral

$$\tau (b)=2\sigma _T\frac{n_{e0}}{f(0)}_b^{\mathrm{}}(1+(\frac{r}{r_{c1}})^2)^{\frac{3}{2}\beta }\frac{f(r)rdr}{\sqrt{r^2b^2}}$$
()
instead of Equation 6. In Equation 7 $`n_{e0}`$ and $`f(0)`$ are the number density and the “fudge” function ($`f(r)`$) at the cluster’s center (MME99). We use for the functional form of $`f(r)`$ an inverse polynomial given as $`f(r>1^{})const0+\frac{const1}{b}+\frac{const2}{b^2}`$, where the best fit values within b$`<`$190 h$`{}_{}{}^{1}{}_{50}{}^{}`$ kpc are $`const0=0.817`$, $`const1=1.05\times 10^{23}cm`$ & $`const2=2.7\times 10^{46}cm^2`$.

The surface brightness parameters used to derive the optical depth (Equation 7) are Te(r=5) = 3.4 keV, Vcirc(r=5)=1600 km s-1 (DB01), $`\beta `$=0.57, $`n_{e0}`$=0.068 cm-3, rc1=7.3 (MME99). In all three plots the solid thick line shows the difference between the redshifted and blueshifted curves. It can be seen that the CMBR intensity variation difference (Figure 2, Top) is maximized at 214 GHz ($`\lambda `$=0.14 cm). This is the optimal frequency to observe the relative kinetic S-Z effect and consequently to determine a velocity map. At that frequency the difference between the relative variations of intensity and temperature are $`\mathrm{\Delta }(\frac{\mathrm{\Delta }I}{I})`$ 5 $`\times `$ 10-4 and $`\mathrm{\Delta }(\frac{\mathrm{\Delta }T}{T})10^4`$ (Figure 2 middle and bottom plots).

The magnitude of the parameters derived above are, naturally, a function of the projected radial distance from the cluster’s center. In order to determine the radial profile of the expected $`\frac{\mathrm{\Delta }T}{T}`$ at 214 GHz we make two assumptions; Firstly, we assume that in the region from 0 to 5 the cooling flow can be approximated by a linear function varying from 3.4 keV to 1 keV at the center. The flattening of the temperature profile at the cluster’s center has been observed for some clusters with Chandra and XMM, and may be a common feature of cooling flows (Fabian 2001). In the case of Centaurus there is evidence that the minimum temperature derived from spectral fittings using cooling flow models is $`>`$ 0.4-1.0 keV (Sanders and Fabian 2001). The central arcmin in Centaurus has an impressive amount of substructures (plumes) and this makes our estimation of the central densities more uncertain. Therefore, we exclude that region from our analysis (Figure 3). Secondly, we assume the gas motion can be approximated by a solid body rotation. This functional form for the gas velocity cannot be supported for larger radii since the gas becomes gravitationally unbound but it is a good observed first order approximation for the inner regions (Dupke & Bregman 2002). Actually, this assumption is supposed to break at $``$ 5-6, when the combined bulk kinetic and thermal energies become greater than the gravitational and the gas becomes unbound (e.g. see Allen & Fabian 1994). Assuming that the velocity “freezes” at 1600 km s-1, we also show in Figure 3 the predicted behavior of $`\frac{\mathrm{\Delta }T}{T}`$ at regions radially greater than that encompassed by ASCA velocity analysis, for illustration purposes.

The kinetic Comptonization of the spectrum that reaches us from a rotating cluster is roughly antisymmetric with respect to the center of the cluster from $``$ 160-270 GHz. This suggests that larger variations of intensity and temperature of the CMBR can be observationally enhanced by comparing the measurements (in the same channel) “within” the cluster, rather than to a background region away from the cluster. The results are maximized at 214 GHz and we show the distribution of CMB temperature variations ($`\frac{\mathrm{\Delta }T}{T}(b)`$) at that frequency in Figure 3. The horizontal errors in Figure 3 show the FWHM of BOLOCAM. The solid lines represent the estimated 1-$`\sigma `$ errors of measurement associated to the X-ray-derived parameters, and are dominated by X-ray velocity measurement errors. S-Z measurements of radial velocities can be performed within these 1-$`\sigma `$ limits with a relatively short exposure time (see below).

Variations of intensity or temperature of that magnitude can be relatively easily observed currently, or in the very near future, with the new generation of bolometers. A total variation of CMB temperature of $``$ 0.3 mK between the “knees”in Figure 3 is equivalent to a flux density of $`\mathrm{\Delta }I_{lim}`$1.2 mJy (FWHM)-2, This could be done with high significance in a relatively short amount of time with, for example, BOLOCAM, given that for that instrument the FWHM is 43′′ and the expected777Samantha Edgington and the Bolocam Team, Personal Comunication noise equivalent flux density (NEFD) is 35 mJy Hz$`^{\frac{1}{2}}`$ at the 1.4mm band888The nominal observation time necessary to detect a limiting flux density $`I_{lim}`$ is $`t_{obs}(\frac{NEFD}{\mathrm{\Delta }I_{lim}})^2`$ sec plus overhead(Glenn et al. 1998).

## 4  Summary

Direct indications of the presence of large bulk motions in clusters of galaxies have been recently obtained through X-ray spatially resolved spectroscopy. However, these observations are limited by detector’s gain fluctuations and this will continue to be the case for the next several years. An alternative way of measuring intracluster gas bulk velocities is through the use of the kinetic S-Z effect. The use of bolometers has improved significantly the precision of radio observations of the S-Z effect, especially within the optimal frequency range for such measurements, i.e., $``$ 214 GHz ($`\lambda `$ = 0.14 cm), so that, it will be possible in the near future to perform such velocity measurements in instruments such as BOLOCAM (Bock et al 1998) with relatively short exposure times. The combination of X-ray and radio measurements can provide us with independent/complementary information crucial to build ICM velocity maps and to determine the degree of evolution and the history of galaxy clusters.

So far, the best candidate for testing the proposed multi-wavelength measurements of velocity gradients is the Centaurus cluster. The symmetry of the gas velocity distribution and the low temperatures measured for this cluster should facilitate the separation of the rotational kinetic component from both thermal and peculiar velocity ones. By combining the analyses of intracluster gas velocity & temperature distributions with other relevant X-ray cluster characteristics we predict that the variations of CMBR temperature between two regions symmetrically located 5 from the clusters center along the line of maximal velocity gradient at $`\lambda `$ = 0.14 cm is $``$0.3 mK.

The magnitude of the S-Z parameters observationally derived in this work, if typical, could significantly enhance possible halo-induced gas rotation contribution to the CMBR angular power spectrum (Cooray & Chen 2002). However, the limiting nature of the search for velocity gradients using current X-ray spectrometers tends to be biased towards clusters with high bulk velocities, so that the current data does not allow us yet to calculate precise velocity distributions and, consequently, to determine the deviations from the assumptions used in statistical analysis of the temperature anisotropies such as those performed by Cooray & Chen (2002). Furthermore, if both “rotational” and “transient” bulk velocities are as common in the central regions of clusters as suggested by a more extended ASCA analysis (Dupke & Bregman 2002), measurements of cluster peculiar velocities using the S-Z effect should be treated with extreme caution, given the potential confusion with an internal gas bulk motion component.

Multiple off-center long exposure velocity measurements of intracluster gas with the Chandra and XMM-Newton satellites coupled with radio observations of the differential kinetic S-Z effect will be crucial to determine precise intracluster velocities in the central and intermediate regions of galaxy clusters, improving significantly the determination of their evolutionary stage, and constraining biases related to the presence of coherent velocity fields in blind S-Z surveys.

We are very grateful to Frits Paerels and Gary Bernstein for the helpful discussions and constructive suggestions. We would like to thank Joe Mohr for kindly providing the tabulated distribution of the ”fudge” function that was used in this paper. We also thank Samantha Edgington and the Caltech BOLOCAM team for helpful suggestions. We thank Hugh Aller & Tim Paglione for helpful discussions. We acknowledge support from NASA Grant NAG 5-3247. This research made use of the HEASARC ASCA database and NED.