text:
# M 33: A Galaxy with No Supermassive Black HoleBased on observations with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. (AURA), under NASA contract NAS5-26555.

## 1. Introduction

The mass of a central black hole (BH) and the stellar orbital structure in a galaxy are direct probes into the formation and evolution history of both the galaxy and the BH. The tight correlation between BH mass $`M_{}`$ and bulge velocity dispersion $`\sigma _e`$ (Gebhardt et al. 2000a; Ferrarese & Merritt 2000) establishes a strong link between BHs and their hosts. Many theories (e. g., Silk & Rees 1998; Haehnelt & Kauffmann 2000; Ostriker 2000; Adams, Graff, & Richstone 2001) predict such a correlation. As data and theory improve, we hope to be able to discriminate between the alternatives.

The $`M_{}`$ – $`\sigma _e`$ correlation is derived using galaxies whose bulge-to-total luminosity ratios $`B/T`$ range from $``$ 0.15 to 1. They even include galaxies with “pseudobulges,” i. e., bulge-like, high-density central components that are believed to have grown by the secular evolution of disks (see Kormendy 1993 for a review). Bulges and pseudobulges appear to satisfy the same $`M_{}`$ – $`\sigma _e`$ correlation (Kormendy & Gebhardt 2001; Kormendy et al. 2001). However, the above references show that BHs do not correlate with the total gravitational potential of the disk. It therefore appears that the main requirement for a supermassive BH is that a galaxy contain some kind of bulge, or at least a mass gradient similar to one. Still, pure disk galaxies have not been studied in nearly the same detail as galaxies that contain bulges.

Observing pure disk galaxies is difficult, because they generally do not have bright centers and because the small expected BH masses require good spectral resolution. The presence of nuclei – compact star clusters that are distinct from the disk – have made it possible to study a few bulgeless galaxies. For example, based on its nuclear velocity dispersion, Filippenko & Ho (2001, see also Moran et al. 1999) conclude that the dwarf Seyfert galaxy NGC 4395 has $`M_{}`$ $`\genfrac{}{}{0pt}{}{_<}{^{}}`$ 80,000 $`M_{}`$. However, a BH of only 100 $`M_{}`$ is needed if the nucleus is radiating at the Eddington luminosity. Therefore NGC 4395 does not put a significant strain on the $`M_{}`$ – $`\sigma _e`$ correlation.

M 33 is the best bulgeless galaxy to use to investigate the relationship between BHs and disks. It is a normal, large galaxy (absolute magnitude $`M_V=19.0`$), but it is only 0.8 Mpc away (van den Bergh 1991). The spatial scale is therefore unusually good, 0.′′26 pc-1. M 33 does not contain a bulge, but it does have a bright nucleus. Kormendy & McClure (1993) measured its velocity dispersion with the Canada-France-Hawaii Telescope and found that $`\sigma =21\pm 3`$ km s-1. This implies an upper limit on the mass of any central black hole of $`M_{}`$ $`\genfrac{}{}{0pt}{}{_<}{^{}}`$$`5\times 10^4`$ $`M_{}`$. M 33 was the first large galaxy in which a dead quasar could strongly be ruled out. Recently, Lauer et al. (1998) obtained WFPC2 photometry and suggested that $`M_{}`$ $`\genfrac{}{}{0pt}{}{_<}{^{}}`$$`2\times 10^4`$ $`M_{}`$, by reducing the upper limits for the core radius of the nucleus. In this paper, we present Hubble Space Telescope (HST) STIS spectroscopy of the nucleus obtained with the 0.′′1 slit. The spatial resolution is an order of magnitude better than that in Kormendy & McClure (1993). This leads to a much stronger limit on $`M_{}`$. Our best three-integral dynamical model has $`M_{}`$ = 0. Even our upper limit of 1500 $`M_{}`$ is significantly below the value expected from the $`M_{}`$ – $`\sigma _e`$ correlation.

It is important to emphasize how much M 33 differs from the galaxies in which BHs have been found. Figure 1 shows its brightness profile and Figure 1 in Kormendy & McClure (1993) shows an image of the central 70′′ $`\times `$ 113′′. There is no bulge; the exponential disk dominates the structure in to 3′′. At the center, there is a compact nuclear star cluster that is dynamically distinct from the disk. At absolute magnitude $`M_B=10.2`$, this is like a large globular cluster; it appears to be more closely related to globular clusters than to bulges in the Fundamental Plane parameter correlations (Kormendy & McClure 1993). The nucleus has a high central stellar density, similar to that of M 31 and M 32 (Kormendy & Richstone 1995; Lauer et al. 1998). But its velocity disperson is like that of a globular cluster, not like that of a bulge. Since mass $`M\sigma ^2`$, the difference in velocity dispersion is the essential reason why we find BHs in M 31 and M 32 but not in M 33.

## 2. Photometry

We use HST WFPC2 photometry from Lauer et al. (1998), who provide details of their reduction procedure. The bandpasses are $`V`$ and $`I`$, and the spatial resolution of the final profiles is 0.′′025. Ground-based data at larger radii come from Kormendy & McClure (1993). There is a color gradient in the nucleus; its center is $`\mathrm{\Delta }(VI)=0.15`$ bluer than its outer parts (Fig. 2). We want to measure the mass profile by using the light profile to represent the stellar potential, so we need to check how much the color gradient affects the results. It would be safest to use a $`K`$-band profile, because this would be least sensitive to the gradient in stellar population. Since we don’t have a $`K`$-band profile at the desired spatial resolution, we construct one by using the $`VI`$ profile to estimate $`VK`$. Giant stars have a monotonic and smooth relationship between $`VI`$ and $`VK`$ (Alves 2000, and references therein). At $`VI1`$ (appropriate for M 33), the relation is approximately linear, $`\mathrm{\Delta }`$$`(VK)=3\mathrm{\Delta }`$$`(VI)`$. This relation provides the approximate $`K`$-band profile shown in Fig. 2. To further check the effect that color gradients have on the results, we have also run models where we change the mass-to-light ratio as a function of radius according to the $`VI`$ profile; for these models we include a variation of 20% in the M/LI starting from the edge of the nucleus decreasing linearly into the center. By either including a decrease in the mass-to-light ratio towards the center or using a flatter surface brightness profile (as suggested by the K-band profile), we are increasing the likelihood of requiring a black hole. Our goal is to place a robust and conservative upper limit on the black hole mass. Dynamical models constructed using the measured $`V`$, measured $`I`$, extrapolated $`K`$ profiles, and $`I`$ with varying M/LI profiles all give very similar results. We conclude that the stellar population gradient does not have a significant affect on our estimates of $`M_{}`$. We discuss the results from each of the profiles below.

## 3. HST STIS Kinematics

We measured the kinematics using the Space Telescope Imaging Spectrograph (STIS; Woodgate et al. 1998) centered at 8561 Å. Six exposures at three different pointings had a total exposure time of 2.05 hours. The 0.′′1 $`\times `$$`52\mathrm{}`$ slit and G750M grating provided pixels of 0.′′05 $`\times `$ 0.55 Å = 19 $`\mathrm{km}\mathrm{s}^1`$. The nucleus of M 33 is slightly elongated (axial ratio $``$ 0.85, Lauer et al. 1998); we put the slit along the major axis. Contemporaneous flats were used to remove pixel-to-pixel sensitivity variations and fringing. We did not use the weekly dark frames provided by STScI, because hot pixels come and go on timescales shorter than a day. Instead, we created a dark frame from the data using the iterative procedure described in Pinkney et al. (2001). The M 33 procedure differed slightly from normal: we have only 3 dither positions instead of 5, and the dither spacing is 4 pixels instead of 20. The steepness of the nuclear brightness profile alleviates the light overlap problem caused by the small dither spacing. It takes 5 iterations to bring the residual galaxy light level below 1 % of the initial value, i. e., acceptably below the noise. This process was used to create the dark frame that was subtracted from the galaxy spectra.

The STIS spectra are well sampled in the spatial direction, so it was possible to combine the dithered spectra without aliasing or the introduction of any artifacts or smoothing that would degrade the spatial resolution (Lauer 1999). We also choose to resample the spectra in the spatial direction to provide the greatest flexibility in sampling the kinematic parameters as a function of radius. Combination of the spectra began with shifting the spectra to a common spatial centroid to enable identification and elimination of cosmic ray events. Cosmic ray events in one spectrum were replaced with an average from the remaining spectra. The “cleaned” spectra were then shifted to a common centroid using sinc function interpolation. This allows well-sampled functions to be reinterpolated without the introduction of smoothing. The final composite spectrum was sampled with 0.′′025 pixels. We also used a cruder combination of the data that set the pixel scale at 0.′′05 pixels; in that case we find an upper limit on the black mass that is a factor of two larger. Thus, the resampling procedure is critical as it allowed superior use of the spatial information present in the data.

The galaxy spectrum was radially binned according to the optimal binning scheme used in the modeling. That is, the observations and models have the same radial bins. We use nearly linear bins at small radii and logarithmic bins at large radii (see Richstone et al. 2001 for details). The reduction procedure discussed in Pinkney et al. (2001) was then applied to the extracted spectra. Each spectrum was divided by an estimate of its continuum. This was based on binned windows from which we measured a robust mean (the biweight) using the highest 1/3 of the points. We then linearly interpolated between the mean values. The continuum division removes the need to flux calibrate the spectra.

