\documentstyle[proceedings,epsf]{crckapb}
\newcommand{\stt}{\small\tt}
\begin{opening}
\title{GRAVOTHERMAL OSCILLATIONS}
\author{Junichiro Makino }
\institute{Department of Graphics and Information Science, \\
College of Arts and Sciences, University of Tokyo\\
3-8-1 Komaba, Meguro-ku, Tokyo 153, Japan}
% If there are more authors at one institute, you should first
% use \author{...} for each author followed by \institute{...}.
\end{opening}
\runningtitle{Gravothermal Oscillations}
\begin{document}
\begin{abstract}
We present the first clear evidence that the gravothermal oscillation
takes place in $N$-body systems. We performed direct $N$-body
simulations of systems of point-mass particles with particle numbers
from 2,048 to 32,768. In the simulation with 32,768 particles, the
central density shows an oscillation with an amplitude of $\sim
10^3$, which is similar to what was observed in more approximate models
such as a conducting gas sphere and one-dimensional Fokker-Planck
calculations. The amplitude is smaller for a smaller number of
particles. The number of particles in the core at the maximum
contraction is $\sim 10$ for all runs, while the number of particles
at the maximum expansion is about $0.01 N$. For 16,384- and
32,768-body runs, the temperature inversion during the expansion phase is
clearly visible.
\end{abstract}
\section{Introduction}
Gravothermal oscillation was first found by
\citeauthor{Sugimoto1983} \shortcite{Sugimoto1983}, who modelled the
post-collapse evolution of globular clusters using a conducting gas
sphere with artificial energy production. In their model the energy
production is expressed as $\epsilon = C\rho^k,$
where $\epsilon$ is the energy production per unit mass per unit time
and $\rho$ is the density. The power-law index they used is 1 or
2. They found similar oscillations in both cases, for a wide range of
the value of coefficient $C$.
Before they found the oscillation, the standard picture of the
evolution of a globular cluster had been the following two-stage
one. The first stage is the gravothermal collapse, or core collapse,
driven by the gravothermal instability. This collapse leads to a
self-similar evolution, in which the core density reaches infinity and
the core mass goes down to zero in a finite time (Hachisu, {\it et al.}
1978, Lynden-Bell and Eggleton 1980).
In a real $N$-body system, contraction is halted by the energy
production by 3-body binaries. It had been believed that the whole
system expands homologously, driven by the energy generation from
binaries after the contraction is halted \cite{Goodman1984}.
What Sugimoto and Bettwieser (1983) found is that this
homologous expansion is unstable in a way similar to the way in which an
isothermal sphere is unstable. A linearized stability analysis by
Goodman (1987) showed that the expansion is unstable against
thermal perturbation if the coefficient for the energy production is
small. If translated back to the total number of particles,
the system is unstable if $N >7,000$. This result is confirmed also with
Fokker-Planck calculations \cite{Cohn1989}.
To determine if such an oscillation actually takes place in real
globular clusters, we have to perform a direct $N$-body simulation with
a sufficiently large number of particles, since gas models and FP models
have many simplifying assumptions that might make the evolution
completely different. For example, both assume spherical symmetry,
while in $N$-body simulation the core is known to wander around (Makino
and Sugimoto 1987, Heggie {\it et al.} 1994). Both assume
that the energy production by binaries is smooth, while in an $N$-body
system binaries are formed stochastically.
An $N$-body simulation with sufficiently large number of particles has
been impossible, simply because the requirement for computer power
has been excessive. The largest simulation which has been tried on
a general-purpose supercomputer is 10,000-body run by Spurzem and
Aarseth (1996), which does not cover a long enough time after the first
collapse.
In the present paper, we describe the result of $N$-body simulations
with number of particles 2,048--32,768. All simulations were performed
on GRAPE-4, a special-purpose computer for collisional
$N$-body simulation. Our main result is the following. First, for
large $N$, the core density and core mass exhibited an oscillation of
large amplitude. The core mass at maximum expansion is almost
independent of the number of particles, and is in good agreement with
FP or gas model results (1-2\% of the total mass). Second, we
confirmed that the observed oscillation is driven by the gravothermal
instability. There were several long expanding periods without any
energy input. The temperature inversion is visible in such expansion
phases. In addition, the behavior of the core density is strikingly
similar to the result of FP calculations with a stochastic heat
source (Takahashi and Inagaki 1991), suggesting that the mechanism is
the same.
The structure of this paper is the following. In section 2, we
describe the initial model, the numerical method and the computer used. In
section 3 we present the results. Section 4 is for discussion.
\section{Model and Numerical method}
\subsection{Initial models and the system of units}
We followed the evolution of isolated systems of point-mass
particles. For all calculations, we used random realization of the
Plummer model as the initial condition. We used the standard
system of units (Heggie and Mathieu 1986), in which $G=1$, $M=1$, and
$E=-1/4$, where $G$ is the gravitational constant, $M$ and $E$ are the total
mass and the total energy of the cluster. All
particles have the same mass $m=1/N$, where $N$ is the total number of
particles. The half-mass crossing time $t_{hc}$ is
$2\sqrt{2}$ in this unit.
\subsection{Numerical method}
For all calculations, we used NBODY4 (Aarseth 1995), modified for
GRAPE-4 (Taiji et al. 1995). The numerical integration scheme adopted
in NBODY4 is the 4th order Hermite scheme (Makino and Aarseth
1992). It implements the hierarchical (block) timestep algorithm
(McMillan 1986, Makino 1991) to use the GRAPE hardware
efficiently. Close two-body encounters and stable binaries are handled
by KS regularization (Kustaanheimo and Stiefel 1965). Special
treatment for compact few-body subsystems is also possible.
The outputs are taken at intervals of a fixed time, which is some
fraction of the crossing time. We recorded the central density, core
radius, number of particles in the core, the radii of Lagrangian
shells, the velocity dispersion within Lagrangian shells, the binaries
and their binding energies. The core parameters are calculated following
Casertano and Hut (1985).
The accuracy of the time integration is adjusted so that the energy
error between two outputs is smaller than a certain prescribed
value. The value we used is $1\times 10^{-5}\sim 1\times 10^{-6}$
depending on $N$. We required higher accuracy for larger $N$, since
the duration of the simulation is longer. When the energy error is
very large, the program automatically reads the output at the previous
checkpoint and restarts with a reduced accuracy parameter.
\subsection{Hardware}
For all calculations, we used GRAPE-4 (Taiji 1995). The GRAPE-4 is a
special-purpose computer designed to accelerate the $N$-body
simulation using the Hermite integrator and hierarchical timestep
algorithm. The fully configured system has a theoretical peak speed of
1.08 Tflops. The simulations reported in the present paper were
performed while the assembly and testing of the GRAPE-4 system were
under way. Thus the number of processors varies during the
calculation. For most of the 32k particle run, we used one quarter of the
machine which has a peak speed of 270 Gflops.
\section{Result}
\subsection{Core parameters}
\begin{figure}
\begin{center}
\leavevmode
\epsfxsize 9cm
\epsffile{iau174figs/osc.ps}
\caption{The logarithm of the central density plotted as a function of
the scaled $N$-body time. Curves for different values of
$N$ are vertically shifted by 3 units.}
\end{center}
\end{figure}
Figure 1 shows the time evolution of the central density for all runs.
The time is scaled so that the thermal timescale is same for all runs.
The scaling factor is $t_r(1000)/t_r(N) = 212.75\log(0.11N)/N$ (Giersz
and Heggie 1994). The core density shows an oscillation with large
amplitude in calculations with large $N$ ($>$16k). No matter what is
the real nature of this oscillation, it is at least clear that the
core density of the $N$-body system shows oscillation with the
amplitude comparable to that observed in gas models or FP
calculations. For runs with small $N$ (2k and 4k), there are some
oscillation-like features but they are hardly distinguishable from
fluctuations. For large values of $N$, however, the oscillation with
large amplitude is clearly visible.
Note that in figure 1 there is no clear transition from stable
expansion to oscillation or from regular oscillation to chaotic
oscillation, which were observed in gas and FP models
(e.g. \citeauthor{Cohn1989} 1989, Heggie and Ramamani 1989). The
reason is that binaries emit energy intermittently at the formation
time and as a result of binary-single-body interactions.
Takahashi and Inagaki \shortcite{Takahashi1991} incorporated this
stochastic nature of the energy source into their FP model and found
that the core density shows chaotic oscillatory behavior even if the
energy production rate is larger than the critical value. They also
found that the amplitude of the oscillation is smaller for larger
energy input (smaller $N$), which is consistent with the present
result. In fact, it would be difficult to distinguish between the
result of $N$-body calculations and their stochastic FP result, except
that their result is smoother while the central density is low.
\begin{figure}
\begin{center}
\leavevmode
\epsfxsize 11cm
\epsffile{iau174figs/nc.ps}
\caption{The number of particles in the core as the function of the
scaled time, for simulations with 2k, 4k, 16k, and 32k particles}
\end{center}
\end{figure}
Figure 2 shows the evolution of the number of particles in the
core. For all runs, the number of particles at maximum contraction
is of the order of 10, while that at maximum expansion is 1--2\%
of the total number of particles. This result is again in good
agreement with gas models and FP calculations.
\begin{figure}
\begin{center}
\leavevmode
\epsfxsize 6.5cm
\epsffile{iau174figs/ncstat.ps}
\caption{Fraction of time for which the number of particles in the
core is smaller than $N_c$, as a function of $N_c/N$.}
\end{center}
\end{figure}
Figure 3 shows the fraction of time for which the number of particles
in the core is smaller than the value $N_c$ as a function of $N_c/N$,
for the post-collapse phase. For $N > 8192$, the median core mass is
around 0.5--0.6\%. This corresponds to $r_c/r_h \sim 0.01$. For
$N=32,768$, the core mass is somewhat smaller than that for 16k or 8k
runs. This difference is because the 32k run is not long enough.
Figure 3 indicates that the typical core mass for a post-collapse
cluster would be less than $1\%$, if gravothermal oscillation
occurs. One interesting question is whether it is possible to
distinguish a core in gravothermal oscillation from a core
dominated by primordial binaries. The theoretical prediction by
Goodman and Hut (1989) gives $r_c/r_h \sim 0.02$, while $N$-body
simulations by \citeauthor{McMillan1990} (1990) give $r_c/r_h \sim
0.2$. \citeauthor{McMillan1990} (1990) used 1136 particles. If their
$N$-body results can be extrapolated to a larger $N$, cores with
primordial binaries and cores in gravothermal oscillation would be
clearly distinguishable. To obtain a definitive answer, we need to
perform simulations of clusters with primordial binaries using a larger
number of particles than employed by \citeauthor{McMillan1990} (1990).
\subsection{Detailed view of the 32k run}
In the previous section, we gave an overview of the post-collapse evolution of
$N$-body point-mass systems with 2k--32k particles. In this section, we
take a closer look of the 32k-particle simulation to see whether we
can find a direct signature of the gravothermal expansion. It is
generally believed that a long expansion phase without significant
energy input and a temperature inversion during the long expansion
phase are the most direct signatures of gravothermal expansion
(Bettwieser and Sugimoto 1984, McMillan and Engle 1995). In
this section we investigate both of them.
\begin{figure}
\begin{center}
\leavevmode
\epsfxsize 10cm
\epsffile{iau174figs/32kosc.ps}
\caption{ The core radius (top) and the binding energy of binaries
(bottom) as a function of time for the 32k-particle run}
\end{center}
\end{figure}
Figure 4 shows an enlarged view of the time variation of the core
radius as compared with the sum of the binding energies of all
binaries. It is clear that most of the energy is generated when the
core is very small. Five peaks account for most of the energy
production.
Continued expansion without energy generation is considered to be
one of the most direct signatures of gravothermal oscillation. In
figure 4, we clearly see two such expansions, for $t=6700-6740$ and
$t=6860-6920$. Both expansions continue for more than 10 half-mass
crossing times, which is hundreds of the core relaxation time. These
expansions cannot be driven simply by binary heating. If the
expansion were driven only by binary heating, it could not
continue without energy input for a timescale longer than the core
relaxation time.
\newlength{\minitwocolumn}
\setlength{\minitwocolumn}{0.5\textwidth}
\addtolength{\minitwocolumn}{- 0.5 cm}
\begin{figure}
\begin{minipage}[b]{\minitwocolumn}
\epsfxsize 5.5cm
\epsffile{iau174figs/tprof1.ps}
\end{minipage}
\begin{minipage}[b]{\minitwocolumn}
\epsfxsize 5.5cm
\epsffile{iau174figs/tprof2.ps}
\end{minipage}
\caption{Velocity dispersion profiles for (a) contracting and (b)
expanding phases. Each profile is obtained by time averaging over 80
snapshots (10 time units). The time interval between curves is 5 time
units. }
\end{figure}
Figure 5 shows temperature profiles for the contracting and
expanding phases. Near the end of the expanding phase a temperature
inversion of the order of 5\% is clearly visible. Since these
profiles are time-averaged over 10 time units (80 snapshots), the
actual inversion might be somewhat stronger. Note that the temperature
inversion is visible only near the end of long expansion phases also
in gas model and FP calculations (Bettwieser and Sugimoto 1984,
\citeauthor{Cohn1989} 1989).
\begin{figure}
\begin{center}
\leavevmode
\epsfxsize 6.5cm
\epsffile{iau174figs/rhosigma.ps}
\caption{The change of the central density and the central velocity
dispersion. Each data points is time-averaged value over 80
snapshots. Arrows indicate the direction of evolution. }
\end{center}
\end{figure}
Figure 6 shows the relation between the central density and the
central velocity dispersion. The trajectory shows a striking resemblance
to what is obtained by gas-model calculation \cite{Goodman1987}. The
fact that the trajectory shows clockwise rotations means that this
is a refrigeration cycle, in which the central region absorbs the heat
when the temperature is low, and release heat when the temperature is
high (Bettwieser and Sugimoto 1984, Bettwieser 1985). In particular,
the later phase of the large expansions (indicated by the arrow marked
``B'') is nearly isothermal. Therefore this phase is driven by the
heat supply. Since the binding energy of binaries is
unchanged during this phase, the heat is supplied from outside the
core. In other words, the expansion is gravothermal.
\section{Discussion}
We have performed direct $N$-body simulation of the post-collapse evolution
of globular clusters. We confirm that gravothermal oscillation
actually takes place in a point-mass $N$-body system.
Whether real globular clusters undergo gravothermal oscillation or not
is a question which requires further research. If clusters contain
many primordial binaries, even after the core collapse they might
still be burning the primordial binaries. In addition, the effect of
two-body binaries on the evolution of the cluster is still unclear. In
fact, the cross section for binary formation by tidal capture and
for merging is not fully understood yet (Mardling 1995a, 1995b).
The most straightforward way to study the effect of primordial
binaries or two-body capture is direct $N$-body simulation. In
principle, we can put primordial binaries and their evolution into FP
calculation by following the distribution of internal binding energy
of binaries. However, the standard one-dimensional FP calculation,
which assumes an isotropic distribution, is not appropriate for
following the binary population, since most binaries are formed in the
core and their orbits are nearly radial. Thus, we have to solve the FP
equation at least in three dimensions ($E$, $J$ and the binary binding
energy $E_b$). This would require prohibitively large computer power.
One could also use the Monte-Carlo approach, but its result must be
compared with $N$-body simulation anyway.
%K I'M SORRY; I JUST DON'T UNDERSTAND ``by introducing the space of
% binary binding energy.'' PLEASE RESTATE IT DIFFERENTLY.
We have demonstrated that $N$-body simulation is now possible with a
number of particles close to that in real globular clusters, thanks to the
extremely powerful special-purpose computer and its full-time
availability. We are now able to use direct $N$-body simulation to
study various aspects of the evolution of globular clusters.
The present 32k-particle calculation took about three months of CPU
time on 1/4 of GRAPE-4. If we tried to do a similar calculation on a
Cray T90 vector supercomputer, it would have taken several years of
CPU time. This is simply impossible to do on present-day supercomputers.
If we want to finish the calculation in, say, one CPU month, we need a
computer 50--100 times faster than a Cray T90, which will be
available 10 years from now.
The CPU time of three months is still very long. However, for many
simulations, we do not need as many as 32k particles. A 16k-particle
calculation was finished in 2-3 weeks, on 1/8 of GRAPE-4. Thus to
run many simulations of 16k-particle systems is now practical.
If we can continue the development of the special-purpose computer, we
will have a system 100--1000 times faster than the present GRAPE-4 in
the next five years. Such a system will make it possible to run
50--100k-particle simulation routinely, while 500k--1M-particle
simulations will still take months.
\section*{Acknowledgements}
It's a pleasure to thank Daiichiro Sugimoto for his effort on GRAPE
project, and also for many useful discussions and comments on the
manuscript. I'm grateful to Sverre Aarseth for his NBODY4 code, Makoto
Taiji, Toshiyuki Fukushige and Toshikazu Ebisuzaki for developping
GRAPE-4 along with myself, Piet Hut and Steve McMillan for many useful
discussion and Ivan King for comments on the manuscript. This work was
supported by the Grant-in-aid for Specially Promoted Research
(04102002) of the Ministry of Education, Science, and Culture.
\begin{thebibliography}{}
\bibitem[\protect\citeauthoryear{Aarseth}{1985}]{Aarseth1985} Aarseth,
S. J. (1985) in {\it Multiple Time Scales}, eds. J. U. Brackhill and
B. I. Cohen, Academic, New York, p. 377.
\bibitem[\protect\citeauthoryear{
Bettwieser}{1985}]{Bettwieser1985}
Bettwieser, E. (1985)
in {\it Dynamics of Star Clusters, IAU Symposium No. 113},
eds. J. Goodman and P. Hut, Reidel, Dordrecht, p. 219.
\bibitem[\protect\citeauthoryear{
Bettwieser and Sugimoto}{1984}]{Bettwieser1984} Bettwieser,
E. and Sugimoto, D. (1984) {\it MNRAS}{\bf ~208}, ~493.
\bibitem[\protect\citeauthoryear{Casertano and Hut}{1986}]{Casertano1985}
Casertano, S. and Hut, P. (1985) {\it Ap. J.}{\bf ~298}, {~80}.
\bibitem[\protect\citeauthoryear{Cohn {\it et al.}}{1989}]{Cohn1989}
Cohn, H., Hut, P., and Wise, M. (1989) {\it Ap. J.}{\bf ~342}, {~814}.
\bibitem[\protect\citeauthoryear{Giersz and Heggie}{1994}]{Giersz1994}
Giersz, M. and Heggie, D. (1994) {\it MNRAS}{\bf ~268}, {~257}.
\bibitem[\protect\citeauthoryear{Goodman}{1984}]{Goodman1984}
Goodman, J. (1984) {\it Ap. J.}{\bf ~280}, {~298}.
\bibitem[\protect\citeauthoryear{Goodman}{1987}]{Goodman1987}
Goodman, J. (1987) {\it Ap. J.}{\bf ~313}, {~576}.
\bibitem[\protect\citeauthoryear{Goodman and Hut}{1989}]{Goodman1989}
Goodman, J. and Hut, P. (1989) {\it Nature}{\bf ~339}, {~40}.
\bibitem[\protect\citeauthoryear{Hachisu, {\it et al.}}{1978}]{Hachisu1978}
Hachisu, I{.}, Nakada, Y., Nomoto, K., and Sugimoto, D. (1978) {\it
Prog. Theor. Phys.} {\bf ~60}, ~393.
\bibitem[\protect\citeauthoryear{Heggie}{1986}]{Heggie1986}
Heggie, D. (1986) in {\it The Use of Supercomputers in Stellar
Dynamics}, eds. S. McMillan and P. Hut, Springer, New York, p. 233.
\bibitem[\protect\citeauthoryear{Heggie and Ramamani }{1989}]{Heggie1989}
Heggie, D. and Ramamani, N. (1989) {\it MNRAS}{\bf ~237}, {~757}.
\bibitem[\protect\citeauthoryear{Heggie \it{et al.} }{1994}]{Heggie1994}
Heggie, D. , Inagaki, S., and McMillan, S.L.W. (1994) {\it MNRAS}{\bf ~271}, {~706}.
\bibitem[\protect\citeauthoryear{Kustaanheimo and
Stiefel}{1965}]{Kustaanheimo1965} Kustaanheimo, P., and Stiefel,
E. (1965) {\it J. Reine. Angew. Math.}{\bf ~218}, {~204}.
\bibitem[\protect\citeauthoryear{Lynden-Bell and
Eggleton}{1980}]{Lynden-Bell1980} Lynden-Bell, D{.} and Eggleton,
P.P. (1980) {\it MNRAS} {\bf ~191}, ~483.
\bibitem[\protect\citeauthoryear{Makino}{1986}]{Makino1986}
Makino, J. (1986) in {\it The Use of Supercomputers in Stellar
Dynamics}, eds. S. McMillan and P. Hut, Springer, New York, p. 151.
\bibitem[\protect\citeauthoryear{Makino}{1991}]{Makino1991}
Makino, J.. (1991) {\it Publ. Astron. Soc. Japan}{\bf ~43}, {~841}.
\bibitem[\protect\citeauthoryear{Makino and Aarseth }{1992}]{Makino1992}
Makino, J. and Aarseth, S. J. (1992) {\it Publ. Astron. Soc. Japan}{\bf ~44}, {~141}.
\bibitem[\protect\citeauthoryear{Makino and Sugimoto }{1987}]{Makino1987}
Makino, J. and Sugimoto, D. (1987) {\it Publ. Astron. Soc. Japan}{\bf ~39}, {~589}.
\bibitem[\protect\citeauthoryear{Mardling}{1995a}]{Mardling1995a}
Mardling, R. A. (1995a)
{\it Ap. J.} {\bf ~450}, {~722}.
\bibitem[\protect\citeauthoryear{Mardling}{1995a}]{Mardling1995a}
Mardling, R. A. (1995b)
{\it Ap. J.} {\bf ~450}, {~732}.
\bibitem[\protect\citeauthoryear{McMillan}{1986}]{McMillan1986}
McMillan, S. L. W. (1986) in {\it The Use of Supercomputers in Stellar
Dynamics}, eds. S. McMillan and P. Hut, Springer, New York, p. 156.
\bibitem[\protect\citeauthoryear{McMillan and Engle}{1995}]{McMillan1995}
McMillan, S. L. W. and Engle, E. A. (1995) in this volume.
\bibitem[\protect\citeauthoryear{McMillan {\it et al.}}{1990}]{McMillan1990}
McMillan, S. L. W., Hut, P. and Makino, J. (1990)
{\it Ap. J.} {\bf ~362}, {~522}.
%\bibitem[\protect\citeauthoryear{Spitzer}{1987}]{Spitzer1987}
%Spitzer, L. (1987) {\it Dynamical Evolution of Globular
%Clusters}, Princeton University Press, New Jersey.
\bibitem[\protect\citeauthoryear{Spurzem and Aarseth}{1996}]{Spurzem1996}
Spurzem, R. and Aarseth, S. J. (1996) preprint.
\bibitem[\protect\citeauthoryear{Sugimoto and
Bettwieser}{1983}]{Sugimoto1983} Sugimoto, D. and Bettwieser,
E. (1983) {\it MNRAS}{\bf ~204}, ~19p.
\bibitem[\protect\citeauthoryear{Takahashi and
Inagaki}{1991}]{Takahashi1991} Takahashi, K. and Inagaki, S. (1991)
{\it Publ. Astron. Soc. Japan}{\bf ~43}, {\
~589}.
\bibitem[\protect\citeauthoryear{Taiji}{1995}]{Taiji1995}
Taiji, M. (1995) in this volume.
\end{thebibliography}
\end{document}
\begin{tabbing}
{\stt $\backslash$cite\{Smith92\}} \>generates: (Smith {\it et al.}, 1992) \\
{\stt $\backslash$shortcite\{Smith92\}} \>generates: (1992) \\
{\stt $\backslash$citeauthor\{Smith92\}} \>generates: Smith {\it et al.}\\
{\stt $\backslash$citeyear\{Smith92\}} \>generates: 1992
\end{tabbing}
memo:
contraction: 6500-6575
expansion : 6670-6700