text:
# The Gravitational Million-Body Problem

## 1. Introduction

The title of this review unfortunately disguises the fact that it is about globular clusters, in the traditional sense. Actually, “a million” is somewhat on the large side in this context. Though the total number of stars in a globular cluster is not known, it may be estimated roughly by converting the luminosity (Harris 1996) to the mass, using the method of Mandushev, Staneva, & Spasova (1991), and adopting a mean stellar mass of $`0.34M_{}`$, which corresponds approximately to a power-law mass function with slope $`1`$ (Kroupa 2001) between turnoff and the hydrogen burning limit. The result is (Fig.1) that the median value is about $`3\times 10^5`$, and roughly 10% of all clusters have more than one million stars. In this review we highlight some of the ways in which $`N`$-body simulations are being used in studying the dynamical theory of globular star clusters.

## 2. Comprehensive Calculations111At Tokyo several people expressed some distaste for the phrase “kitchen-sink calculations”, which I used in my talk.

No simulation of a globular cluster will ever be complete, but rapid strides are being made to extend the traditional techniques of stellar dynamics so as to include stellar evolution. As a result there is a growing body of simulations with the following typical characteristics: (i) 1 star is represented by 1 particle (“star-by-star calculation”); (ii) forces include star-star gravity, a galactic tide (for a cluster on either a circular or elliptic galactic orbit), and disk shocking; (iii) treatment of stellar collisions; (iv) inclusion of stellar properties (luminosity, colour, metallicity) and stellar evolution; and (v) inclusion of (primordial) binaries and binary evolution.

There is a strong synergy between the development of such codes and the development of special-purpose hardware (see the papers by Aarseth, Makino and the “starlab” group in this volume.) Table 1 lists some recent comprehensive simulations, where $`f_{bin}`$ is the initial fraction of binaries. The scientific value of such work is well illustrated by Hurley et al (2001), who quantified the role of dynamical influences on the population of blue stragglers in an old open cluster. Such a result has been a long standing goal of this programme of research.

Despite such successes, it is disappointing that not one of these is a simulation of a globular cluster, in the traditional sense of the term. Even with GRAPE-4 it was possible to study a model with $`N=32768`$ well past core collapse (Makino 1996), which includes at least the smallest 10% of the existing galactic globular clusters.

There are two reasons why simulation of any globular cluster is much harder than is suggested by the particle number alone. One is the requirement for a substantial population of primordial binaries. Though it has been argued that the computational load associated with binaries becomes negligible when $`N`$ is sufficiently large (Makino & Hut 1990), this asymptotic limit is not yet in sight. The second reason is that the present particle number in a globular cluster may be a poor guide to what is required. Hurley et al (2001) estimated that even M67 had 7 times as many stars at birth as at present, and were therefore unable to simulate its early evolution. The present day mass function of the globular cluster NGC 6712 can be understood if it has lost about 99% of its original mass (Takahashi & Portegies Zwart 2000), and since its current mass is estimated to be of order $`10^5M_{}`$ (Pryor & Meylan 1993), its likely original particle number puts it far beyond reach.

One conclusion from this discussion is that, if a comprehensive simulation of a globular cluster is carried out within the next few years, it will tell us only about the evolution of the few clusters (if there are any) which are of low mass now and always have been. Study of the dynamical evolution of globular clusters cannot wait until such simulations can cope with all clusters.

One of the aims of comprehensive simulations, as we have seen, is the study of the interaction between stellar and dynamical evolution, and the effects of collisions. While the dynamics of point masses can be computed with precision, the effect of collisions is very uncertain. As Sills et al (2001) make clear, the effects of rotation on a merger remnant are dramatic. Since this is likely to be the largest source of uncertainty in studying the effects of collisions in globular clusters, it may be perfectly acceptable to adopt a model in which the point-mass dynamics is treated approximately, such as a Monte Carlo model (Sec.5). These are much faster than direct $`N`$-body methods, and can reach much larger “particle” numbers, ensuring that rare (but observable) species are not neglected.

## 3.  Restricted Simulations

Considerable progress has been made in recent years in understanding important aspects of the dynamics of globular clusters without attempting a comprehensive simulation. The example we shall concentrate on here is the simulation of tidal tails. It is a theoretical problem that has been stimulated by observations (Grillmair et al 1995, 1996; Kharchenko, Scholz, & Lehmann 1997; Leon, Meylan, & Combes 2000; Testa et al 2000; Odenkirchen et al 2001). For this purpose one need include only (i) point mass dynamics, (ii) relaxation, (iii) the galactic tide, and (iv) the galactic orbit. Binaries and stellar evolution can be ignored. This is a pure stellar dynamics problem, and several sets of simulations have been carried out: Moore (1996); Combes, Leon & Meylan (1999), Murali & Dubinski (1999), Johnston, Sigurdsson & Hernquist (1999), Dehnen (pers. comm., and Fig.2), and others. The list would be much longer if one added comparable work on satellite galaxies.

What is not always clear from simulations is the mechanism by which stars escape from a cluster into the tidal tail. The most obvious possibilities are time-dependent tides (disk and bulge shocking) and relaxation. But such simulations usually cover only one or two Gyr of evolution, and it is possible that the choice of initial conditions plays a role, just as in the case of circular orbits (Fukushige & Heggie 2000). The other hard thing to get right is relaxation, though the simulations of Combes et al (1999) show that this can be done.

The role of these various mechanisms can be considerably clarified by a judicious choice of initial conditions and parameters. Baumgardt (this volume) begins at apogalacticon, with a cluster whose radius corresponds to the Roche radius at perigalacticon. The rationale for this is not the assumption that, over the course of many galactic orbits, the radius of a cluster is set by its Roche radius at perigalacticon; such an assumption would be in contradiction with observational evidence (Ninkovic 1985, Allen & Martos 1988, Meziane & Colin 1996; Brosche, Odenkirchen & Geffert 1999). Rather, if the initial radius of the cluster is either much greater or much smaller than the perigalactic tidal radius, then it is clear that initial conditions have much to do with the subsequent mass loss.

Baumgardt’s finding from his $`N`$-body simulations is that, if one estimates the lifetime on the basis of mass lost in the first few orbits, as is often done, one finds that the lifetime increases with $`N`$ roughly as the relaxation time, until sufficiently large $`N`$ is reached, and then the lifetime is approximately independent of $`N`$. This is what would be expected if relaxation dominates for small $`N`$ and tidal heating dominates for large $`N`$. If, however, one waits for as many galactic orbits as are required to lose a certain fraction (such as 25%) of the mass, one finds that the mass loss is dominated by relaxation, for all $`N`$ up to the largest values he was able to study. Though this $`N`$ is still much smaller than the number of stars in the larger globular clusters, Baumgardt’s $`N`$-body models used particles of equal mass, and their relaxation time actually corresponds to the relaxation time of a much larger cluster with unequal masses.

One draws two conclusions from this study. One is that the initial conditions matter; it is only after the lapse of a sufficient number of galactic orbits that the mass loss rate reaches its asymptotic size. Second, the escape rate is dominated by relaxation, and not by bulge shocking, even on a very eccentric galactic orbit. It remains to be seen, however, whether disk shocking is ever important.

## 4.  New Stellar Dynamics

While it is clear that there are quantitative questions which can only be answered from $`N`$-body simulations, has anything qualitatively new come out of them? Indeed, yes, but the list is not a long one. Most recently, we have learned that the lifetime of clusters does not scale with $`N`$ in the expected manner, and the present section summarises this discovery.

The aim of a collaborative experiment carried out in 1997 (Heggie et al 1998; Heggie, this vol.) was to compare results of various kinds of codes (and not only $`N`$-body codes) in studying the dynamical evolution of a specific object. Since the number of particles was of order $`250000`$, it was necessary to scale $`N`$-body simulations by $`N`$. When this was done by the relaxation time, it was found that the lifetime of the system depended on $`N`$, though its internal evolution (e.g. time to core collapse) was fairly insensitive to $`N`$.

Fukushige suggested (pers. comm.) that the problem lay in the time taken for stars to escape, once they had enough energy to do so. If, however, it was assumed that the escape time scale is proportional to the crossing time, $`t_{cr}`$, then it seemed hard to understand how the problem persisted for the large values of $`N`$ (up to 64K) for which results were available. The lifetime of a cluster should vary nearly with the relaxation time, except for $`N`$ up to a few hundred, as in the escape theory of King (1959). Fukushige & Heggie (2000), however, found that the time taken to escape is a sensitive function of the energy, $`E`$, of a star: $`t_{esc}(EE_{esc})^2`$, where $`E_{esc}`$ is the energy needed to escape (i.e. the potential at the Lagrange point, where the potential well of the cluster opens out into the rest of the galaxy).