Obtaining the internal kinematic information requires a deconvolution of the observed galaxy spectrum with a representative set of template stellar spectra. Both the deconvolution procedure and the choice of templates stars affect the derived velocity profiles. For the template, we use the K giant star HR7615. M 33 has younger stars near its center than farther out in the nucleus, so using using only a K III standard star may not be appropriate. In fact, M 33’s nucleus clearly contains A stars (Gordon et al. 1999, and references therein). Fortunately, the Ca triplet region is not as sensitive to template variations as is the more traditional Mg spectral region (Dressler 1984). M 33 shows a small radial variation in Ca triplet equivalent width from 10 Å at the center to 8 Å at 0.′′5. Our template star, HR7615, has a Ca triplet equivalent width of 8 Å. García-Vargas et al. (1998) show that such a small difference in equivalent width has a small effect on age and metalicity. In Section 5, we discuss the possibility of this equivalent width difference affecting the deduced kinematics.

For the deconvolution, we reduce each spectrum using a maximum-penalized-likelihood (MPL) technique that produces a non-parametric line-of-sight velocity distribution (LOSVD). First, we guess initial velocity profiles in bins. We convolve these profiles with the template and calculate the residuals from the galaxy spectra. The program then varies the velocity profile parameters – i. e., the bin heights – to provide the best match to the galaxy spectrum. The MPL technique is similar to that used in Saha & Williams (1994) and in Gebhardt et al. (2000b). The complete set of spectra and velocity profiles are shown in Figure 3.

Our models are axisymmetric, so at any absolute radius $`|r|`$, they have the same LOSVDs on opposite sides of the center, except for the sign of the velocity. Thus, in order to maximize the signal-to-noise ratios in our observations, we fit the same velocity profile shape to both sides of the galaxy. In Figure 3, we plot these velocity profiles together with the two spectra at $`\pm |r|`$. We have also measured the kinematics without symmetrizing. As expected, we obtain similar results but with larger uncertainties. Since we are trying to extract maximum information given low signal-to-noise levels in some spectra, the symmetrized versions are preferred.

Monte Carlo simulations are used to measure the uncertainties in the velocity profiles. We convolve the template star with the measured LOSVD to fit to the galaxy spectrum as illustrated by the red line in Figure 3. From this synthetic galaxy spectrum, we then generate 100 simulated observed spectra with noise, and we determine their velocity profiles. Each simulated spectrum contains random Gaussian noise with the standard deviation given by the RMS of the initial fit. The 100 simulated spectra and reductions provide a distribution of LOSVDs from which we estimate the confidence bands. The dashed lines in the right panels of Figure 3 represent the 68 % confidence bands of the distributions. This bootstrap techniques estimates both the random error and the bias of the estimator; however for these profiles the bias is negligible. We do not test for the bias due to template mismatch, but only estimator bias. The modeling described below uses these velocity profiles directly; we do not use a parametric representation of them. However, for convenience, Table 1 provides parameters and uncertainties of Gauss-Hermite fits to the velocity profiles.

Gebhardt et al. (2000a) use the luminosity-weighted velocity dispersion integrated out to the half-light radius of the bulge to define the “effective dispersion”. Here, we follow the same procedure for the nucleus. Integrating the STIS spectra from the center to the edge of the nucleus, we find an effective velocity dispersion of $`\sigma _e=24.0\pm 1.2`$ $`\mathrm{km}\mathrm{s}^1`$.

## 4. Three-Integral Models

The 3-integral models are the same as in Gebhardt et al. (2000b) and are described in detail in Richstone et al. (2001). The models are orbit-based (Schwarzschild 1979) and constructed using maximum entropy (Richstone & Tremaine 1988). We run a representative set of orbits in a specified potential and determine the non-negative weights of the orbits to best fit the available data. The library of orbits is convolved with the STIS PSF (as in Bower et al. 2001) and binned according to the STIS pixels. Maximizing the entropy helps to provide a smooth phase-space distribution function. Rix et al. (1997) present a similar code which has been applied to the spherical case; Cretton et al. (1999) and van der Marel et al. (1998) developed a fully-general axisymmetric code that is also orbit-based.

Because of M 33’s small dispersion and size we did not use the standard setup parameters as in Pinkney et al. (2001). Instead, we designed the model bins to include the high-resolution sampling that we used for the kinematic measurements. The central bins are 0.′′025 wide.

The model consists of 15 radial, 5 angular, and 13 velocity bins. However, the force calculations use 4 times finer spatial sampling in order to achieve energy conservation to better than 1% for most stars. On average, each model contains about 3500 stellar orbits. We tried a variety of inclinations for the nucleus from edge-on to nearly face-on. In all cases we find a similar upper limit for the black hole mass. Zaritsky, Elston, & Hill (1989) measure an inclination of 56° for the disk.

The data consist of six HST radial bins and two ground-based bins with 13 velocity bins each. However, the actual number of independent measurements is less than 13 $`\times `$ 8, for several reasons. First, the ground-based measurements add only six independent values: a velocity and a dispersion at each of two positions. Second, each STIS pixel is 19 $`\mathrm{km}\mathrm{s}^1`$. Given the 120 $`\mathrm{km}\mathrm{s}^1`$ region where the velocity profile is non-zero, this implies about 6 independent measurements. Third, the smoothing length in the LOSVD estimate correlates some of the points. Thus, we have roughly 5 independent measurements per the six HST velocity profiles combined with the 4 ground-based measurements, giving about 34 degrees of freedom.

## 5. Results

For each of the three light profiles, we ran nine models with BH masses ranging from zero to $`1\times 10^4M_{}`$. The total $`\chi ^2`$ values are shown in Figure 4 as a function of BH mass using the $`I`$-band light profile. For each BH mass, the best-fitting constant stellar $`M/L`$ ratio is 0.70 in $`I`$ band. Figure 4 only shows the one-dimensional $`\chi ^2`$ values, and we do not plot these as a function of $`M/L`$ since there is little information there. The minimum $`\chi ^2`$ value is around 27, which makes the reduced $`\chi ^2`$ near unity.

Zero BH mass provides the best fit. Figure 4 only shows the results for the $`I`$ band where we have allowed the M/LI to vary by 20%, but using any of the other light profiles also results in a minimum at zero BH mass. The difference between the various light profiles is in the upper limit for the black hole mass. For the $`I`$ band profile with varying M/LI, $`\mathrm{\Delta }\chi ^2=1`$ at 1400 $`M_{}`$; for $`I`$ with constant M/LI, the upper limit is 600 $`M_{}`$; using the $`K`$ profile, the upper limit occurs at 1200 $`M_{}`$. We generally use a $`\mathrm{\Delta }\chi ^2`$ equal to one to signify the extent of the viable models, but, in this case, given the uncertainty of the light profile we choose a conservative upper limit of 1500 $`M_{}`$.

Figure 5 plots the first four Gauss-Hermite moments ($`h_3`$ and $`h_4`$ are measures of skewness and tail weight) of the velocity distributions for both model and data. There are two different BH models plotted in Figure 5, yet one would be very hard pressed to discern any difference between them. In fact, the difference in $`\chi ^2`$ between the two models shown in Fig. 5 is 20, which significantly excludes the large black hole mass model. However, we do not use these Gauss-Hermite values during the fitting since the full LOSVD determines the $`\chi ^2`$ minimization. Thus, it is not optimal to use these moments as a figure of merit to judge the quality of the fit. Instead, one should inspect the comparison with the velocity profiles. Figure 6 plots the LOSVD for three STIS bins with the same three BH models. The difference between the various models can now be seen. Furthermore, inspection of the other STIS velocity profiles reveal similar results—i.e., the model with no BH consistently provides the best fit. Only by looking at the total $`\chi ^2`$ (Figure 4) can one judge the model fits.

The 3-Integral models also provide the stellar orbital distribution. In order to constrain adequately the stellar orbits, we must have two-dimensional spatial kinematic information. Major axis data allow for too large a parameter space for the models. Based on comparison with analytic models, the maximum entropy solution using only major kinematic axis data does not recover the input orbital distribution (Richstone et al. 2001). Having only one additional position angle helps considerably. Unfortunately, we only have major axis data for M 33 and so do not report the orbital structure. However, with only the major axis data, the estimate for the BH mass is unbiased.

### 5.1. Possible Problems

These 3-Integral models are limited by various assumptions. First, we assume that the nucleus is spheroidal (elliptical isophotes with ellipticity constant with radius). Second, the models are axisymmetric, so we do not consider triaxial shapes. Third, we assume specific radial forms for the $`M/L`$ variation. Any of these assumptions will reflect themselves as a bias in the models and restrict the measured uncertainties. The most important is likely the assumption of the $`M/L`$ variation. However, the color gradient is rather small ($`\mathrm{\Delta }VI=0.15`$) and when we include a $`M/L`$ gradient, the results do not change significantly. This is apparent since we find similar results when we use the extrapolated $`K`$-band light profile.

No model can reproduce the $`h_4`$ values seen in the central few bins. We have thoroughly checked the negative $`h_4`$ values in these bins using extensive simulations to confirm that they are not an artifact of either the data analysis or velocity profile estimation. They appear to be a real feature in the data that is not reproducible in the models. The absolute difference between the model $`h_4`$ of $`0.02`$ and the data $`h_4`$ of $`0.04`$ is small but significant given the uncertainty. This difference may either be due to use of an improper template or to our lack of knowledge of the light distribution interior to our spatial resolution. Since M 33’s light profile is bluer in the center, template mismatch is likely since we are only using a K giant star. An A star will likely have a lower equivalent width for the Ca triplet, depending on the metalicity (García-Vargas et al. 1998); thus, if the center of M 33 consists of A stars, a fit using a higher equivalent width may result in an LOSVD with truncated wings, such as is seen. However, the absolute difference is so small (0.02) that it will have a negligible effect on our results.

Gordon et al. (1999) advocate a strong foreground of absorbing dust (V-band optical depth of 2), and we are not using an observed K-band light profile. Instead, we have estimated the K light using the V and I profiles. If the dust is distributed with a steep radial profile, then our K-band estimation will be in error. However, the effect is to make the upper black hole mass limit even tighter. As one includes more dust near the center, the result is to steepen the actual K-band profile and thereby increasing the amount of mass in the stellar component. Thus, unless the dust is distributed in a clever way, the upper limit that we derive on the black hole mass is robust.

We have assumed that the dynamical center of the galaxy is at the nucleus. Photometrically, the center of the galaxy is difficult to identify given the amount of dust; high spatial resolution IR imaging is required to determine the photometric center. X-ray observations (Colbert & Mushotzky 1999) show a compact source within 2″ of the nucleus, which is well within the ROSAT HRI uncertainty; but given that M33 does not appear to have a supermassive black hole, x-ray observations are unlikely to point to the galaxy center. Kinematically, the rotation curve is symmetric about the nucleus giving weight that we are looking at the galaxy center. However, the nature of the nucleus is unclear; it is similar to a large core-collapse globular cluster in both photometric and kinematic properties. Given the shallow potential well in the center of M33, such a dense system would be free to wander. Thus, even though there is no evidence to suggest a problem, a two-dimensional kinematic map of the nuclear region would be very useful to pinpoint the dynamical center.

## 6. Conclusions

The nuclear effective velocity dispersion of M 33 is $``$ 24 $`\mathrm{km}\mathrm{s}^1`$. If the mass of a central BH were related to this dispersion of a nucleus like the known BHs are related to the dispersions of their bulges; i. e., if $`M_{}=1.3\times 10^8(\sigma _e/200)^{3.65}`$ (Gebhardt et al. 2001), then we would expect a BH of $`M_{}=5.6\times 10^4`$ $`M_{}`$. The corresponding relationship from Ferrarese & Merritt (2000) implies a mass of 7000 $`M_{}`$. The upper limit derived from the HST data is 1500 $`M_{}`$ and the best-fit value is zero. Thus, M 33 appears to lie below the $`M_{}`$ – $`\sigma _e`$ correlation for bulges. This suggests that pure disk galaxies—even ones with nuclei—are different from bulges and ellipticals in their BH properties. Figure 7 plots M33’s upper limit in the $`M_{}`$ – $`\sigma _e`$ correlation. The line is given by the above relation and the 68% confidence bands are determined through a bootstrap procedure (Gebhardt et al. 2001).

There is a caveat. Because there are no firm dynamical BH detections in the mass range $`10^3`$ to $`10^6`$ $`M_{}`$, there is a big gap in the $`M_{}`$ – $`\sigma _e`$ correlation between the smallest known supermassive BHs and the mass limit in M 33. That is, we do not know whether the correlation extends down to the small masses of interest here; this is not a problem: our observation may be a sign that it does not. But given the large gap, we also do not know that the correlation—if it exists—is linear between $`10^6`$ and $`10^3`$ $`M_{}`$. If the correlation exists but is not linear, then our limit of $`M_{}`$ $`\genfrac{}{}{0pt}{}{_<}{^{}}`$ 1500 $`M_{}`$ may not be below the correlation. Furthermore, the confidence band in Figure 7 assumes a linear relation; thus, deviations from a linear relation will have a significant affect on the estimated uncertainties.

On the other hand, Kormendy et al. (2001) show that our $`M_{}`$ limit in M 33—and published masses and mass limits on BHs in other disk galaxies—are far below the correlation of $`M_{}`$ with the total luminosities of their host galaxies. M 33 is by far the most extreme case. There is a corresponding argument for the $`M_{}`$ – $`\sigma _e`$ correlation. The H I rotation curve of the outer disk of M 33 rises slowly from about 100 $`\mathrm{km}\mathrm{s}^1`$ at $`r=15^{}`$ = 3.5 kpc to about 135 $`\mathrm{km}\mathrm{s}^1`$ at $`r=77^{}`$ = 18 kpc (Corbelli & Salucci 2000). If these are circular velocities in the gravitational potential of an almost isothermal dark halo, then the corresponding one-dimensional velocity dispersion is $`\sigma 85\pm 10`$ $`\mathrm{km}\mathrm{s}^1`$. If a BH in M 33 were related to the dark matter potential well like known BHs are related to their bulge potential wells (both represented by $`\sigma _e`$), then M 33 should contain a BH of mass $`M_{}6\times 10^6`$ $`M_{}`$. Clearly it does not. If M 33 is typical of disk galaxies, then BH masses are not related in the same way to every gravitational potential well (as measured by $`\sigma _e`$). In particular, they are not related to the potential well of the disk and dark matter. Taken together, our results and those of Kormendy et al. (2001) strongly suggest that BH masses are causally connected with the properties—luminosity, mass, and especially velocity dispersion—of their host bulges and not with their disks.

These results are still based on only a few, well-observed galaxies. We therefore need to study the $`M_{}`$ – $`\sigma _e`$ correlation more fully at the low-dispersion end. This should provide further insight into the fundamental differences between disk and bulge galaxies, and it may help to tell us why such a correlation exists in the first place.

K.G. was partially supported by NASA through Hubble Fellowship grant HF-01090.01-97A awarded by the Space Telescope Science Institute. This work was also supported by NASA Guaranteed Time Observer funding to the STIS Science Team.


format: text

text:
# Scaling in the growth of geographically subdivided populations: invariant patterns from a continent-wide biological surveyIn press, Philosophical Transactions of the Royal Society of London, Series B.

## 1 Introduction

Perhaps one of the most intriguing patterns in ecology is Taylor’s law (Anderson et al., 1982; Soberón and Loevinsohn, 1987; Routledge and Swartz, 1991; Leps̆, 1993; Curnutt et al., 1996; Maurer, 1999). Taylor (1961) was the first to notice that when the mean $`S`$ of a population survey is plotted versus its variance $`\sigma ^2(S)`$, either in space or time, the relationship is typically a power-law with a fractional exponent

$$\sigma ^2(S)S^\gamma .$$
(1)
Taylor was originally interested in the slope of the power-law relationship as a scale-free measure of spatial contagion or dispersion — values greater (less) than one indicate spatial clustering (over-dispersion). Later, Taylor used both spatial and temporal scaling as a basis for comparative studies of, in his words, “synoptic population dynamics” across taxonomic groups (Taylor and Woiwod, 1982; Taylor, 1984).

Taylor’s synoptic approach is, in many respects, a precursor to the recent development of “macroecology,” (Brown, 1995), a sub-discipline of ecology and biogeography (MacArthur and Wilson, 1967) that seeks to identify broad patterns in species’ abundance and distribution. Macroecology has largely focussed on static patterns, such as spatial relationships between abundance and environmental factors (Brown et al., 1995) and relationships between metabolic energy use and geographic distribution (Brown and Maurer, 1987). Thus, relatively few continental-scale macroecological studies (Maurer, 1994, 1999) have explored interactions between spatial distribution and population variability through time.

In this paper, we adopt Taylor’s synoptic approach and analyze one of the most comprehensive macroecological data sets available, the North American Breeding Bird Survey (Peterjohn, 1994). The data are estimates of local abundance (counts) for over 600 bird species recorded at 2–3 thousand sites (routes) across North America for the years 1966 to 1998. Unlike many previous macroecological studies, our focus is on linking geographic distribution to population dynamics. We report two new scaling laws, closely related to Taylor’s power-law, for these data: one relating variability of population time series to their mean, and another relating number of sites occupied to total population size. In addition, we show how these patterns may be related, and discuss a simple stochastic growth model for a geographically subdivided population that formalizes the relationship.

## 2 Scaling of species growth rates

Our goal is to understand temporal variation in abundance at the entire population level. We therefore compute time series of total counts for each species by summing over all routes surveyed in a given year (Fig, 1a). (We have previously analyzed these data at the individual route level, see Keitt and Stanley, 1998.) The aggregated time series should be relatively robust to observer errors inherent in the route-level data (Kendall et al., 1996) since random under-counts and over-counts will cancel in the summation. After performing a bias correction (see Appendix for details), the resulting time series are relatively free of systematic trends, particularly in the early years of the survey when the number of routes was increasing rapidly (compare dashed to solid lines in Fig. 1a).

We choose as our measure of the magnitude of time-series fluctuations, the logarithm of the ratio of successive counts,

$$g(t)\mathrm{log}\left(\frac{S(t+1)}{S(t)}\right),$$
(2)
where $`S(t)`$ and $`S(t+1)`$ are the total number of birds of a given species counted in years $`t`$ and $`t+1`$, respectively. This measure has several nice properties. First, any multiplicative, time-independent sample bias cancels in the ratio. Second, the measure has a natural interpretation in terms of population demography since, in a closed population, $`\mathrm{exp}[g(t)]1`$ \+ (per capita birth rate - per capita death rate).

As shown in Fig. 1b, the standard deviation $`\sigma (g)`$ of population growth rates is strongly related to the average total population size. The relationship follows a power-law