The next question was to find the appropriate energy to use in this formula. The answer is provided in Baumgardt (2001a). (His argument is more detailed than the following, but leads to the same scalings.) Unescaped stars with energies just above $`E_{esc}`$ are relaxing during the time they take to escape, and so their excess energy is $`EE_{esc}(t_{esc}/t_{rlx})^{1/2}`$, since relaxation is a diffusion process. (Here, $`t_{rlx}`$ is the relaxation time.) By solving these two relations we find the following conclusions: (i) $`t_{esc}\left(t_{cr}t_{rlx}\right)^{1/2}`$, much longer than the crossing time, and quite different from previous expectations; this helps to explain why the effects of the escape time persist for much larger $`N`$ than expected; (ii) $`EE_{esc}t_{rlx}^{1/4}`$; this is also the fraction of stars inside the cluster with energy above $`E_{esc}`$; this fraction is still a few percent, even for systems of the size of a typical globular cluster; and (iii) from the previous two conclusions, the time to lose some fraction (say, half) of the mass of the cluster is $`t_Mt_{rlx}^{3/4}`$.

The last is the most far-reaching conclusion. Though the theory applies only for circular galactic orbits, Baumgardt has shown empirically (these proceedings) that the same conclusion appears to hold for clusters on elliptic galactic orbits. And yet Vesperini & Heggie (1997) studied the evolution of the mass function of clusters assuming that $`t_Mt_{rlx}`$, which is now known to be wrong. Indeed a great deal of recent modelling of globular cluster systems (e.g. Capuzzo-Dolcetta & Tesseri 1997, Gnedin & Ostriker 1997, Murali & Weinberg 1997, and other work of authors already quoted) may require revision. In effect, the new scaling of $`t_M`$ means that the lifetime of a globular cluster through dynamical processes is smaller (for large $`N`$) than one would have thought. Baumgardt (2001b) shows that one consequence of this is that the initial conditions of the cluster system are more completely erased.

## 5.  Calibration of Faster Methods

### 5.1. Faster methods for star cluster dynamics

Besides $`N`$-body simulations, there are several other methods for studying the dynamical evolution of globular clusters (Fig.3). There are scaling methods, in which the cluster is characterised by nothing more than its mass and length scale (for example), Fokker-Planck methods, gas methods, and various hybrids. Even tree methods may function well for some important problems (McMillan & Aarseth 1993, Arabadjis & Richstone 1998, Miocchi & Capuzzo-Dolcetta 2001). All these methods are faster and involve more approximations, and $`N`$-body simulations are required in order to calibrate them.

We have already seen (Sec.2) that such methods may be quite appropriate for problems where the accuracy of the stellar dynamics is not the main uncertainty. Fast methods are also appropriate for the following kinds of problems:

1. Evolution of cluster systems: For references see Sec.4. For this purpose it must be possible to study the dynamical evolution of many globular cluster models with the full range of realistic $`N`$. At present this is possible only with fast approximate methods.
2. Modelling individual objects: Fokker-Planck and (occasionally) fluid methods have been used to construct evolutionary models for comparison with observations of specific clusters, including M3 (Angeletti, Dolcetta & Giannone 1980), M71 (Drukier, Fahlman & Richer 1992), NGC 6624 (Grabhorn et al 1992), NGC 6397 (Drukier 1993, 1995), M15 (Grabhorn et al 1992, Phinney 1993, Dull et al 1997) and M30 (Howell, Guhathakurta & Tan 2000). This requires extensive trial and error, with variations of the initial conditions of the model, and fast methods are essential. One motivation for such work has been the need to select suitable observational fields in clusters such that the local mass function requires least correction to the global mass function.
3. Predicting initial conditions: This is the theorists’ counterpart of the previous point. If one wants to construct an $`N`$-body model of a given object, considerable trial and error is needed in order to find appropriate initial conditions. Faster methods should be able to give a rapid but approximate answer, which can then be refined with a minimal number of $`N`$-body simulations.

### 5.2. Recent examples

#### Scaling models

As part of a study of the dynamical evolution of the mass function of a star cluster, Vesperini & Heggie (1997) gave a simple formula for the evolution of the mass, $`M(t)`$, i.e. $`M(t)=M(0)\mathrm{\Delta }M_{s.e.}0.828M(0)t/F_{CW}`$, where $`\mathrm{\Delta }M_{s.e.}`$ is the contribution from stellar evolution (which also depends on the mass function), and $`F_{CW}`$ is the “family” parameter of Chernoff & Weinberg (1990). For a given mass function, it is proportional to the relaxation time for average conditions inside the tidal radius, and is defined by $`F_{CW}=(M/M_{})(R_g/\mathrm{kpc})(220kms^1/v_g)\gamma ^{3/2}/\mathrm{ln}N`$, where $`M`$ is the (original) mass of the cluster, moving on an orbit of radius $`R_g`$ at speed $`v_g`$. $`\gamma `$ depends on the distribution of mass in the galaxy, is unity for a point-mass galaxy, and $`(3/2)^{1/3}`$ for an isothermal potential.

The above formula for $`M(t)`$ has been used by Vesperini (1998, 2000, 2001) in studies of the evolution of cluster systems. Its relevance to the present review is that it was based on the simplest scaling considerations, but the numerical coefficient was obtained by examination of $`N`$-body simulations. In other words, $`N`$-body results were used to calibrate a simpler theoretical model.

We now briefly consider the lifetime of a cluster. In the absence of stellar evolution, the above formula gives a cluster lifetime $`t_{VH}=F_{CW}/0.828`$ in Myr. The coefficient depends weakly on the structure of the initial model. A similar formula for $`M(t)`$ has been given by Portegies Zwart et al (2001b), which leads to a lifetime $`t_{PZMMH}=0.29t_{rxt}`$. Here $`t_{rxt}`$ is another measure of the relaxation time within the tidal radius, defined by $`t_{rxt}=0.138M^{1/2}r_t^{3/2}/(\overline{m}G^{1/2}\mathrm{log}_{10}\mathrm{\Lambda })`$ (cf. Spitzer 1987, eq.(2-63), except that Spitzer uses the natural logarithm in the Coulomb factor $`\mathrm{log}\mathrm{\Lambda }`$). In this formula, $`\overline{m}`$ is the mean stellar mass.

For the galactic mass law used by Portegies Zwart et al, $`\gamma =(3/1.8)^{1/3}`$, and so it is easy to show, if we equate the two forms of the Coulomb logarithm, that the ratio of the lifetimes is $`t_{VH}/t_{PZMMH}0.34M_{}/\overline{m}`$. The dependence on $`\overline{m}`$ arises because $`F_{CW}/t_{rxt}`$ also depends on $`\overline{m}`$. If it is assumed that the lifetime is proportional to $`t_{rxt}`$ (independent of the mass function, cf. de la Fuente Marcos 1995), then it is necessary to modify $`t_{VH}`$ to $`t_{VH}=(F_{CW}/0.828)(\overline{m}_{VH}/\overline{m})`$, where $`\overline{m}_{VH}`$ is the value in the simulations by Vesperini & Heggie (1997). In fact for most of their work, which used a mass function $`f(m)m^{2.5}`$ for $`0.1M_{}<m<15M_{}`$ we have $`\overline{m}_{VH}0.28M_{}`$. The effect of this change (which is also necessary in the formula for $`M(t)`$) is to bring the two independent estimates of cluster lifetime into fair agreement. It must be borne in mind, however, that the assumption that the lifetime is proportional to the relaxation time is incorrect (Sec.4), and that $`\overline{m}`$ is time-dependent.

Despite the fact that lifetime (in the absence of stellar evolution) does not depend much on the initial structure (if the cluster initially fills its Roche lobe), it would be useful to have an equally simple way of estimating the time-dependence of the concentration parameter.

#### Fokker-Planck models

The comparison between $`N`$-body and Fokker-Planck models (computed with finite differences) has a considerable history, but recent work appears to have brought the two methods into reasonable agreement, with one caveat, discussed below. The central issue is the escape criterion, as shown by Takahashi & Portegies Zwart (1998, 2000). These papers show how a free parameter in the Fokker-Planck specification can be determined by calibration using the results of $`N`$-body models. Incidentally, the second of these is also an extensive survey of globular cluster models, which essentially supersedes the comparable extensive Fokker-Planck work of Chernoff & Weinberg (1990).

As Takahashi & Portegies Zwart themselves point out, the caveat in this is that the $`N`$-body models used in this survey adopt a tidal cutoff. Since it is known that the lifetime varies with $`N`$ in a different way when a tidal field is used (Sec.4), it remains to be seen how well the agreement between Fokker-Planck and $`N`$-body models survives in the more realistic case.

#### Monte Carlo models

Though these can be regarded also as Fokker-Planck models, they offer some advantages over finite-difference models, and have received increasing attention in recent years. While part of this effort involves calibration with $`N`$-body models, the real purpose of this little subsection is to say something about the present state of play.