$$\sigma (g)S^\beta $$
(3)
for over four orders of magnitude in $`S`$, the total count averaged across years. For these analyses, we are not interested in predicting a “dependent” variable from an “independent” variable. Rather, we are interested in modeling the functional form of the interdependence among variables. We therefore use major-axis regression (Sokal and Rohlf, 1995) to estimate model parameters. Major-axis regression is based on computing the leading eigenvector of the covariance matrix, and minimizes squared errors measured perpendicular to the trend line. We also restrict our analysis to non-zero time series of at least 25 years in length and with a minimum average total count of no fewer than 5 birds per year. Using major-axis regression with bootstrap precision estimates, we find $`\beta =0.36\pm 0.02`$. Taylor’s exponent (replicated across species) is simply $`\gamma =2(1\beta )=1.44\pm 0.04`$.

We also study the scaling properties of population growth rates by examining changes in the distribution of growth rates with increasing population size, a technique familiar to statistical physicists. First, we separate the observed growth rates into three bins according to the total count $`S(t)`$ and then construct a histogram to estimate the conditional probabilities $`p(g|S)`$. The resulting distributions are roughly triangular in shape with the width depending on $`S`$ (Fig. 1a). (The triangular shape may result from summing over a large number of time series with different local variances, see Amaral et al., 1998). If the distributions are “self-similar” (i.e., exhibit scaling), then we should be able to identify a function $`f`$ that rescales the distributions so that they “collapse” onto each other (Stanley, 1971). We plot the scaled quantities:

$$\sigma (S)p\left(\frac{g}{\sigma (S)}|S\right)vs.\frac{g}{\sigma (S)},$$
(4)
(Fig. 2a) and find that indeed the three curves do collapse onto each other (Fig. 2b), suggesting that $`p(g|S)`$ follows a universal scaling form

$$p(g|S)\frac{1}{\sigma (S)}f\left(\frac{g}{\sigma (S)}\right).$$
(5)
These results are interesting for a number of reasons. In statistical physics, the presence of non-trivial scaling is usually taken to mean that the dynamics are largely governed by simple geometric properties of the system and do not depend strongly on detailed properties of the system subcomponents (Wilson, 1983). Thus, it is remarkable that we should find strong evidence for scaling across such a taxonomically and ecologically diverse set of species as found in the North American Breeding Bird Survey. These results suggest that the dynamics of North American Breeding Birds are unexpectedly “simple” and depend primarily on common patterns of internal population structure across species ranges, rather than details of individual species life-histories.

## 3 Spatial structure of subpopulations

Another reason the growth-scaling results are interesting is the large variability of highly abundant species. Imagine the null model that each population is subdivided into $`n`$ equally sized, independent subpopulations, and that the number of these subpopulations depends linearly on $`S`$. The expectation, according to the central limit theorem, is that the standard deviation in growth rates should decay as the -1/2 power of $`S`$. The observed decrease in fluctuations is considerably slower (i.e., $`\beta <1/2`$), such that highly abundant species are considerably more variable than expected under the null model.

To account for the increased variability for highly abundant species, we require that the number of subpopulations does not scale in a simple, linear fashion with increasing $`S`$, but instead takes the form:

$$nS^{1\alpha },$$
(6)
with $`\alpha 0`$. Values of $`\alpha >0`$ will be found, for example, when the “typical” size of the subpopulations also scales with total abundance, i.e., there is a positive relationship between regional and local abundance, a well documented pattern in macroecology (Gaston and Lawton, 1988; Gaston, 1996). The positive correlation between local and regional abundance results in fewer subpopulations for a given total population size as each subpopulation accounts for more individuals. Again appealing to our observation that, under the central limit theorem, $`\sigma (g)n^{1/2}`$, and in combination with equations (3) and (6), it is strait forward to show that for roughly equal-sized subpopulations, the estimated exponents must obey

$$\beta =\frac{1\alpha }{2}.$$
(7)
We do not have access to a precise estimate of the number of subpopulations for each species in the survey. However, we can use as a proxy the number of survey routes where a species had a non-zero count in a given year. To test the assertion in Eq. 7, we plot number of survey routes with non-zero counts $`\stackrel{~}{n}(t)`$ versus the (uncorrected) total count $`S(t)`$ for all bird species recorded in the survey in 1997, excluding species seen at fewer than 5 routes or with fewer than 5 total individuals counted. The data follow closely the power law dependence predicted by Eq. (6) with an exponent $`\alpha =0.25\pm 0.03`$, again using major-axis regression with bootstrap precision estimates (Fig. 3a). Remarkably, the estimate of $`\alpha `$ predicts a value of $`\beta =0.38\pm 0.02`$, very close to the estimate ($`\beta =0.36\pm 0.02`$) obtained by measuring the standard deviation in growth rates directly (Fig. 1b). Even more striking is the consistency of our estimate of $`\alpha `$ across years, despite large changes in the number and spatial distribution of sampling locations through time (Fig. 3b). These results directly imply that average local abundance $`s=S(t)/\stackrel{~}{n}(t)`$, scales with total (regional) abundance according to $`sS(t)^\alpha `$.

We can gain further insight into the organization of a species population in different routes by considering how the distribution of number of routes with non-zero counts depends on total counts. That is, we may quantify the organization of the subpopulations through the conditional probability density $`\rho (\stackrel{~}{n}|S)`$, which measures the probability to find a bird species with S total counts having non-zero counts in $`\stackrel{~}{n}`$ distinct routes. Fig. 3 suggests that $`\rho (\stackrel{~}{n}|S)`$ will have a peak that increases as a power law with $`S`$. As shown in Fig. 4a this is indeed the case. If the data exhibit scaling, we should be able to identify a universal scaling function $`h`$ such that

$$\rho (\stackrel{~}{n}|S)\frac{1}{S^\alpha }h\left(\frac{\stackrel{~}{n}}{S^\alpha }\right),$$
(8)
We test the scaling hypothesis in Eq. (8) by plotting the scaled variables:

$$S^\alpha \rho \left(\frac{\stackrel{~}{n}}{S^\alpha }|S\right)vs.\frac{\stackrel{~}{n}}{S^\alpha }.$$
(9)
Fig. 4b shows that all curves collapse onto a single curve, which yields the scaling function $`h(u)`$.

## 4 Discussion

Our analysis differs from Taylor’s original studies (Taylor, 1961; Taylor and Woiwod, 1982; Taylor, 1984) in an important way. Taylor was interested in comparative analysis and so calculated a separate exponent for each species. He did this by analyzing multiple samples, replicated across time or space, for each species. Here, we calculate exponents replicating *across* species. One advantage of this approach is that we analyze the time series of total counts, summed over the entire survey. These total counts are considerably more robust estimates of abundance than local counts taken at individual routes.

Another advantage of analyzing scaling across species is that is allows us separate general patterns or “laws” that are invariant across taxonomic groups from general rules that may explain *deviations* from these laws. Our reasoning is that when the physical dimensions of a problem, such as energy or material flows, or spatial population structure, predominate, we should observe scaling laws that do not depend strongly on the biological differences among species, but that species-specific differences should appear as residual variation after the common scaling laws are factored out. That so many species fall along a single scaling relationship describing variability as a function of population size (Fig. 12) suggests that there may be universal features to the way in which North American breeding bird populations are subdivided spatially. We find exactly these features in the invariant, 3/4-power scaling of number of occupied survey routes versus total population size (Fig. 34).

However, not all of the variability in the data is accounted for by these scaling laws. For example, species with average total counts of approximately 250 individuals exhibit more than two orders of magnitude range in growth rate standard deviation (Fig. 1b). We believe that this residual variation does reflect important aspects of the ecology of individual species. There is a strong correlation between the residuals in Fig. 1b and the area of the corresponding species ranges, measure in terms of the average number of non-zero routes (T. Keitt, unpublished results). A likely explanation for this pattern is that fluctuations in the abundances of broadly distributed species will tend to average out spatially because different regions are influenced by geographically distinct climate regimes. Thus, it appears that species whose life histories tend to produce strongly aggregated distributions (i.e., species that are locally common, but regionally rare) are the ones that fluctuate the most relative to their total abundance. Species that have broad spatial distributions (i.e., locally rare, but regionally common) are therefore expected to fluctuate less than similar species with more restricted geographic ranges.

Ranking species in terms of their residual deviation from the growth-scale law (Fig. 1b) supports our hypothesis that locally common, but regionally rare species fluctuate more than expected, and vise versa. Large positive residuals correspond to species with restricted geographic ranges, such as Golden-cheeked Warbler (Dendroica chrysoparia; 2.5 times more variable), species that are habitat specialists and nest in large colonies, such as Tricolored Blackbird (Agelaius tricolor; 13.7 times more variable), species that breed in large groups called “leks”, such as Greater Prairie Chicken (Tympanuchus cupido; 3 times more variable), and species that show strong local migration patterns in response to changes in resource availability, such as White-winged Crossbill (Loxia leucoptera, 3.5 times more variable) and Red Crossbill (Loxia pytyopsittacus; 3.7 times more variable). Species that show low variability in relation to the scaling-law are typically solitary, territorial breeders such as Yellow-throated Warbler (Dendroica dominica; 2.8 times less variable), Prairie Falcon (Falco mexicanus; 2.6 times less variable), Swamp Sparrow (Melospiza georgiana; 2.5 times less variable), Kentucky Warbler (Oporornis formosus; 2.4 times less variable), and Chuck-will’s-widow (Caprimulgus carolinensis; 2.3 times less variable). The important point is that had we started from a purely autecological standpoint and ignored the important physical dimensions of the problem (e.g., structure of geographic ranges), we could easily have missed key pattern in terms of deviations from general scaling laws.

We should however mention several caveats. We do not as of yet know whether our results can be generalized to include other, non-avain taxonomic groups, or to other continents and climate regimes. Also, despite our use of highly aggregated, and therefore more robust time series, we suspect that there remain sources of variation in our analysis unrelated to actual population fluctuations. One vexing problem is repeated local migration between sampled and unsampled locations (we call this “sloshing”). Even if there is no variation in the true abundance across years, sloshing will lead to a given individual being counted in some years and not others, leading to measurement errors in the time series. We suspect this effect is not a significant component of variation in most of our time series, but may be substantial in a few cases. (Sloshing may contribute to the extreme variability of Tricolored Blackbird, for example.) Additional data, such as mark-recapture, may be need to resolve this issue.

There are two additional mechanisms related to our model for geographically subdivided populations that we have not discussed. We have shown how a non-linear dependency of the number of subpopulations versus total population size may explain the observed deviation from 1/2-power scaling of population fluctuations. Our basic hypothesis depends on the average local abundance scaling with total abundance in independently fluctuating subpopulations of roughly equal size. However, there are other patterns that may influence the “effective” number of independently fluctuating subpopulations and thus partially amount for the observed exponent in the growth-scaling law. First, large spatial variation in local abundance (Brown et al., 1995) could cause wide-spread species to fluctuate with greater magnitude than if all subpopulations have the same local abundance, since most of the variation would be driven by a few, high abundance sites. Second, strong spatial autocorrelation in population growth increments or “spatial synchrony” among fluctuating subpopulations (Grenfell et al., 1998; Bjørnstad et al., 1999; Kendall et al., 2000; Lundberg et al., 2000) may also cause a reduction in the effective number of independent subpopulations, and thus account for the increased magnitude of fluctuation in broadly distributed species. Temporal autocorrelation may act similarly to increase or decrease variability relative to our model. The consequences of these mechanisms need further exploration.

A surprising result of our analysis is the, to our knowledge, previously unreported 3/4-power scaling of spatial distribution as a function of total population size (Fig. 3). This result is closely related to, but not the same as, the “Distribution-Abundance” curve of Hanski and Gyllenberg (1997) that describes the fraction of regional habitats occupied as a function of average local abundance. We do not as of yet have an explanation for why the exponent should take this particular value, nor why it is so consistent through years. Recently, there has been considerable interest in explanations for the apparent 3/4-power scaling law relating body mass to metabolic output (Enquist et al., 1998; West et al., 1999; Dodds et al., 2001; Niklas and Enquist, 2001). One explanation posited to explain 3/4-power scaling is optimal structuring of a fractal transport network, such as the vascular system of plants and animals (West et al., 1999). This suggests an interesting hypothesis to explain 3/4-power scaling in our analysis: if the geographic ranges of species are subdivided in according to a particular fractal pattern, perhaps because of the fractal nature of the physical environment (e.g., Rinaldo et al., 1995), then it might lead to our observed scaling laws. Testing this hypothesis will require additional study.

It is interesting to note that our results are in striking qualitative agreement with similar studies from a broad range of social systems, ranging from growth of companies in the U.S. economy to the GDP of countries (Stanley et al., 1996; Lee et al., 1998; Plerou et al., 1999), suggesting that our simple model of growth may apply quite broadly (Amaral et al., 1998). Our observation that more “specialized” birds (in terms of smaller number of subpopulations) fluctuate more in total number than those that average fluctuations over many subpopulations may have an interesting parallel in social organizations: those that specialize on a few economic activities, e.g., countries with a single export product, may fluctuate considerably more than similarly sized organization with diverse economic activities, e.g., countries that produce a range of products. Putting all of one’s eggs in a single basket, as the saying goes, some times leads to catastrophes, and, it appears, greater variability as well.

## 5 Acknowledgements

T. Keitt thanks the Santa Fe Institute and the National Center for Ecological Analysis and Synthesis for support during the initial phase of this research. This research was made possible by the efforts by thousands of U.S. and Canadian BBS participants in the field, as well as, USGS and CWS researchers and managers.

## Appendix: Corrections applied to time series

Let $`s_{tuv}`$ be the number of birds of species $`u`$ counted at route $`v`$ in year $`t`$. The raw total counts $`R_{tu}=_v^{N_t}s_{tuv}`$ contain information about the abundance of species $`u`$ in year $`t`$ as well as information about the number $`N_t`$ and distribution of routes surveyed in year $`t`$. The goal is to remove the bias in the counts $`R_{tu}`$ introduced by variation in the number and distribution of survey routes through time. We do this by replacing each count $`s_{tuv}`$ for a given species at a given route with the time average $`\mu _{uv}=T_v^1_ts_{tuv}`$ for that route and species, where $`T_v`$ is the number of years that route $`v`$ was surveyed. We then construct new, surrogate time series $`M_{tu}=_v^{N_t}\mu _{uv}`$ whose variation only reflects changes in the number and distribution of survey routes through time (because the same $`\mu _{uv}`$ is used in each year), and not any real change abundance. We can then generate a bias corrected time series by subtracting these new time series from the raw totals:

$$S_{tu}=R_{tu}M_{tu}+\overline{M}_u$$
(10)
where $`\overline{M}_u`$ is the time average of $`M_{tu}`$ for species $`u`$. The advantage of this approach is that survey routes added or removed outside a species range will not influence the corrected total, because these routes will have $`\mu _{uv}=0`$.


text:
# Thermal Treatment of the Minority Game

## I Introduction

The Bar Attendance Model (BAM)  and the Minority Game (MG) (see Refs. \- ) have recently became regular testing grounds to investigate how the individual actions of a system of independent agents give rise to some kind of macroscopic ordering. In the MG, the agents have to make a binary decision which for the sake of concreteness, it is usually taken to be associated to going or not going to a bar. The winning option is that of the minority. The MG is a particular case of the BAM which has in turn been introduced to show how an ensemble of agents that perform inductive reasoning can self organize to match some condition that is generally accepted to be the most adequate. In the case of the BAM this corresponds to the largest acceptable attendance without incurring in some discomfort.

Both models have been compared with each other in Refs. and  working out a generalized version of the MG (the GMG)in order to consider situations in which the minority is replaced by an arbitrary fraction $`\mu `$ of the ensemble of players. This is fixed externally as a control parameter. In all these models the players update their attendance probabilities with a random correction, depending upon the past record of successes and failures. Asymptotic stable configurations are always reached. These are, however, of quite different nature depending upon the values of the control parameters, of the initial conditions and on the updating rules involved in each model.

In the present work we are interested in the cases in which the asymptotic stable distribution can be assimilated to a kind of thermodynamic equilibrium. In these situations the agents continue to update their attendance probabilities but the corresponding probability density distribution remains stationary. The stochastic dynamics that has been developed for the BAM in ref. always leads the system to these type of configurations while in the cases studied for the GMG, when $`\mu `$ is significantly larger (or smaller) than 1/2, the system gets stuck in quenched configurations that strongly depend upon the initial conditions. Updating stops because agents have accumulated a great number of successes. However, these “glassy” states can nevertheless be “melted” into equilibrium if the memory of past successes is repeatedly eliminated in an iterative process that can be assimilated to an annealing procedure.

A remarkable result that has been obtained in all numerical simulations is that the equilibrium configuration entails a diversity in the individual actions. The population is drastically partitioned into two subsets, one that always goes to the bar and the other that never goes. It therefore seems that in spite of the fact that the agents do not exchange information, they manage to coordinate their actions to proceed in two opposite ways. The number of agents in both subsets are in a ratio that is equal to $`\mu /(1\mu )`$. Such polarization is not an intuitive result. A naïve guess is to assume that all agents should choose the same probability of attendance and this should be equal to $`\mu `$. However this turns out to be not a stable distribution because parties that are larger or smaller than the accepted crowding occur with a great chance.

The fact that all agents adjust their attendance probabilities in order to minimize their failures (i.e. to go when the bar is crowded or not go when the bar is empty) leads to an aggregate behavior that minimizes a global cost associated to inadequate attendances. We propose to express such cost by the second moment of the attendance with respect to the acceptable level $`\mu `$.

The purpose of the present paper is to investigate the effects of introducing that cost function in the relaxation dynamics of the system. We show that this is a Lyapunov function for the many agent system, i.e. it is possible to derive a deterministic dynamics as the descent along its gradient, that monotonically reduces its value. This corresponds to a heavily coordinated, synchronous evolution.

We prove that the cost function meets the requirements of an internal energy of the many agent system. We also introduce a temperature parameter through a Langevin-like approach that can be defined in terms of the fluctuations of the attendance strategies. Except for finite size effects this can be proven to be an intensive parameter. We also superimpose thermal fluctuations to the deterministic dynamics mentioned above. Depending upon the amplitude of these fluctuations, the polarization is gradually smeared until a point in which completely disappears.

The thermally modified, relaxation process that we define here is completely different from those involved in the GMG or BAM approaches that involve the independent and uncoordinated actions of all the agents. The latter involves a random updating of individual attendance strategies governed by a (small) uncertainty amplitude that is interpreted as the precision of such updating. We prove that in the limit of low temperature, and small uncertainty amplitude both dynamics lead to entirely equivalent asymptotic equilibrium configurations. The thermal interpretation of the uncertainty amplitude also allows to cast the annealing process presented in Refs. and  into a thermal framework as the well known case of simulated annealing .

In section II we derive the cost function, and in section III we investigate the dynamics that corresponds to the descent along its gradient. In section IV we present a Langevin approach to define the temperature in terms of the fluctuations that are present in the asymptotic equilibrium configuration. In Sec. V we compare this with more traditional approaches for the relaxation process. In section VI we draw the conclusions.

## II The cost function