The largest models to date are those of Freitag & Benz (2001), who present results of Monte Carlo simulations with up to $`2\times 10^6`$ shells. At such levels it becomes possible to represent many globular clusters with one shell per star. These models, which are motivated by problems of galactic nuclei, do not include binaries, but have been extended to include stellar collisions (Fig.4). For models including the effects of dynamically formed binaries the record is held, at time of writing, by Giersz (this volume). Reaching well past core collapse with $`N=10^6`$ shells, as he does, is a notable landmark. Primordial binaries are excluded from these models, but a 10% population of primordial binaries are present in computations with up to $`3\times 10^5`$ shells by Rasio’s group (Rasio, Fregeau & Joshi, 2001; Fig.5). Soon they will compute few-body interactions “on the fly”, instead of using cross sections. All these developments give great promise for the Monte Carlo method in the coming years.

## 6. Conclusions

This rather selective review has concentrated on what $`N`$-body models are telling us about the dynamics of globular star clusters. Comprehensive simulations, which include details of evolution of single and binary stars, are at an exciting stage, but are still limited to rich open clusters. It seems unlikely that comprehensive simulations of globular clusters will be possible in the near future. Much progress can be made on some interesting kinds of dynamical problems about globular clusters (e.g. tidal tails) using more restricted kinds of simulation, which are perfectly feasible now. One of the dynamical lessons that has been learned in recent years is that it is vital to implement a correct treatment of tidal effects, even in the case of a steady tide (i.e. a circular galactic orbit). Monte Carlo methods now have the dynamical aspects of globular clusters within reach. When a detailed treatment of stellar evolution has been incorporated (at the level which has already been successful in the best $`N`$-body codes), these codes may well prove optimal for the comprehensive study of globular star clusters, until the time when real $`N`$-body models are fast enough.

##### Acknowledgments.