Consider a set of $`N`$ agents that have a probability $`p_i(i=1,2,\mathrm{},N)`$ to go to the bar. The distribution of the $`p_i`$’s is given by the probability density function $`P(p)`$. As we shall shortly explain the $`p_i`$ are updated in time according to some dynamics and therefore the function $`P(p)`$ also changes in time.

In the ordinary rules of the GMG when a player goes to the bar and finds it is crowded or when she does not go and the bar is empty, loses a point. If the opposite happens she gains a point. The level of crowding is specified by the value of the control parameter $`\mu `$. When her account of points falls below zero she updates her attendance probability choosing at random a different value within the interval $`(p_i\delta p/2,p_i+\delta p/2)`$. When equilibrium is reached, the resulting distribution $`P(p)`$ concentrates the population in the immediate neighborhood of $`p0`$ and $`p1,`$ plus an almost vanishing contribution from intermediate values. The ratio of the areas below these two peaks is close to $`\mu /(1\mu )`$.

The aggregate behavior is associated to the density distribution $`𝒫(A)`$ that gives the probability of occurrence of a party of $`A`$ customers attending the bar. The function $`𝒫(A)`$ is of course completely determined by $`P(p)`$. In order to calculate it let us assume without loss of generality that all the agents distribute themselves into $`D+1`$ different bins of $`n_d(d=0,1,\mathrm{},D)`$ agents each, with strategies $`p_d=d/D`$. The density distribution $`P(p)`$ can then be written as:

$$P(p)=\underset{d=0}{\overset{D}{}}\frac{n_d}{N}\delta (pp_d).$$
(1)
With this assumption, the distribution $`𝒫(A)`$ can be written as:

$`𝒫(A)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}_0=0}{\overset{n_0}{}}}\mathrm{}{\displaystyle \underset{\mathrm{}_D=0}{\overset{n_D}{}}}{\displaystyle \underset{d=0}{\overset{D}{}}}[\left({\displaystyle \genfrac{}{}{0pt}{}{n_d}{\mathrm{}_d}}\right)p_d^\mathrm{}_d(1p_d)^{n_d\mathrm{}_d}]`$ (3)
$`\times \delta (A{\displaystyle \underset{d=0}{\overset{D}{}}}\mathrm{}_d).`$
We define the cost function for the whole ensemble of agents as in ref., namely as the second moment  with respect to the tolerated crowding level $`\mu `$:

$$𝒞=\underset{A=0}{\overset{N}{}}(AN\mu )^2𝒫(A)$$
(4)
In order to calculate it, we introduce Eq. (3) into the definition of Eq.(4) and perform first the summation over $`A`$ taking advantage of the $`\delta (A_d^D\mathrm{}_d)`$. Once this is done, one can perform the summations involved in each of the terms in which $`(N\mu \mathrm{}_d)^2`$ splits down. The summations over different $`\mathrm{}`$’s decouple from each other and result either in a 1; or in $`n_dp_d`$; in $`n_d^2p_d^2+n_dp_d(1p_d)`$ or in $`(n_dp_d)(n_d^{}p_d^{})`$. These terms can be gathered again to yield:

$`𝒞`$ $`=`$ $`(N\mu {\displaystyle \underset{d=0}{\overset{D}{}}}n_dp_d)^2+{\displaystyle \underset{d=0}{\overset{D}{}}}n_dp_d(1p_d)`$ (5)
$`=`$ $`N^2(\mu <p>)^2+N(<p><p^2>)`$ (6)
where $`<p^m>`$ stands for $`_pp^mP(p)=_dp_d^mn_d/N`$ for $`m=1,2`$. The expression of $`𝒞`$ given in Eq. (6) contains no assumption about the system being in equilibrium. This is the reason why $`𝒞`$ is proportional to $`N^2`$ instead of being proportional to the size $`N`$ of the system, as befits to an extensive magnitude. The numerical simulations however indicate that in equilibrium $`<p>=\mu `$ and therefore this term cancels except for possible fluctuations. Actually the $`O(N^2)`$ term is eliminated by any distribution $`P(p)`$ whose mean has the required value $`\mu `$. For an initial condition with uniformly distributed $`p_i`$’s and $`P_o(p)=1/N`$, as it is used for most simulations, the cost is $`𝒞=N^2(\mu 1/2)^2+N/6`$. Such initial condition is a good guess for the final distribution when $`\mu 1/2`$ (as for the most traditional settings of the MG), but it is indeed very poor for the GMG when $`\mu 1/2`$. In the next sections we discuss in greater detail the value of $`𝒞`$ in equilibrium.

The naïve guess $`P(p)=\delta (p\mu )`$ is also seen to cancel the $`O(N^2)`$ terms in $`𝒞`$. However such distribution causes that parties with $`A`$ close to, but different from $`N\mu `$ occur with a sizable probability. The $`O(N)`$ in $`𝒞`$ are minimized precisely when the probability of occurrence of such parties tends to zero by polarizing the population into two subsets with opposite attendance strategies. To see this we approximate the two peaked equilibrium distribution that is usually obtained in numerical simulations by

$$P(p)=\frac{n_1}{N}\delta (pp_1)+\frac{n_2}{N}\delta (pp_2)$$
(7)
One readily sees that the $`O(N^2)`$ terms are eliminated when $`n_1p_1+n_2p_2=\mu N`$ and the $`O(N)`$ terms are also eliminated if the two peaks are $`p_1=0;n_1=N(1\mu )`$ and $`p_2=1;n_2=\mu N`$. The relaxation dynamics that tends to minimize individual losses is therefore seen to also optimize the global cost function defined in Eqs. (4) and (6)

## III A deterministic dynamics for the GMG

All the agents of the system, through uncoordinated actions minimize the total cost $`𝒞`$ that is an aggregate function defined for the whole system. This fact suggests an alternative representation of the actions of the agents as a synchronous, deterministic dynamics associated to the descent along the gradient of $`𝒞`$. This is described by the following set of coupled differential equations for the $`p_i`$’s:

$$\frac{dp_i}{dt}=\eta \frac{𝒞}{p_i}=\eta [2N(\mu <p>)(12p_i)]$$
(8)
In Eq. (8) $`\eta `$ stands for a positive free parameter that - as we shall shortly see - provides the scale for the time evolution of the system. The $`O(N^2)`$ and $`O(N)`$ terms in Eq. (6) are translated into a fast and a slow dynamics that involve corrections of the $`p_i`$ that are respectively $`O(N)`$ and $`O(1)`$. To see this we first derive the dynamics followed by $`<p>`$ by calculating the average over $`i`$ in both sides of Eq 8. We thus obtain:

$$\frac{dW(t)}{dt}=2\eta (N1)W(t)2\eta (\frac{1}{2}\mu )$$
(9)
where we have set $`W(t)=(<p>\mu )`$. This can explicitly be integrated. The solution is:

$$W(t)=\frac{\mu 1/2}{N1}+W_o\text{e}^{2\eta (N1)t}$$
(10)
with $`W_o`$ standing for the initial value of $`W(t)`$. This expression allows in turn to find an approximate solution of the equations of motion for the individual $`p_i`$’s. To this end we write an asymptotic approximation of Eq. (8) in which we assume that a long enough time has elapsed so that $`<p>\mu `$ can be approximated by the constant term of $`O(1/N)`$ in Eq. (10). By keeping only the leading order in $`N`$ we obtain:

$$\frac{dp_i}{dt}=2\eta (p_i\mu ).$$
(11)
Note that dependence of $`p_i(t)`$ involves a positive exponential. However, this equation is not valid for $`t\mathrm{}`$ because the fact that the $`p_i`$’s are probabilities, and are therefore bounded between 0 and 1, it is not included in the equations but rather in the boundary conditions of Eqs. (8).

Eqs.(9) and (11) correspond respectively to the fast and slow dynamics that have been mentioned above. In the first place we see that except for terms that are $`O(1/N)`$, $`<p>`$ approaches $`\mu `$ exponentially with the very short time constant $`\lambda =1/(2\eta N)`$ that tends to zero as the system involves a larger number of individuals. On the other hand, the differences $`p_i(t)\mu `$ instead grow exponentially for all $`i`$ indicating that the $`p_i`$’s depart exponentially from the average $`\mu `$ and eventually saturate at its largest or smallest possible values: 1 or 0, thus polarizing the population of agents. This process however takes place with a time constant $`1/(2\eta )`$, that is $`O(N)`$ longer than the one involved in the evolution of $`<p>`$ and is independent of the size of the system. While the average $`<p>`$ approaches very fast to the value $`\mu `$, the individual $`p_i`$’s depart slowly from the same value.

Eqs.(8) can be tested numerically by approximating them by finite differences. The individual attendance probabilities $`p_i`$ are thus taken to be updated as $`p_i(t+1)=p_i(t)+\mathrm{\Delta }(p_i)`$ where :

$$\mathrm{\Delta }(p_i)=\eta [2N(\mu <p>)(12p_i)]$$
(12)
The resulting density distributions $`P(p)`$ that are obtained with this dynamics are shown in Fig. 1. The value of $`\eta `$ and therefore that of the time constant $`\lambda `$ is in principle arbitrary. However if $`\lambda 1`$ the only effects that are noticeable are those of the fast dynamics while if $`\lambda 1`$ the descent towards the minimum keeps bouncing at opposite sides of the quadratic well and never reaches its bottom. When $`1/2\stackrel{<}{}\lambda \stackrel{<}{}2`$ the descent is gradual enough so that the interplay of both terms in $`\mathrm{\Delta }(p_i)`$ leads the system to a minimum of $`𝒞`$.

The intermediate stages in the gradient descent are also shown in Fig. 1. In the first few steps the (fast) uniform correction of $`O(N)`$ is seen to shift rigidly the initial distribution to one side with the aim of adjusting the value of $`<p>`$ to that of $`\mu `$. As a consequence, agents are piled up in one end while the other is completely cleared. Once the leading term in $`𝒞`$ is nearly canceled, the slow dynamics gradually gathers agents at both ends of the distribution producing minor fluctuations in the value of $`<p>`$. The density distribution $`P(p)`$ that is finally obtained is seen to correspond to a strongly polarized population thus reproducing the main feature of the equilibrium distributions obtained with the rules traditionally used in the GMG or the BAM.

The present approach yields a density distribution that displays the same polarization that is found in the GMG or in the BAM. It is remarkable that such a general qualitative agreement is found, although those frameworks differ deeply from the deterministic formulation. The conceptual difference between the two approaches lies in the special role played by the record of successes and failures that is kept in the BAM or GMG and that is completely absent in the present treatment. The usual rules of the GMG can thus be considered to correspond to a dynamics constrained by the (positive) balance of points that have been accumulated in the past instances of the game. There are other differences that deserve further discussion. These are related to the stochastic elements of the dynamics used in that framework which are absent from the present one. Within this approach, these can be assimilated to the effects of a finite temperature. We turn to this point in the next section.

## IV Thermal fluctuations

The usual rules of the BAM or the GMG involve a stochastic updating of the attendance probabilities of each customer. When the account of points of the $`i`$th player falls below zero a new value of $`p_i`$ is chosen at random from the interval $`(p_i\delta p/2,p_i+\delta p/2)`$. This can be interpreted as a kind of thermal fluctuation in which $`\delta p`$ can be related to the temperature.

A few qualitative features support this. In equilibrium, the population is drastically polarized into those that consistently go to the bar (and therefore $`p_i=1`$) and those that do not go ($`p_i=0`$). A small fraction having $`p_i`$’s with intermediate values continuously migrate between both extreme strategies. This migration causes that the value of $`<p>`$ fluctuates around $`\mu `$. These random values of $`<p>`$ have a distribution that is sharply peaked at that value and has a width that is regulated by $`\delta p`$. In what regards the density distribution $`P(p)`$, a small value of $`\delta p`$ produces sharp peaks at $`p=0`$ and $`p=1`$ and $`P(p)0`$ for intermediate values. For larger values of $`\delta p`$ there is a larger fraction of players that migrate between $`p=0`$ and $`p=1`$ thus producing a rising in the “bottom” of the distribution $`P(p)`$.

The above qualitative arguments provide hints to introduce thermal fluctuations in the deterministic dynamics presented in the preceeding section and also about their relationship with $`\delta p`$ for the case of the GMG. However a singular situation occurs for $`\delta p0`$ that is associated to an infinitely long relaxation process or when $`\delta p>1`$ in which this parameter loses its physical meaning of a being a probability.

Thermal-like fluctuations can formally be introduced following the same steps as the Langevin approach to describe a Brownian particle.In the present situation we start with the Eq. (9) for the motion of the average value $`<p>`$, and we add a stochastic term $`L(t)`$ that accounts for the random fluctuations

$$\frac{dW_s(t)}{dt}=2\eta (N1)W_s(t)2\eta (\frac{1}{2}\mu )+L(t)$$
(13)
We have added an index $`s`$ to $`W(t)`$ in Eq. (9) to stress the fact that this is the value of $`W(t)`$ in the presence of stochastic external fluctuations. The source of noise $`L(t)`$ can be taken to be the average of $`N`$ uncorrelated sources of random fluctuations affecting all the independent agents. One still has to specify a parameter related to the statistical properties of the distribution of the stochastic function $`L(t)`$. We will shortly prove that this is closely related to the temperature. As usual we assume:

$`\overline{L(t)}`$ $`=`$ $`0`$ (14)
$`\overline{L(t)L(t^{})}`$ $`=`$ $`\mathrm{\Gamma }\delta (tt^{})`$ (15)
In Eqs. (14) and (15) and in all what follows $`\overline{(\mathrm{})}`$ denotes an average over a suitable ensemble of replicas of the $`N`$ agent system. The parameter $`\mathrm{\Gamma }`$ is a constant that represents the mean square amplitude of instantaneous, uncorrelated perturbations. The stochastic differential equation (13) can explicitly be integrated. The result is