I thank the organisers of the meeting for their support, and I am grateful for comments from H. Baumgardt, S.P.F. Portegies Zwart and E. Vesperini on scientific matters. I also thank W. Dehnen, M. Giersz, M. Freitag and F. Rasio for the provision of unpublished results. The picture of M15 in Fig.3 was obtained using NASA’s SkyView facility (http://skyview.gsfc.nasa.gov) located at NASA Goddard Space Flight Center.

## References

Allen, C. & Martos, M. A. 1988, Revista Mexicana de Astronomia y Astrofisica, 16, 25

Angeletti, L., Dolcetta, R., & Giannone, P. 1980, Ap&SS, 69, 45

Arabadjis, J.S. & Richstone D.O. 1998, astro-ph/9810192

Baumgardt, H. 2001a, MNRAS, 325, 1323

Baumgardt, H. 2001b, in ‘Modes of Star Formation’, eds. E. Grebel & W. Brandner, ASP Conf. Ser (ASP, San Francisco), in press

Brosche, P., Odenkirchen, M., & Geffert, M. 1999, New Astronomy, 4, 133

Capuzzo-Dolcetta, R. & Tesseri, A. 1997, MNRAS, 292, 808

Chernoff, D. F. & Weinberg, M. D. 1990, ApJ, 351, 121

Combes, F. , Leon, S. , & Meylan, G. 1999, A&A, 352, 149

Dehnen, W. & Binney, J. 1998, MNRAS, 294, 429

de la Fuente Marcos, R. 1995, A&A, 301, 407

Drukier, G. A. 1993, MNRAS, 265, 773

Drukier, G. A. 1995, ApJS, 100, 347

Drukier, G. A., Fahlman, G. G., & Richer, H. B. 1992, ApJ, 386, 106

Dull, J. D., Cohn, H. N., Lugger, P. M., Murphy, B. W., Seitzer, P. O., Callanan, P. J., Rutten, R. G. M., & Charles, P. A. 1997, ApJ, 481, 267

Freitag, M. 2000, PhD thesis, Université de Genève;
http://obswww.unige.ch/$``$freitag/thesis/

Freitag, M., & Benz, W. 2001, A&A, 375, 711

Fukushige, T. & Heggie, D. C. 2000, MNRAS, 318, 753

Gnedin, O. Y. & Ostriker, J. P. 1997, ApJ, 474, 223

Grabhorn, R.P., Cohn, H.N., Lugger, P.M., & Murphy, B.W. 1992, ApJ, 392, 86

Grillmair, C. J., Freeman, K. C., Irwin, M., & Quinn, P. J. 1995, AJ, 109, 2553

Grillmair, C. J., Ajhar, E. A., Faber, S. M., Baum, W. A., Holtzman, J. A., Lauer, T. R., Lynds, C. R., & O’Neil, E. J. 1996, AJ, 111, 2293

Harris, W.E. 1996, AJ, 112, 1487

Heggie, D. C., Giersz, M., Spurzem, R., & Takahashi, K. 1998, Highlights in Astronomy, 11, 591

Howell, J. H., Guhathakurta, P., & Tan, A. 2000, AJ, 119, 1259

Hurley, J.R., Tout, C.A., Aarseth, S.J., Pols, O.R. 2001, MNRAS, 323, 630

Johnston, K. V., Sigurdsson, S., & Hernquist, L. 1999, MNRAS, 302, 771

Kharchenko, N., Scholz, R.-D., & Lehmann, I. 1997, A&AS, 121, 439

King, I. 1959, AJ, 64, 351

Kroupa, P. 2001, in S. Deiters, R. Spurzem, et al., eds, STAR2000: Dynamics of Star Clusters and the Milky Way, ASP Conf. Ser. (ASP, San Francisco), in press; astro-ph/0011328

Kroupa, P., Aarseth, S., Hurley, J. 2001, MNRAS, 321, 699

Leon, S., Meylan, G., & Combes, F. 2000, A&A, 359, 907

Makino, J. 1996, ApJ, 471, 796

Makino, J., Hut, P. 1990, ApJ, 365, 208

Mandushev, G., Staneva, A., & Spasova, N. 1991, A&A, 252, 94

McMillan, S. L. W. & Aarseth, S. J. 1993, ApJ, 414, 200

Meziane, K. & Colin, J. 1996, A&A, 306, 747

Miocchi, P. & Capuzzo-Dolcetta, R. 2001, A&A, submitted; astro-ph/0104152

Moore, B. 1996, ApJ, 461, L13

Murali, C. & Dubinski, J. 1999, AJ, 118, 911

Murali, C. & Weinberg, M. D. 1997, MNRAS, 291, 717

Ninkovic, S. 1985, AN, 306, 237

Odenkirchen, M. et al. 2001, ApJ, 548, L165

Phinney, E.S. 1993, in S.G. Djorgovski and G. Meylan, eds, Structure and Dynamics of Globular Clusters, ASP Conf. Ser. 50 (ASP, San Francisco)

Portegies Zwart, S.F., McMillan, S.L.W., Hut, P., Makino, J. 2001a, MNRAS, 321, 199

Portegies Zwart, S.F., Makino, J., McMillan, S.L.W., Hut, P. 2001b, ApJ, 546L, 101

Pryor, C. & Meylan, G. 1993, in S.G. Djorgovski and G. Meylan, eds, Structure and Dynamics of Globular Clusters, ASP Conf. Ser. 50 (ASP, San Francisco)

Rasio, F.A., Fregeau, J.M., Joshi, K.J. 2001, in ‘The Influence of Binaries on Stellar Population Studies’, ed. D. Vanbeveren (Kluwer, Dordrecht), in press; astro-ph/0103001

Sills, A., Faber, J.A., Lombardi, J.C., Rasio, F.A., & Warren, A.R. 2001, ApJ, 548, 323

Spitzer, L., Jr., 1987, Dynamical Evolution of Globular Clusters (Princeton, Princeton University Press)

Takahashi, K. & Portegies Zwart, S.F. 1998, ApJ, 503, L49

Takahashi, K. & Portegies Zwart, S.F. 2000, ApJ, 535, 759

Testa, V., Zaggia, S. R., Andreon, S., Longo, G., Scaramella, R., Djorgovski, S. G., & de Carvalho, R. 2000, A&A, 356, 127

Vesperini, E. 1998, MNRAS, 299, 1019

Vesperini, E. 2000, MNRAS, 318, 841

Vesperini, E. 2001, MNRAS, 322, 247

Vesperini, E. & Heggie, D. C. 1997, MNRAS, 289, 898


format: text

text:
# Particle and Astroparticle Searches for Supersymmetry

(November 22, 2001)

## Abstract

Supersymmetry may be discovered at high energy colliders, through low energy precision measurements, and by dark matter searches. We present a comprehensive analysis of all available probes in minimal supergravity. This work extends previous analyses by including the focus point branch of parameter space and the full array of promising indirect dark matter searches. We find that particle and astrophysical searches underway are highly complementary: each separately provides only partial coverage of the available parameter space, but together they probe almost all models. Cosmology does not provide upper bounds on superpartner masses useful for future colliders. At the same time, in the cosmologically preferred region, if supersymmetry is to be observable at a 500 GeV linear collider, some signature of supersymmetry must appear before the LHC.

preprint: MIT–CTP–3199, UCI–TR–2001–32, CERN–TH/2001-297, hep-ph/0111295, Snowmass P3-09

Weak-scale supersymmetry has consequences for a wide variety of particle physics experiments. At the same time, it is now well-established that luminous matter makes up only a small fraction of the mass of the observed universe Primack:2000dd . Neutralinos in supersymmetry are well-motivated candidates to provide much or all of the non-baryonic dark matter. While it is conceivable that all standard model superpartners are unstable, decaying, for example, into gravitinos or axinos, in many model frameworks, a neutralino is the lightest supersymmetric particle (LSP) and stable. Further, neutralino thermal relic densities of the desired amount are remarkably generic in supergravity models neutralinos . If neutralinos make up a significant portion of the halo dark matter, astrophysical searches for dark matter open up many additional avenues for the discovery of supersymmetry.

Here we summarize results of two studies of current and near future prospects for supersymmetry searches Feng:2000gh ; Feng:2001zu . We consider minimal supergravity, a simple framework that is conservative in the sense that many low energy signals generic in supersymmetry are suppressed. This work extends previous analyses in two ways. First, we include the focus point branch of supersymmetry parameter space FP , previously neglected. As we will see, neutralino dark matter in the focus point region has novel properties, with strong, positive implications for dark matter searches. Second, we consider not only collider experiments, low energy probes, and direct dark matter searches, but also the full array of indirect dark matter searches. We find that the various searches are highly complementary, and this comprehensive approach leads to strong conclusions that may be reached only by considering all available experiments.

The signals we consider, the projected sensitivities, and the experiments likely to achieve them, are listed in Table 1. Conventional particle physics probes include searches at LEP and the Tevatron, the improved measurement of the the $`BX_s\gamma `$ branching ratio at $`B`$ factories, as well as the projected final sensitivity of the Brookhaven $`g_\mu 2`$ experiment. On the astrophysics side, we consider the projected reach of detectors searching for direct dark matter interactions, including the CDMS cdms , CRESST cresst  and GENIUS genius  experiments.

We also include the most promising indirect dark matter searches. In the next five years, an astounding array of experiments will be sensitive to potential neutralino annihilation products. These include under-ice and underwater neutrino telescopes (AMANDA amanda , NESTOR nestor , ANTARES antares ), which will be sensitive to muons produced by neutrinos from relic annihilations in the cores of the Sun and Earth neutrinos ; atmospheric Cherenkov telescopes (STACEE stacee , CELESTE celeste , ARGO-YBJ argo , MAGIC magic , HESS hess , CANGAROO cangaroo , VERITAS veritas ), and space-based $`\gamma `$ ray detectors (AGILE agile , AMS/$`\gamma `$ amsgamma , GLAST glast ), which may detect gamma rays from the galactic center photons ; and anti-matter/anti-particle experiments (PAMELA pamela , AMS ams ), which may detect positrons positrons , antiprotons antiprotons  and antideuterium antideuterium  from the galactic halo. In many cases, these experiments will improve current sensitivities by several orders of magnitude. The search for supersymmetric dark matter is reviewed in a number of excellent articles indirectreviews ; see also Refs. RDrecent ; direct ; recent  and the bibliographies of Refs. Feng:2000gh ; Feng:2001zu .

The main results from our combined study Feng:2001zu  of dark matter and collider signals in supersymmetry prior to the start of the LHC are shown in Fig. 1. The regions of cosmologically interesting relic densities are shown; see Ref. Feng:2000gh  for details. A generous range is $`0.025\mathrm{\Omega }_\chi h^21`$, where the lower bound is the requirement that neutralino dark matter explain galactic rotation curves and the upper bound follows from the lifetime of the universe. Above this shaded region, $`\mathrm{\Omega }h^2>1`$; below, $`\mathrm{\Omega }h^2<0.025`$. The range $`0.1\mathrm{\Omega }_\chi h^20.3`$, preferred by current limits, is also shown.

Several striking conclusions follow from Fig. 1:

No Useful Upper Bounds from Cosmology. For all $`\mathrm{tan}\beta `$, cosmologically interesting densities are possible for either relatively small $`m_0`$, where the dark matter is Bino-like and coannihilation effects may be relevant coann , or in the focus point region at large $`m_01\mathrm{TeV}`$, where the dark matter contains a small but phenomenologically significant Higgsino component. As is evident from Fig. 1, for all values of $`\mathrm{tan}\beta `$, the cosmologically preferred region extends to very large $`m_0`$ and $`M_{1/2}`$; discovery at near future colliders cannot be guaranteed by cosmological arguments. The cosmologically preferred region is ultimately cutoff around $`M_{1/2}6\mathrm{TeV}`$, but this is far beyond the reach of the LHC and planned linear colliders.

Excellent Prospects for Combined Searches. At the same time, by combining all available searches, nearly all of the cosmologically preferred models predict signals in at least one experiment. This is strictly true for $`\mathrm{tan}\beta =10`$. For $`\mathrm{tan}\beta =50`$, some of the preferred region escapes all probes, but this requires $`M_{1/2}450\mathrm{GeV}`$ and $`m_01.5\mathrm{TeV}`$, and requires significant fine-tuning of the electroweak scale. In the most natural regions, all models with significant neutralino dark matter will yield some signal before the LHC begins operation.

Complementarity of Particle and Astrophysical Searches. Also noteworthy is the complementarity of traditional particle physics searches and indirect dark matter searches. Collider searches require, of course, light superpartners. High precision probes at low energy also require light superpartners, as the virtual effects of superpartners quickly decouple as they become heavy. Thus, the LEP and Tevatron reaches are confined to the lower left-hand corner, as are, to a lesser extent, the searches for deviations in $`BX_s\gamma `$ and $`a_\mu `$. These bounds, and all others of this type, are easily satisfied in the focus point models with large $`m_0`$, and indeed this is one of the virtues of these models. However, in the focus point models, all of the indirect searches are maximally sensitive, as the dark matter contains a significant Higgsino component. Direct dark matter probes share features with both traditional and indirect searches, and have sensitivity in both regions. It is only by combining all of these experiments, that the preferred region may be explored completely.

Implications for Linear Colliders. Finally, these results have implications for future colliders. In the cosmologically preferred regions of parameter space with $`0.1<\mathrm{\Omega }_\chi h^2<0.3`$, all models with charginos or sleptons lighter than 300 GeV will produce observable signals in at least one experiment. This is evident for $`\mathrm{tan}\beta =10`$ and 50 in Fig. 1. In Fig. 2, we vary $`\mathrm{tan}\beta `$, fixing $`M_{1/2}`$ to 400 GeV, which roughly corresponds to 300 GeV charginos. We see that the preferred region is probed for any choice of $`\mathrm{tan}\beta `$. (For extremely low $`\mathrm{tan}\beta `$ and $`m_0`$, there appears to be a region that is not probed. However, this is excluded by current Higgs mass limits for $`A_0=0`$. These limits might be evaded if $`A_0`$ is also tuned to some extreme value, but in this case, top squark searches in Run II of the Tevatron stop  will provide an additional constraint.) These results imply that if any superpartners are to be within reach of a 500 GeV lepton collider, some hint of supersymmetry must be seen before the LHC begins collecting data. This conclusion is independent of naturalness considerations.

While our quantitative analysis is confined to minimal supergravity, we expect this result to be valid more generally. For moderate values of $`\mathrm{tan}\beta `$, if the dark matter is made up of neutralinos, they must either be light, Bino-like, or a gaugino-Higgsino mixture. If they are light, charginos are typically also light, and can be discovered at colliders. If the dark matter is Bino-like, light sfermions are required to mediate their annihilation, and there will be anomalies in low energy precision measurements. And if they are a gaugino-Higgsino mixture, at least one indirect dark matter search will see a signal. For large $`\mathrm{tan}\beta `$, low energy probes become much more effective and again there is sensitivity to all superpartner spectra with light superpartners. These conclusions are not airtight — they may be evaded by highly split charginos and neutralinos, highly split sfermions, or enhanced annihilation from Higgs poles. For example, in the case of highly split sfermions, one could arrange for light staus to mediate neutralino annihilation, but very heavy smuons to suppress $`g_\mu 2`$. Barring such possibilities, however, if supersymmetry exists at the weak scale and neutralinos provide the bulk of dark matter, there are excellent prospects for supersymmetry discovery prior to LHC operation.

###### Acknowledgements.

The work of J.L.F. and F.W. was supported in part by the US Department of Energy under cooperative research agreement DF–FC02–94ER40818.


text:
# Tachyon Condensation and Black Hole Entropy

(October 2001)

## Abstract

String propagation on a cone with deficit angle $`2\pi (1\frac{1}{N})`$ is considered for the purpose of computing the entropy of a large mass black hole. The entropy computed using the recent results on condensation of twisted-sector tachyons in this theory is found to be in precise agreement with the Bekenstein-Hawking entropy.

preprint: HUPT-01/A051preprint: TIFR/TH/01-42preprint: hep-th/0111004

Introduction

The statistical interpretation of Bekenstein-Hawking entropy of a black hole remains an outstanding problem in quantum gravity. For a black hole in $`d`$ spacetime dimensions, the entropy $`S`$ is given by a universal formula

$$S=\frac{A}{4G}$$
(1)
that depends only on the area $`A`$ of the event horizon and the $`d`$-dimensional Newton’s constant $`G`$. Thermodynamically, this entropy behaves in every respect like ordinary entropy and unites the second law of thermodynamics with the area theorems of classical general relativity into an elegant generalized second law Beke ; Hawk . These beautiful results demand that, like any other entropy, the black hole entropy also must have a statistical interpretation in terms of underlying microstates.

Recent progress in string theory has shown that this is indeed true for a large class of supersymmetric black holes. The microstates of these special black holes are precisely countable and can completely account for their entropy Sen ; StVa . These striking results can be extended to near-extremal black holes as well as to certain nonsupersymmetric charged black holes. At present, however, these methods cannot be applied to the more general case of a nonsupersymmetric neutral black hole.

We will address the problem of microstates of a Schwarzschild black hole from a very different perspective based on earlier ideas of ’t Hooft Hoof , Susskind Suss , and others ZuTh ; BKLS ; Sred ; CaWi . ’t Hooft Hoof  has advocated that the thermal entropy of the heat bath seen by a Schwarzschild observer should account for the Bekenstein-Hawking entropy. This would then offer the desired statistical interpretation of the entropy in terms of the near-horizon microstates of the heat bath. In field theory, the leading contribution to the thermal entropy is quadratically divergent in the ultraviolet Hoof  and is proportional to the area of the horizon. If the cut-off is of the order of the Planck length, then the thermal entropy is of the right order of magnitude to be identified with the Bekenstein-Hawking entropy.

In the context of field theory, there are several difficulties with this appealing idea. For example, the thermal entropy depends on an arbitrary cutoff and its precise identification with the Planck length is not clear. The thermal entropy depends on the species and the couplings of the various particles in the theory whereas the black hole entropy is species-independent. Finally, since the thermal entropy always starts at one loop, it is difficult to see how it can possibly account for the tree level black hole entropy which is inversely proportional to the coupling constant. ’t Hooft has argued that it is necessary to understand the ultraviolet structure of the theory to address these questions and therefore these difficulties will be resolved in the correct short-distance theory of quantum gravity.

In string theory, ultra-violet divergences are expected to be appropriately controlled and Susskind, in particular, has argued that string theory offers a suitable framework for realizing this proposal Suss . We will pursue these ideas further by considering string theory on the near horizon geometry in Euclidean formalism following earlier work in Dabh ; DabhI .

Strings on a Cone

Consider a Schwarzschild black hole in four dimensions with a large mass $`M`$. The metric is given by

$$ds^2=(1\frac{2GM}{r})dt^2+(1\frac{2GM}{r})^1dr^2+r^2d\mathrm{\Omega }^2$$
(2)
We are interested in the thermal entropy of string modes in this background near the horizon as seen by a Schwarzschild observer. To focus on the near-horizon region we choose coordinates $`\rho =\sqrt{8GM(r2GM)}`$ and $`\eta =\frac{lt}{2GM}`$ in which the metric becomes

$$ds^2=(\frac{\rho }{l})^2d\eta ^2+d\rho ^2+(2GM)^2d\mathrm{\Omega }^2.$$
(3)
Here $`l`$ is an arbitrary parameter that sets the length scale which in our context will be taken to be the string length. We see from the form of the metric that $`\rho `$ measures the proper distance from the horizon at $`\rho =0`$ and $`\eta `$ is the proper time measured by a Schwarzschild observer located at proper distance $`\rho =l`$. We will refer to this location as the ‘stretched horizon’. Energy at the stretched horizon is related to asymptotic energy by the red-shift factor $`l/4GM`$.

In the limit $`GM>>l`$, the two-sphere at the horizon can be approximated by a flat 2-dimensional transverse space and then the space has flat Minkowski geometry. Thus, in this limit, the Schwarzschild observers are exactly like Rindler observers in uniform acceleration in Minkowski space. More generally, for a d-dimensional black hole the near horizon geometry is given by the $`d`$-dimensional Rindler metric with topology $`𝐑^2\times 𝐑^{d2}`$

$$ds^2=\rho ^2d\eta ^2+d\rho ^2+\underset{i=1}{\overset{d2}{}}(dx^i)^2,$$
(4)
where $`x^i`$ are the coordinates of the $`(d2)`$-dimensional transverse space and we have chosen units to set $`l=1`$.

The observer at the stretched horizon sees a thermal bath at Rindler temperature $`T=\frac{1}{2\pi }`$ which, as usual, is inversely proportional to the periodicity of Euclidean time. Under Euclidean continuation, $`\eta =i\theta `$, the metric in the $`𝐑^2`$ factor is given by the flat metric in polar coordinates

$$ds^2=\rho ^2d\tau ^2+d\rho ^2$$
(5)
The Euclidean time $`\tau `$ is thus an angular variable and when its periodicity is $`2\pi `$ corresponding to the Rindler temperature, the metric is smooth at the origin. The Rindler temperature at the stretched horizon is red-shifted to the Hawking temperature $`T_H=\frac{1}{8\pi GM}`$ seen by the asymptotic observer.

The partition function $`Z`$ of this thermal bath in string theory would be given by a functional integral of all string modes on the Euclidean Rindler space. The thermal entropy of this heat bath $`S`$ would be given, as usual, by $`S=\frac{F}{T}`$ as a derivative of the free energy, $`FT\mathrm{ln}Z`$. Therefore, to calculate the entropy, we need to find the variation in the free energy of a string gas to order $`\delta `$ as we vary the temperature seen by the observer at the stretched horizon to $`\frac{1}{2\pi }+\delta `$. Changing the temperature changes the periodicity of Euclidean time in (5) to $`2\pi \delta `$ to leading order. The Euclidean geometry is now conical with deficit angle $`\delta `$ and with a curvature singularity at the tip of the cone.

One of the difficulties in evaluating the partition function for strings on a cone with arbitrary deficit angle is that the cone is not a solution of the string equations of motion. The Einstein equation is not satisfied unless there is an explicit source at the tip to account for the curvature. One would thus require an off-shell formulation of string theory to evaluate this partition function.

We will instead proceed differently. We will consider string theory on a cone obtained as an orbifold $`𝐂/𝐙_N`$ in Dabh ; DabhI ; LoSt . This corresponds to the temperature $`T=\frac{N}{2\pi }`$ at the stretched horizon and the deficit angle is $`\delta =2\pi (1\frac{1}{N})`$. Because the orbifold is a conformal field theory, for these special values of the deficit angle labeled by an integer, the tree-level string equations of motion are indeed satisfied. In addition to the bulk modes, there are states from $`N1`$ twisted sectors labeled by $`k=1,\mathrm{},N1`$. We will take both $`k`$ and $`N`$ to be odd to simplify the subsequent discussion. The ground state in each sector is tachyonic and its mass is given by $`\alpha ^{}m^2=2(1k/N)`$. In the next section we evaluate the free energy $`F(N)`$ as a function of $`N`$ by properly taking into account infrared divergences due to the tachyons and then compute the entropy analytically continuing in $`N`$.

Tachyon Condensation and Black Hole Entropy

The one-loop partition function for strings on the orbifold $`𝐂/𝐙_N`$ can be easily written down. It is modular invariant and therefore ultraviolet finite as expected for a string partition function. This is already an improvement over the field theory calculation which was UV divergent. However the partition function has severe infrared divergences because of the tachyons.

Infrared divergences, unlike ultraviolet divergences, are not a matter of renormalization but signify important physics. In this context, the existence of tachyons in the spectrum suggests that the thermal ensemble is unstable. It was speculated in Dabh ; DabhI  that these tachyons can condense causing a phase transition and the latent heat of this transition could account for the tree level black hole entropy. However, in the absence of a good candidate for the endpoint of tachyon condensation, it was not possible to pursue this further.

Recent work of Adams, Polchinski, and Silverstein APS  has provided new insights about condensation of these tachyons. They have argued that tachyon condensation relaxes the cone to flat space. The most convincing evidence for this claim comes from the geometry seen by a D-brane probe in the sub-stringy regime. In the probe theory, one can identify operators with the right quantum numbers under the quantum $`𝐙_N`$ symmetry of the orbifold that correspond to turning on tachyonic vevs. By selectively turning on specific tachyons, the quiver theory of the probe can be ‘annealed’ to successively go from the $`𝐙_N`$ orbifold to lower $`𝐙_{N2}`$ orbifold all the way to flat space. The deficit angle seen by the probe in this case changes appropriately from $`2\pi (1\frac{1}{N})`$ to $`2\pi (1\frac{1}{N2})`$.

These results are consistent with the assumption that in the field space of tachyons there is a potential $`V(T)`$ where we collectively denote all tachyons by $`T`$. The $`𝐙_N`$ orbifold sits at the top of this potential, the various $`𝐙_M`$ orbifolds with $`M<N`$ are the other extrema of this tachyonic potential, and flat space is at the bottom of this potential. Such a potential can also explain why a conformal field theory exists only for special values of deficit angles. We will be concerned here with the static properties such as the endpoint of tachyon condensation and the effective height of the tachyon potential and not so much with dynamical details of the process of condensation.

Let us return now to the computation of the black hole entropy. We would like to evaluate the free energy of a $`𝐙_N`$ orbifold as a function of $`N`$. Now, the existence of tachyons in the orbifold implies that we are expanding the string field theory functional integral around a maximum. The asymptotic expansion provided by string perturbation theory around this point is as a result IR divergent and essentially useless. To correctly evaluate the partition function, we must expand around the stable saddle point at the minimum of the potential. The leading semiclassical contribution to the partition function will be given by $`Z\mathrm{exp}S_E`$ where $`S_E`$ is the classical Euclidean action after condensation. To elaborate this point let us consider a toy model in field theory of a single scalar field $`\varphi `$ with double well potential $`m^2\varphi ^2+g^2\varphi ^4`$. The perturbative expansion for the partition function around $`\varphi =0`$ is IR divergent which signifies that we have expanded around the wrong saddle point. The stable saddle point is at $`\varphi ^2=\frac{m^2}{2g^2}`$ and the leading semiclassical contribution to the partition function will be given by $`Z\mathrm{exp}S_E`$ where $`S_E=\frac{m^4}{4g^2}`$ is the change in the classical Euclidean action after condensation.

It may seem difficult to evaluate the change in the classical action between the $`𝐙_N`$ orbifold and flat space but we are helped by the fact that, for the orbifold conformal field theory, the equations of motion for the dilaton and the graviton are satisfied exactly. To extract this information, let us consider the Lorentzian string effective action, for concreteness, first to leading order in $`\alpha ^{}`$

$`S={\displaystyle \frac{1}{16\pi G}}{\displaystyle _M}\sqrt{g}e^{2\varphi }[R+4(\varphi )^2\delta ^2(x)V(T)]`$
$`+{\displaystyle \frac{1}{8\pi G}}{\displaystyle _M}\sqrt{g}e^{2\varphi }K,`$ (6)
where $`K`$ is the extrinsic curvature and $`\delta ^2(x)V(T)`$ denotes the tachyon potential localized at the tip of the cone. The extrinsic curvature term is as usual necessary to ensure that the effective action reproduces the string equations of motion for variations $`\delta \varphi `$ and $`\delta g`$ that vanish at the boundary.

The action is very similar to the one for a cosmic string in four dimensions. The tachyon potential supplies an 8-brane source term for gravity. Einstein equations imply $`R=\delta ^2(x)V(T)`$ and therefore a conical curvature singularity at $`x=0`$. Because of this equality, there is no source term for the dilaton and as a result the dilaton equations are satisfied with a constant dilaton. We see that the bulk contribution to the action is zero for the solution. The boundary has topology $`𝐑^8\times 𝐒^1`$. For a cone, the circle $`𝐒^1`$ has radius $`r`$ but the angular variable will go from $`0`$ to $`\frac{2\pi }{N}`$. The extrinsic curvature for the circle equals $`1/r`$ and thus the contribution to the action from the boundary term equals $`\frac{A}{4GN}`$. We have to remember a factor of $`i`$ in Euclidean continuation of $`\sqrt{g}`$.

In the conformal field theory, we should worry about the higher order $`\alpha ^{}`$ corrections to the effective action. These corrections are dependent on field redefinitions or equivalently on the renormalization scheme of the world-sheet sigma model. However, the total contribution of these corrections to bulk action must nevertheless vanish for the orbifold because we know that the equations of motion of the dilaton are satisfied with a constant dilaton which implies no source terms for the dilaton in the bulk. Thus, the entire contribution to the action comes from the boundary term even when the $`\alpha ^{}`$ corrections are taken into account and we can reliably calculate it in a scheme independent way using the conical geometry of the exact solution at the boundary.

From the change in the classical action between the $`𝐙_N`$ orbifold and flat space, we can calculate the leading semi-classical contribution to the free energy up to an $`N`$-independent additive constant:

$$F(N)=\frac{A}{4G}\frac{(N1)}{2\pi }$$
(7)
The resulting entropy is

$$S=2\pi \frac{F(N)}{N}=\frac{A}{4G},$$
(8)
which is in precise agreement with the Bekenstein-Hawking entropy. Note that the entropy is independent of $`N`$ even though the latent heat scales with $`N`$.

Unlike the change in the classical action, the height of the tachyon potential is dependent on the renormalization scheme. It seems likely, however, that there exists a particular renormalization scheme of the $`N=2`$ supersymmetric sigma model of the cone in which the beta function equations are reproduced exactly by (6) to all orders in $`\alpha ^{}`$ and not just to leading order. This would be analogous to Calabi-Yau compactifications where one can argue that there always exists a scheme in which the metric is Ricci flat NeSe . In such a renormalization scheme, the dilaton would be constant, the metric would be precisely conical as for the exact solution of the orbifold. In this scheme, the height of the tachyon potential $`\mu _N`$ for the $`𝐙_N`$ would then be directly related to the deficit angle $`\delta `$ of the solution using Einstein equations from (6) would be given by

$$\mu _N=2\delta =4\pi (1\frac{1}{N}).$$
(9)
It would be interesting to verify this prediction in some form of closed string field theory if the scheme-dependence can somehow be taken into account.

Conclusions and Discussion

We end with a few comments and open problems.

The calculation here may seem similar to the calculation of Gibbons and Hawking or Teitelboim GiHa ; Teit  where the entropy comes from the classical gravitational action. However, conceptually, it is fundamentally different. The Gibbons and Hawking calculation in canonical gravity is in an effective low energy theory and is insensitive to the high energy modes. In canonical gravity, a cone is not a saddle point. One can introduce energy density $`\mu `$ as a chemical potential for holding the area of the horizon fixed and then the cone of arbitrary deficit angle can be a saddle point Teit ; Suss . But this chemical potential is not related to the potential of any dynamical field. There are no IR divergences in these calculations and the UV divergences simply renormalize Newton’s constant in the effective theory.

Somehow, string theory, as a more complete theory of gravity, automatically supplies additional tachyonic fields at the Euclidean horizon. The string partition function is a functional integral over the bulk fields such as the metric as well as the twisted fields including the tachyons and thus is clearly different from the Gibbons and Hawking path integral. The cone is now a saddle point of the full string equations of motion only for integer $`N`$ but without having to introduce any chemical potential by hand. The UV finite partition function is now IR divergent, which implies a dynamical instability. One of the main motivations of this paper is to take seriously the dynamics of the new fields supplied by string theory. This dynamics at the tip seems to capture the physics of the heat bath close to the horizon. The entropy computed this way is universal in that it does not depend on what bulk theory we started with but is determined entirely by the potential of the tachyonic degrees of freedom near the horizon.

The UV divergence of entropy in field theory is intimately related to the puzzle of loss of information in black hole evaporation. If the entropy does have a statistical interpretation in terms of counting of states, then its divergence would suggest an infinite number of states associated with a finite mass black hole. As long as the black hole has an event horizon, it can apparently store an arbitrary amount of information in terms of correlations between the outgoing radiation and the high energy modes near the horizon. When the horizon eventually disappears, the information in these correlations is irretrievably lost. Finiteness of the entropy as we have found both in UV and IR, on the other hand, implies unitary evolution without information loss.

The tachyons at the tip of the cone are very similar to the thermal tachyons in string theory at finite temperature that come from the winding modes around the Euclidean time direction Sath ; Koga . They cannot be interpreted as states in the spectrum but are rather new order-parameter fields. In the Lorentzian continuation, the tachyon condensate near the tip of the cone would seem to imply a Hagedorn-like phase near the Lorentzian horizon reminiscent of the ‘membrane paradigm’ TPM . Apparently, the entropy of this phase accounts for the black hole entropy. A useful analogy is the deconfinement transition in large $`N`$ QCD. It is as though we are able to identify the order parameter for the deconfinement transition in the form of tachyons and compute the latent heat of the transition precisely from the potential to find a tree level contribution of order $`N^2`$ that is suggestive of gluon degrees of freedom. It would be desirable to figure out in a Hamiltonian formalism, the analog of the gluon degrees of freedom. Since the thermal entropy accounts for the entire tree-level entropy, it suggests that the Einstein action is wholly induced by the interactions of the UV degrees of freedom. Matrix theory BFSS  or a gauge theory dual AGMO  may be the right framework to address this question. It is interesting to explore the tachyon landscape irrespective of the problem of black hole entropy and to verify the form of the potential in a closed string field theory. It would also be interesting to pursue the implications of these results for the Hagedorn transition in string theory AtWi  as well as for holography Hoof1 ; Suss1 ; AGMO .

Acknowledgments

I am grateful to Rajesh Gopakumar, Kentaro Hori, and Sandip Trivedi for many valuable discussions.


text:
# Local Simulation Algorithms for Coulomb Interactions

## Abstract

Long ranged electrostatic interactions are time consuming to calculate in molecular dynamics and Monte-Carlo simulations. We introduce an algorithmic framework for simulating charged particles which modifies the dynamics so as to allow equilibration using a local Hamiltonian. The method introduces an auxiliary field with constrained dynamics so that the equilibrium distribution is determined by the Coulomb interaction. We demonstrate the efficiency of the method by simulating a simple, charged lattice gas.

The electrostatic interaction between two point charges in a medium with uniform dielectric constant $`ϵ_0`$ varies as $`e_1e_2/4\pi ϵ_0r`$. The large numerical value of this energy together with its long range are such that it is very often the most costly component in the simulation of charged condensed matter systems. Naive evaluation of the electrostatic energies in molecular dynamics and Monte-Carlo algorithms leads to inner loops where the summation over all pairs takes a time which scales as $`O(N^2)`$ for a single step in which all $`N`$ particles are updated.

Many methods are used to improve this poor scaling: The optimized Ewald algorithm splits the summation between real and Fourier space and has a complexity of $`O(N^{3/2})`$ Beckers et al. (1988); Perram et al. (1988). By interpolation of the charge distribution onto a grid, fast Fourier transform methods allow a scaling in $`O(N\mathrm{log}(N))`$ Essmann et al. (1995). Finally, a popular method in very large simulations is an expansion of the charge distribution using hierarchical multipoles Barnes and Hut (1986); Greengard and Rokhlin (1987). The asymptotic improvement in efficiency comes, however, with great increases in the complexity of the coding, especially when distributed on multiprocessor computers. The numerical prefactors in these scaling laws are uncomfortably high: Despite the great effort put into optimizing the electrostatic loop, it is found that in the simulation of a large biomolecule (with $`N10^5`$) the great majority of the CPU time is still used in the Coulomb loop Schlick et al. (1999) in even the most sophisticated numerical codes. Most of these “fast” methods can only be used efficiently in molecular dynamics simulations; there are many occasions where one would like to perform efficient Monte-Carlo simulations due to the stability and simplicity of the method.

The classical methods for treating charged systems have another disadvantage, their inability to treat systems with inhomogeneous dielectric constants. Dielectric inhomogeneities have drastic effects on material properties. For instance the dielectric contrast between water and the core of proteins leads to expulsion of counter-ions from a $`3A`$ thick hydration layer Perutz (1979). To treat charging effects in proteins quite arbitrary, uncontrolled approximations are made Solmajer and Mehler (1991) on effective electrostatic interactions in the vicinity of a protein in order to reduce interactions to effective pair-wise additive potentials. Similarly much work Kantor et al. (1992) has been performed on the phase structure of charged synthetic polymers while neglecting the large dielectric contrasts between water based solvents and oily backbone structures which are surely important in the discussion of the stability of the necklace structures predicted in these systems. At present the most promising algorithms are based on the non-local Marcus energy functional Marcus (1956); Marchi et al. (2001).

This letter introduces a local algorithm with a propagating field $`𝐄`$ with purely local dynamics on an interpolating grid; it has a complexity in $`O(N)`$ and is elementary to implement. In contrast to conventional grid methods we do not solve for all the field variables at each integration step; we let the field evolve with its own intrinsic dynamics. We were motivated by the observation that Maxwell’s equations, which are local, produce Coulomb interactions due to the propagation of a vectorial field. These dynamic equations are, as we shall see, not the only dynamic way of generating the Coulomb interaction. Our method allows a direct, local implementation of dielectric inhomogeneities. As a demonstration of the method we present an explicit implementation of a local Monte-Carlo algorithm for a charged lattice gas. We note that techniques which interpolate charge degrees of freedom onto a lattice are already very well understood; they form part of standard packages such as Amber Essmann et al. (1995).

We proceed by showing that the Coulomb interaction can be derived from a constrained variational problem. We then show that the constraint equations are solved locally if we allow electric fields which have both gradient and rotational degrees of freedom. This freedom can be used to produced a local Monte-Carlo algorithm. Finally we present a numerical verification of the method.

The energy of a system of charged particles in a uniform dielectric background is expressed as a function of the electric field $`𝐄`$

$$𝒰=ϵ_0\frac{𝐄^2}{2}\mathrm{d}^3𝐫$$
(1)
where the electric field is constrained by Gauss’s law

$$\mathrm{div}𝐄\rho /ϵ_0=0$$
(2)
It is known from classical electrostatics that one solution of equation (2) is given by $`𝐄=\mathrm{grad}\varphi `$ so that $`^2\varphi =\rho /ϵ_0.`$ The general solution to the constraint eq. (2) is thus

$$𝐄=\mathrm{grad}\varphi +\mathrm{curl}𝐐$$
(3)
where $`\varphi `$ is unique to within an additive constant and $`𝐐`$ arbitrary. In Fourier space the electric field can be expressed as $`𝐄(𝐤)=i𝐤\varphi +i𝐤𝐐`$. The second term of this expression is perpendicular to $`𝐤`$ so that there are two physical degrees of freedom in the $`𝐐`$ field, corresponding to two independent polarization states. We can consider that the field is due to a static potential plus transverse photons.

Let us study the stationary states to the variational problem posed by eq. (1) and eq. (2) by using a Lagrange multiplier, with the functional

$$𝒜=\left[ϵ_0\frac{𝐄^2}{2}\lambda (𝐫)(ϵ_0\mathrm{div}𝐄\rho )\right]\mathrm{d}^3𝐫,$$
(4)
implying that $`𝐄+\mathrm{grad}\lambda =0`$. The Lagrange multiplier is identical to the static electrostatic potential, $`\varphi `$ and the minimum energy is $`𝒰_{Coulomb}=\frac{ϵ_0}{2}(\mathrm{grad}\varphi )^2\mathrm{d}^3𝐫`$. Consider the energy eq. (1) for an arbitrary $`𝐄`$ satisfying the constraints then

$`𝒰={\displaystyle \frac{ϵ_0}{2}}{\displaystyle \left[(\mathrm{grad}\varphi )^2+(\mathrm{curl}𝐐)^2\right]\mathrm{d}^3𝐫}`$ (5)
Cross terms vanish, as is shown by integrating by parts.

We now turn to the statistical mechanics of a field with the energy of eq. (1) constrained by Gauss’s law. We do not impose that the electric field is calculated from a potential. The partition function of a fixed set of charges in presence of a fluctuating field $`𝐄`$ is given by

$$𝒵(\{𝐫\})=𝒟𝐄e^{{\scriptscriptstyle {\scriptscriptstyle \frac{\beta ϵ_0𝐄^2}{2}}\mathrm{d}^3𝐫}}\underset{𝐫}{}\delta (\mathrm{div}𝐄\rho (𝐫)/ϵ_0)$$
(6)
The argument $`\{𝐫\}`$ denotes the fact that the integral over the particle positions has not yet been performed. There are two ways of treating this equation. Either one introduces an integral representation of the delta function or one notes that the integral over the field $`𝐄`$ decomposes into a (unique) gradient term and a (non unique) rotation so that

$$𝒵(\{𝐫\})=e^{\frac{\beta ϵ_0}{2}{\scriptscriptstyle (\mathrm{grad}\varphi )^2\mathrm{d}^3𝐫}}𝒟𝐄_te^{\frac{\beta ϵ_0}{2}{\scriptscriptstyle 𝐄^2\mathrm{d}^3𝐫}}$$
(7)
where $`𝒟𝐄_t=_𝐫\delta (\mathrm{div}𝐄)𝒟𝐄`$ performs the summation over all the rotational degrees of freedom of the field described by the potential $`𝐐`$. All the dependence on the particle positions is in the prefactor characterized by the electrostatic potential $`\varphi `$ found by solving Poisson’s equation. This prefactor gives the Coulomb interaction between the particles. The remaining integral is independent of the positions of the charges; rather remarkably integration over the full set of fields allowed by the constraint multiplies the standard partition function by a simple constant. This extra factor in the partition function can be ignored.

In the presence of non-uniform dielectric media one proceeds in a similar manner with the energy

$$𝒰=\frac{𝐃^2}{2ϵ(𝐫)}\mathrm{d}^3𝐫$$
(8)
and the constraint $`\mathrm{div}𝐃\rho =0`$. We deduce that the displacement is given by $`𝐃=ϵ\mathrm{grad}\varphi +\mathrm{curl}𝐐`$ with $`\mathrm{div}(ϵ\mathrm{grad}\varphi )=\rho `$ and

$$𝒵(\{𝐫\})=e^{\frac{\beta }{2}{\scriptscriptstyle ϵ(𝐫)(\mathrm{grad}\varphi )^2\mathrm{d}^3𝐫}}𝒟𝐄_te^{{\scriptscriptstyle {\scriptscriptstyle \frac{\beta }{2ϵ(𝐫)}}𝐄^2\mathrm{d}^3𝐫}}$$
(9)
so that

$$𝒵(\{𝐫\})=𝒵_{Coulomb}(\{𝐫\})𝒵_{fluct}(\{𝐫\})$$
(10)
This time the normalization is a function of the distribution of dielectric inhomogeneities. When implemented in inhomogeneous media our treatment leads to potentials which are the sum of the Coulomb and a fluctuation potential Landau and Lifshitz (1999) which varies as $`1/r^6`$ for two widely separated particles. This term comes from thermally driven dipole-dipole interactions: Fluctuations in the field produce an inhomogeneous polarization, $`𝐏`$, of the dielectric background. This produces an equivalent charge density of $`\mathrm{div}𝐏`$ which interacts via Coulomb’s law. Such fluctuation potentials are to be expected from the Lifshitz theory of dielectrics.

We now propose a lattice version of the above equations suitable for numerically studying the thermodynamics of charged systems. The trick is to use the arbitrary vector potential $`𝐐`$ to simplify the calculation of the updated fields after the motion of a charged particle. We need, also, to sum over all rotational degrees of freedom of the field in order to calculate statistical weights from the partition function eq. (6).

The observation that we use in order to implement a local algorithm is that the constraint equation, eq. (2) can be updated locally in a system in which charge is conserved. We firstly reinterpret the constraint in terms of Faraday’s concept of conserved electric flux: Consider, fig. 1 (top), a network where the charges are confined to the vertices $`\{i\}`$ and the field is associated with links between two sites, $`\{i,j\}`$. Around each lattice point we imagine a cube and write the constraint in integral form $`𝐄.d𝐒=e_i/ϵ_0`$. The integral is over the surface of the cube, $`e_i`$ is the enclosed charged at the site. We use the notation $`E_{1,2}`$ to denote the total flux leaving $`1`$ towards $`2`$; clearly $`E_{1,2}=E_{2,1}`$.

The discretized version of the integral constraint is $`_jE_{i,j}=e_i/ϵ_0`$ The discretized energy is given by

$$𝒰=\frac{ϵ_0}{2}\underset{links}{}E_{i,j}^2$$
(11)
Here and in what follows we assume that the lattice spacing is unity.

Start with a system where the constraint is satisfied, fig. 1 (top), and displace a charge, $`e`$, situated on the leftmost lattice site, $`1`$, to the rightmost site, $`2`$. The constraint is again satisfied at both sites if the field associated with the connecting link is updated according to the rule $`E_{1,2}E_{1,2}e/ϵ_0`$. This is our Monte-Carlo move for the particles, involving a correlated update of a single charge and the field on the link connecting two sites. To update the field configurations fig. 1 (bottom) we update all the field values of a single plaquette while conserving the constraint at each vertex. In fig. 1 (bottom) $`E_{1,2}`$ and $`E_{4,1}`$ increase by an increment $`\mathrm{\Delta }`$ whereas $`E_{4,3}`$ and $`E_{3,2}`$ decrease by $`\mathrm{\Delta }`$ so that at each vertex the sum of the entering and leaving fields is constant. It is this last update that performs the integration over all the rotational degrees of freedom in the $`𝐄`$ field.

The two moves are not quite sufficient to equilibrate a system with periodic boundary conditions in all situations. This problem is linked with the solution $`\varphi =\overline{𝐄}.𝐫`$ or $`𝐄=\overline{𝐄}`$ of the Laplace equation on a torus where $`\overline{𝐄}`$ is an arbitrary constant vector. Motion of the charges generates fluctuations in $`\overline{𝐄}`$, while updates such as those in fig. 1 (bottom) preserve $`\overline{𝐄}`$; similar phenomena also occur with Maxwell’s equations. In order to be absolutely sure that the algorithm is ergodic we introduced a third possible Monte-Carlo step which consists of a shift in $`\overline{𝐄}`$. By keeping track of the evolution of $`\overline{𝐄}`$ as the particles move this last update can be efficiently implemented without destroying the $`O(N)`$ scaling of the algorithm. In the largest systems fluctuations in this single mode should give a small contribution to the thermodynamics if the initial condition is typical; in such cases this update can be eliminated.

We have performed two initial verifications of the algorithm. Firstly we randomly placed four positive and four negative charged particles on a $`4\times 4\times 4`$ lattice with periodic boundary conditions. We performed field updates using the Metropolis algorithm at zero temperature in order to anneal the field $`𝐄`$. We then solved for the electric fields using a standard linear algebra package. The results were identical to within numerical errors. Annealing $`\overline{𝐄}`$ was crucial in order to get agreement between the two methods with frozen charges. A second check was then performed with two different implementations of the Metropolis algorithm now run at several finite temperatures: 36 mutually avoiding, charged particles were distributed in a $`6\times 6\times 6`$ cube with periodic boundary conditions. The dielectric constant was uniform. In the first simulation the linear solver was used to calculate the exact interaction energy of the charged particles at each Monte-Carlo step. In the second simulation we implemented the above local algorithm. For the two methods we then compared the static structure factor for charge-charge correlations, finding good agreement.

We verified the efficiency of the algorithm by determining the autocorrelation time of the slowest density and charge modes in a cube of dimension $`L`$ using the method described in Maggs (2001), fig. (2). We plot the relaxation time in Monte-Carlo sweeps; during each sweep an average of one Monte-Carlo trial is performed for each degree of freedom in the simulation. As expected in a Monte-Carlo algorithm dominated by diffusive motion the slowest mode relaxes in a time which varies as $`L^2`$. There is no anomalous or critical slowing down due to the coupling between the particles and the electric field. The saturation of the relaxation time for charge fluctuations at large $`L`$ is a consequence of screening; charge fluctuations relax by diffusion over the Debye length.

We have shown that the thermodynamics of charged systems can be simulated locally by introducing a propagating vector field so that particles interact via retarded, diffusing fields. The dynamic properties of the system are strongly modified but by construction thermodynamics is an invariant of the propagation dynamics. Our treatment of the dynamics of the field $`𝐄`$ is similar to the Coulomb or radiation gauge in classical electrodynamics: Normally one writes that

$$𝐄=\mathrm{grad}\varphi \frac{𝐀}{t}$$
(12)
where in the Coulomb gauge $`\mathrm{div}𝐀=0`$. There is a rather close analogy between our dynamic scheme of solving for the constraint equations with certain methods of quantization in the Coulomb gauge: In Dirac’s quantization of the electrodynamic field Gauss’s law is a weak identity Dirac (2000) dependent on the choice of the initial wave function. In our simulation Gauss’s law is the result of a restricted choice of possible moves in the Monte-Carlo algorithm together with a special initial condition.

What are the advantages of the present method over direct integration of Maxwell’s equations which are also an example of an $`O(N)`$ algorithm 111We thank the reviewer for this remark? Monte-Carlo algorithms are particularly easy to implement and have good stability with large step sizes. In addition, we have checked that the fast equilibration of the electric degrees of freedom in fig. (2) allows one to perform field updates far more rarely than particle updates leading to additional acceleration of the algorithm. Such multiple time step ideas have been applied to conventional electrostatic solvers but in molecular dynamics are sometimes prone to numerical instabilities Schlick et al. (1999).

In the implementation of the algorithm we were inspired by recent work on hydrodynamic interactions (varying as $`1/r`$) via a Lattice-Boltzmann algorithm in simulations of polymer solutions Ahlrichs et al. (2001). Another analogous problem to ours has been treated by Car and Parrinello who have shown Car and Parrinello (1985) that introduction of a fictitious dynamics leads to much improved efficiencies in solving for constraints; our work differs in that the constraint of Gauss’s law is solved exactly at each simulation step whereas the Car-Parrinello algorithm leads to a simulated annealing solution of the constraint equations (2).

Finally let us note that other functionals do exist for the electric potential in the presence of sources Jackson (1999). In particular the functional

$$\left[\frac{ϵ_0}{2}(\varphi )^2\rho \varphi \right]\mathrm{d}^3𝐫$$
(13)
seems, at first sight particularly simple and attractive. Unfortunately, the minimum of this functional is minus the correct electrostatic energy. It can not be used as a functional for both the field evolution and the particle motion. Application of the algorithm to large atomistic systems remain to be tested, but our method provides an alternative to existing treatments of Coulomb interactions.

We would like to thank R. Everaers for many crucial discussions in the formulation of this work, in particular for his remarks on the importance of dielectric effects.