$$W_s(t)=W(t)+\text{e}^{2\eta (N1)t}_0^t\text{e}^{2\eta (N1)\omega }L(\omega )𝑑\omega $$
(16)
where $`W(t)`$ is the solution given in Eq(10) in which no fluctuations are present. If an average is made on both sides of Eq. (16), over a sub-ensemble of systems having the same initial conditions $`W_o`$ appearing in Eq. (10), one can immediately see that Eq. (14) implies that $`\overline{W_s(t)}=W(t)`$ and therefore the convergence of $`<p>`$ to $`\mu `$ (up to terms $`O(1/N)`$ is also insured within the stochastic dynamics. If the mean square fluctuations of $`W_s(t)`$ are calculated with the aid of Eq. (15), we get:

$$\overline{W_s^2(t)}=W^2(t)+\frac{\mathrm{\Gamma }}{4\eta N}\left[1\text{e}^{4\eta Nt}\right]$$
(17)
The effect of the stochastic term in $`W_s(t)`$ produces a non vanishing value $`\overline{W_s^2(\mathrm{})}`$. In ordinary statistical mechanics, the mean square fluctuations of the stationary solution of the velocity of Brownian particles is directly related to its average kinetic energy and can be set equal to $`kT`$. By analogy we formally define a temperature parameter $`T`$ that is independent from the size of the system, as the mean square fluctuations of $`<p>`$ in an equilibrium configuration, scaled by the number of agents of the system. Neglecting terms $`O(1/N^2)`$ we obtain:

$$TN\overline{(<p>\mu )^2}=\frac{\mathrm{\Gamma }}{4\eta }$$
(18)
The parameter $`\eta `$ is a factor relating $`T`$ with the amplitude of the random fluctuations and plays a similar role than the Boltzmann constant.

Eq. (18) allows to write the ensemble average of the cost $`\overline{𝒞}`$ for an equilibrium configuration and for finite temperature. Up to the leading order in $`N`$ we obtain:

$`\overline{𝒞}`$ $`=`$ $`N^2\overline{(\mu <p>)^2}+N\overline{(<p><p^2>)}`$ (19)
$`=`$ $`N[T+\mu \overline{<p^2>}]`$ (20)
$`\overline{𝒞}`$ is a positive, extensive magnitude which, in equilibrium, grows linearly with the size of the system and can therefore be taken to play the role of an internal energy.

The linear dependence of $`\overline{𝒞}`$ with the size of the system can be checked for the GMG. To do so we have calculated the cost using the definition of Eq.(4), with different number of agents. We first allowed the system to relax to the asymptotic equilibrium configuration and performed a suitable ensemble average over several replicas of the system. The linear dependence is shown in Fig. 2. The last iteration steps are used to estimate the dispersion of the numerical result and is shown with a pair of dotted lines. The slope of these lines change slightly with the parameter $`\delta p`$ of the GMG. This is due the relation between $`T`$ and $`\delta p`$ that we discuss later.

## V Thermal relaxation

To include thermal fluctuations into a numerical treatment of the deterministic dynamics amounts only to introduce a random additive term in Eq.(12), namely:

$$p_i(t+1)=p_i(t)+\mathrm{\Delta }(p_i)+L_\tau ^{(i)}(t)$$
(21)
where $`L_\tau =\tau (1/2r)`$ and $`r`$ is a random number uniformly distributed in the interval $`[0,1]`$. This function represents the fluctuations produced on the $`i`$-th agent by a thermal bath. The temperature is defined by the second moment $`\mathrm{\Gamma }`$ of the probability density of the $`L_\tau ^{(i)}(t)`$.

The limit in which $`L_\tau ^{(i)}(t)`$ has zero width (and therefore $`\tau =0`$) corresponds to the deterministic dynamics discussed in Sec. II. Larger values of $`\tau `$ are associated to fluctuations that may eventually override the updating amplitude $`\mathrm{\Delta }(p_i)`$ and tend to smear the distribution with two sharp $`\delta `$-functions, increasing the fraction of the population that have strategies $`p_i0`$ or 1. (see Fig. 3$`(a)`$). If $`\tau `$ is further increased the polarization is progressively destroyed because the drift of the $`p_i`$’s towards 0 or 1 has to equilibrate against random shocks that prevent them to reach those limiting values.

Given the stochastic dynamics of Eq.(21) together with the definition in Eq.(18) it is possible to calculate the value of $`T`$ in an equilibrium configuration, and relate $`T`$ with $`\tau `$. The parameter $`\eta `$ has to be chosen such that the relaxation of the deterministic dynamics is guaranteed i.e. when the time constant $`\lambda =1/(2\eta N)`$ introduced in Sect. 3 is $`\lambda 1`$. In Fig. 4 we show that, as expected, $`T\tau ^2`$.

Eq.(18) allows also to calculate $`T`$ in any configuration reached through the stochastic dynamics of the GMG or the BAM. With this we can check two important features. The first is an estimation of the finite size corrections in the definition of $`T`$ given in Eq.(18), i.e. the regime in which $`T`$ is independent of the size of the system. The second outcome is to establish a quantitative relationship between $`T`$ and $`\delta p`$ that goes into the relaxation dynamics of the GMG.

We have calculated $`W_s^2(t)`$ for the GMG using several values of $`\delta p`$ and $`N`$. We have allowed $`t`$ to be large enough to reach equilibrium. We have then performed an ensemble average over several replicas of the system. The last steps have been used to gauge the dispersion of the numerical values. The results are shown in Fig.4$`(b)`$ where we plot $`N\overline{W_s^2(\mathrm{})}`$ as a function of $`\delta p`$.

All the above mentioned features can be extracted from Fig. 4. Firstly finite size effects are clearly seen to affect only the smallest systems up to $`N500`$. Second the independence of $`N\overline{W_s^2(\mathrm{})}`$ from the size of the system as assumed in the definition of Eq.(18) follows from the fact that the curves for $`N500`$ lump tightly together. In the third place a linear regression of all the curves establishes that $`\delta p`$ and $`T`$ have the same physical interpretations, and within the interval considered are nearly proportional to each other, namely $`T=K\delta p`$, with $`K=(320\pm 20)10^4`$.

The fact that $`T`$ and $`\delta p`$ are conceptually equivalent leads to extend the GMG simulations to higher values of $`\delta p`$. These values have seldom been explored  in the literature because this parameter measures the minor adjustments performed by the agents that try to find the “best” attendance probability. Large values of $`\delta p`$ could for instance correspond to irresolute or hesitating agents.

There are however important points that have to be considered. In the first place the value of $`\delta p`$ can not be taken arbitrarily large. This is so because it measures the uncertainty of the value of a probability. Values of $`\delta p\stackrel{>}{}1`$ have therefore little physical meaning. In addition, if $`\delta p`$ is nevertheless extended to values higher than 1 by any plausible analytical extension (for instance using periodic or reflective boundary conditions), the fluctuations $`\overline{W_s^2}`$ for $`\delta p>1`$ are seen to saturate at an approximately constant value (see inset in Fig. 4$`(b)`$). These facts cause that the correspondence between $`\delta p`$ and $`T`$ necessarily breaks down.

A comparison of the probability density distributions $`P(p)`$ obtained with both approaches further supports this departure. In Fig. 3$`(b)`$ we show the equilibrium density distributions that are obtained with the stochastic, asynchronous updating rules of the GMG for two values of $`\delta p`$ (and $`\mu =0.6`$). It is seen that these diverge from those of Fig. 3$`(a)`$ that are obtained with the dynamics given in Eq.(21). Note however that there are noticeable ressemblances for small amplitude fluctuations. See for instance the distributions plotted in full line in Fig.3(b) and the one for $`\tau =0.003`$ in Fig.3(a).

As mentioned before, the origin of the departure between both dynamics can be found in the scoring of successes and failures that is used in the GMG, that is absent in the present approach. Some customers can be considered to be excluded from the updating dynamics as a consequence of their great accumulation of points. This, for instance, produces the large value of $`P(p=1)`$: many players that have accumulated a large positive account attending the bar do not change strategy. The scoring of each player works as a kind of “Maxwell Deamon” that classifies agents into different groups, endowing each one with a different updating rate.

The equilibrium configuration that is reached in the GMG therefore entails a distribution of updating rates in which some players are essentially frozen while others modify their attendance strategies frequently. This situation is completely different to the one obtained with the dynamics of Eq.(21) in which all agents undergo stochastic perturbations in every time step.

In order to show this we present in Fig. 5 some results of the GMG, in which we have used a large value of $`\delta p`$ ($`\delta p=0.8`$) and we have arbitrarily partitioned the ensemble of 1001 players into two sets. One of the sets gathers all players having at most 10 points the other contains all the rest. We have plotted their respective density distributions $`P(p)`$. The agents having less that 11 points are the ones that participate more strongly in the dynamics because undergo more frequent updatings.

The above comparison indicates that the GMG and the thermal relaxation dynamics of Eq.(21) strictly coincide only in the limit of $`T0`$. However the strong qualitative resemblance of the results for $`\delta p.6`$ allows to interpret $`\delta p`$, with these limitations, as equivalent to a thermal fluctuation.

The thermal interpretation of $`\delta p`$ has one interesting consequence. The most remarkable feature of the relaxation processes of the GMG performed with large $`\delta p`$ is that the high fluctuations prevents quenching (see Fig. 6). This allows to provide a new framework to the annealing procedure presented in Refs. and  that resembles more closely the traditional protocol of Ref..

The method presented in Ref. requires an iterative procedure which involves a short evolution of the $`N`$ agent system and the subsequent elimination of all points accumulated in the system. This is repeated until a moment in which the distribution $`P(p)`$ remains stationary. With the present interpretation of $`\delta p`$, a thermal annealing relaxation for the GMG can be performed for the cases in which $`\mu `$ is significantly different from 1/2. This new protocol can be assumed to take place in episodes. In the first episode, relaxation is allowed using a value of $`\delta p`$ that is large enough to insure that equilibrium is reached and quenching is prevented. The following episodes start from the equilibrium reached in the preceeding one, and a new relaxation process is allowed with a smaller value of $`\delta p`$ that is still large enough to avoid the appearance of quenching. The process continues until a lower bound of $`\delta p`$ is reached. Following this “cooling” protocol quenching never occurs, an absolute minimum of $`𝒞`$ is obtained and the population remains strongly polarized.

## VI Conclusions

In the present paper we provide and alternative description of the dynamics of a system composed by many agents that play at the GMG. This is given in terms of the optimization of a single global magnitude, instead of doing it in terms independent actions of the $`N`$ agents. We do this by studying the effect of introducing a cost function $`𝒞`$ that is associated to the second moment of the probability distribution of the size of the attending parties.

We have proven that $`𝒞`$ has the relevant properties of an internal energy. In equilibrium, it is a positive extensive quantity that scales linearly with the number of agents $`N`$ and its minima correspond to equilibrium configurations with a highly polarized population, as found in the BAM or the GMG without quenching.

In addition, the deterministic dynamics that is derived from the descent along the gradient of $`𝒞`$ leads the system to configurations that have an equivalent polarization as that found with the traditional stochastic updating of the BAM or the GMG. This is a non trivial equivalence between two completely different organization schemes of the N-agent system. On the one hand the gradient descent gives rise to a set of coupled differential equations that represents a coordinated evolution of all the agents as would be the result of the action of a “central planner” of the whole system. On the other hand, within the GMG all the agents act independently from each other adjusting their attendance strategies with the purpose of optimizing their individual utilities. Even though the two relaxation mechanisms are very different, the final configurations of the system turn out to have equivalent features.

The definition of $`𝒞`$ in terms of the second moment of the probability distribution of attending parties is reminiscent of the many body Hamiltonian introduced in Ref. to cast a version of the MG into the spin glass formalism. In the present case $`𝒞`$ can also be considered as a many body Hamiltonian with one- and two-body interactions in which the $`N`$ dynamic variables are the attendance probabilities $`p_i`$’s, with $`i=1,2\mathrm{}N`$.

The introduction of $`𝒞`$ and the associated relaxation process allows to define a temperature parameter through a Langevin-like approach. The value of $`T`$ remains associated to the ensemble average of the square of the fluctuations of the attendance, scaled by the number of agents. Its introduction in $`𝒞`$ provides the proof that this quantity, in thermal equilibrium, scales linearly with the size $`N`$ of the system and therefore qualifies as an extensive parameter.

On the other hand, in order to be an intensive parameter, $`T`$ should be independent of the size of the system. This has been checked numerically for the case of the GMG. However finite size effects in the definition of $`T`$ become negligible only for systems that are significantly larger than the minimal ones that already display the self organization features and that have spurred the popularity of the Minority Game.

Thermal fluctuations can be included in the dynamics that corresponds to the descent along the gradient of $`𝒞`$. The corresponding distributions $`P(p)`$ can readily be found and a comparison can be made of $`T`$ with $`\delta p`$ involved in the relaxation of the GMG or the BAM. A direct relationship can be established between both parameters but only in the limit of $`\delta p0`$. We have also considered the dynamics of the GMG with moderately large values of $`\delta p`$ when still the divergence between the GMG and the thermal dynamics is not important. A stochastic updating that involves large values of $`\delta p`$ could be thought to be associated to irresolute or badly informed agents that correct their attendance probabilities performing significant changes in each correction.

The GMG relaxation for large values of $`\delta p`$ avoids quenching even for $`\mu `$ significantly different from 1/2. This fact, together with the thermal interpretation of $`\delta p`$ allows to cast the annealing procedure presented in Ref. into the more traditional framework in which $`T`$ is progressively reduced in successive epochs. This “cooling” protocol could well be assimilated to a succession of learning episodes of the many agent system. In the first episodes in which agents have little “experience” and the information about the past is scarce, all agents perform large amplitude - even random - corrections. In the last episodes of the relaxation process, as there is a richer information about the past history of the system the agents perform finer corrections, the fluctuations are smaller and the cost paid by a wrong attendance are also smaller.

The fact that on the one hand an extensive magnitude can be defined playing the role of an internal energy, and that on the other, a microscopic definition of the temperature can be made, opens the way to a the full thermodynamic description of a system of N-agent performing a GMG. This amount to introduce a Gibbs distribution defined as $`\mathrm{\Phi }(𝒞)=e^{𝒞/T}/Z`$, where $`Z`$ stands for the partition function. All thermodynamic functions should follow from this.

E.B. has been partially supported by CONICET of Argentina, PICT-PMT0051;H.C. and R.P. were partially supported by EC grant ARG/b7-3011/94/27, Contract 931005 AR.