# A Collaborative Experiment in Cluster Evolution

I have asked a number of colleagues if they would be interested in participating in a comparative study of the dynamical evolution of a model star cluster. Our aim is to study the evolution of a cluster specified by given initial and boundary conditions using the widest possible variety of methods and then to compare the results. The motivation for doing so at this point is that I have to provide a comparative discussion of different methods in a presentation at the IAU General Assembly in Kyoto in August 1997. Another stimulus was M. Lecar's study, about 30 years ago, of 11 numerical integrations of a 25-body problem. I think it would be a good opportunity to show to our observational colleagues how much can be said about the evolution of simple model star clusters.

## SPECIFICATION OF THE PROBLEM

The initial conditions I propose are as follows:

• King Model, W_0 = 3
• non-rotating
• r_t = 30pc
• M = 6x10^4 solar masses
• IMF n(m)dm proportional to m^{-2.35} dm for 0.1 < m < 1.5 (m in solar masses)
• no (primordial) binaries
• no initial mass segregation.

The boundary conditions I propose are as follows:

• Tidal boundary conditions corresponding to circular motion round a point-mass galaxy
• heating by 3-body binaries, modelled as appropriate to your code
• no stellar evolution and no collisions.

To specify the output required, adopt the following notation:

```M	mass within the current tidal radius r_t
t_h 	time at which M falls to half of its initial value
t_cc	time of core collapse
t_end	time at which M = 0
M_cc	value of M at time t_cc
mbar	mean stellar mass (within r_t)
mbar_h	mean stellar mass (within r_h)
mbar_10	mean stellar mass (within r_10)
beta	anisotropy (within r_t)
beta_h	anisotropy (within r_h)
beta_10 anisotropy (within r_10)
```
(Here beta = 1 - <vt**2>/(2<vr**2>), where vt and vr are the radial and transverse components of the velocity of an individual star, and the average is not mass-weighted.)

Then the minimal output required is

```t_h	t_cc	t_end	M_cc
```

and, at both times t_h and t_cc, the quantities

```r_h	r_10
mbar	mbar_h	mbar_10
beta	beta_h	beta_10
```

You may also wish to provide tables of the values of M, r_h, r_10, mbar, mbar_h, mbar_10, beta, beta_h, beta_10 at many times. Also please provide a brief note of the code, hardware, and approximate CPU time, and indicate to me whether or not you wish your output to be posted on the web page.

If you agree to participate in this project please feel free to modify these specifications if this makes them more appropriate to your code. For example, a Fokker-Planck code may adopt a tidal cutoff rather than a tidal field. A Fokker-Planck code will also need to discretise the mass function, and this is something for the experimenter to decide on, as a compromise between accuracy and economy (cf. Chernoff & Weinberg, 1990, ApJ, 351, 121.) An N-body code will probably want to refer radii to the density centre, and you will have to decide how to scale your calculation to a cluster of the above size. In an isotropic code, values of beta may of course be omitted. Another issue for continuum models is the heating term; that adopted by Lee, Fahlman & Richer (1991, ApJ, 366, 455) is a possibility. If you do not wish to include any form of heating, you might just stop at core collapse.

The proposed initial and boundary conditions are modelled in part on the style of Chernoff & Weinberg (1990), with modifications to the IMF and the initial cluster mass, and no stellar evolution. (Preliminary study suggests that the dynamical evolution of the above cluster is far advanced after 20Gyr.) The definition of beta is modelled on Binney & Tremaine, p.204. The quantity mbar is required as a measure of the mass function and mass segregation.

Here is a gzipped file (0.5Mb) of coordinates suitable (after some treatment) for N-body experiments. Each line contains a mass, x,y,z and u,v,w for a single star in an equal-mass King model with W_0 = 3, in the units used by King (1966, AJ, 71, 64-75). Data are given for N = 32768 particles, and the last line contains the tidal radius in King units. With suitable assignment of masses and scaling these data can be used as initial conditions for any N-body integration with N<=32768. Note that this file is not part of the specification of the collaborative experiment; for the N-body experiments reported below the initial conditions were generated independently by each group.

## TIMETABLE

I suggest the following:

Up to 16 May: Please send me any comment on the above proposal. This may indicate whether or not you wish to participate, and you might wish to suggest modifications of this proposal, or other participants. Those of you involved in collaborations may like to agree on a single participant.

By 23 May: I circulate any suggested and agreed alterations of the above proposal.

By 18 July: send me your data.

25 August: results presented in Kyoto at JD15

A version of my presentation, in which the cooperation of all those who send me data will be acknowledged, will be published in the IAU Transactions, and we may wish to consider a more definitive publication elsewhere in due course. I hope all this sounds like worthwhile fun!

The mailing list for this experiment is the following:

```x  	gquinlan@physics.rutgers.edu,
x	takahasi@vega.ess.sci.osaka-u.ac.jp,
x	inagaki@kusastro.kyoto-u.ac.jp,
x	hmlee@uju.es.pusan.ac.kr
x	drukier@pegasus2.astro.indiana.edu,
x	cohn@pegasus2.astro.indiana.edu,
x	spurzem@auriga.ari.uni-heidelberg.de,
x	makino@chianti.c.u-tokyo.ac.jp,
x	steve@zonker.drexel.edu,
x	piet@ias.edu,
x	mirek@auriga.ari.uni-heidelberg.de,
x	mirek_giersz@camk.edu.pl,
x	crwe@auriga.ari.uni-heidelberg.de,
x	chernoff@spacenet.tn.cornell.edu,
x	weinberg@phast.umass.edu,
x	s.f.portegieszwart@fys.ruu.nl,
x	fukushig@chianti.c.u-tokyo.ac.jp
x	vesperin@falcon.phast.umass.edu,
x	sverre@ast.cam.ac.uk
x	teuben@astro.umd.edu
```
Interested replies have been received from those with a cross (x)

# RESULTS

## D.C. Heggie & S.J. Aarseth

```N = 4096, NBODY4
Dec Alpha 3000/700 + 70 GRAPE-4 pipelines, 4.5(3)hrs, 8 runs,
starting 28/5/1997.

Times in Gyr, masses in solar masses, radii in pc; standard
deviation in brackets in units of last decimal place given.

t_h	 t_cc	  t_end	  M_cc
13.5(4)	 8.9(6)	  22.1(5) 43700(1500)

Conditions at t_h:		Conditions at t_cc

r_h	 r_10			r_h      r_10
6.9(4)	 1.4(1)			7.1(3)   1.5(2)

mbar	 mbar_h	  mbar_10	mbar     mbar_h   mbar_10
0.309(2) 0.44(2)  0.82(6)	0.266(4) 0.35(1)  0.72(7)

beta	 beta_h	  beta_10	beta     beta_h   beta_10
-0.10(5) -0.02(6) 0.0(3)	-0.07(3) -0.01(4) 0.2(2)

```
```N = 8192, NBODY4
Dec Alpha 3000/700 + 70 GRAPE-4 pipelines, 16(1)hrs, 4 runs,
from 3/6/1997.

Times in Gyr, masses in solar masses, radii in pc.

t_h	 t_cc	  t_end	  M_cc
12.5(2)	 10.3(7)  20.7(4) 37700(1900)

Conditions at t_h:		Conditions at t_cc

r_h	 r_10			r_h      r_10
6.4(2)	 1.33(6)		6.4(2)   1.2(1)

mbar	 mbar_h	  mbar_10	mbar     mbar_h   mbar_10
0.306(3) 0.428(7) 0.82(4)	0.280(4) 0.38(2)  0.80(6)

beta	 beta_h	  beta_10	beta     beta_h   beta_10
-0.11(2) 0.01(5)  0.0(1)	-0.16(2) -0.08(3) -0.1(1)

```
```N = 16384, NBODY4
Dec Alpha 3000/700 + 70 GRAPE-4 pipelines, 62(2)hrs, 2 runs,
from 9/6/1997.

Times in Gyr, masses in solar masses, radii in pc.

t_h	  t_cc	   t_end	M_cc
11.0(5)	  10.4(10) 18.6(2)	32000(2000)

Conditions at t_h:		Conditions at t_cc

r_h	  r_10			r_h      r_10
6.4(2)	  1.1(1)		6.2(2)   1.02(5)

mbar	  mbar_h   mbar_10	mbar     mbar_h   mbar_10
0.303(5)  0.41(2)  0.85(2)	0.29(1)  0.40(3)  0.822(7)

beta	  beta_h   beta_10	beta     beta_h   beta_10
-0.159(4) -0.07(2) 0.0(1)	-0.18(1) -0.09(3) 0.0(1)

```
```N = 32768, NBODY4
Dec Alpha 3000/700 + 70 GRAPE-4 pipelines, 286hrs, 1 run,
finished 2/7/1997.

Times in Gyr, masses in solar masses, radii in pc.

t_h	t_cc	t_end	M_cc
10.3	10.4	17.1	29000

Conditions at t_h:		Conditions at t_cc

r_h	r_10			r_h     r_10
5.5	1.0			5.6	0.93

mbar	mbar_h	mbar_10		mbar    mbar_h  mbar_10
0.305	0.42	0.85		0.307	0.43	0.88

beta	 beta_h	beta_10		beta    beta_h  beta_10
-0.182	-0.08	0.0		-0.19	-0.03	0.0

```
Graphs of results for N = 32768 Abscissa is time in Myr. Clockwise from top left these show
• Mean masses mbar_10,mbar_h, mbar (from top to bottom)
• Anisotropy beta and beta_h (above)
• Lagrangian radii r_h, r_10 and r_1 (i.e. mass fraction 1%)
• Total mass
• ## G. Quinlan & D.C. Heggie

```G. Quinlan's Isotropic Fokker-Planck code, 8 mass groups,
no energy generation
70 MHz Sparc 5, 10.5hrs, 8/6/1997

Times in Gyr, masses in solar masses, radii in pc

t_h	t_cc	t_end	M_cc
-	12.9	-	31000

Conditions at t_h:		Conditions at t_cc

r_h	r_10			r_h     r_10
-	-			4.9     0.89

mbar	mbar_h	mbar_10		mbar    mbar_h  mbar_10
-	-	-		0.314   0.459    0.869

beta	beta_h	beta_10		beta    beta_h  beta_10
-	-	-		-	-	-
```
```G. Quinlan's Isotropic Fokker-Planck code, 16 mass groups,
no energy generation
70 MHz Sparc 5, 14.5hrs, 10/6/1997

Times in Gyr, masses in solar masses, radii in pc

t_h	t_cc	t_end	M_cc
-	12.8	-	31000

Conditions at t_h:		Conditions at t_cc

r_h	r_10			r_h     r_10
-	-			4.3     0.94

mbar	mbar_h	mbar_10		mbar    mbar_h  mbar_10
-	-	-		0.306   0.547    0.838

beta	beta_h	beta_10		beta    beta_h  beta_10
-	-	-		-	-	-

```

## K. Takahashi

Notes on the code and data are given at the end of this section of data. These data were received 18/7/97
```Isotropic FP code, 8 mass groups
100 MHz HP 9000/715, 0.7 hrs

t_h      11.0
t_cc     10.5
t_end    15.8
M_cc     32800

Conditions at t_h

r_h      5.36
r_10     1.12
mbar     0.313
mbar_h   0.424
mbar_10  0.803

Conditions at t_cc

r_h      5.42
r_10     0.95
mbar     0.299
mbar_h   0.405
mbar_10  0.836

-----------------------------------------------

isotropic FP code, 16 mass groups
100 MHz HP 9000/715, 1.0 hrs

t_h      11.0
t_cc     10.3
t_end    15.7
M_cc     33800

Conditions at t_h

r_h      5.43
r_10     1.16
mbar     0.313
mbar_h   0.423
mbar_10  0.805

Conditions at t_cc

r_h      5.54
r_10     1.00
mbar     0.296
mbar_h   0.398
mbar_10  0.829

-----------------------------------------------

Anisotropic FP code, 8 mass groups, apocenter cutoff condition
HP Exemplar, 26.5 hrs

t_h      14.5
t_cc     10.9
t_end
M_cc     42100

Conditions at t_h

r_h      6.80
r_10     1.38
mbar     0.309
mbar_h   0.436
mbar_10  0.828
beta    -0.34
beta_h  -0.080
beta_10  0.0011

Conditions at t_cc

r_h      6.83
r_10     1.25
mbar     0.271
mbar_h   0.361
mbar_10  0.769
beta    -0.16
beta_h  -0.054
beta_10  0.019

-----------------------------------------------

Anisotropic FP code, 8 mass groups, energy cutoff condition
HP Exemplar, 20.0 hrs

t_h      11.6
t_cc     10.5
t_end
M_cc     34300

Conditions at t_h

r_h      5.81
r_10     1.20
mbar     0.309
mbar_h   0.426
mbar_10  0.818
beta    -0.17
beta_h  -0.056
beta_10  0.0015

Conditions at t_cc

r_h      5.82
r_10     1.02
mbar     0.292
mbar_h   0.396
mbar_10  0.830
beta    -0.13
beta_h  -0.050
beta_10  0.021

-------------------------------------------------

The original version of the isotropic FP code was written by Shogo
Inagaki. I modified it to include energy generation and tidal cutoff.
The anisotropic FP code was written by myself.

I tried two kinds of tidal-cutoff condition in the anisotropic FP
simulations (Takahashi, Lee & Inagaki 1997, astro-ph/9705012). These are

Apocenter condition: stars with r_apocenter > r_t are assumed to escape,
Energy condition:    stars with E > E_t = -GM/r_t are assumed to escape.

In the isotropic FP simulations, only the energy condition was used.

Concerning the Coulomb logarithm ln(gamma*N), I assumed gamma=0.11 in
my simulations to include the binary heating effect.
As you know, the value of gamma also affects the conversion of time unit
from FP simulation unit to real unit. For example, if we adopt gamma=0.015
(Giersz & Heggie, 1996) simply for the time-unit conversion, times t_cc,
t_h, and t_end listed below increase by a factor of 1.24.

Times in Gyr, masses in solar masses, radii in pc

H.-M. Lee

-----------------------------------------------------------------------

1. The code: Isotropic Fokker-Planck code with
tidal cut-off,
three-body binary heating.
time step chosen to suppress gravothermal oscillations

2. Binning of the mass: 14 mass components with 0.1 solar mass interval
from 0.15 to 1.45 solar mass.

3. Computer: Silicon Graphics Indy R4000PC (100 MHz)

4. CPU time: 90 minutes to complete (up to M=0.04 M(0)).
1875 poisson steps,
7500 Fokker-Planck Steps

5. Coulomb logarithm taken as 0.4N

-------------------------------------------------------------------------

Summary Table: (time in Gyr, mass in solar mass, radius in pc)
-------------------------------------------------------------------------

(1) core collapse:
t_cc =9.09, M_tot=3.56E+04
m_b(tot)=0.31  m_b(rh)=0.41  m_b(10%)=0.83  (solar mass)
r_c=0.0096     r_10= 1.03    r_h=5.67      r_t=25.2 (pc)

(2) 0.5 x M(t=0), M_tot=3.00E+04
t_h =  10.02
m_b(tot)=0.34  m_b(rh)=0.45  m_b(10%)=0.82
r_c=0.034      r_10=1.19     r_h=5.55        r_t=23.4

(3) M--> 0
t_end= 14.8  (Gyr)
------------------------------------------------------------------------------

Graphs of results
Abscissa is time in yr.  Clockwise from top left these show
Mean masses mbar_10,mbar_h, mbar (from top to bottom)
Total mass

C.R.W. Einsel

These data were sent 8 August

code:            	isotropic FP-code with allowance for rotation (E,J_z)
without binary heating terms
number of mass bins:  	8
machine:         	Sparc 20 (hyper Sparc) with 100 Mhz

data derived for / at collapse time

t_cc = 9.55 Gyrs (using ln 0.4N)
=10.77 Gyrs (using ln 0.11N)
M_cc = 40000 M_o
r_h       = 5.91 pc
r_10      = 1.098 pc
m_bar     = 0.287 M_o
m_bar_h   = 0.323 M_o
m_bar_10  = 0.680 M_o

T. Fukushige and J. Makino

Results sent on 11 August

Simulation Model:

Scaling and implimentation of tidal field are the same as that in
Fukushige and Heggie (1995). I used softened potential (eps=1/2048)
for all simulations.  Therefore, I stopped simulations at core
collapse. I performed simulations of N=4k,8k,16k,32k,and 64k.

I scaled the simulation results to real astronomical units using
relation t_rel \propto N/ln[min(1/2eps,N/8)]. This relation is derived
from mbar_h of N=16k simulations with different softening parameters
(eps=1/32, 1/256, 1/2048, and 0).

Code:

The 4-th order Hermite integration scheme
with an individual (hierarchical) time step algorithm.

Hardware:

DEC Alpha 3000/700 + GRAPE-4(9 board, 423 chips)

CPU time:

N=4k: 2.5 hr.
N=8k: 6.5 hr.
N=16k: 21.5 hr.
N=32k: 94 hr.
N=64k: 459 hr.

DATA AT CORE COLLAPSE (masses in solar masses, time in Gyr, radii in pc)

N	4K		8K		16K		32K		64K
time 	16.0465		14.1672		12.5406		11.75		10.7231
M 	36926.3 	35281 		30500.5 	25520.8 	22027.7
r_h 	6.49641 	5.89375 	5.7151 		4.93929 	4.62576
r_10	1.47772 	1.18599 	0.978642 	0.727103 	0.841263
r_1 	0.00888032 	0.151417 	0.0124362 	0.0505861 	0.0146466
mbar 	0.288844 	0.294076 	0.304509 	0.326961 	0.343717
mbar_h 	0.407689 	0.40828 	0.414922 	0.468283 	0.487696
mbar_10 0.728987 	0.866763 	0.857271 	0.948375 	0.93964
beta	-0.189298 	-0.171871 	-0.16977 	-0.171304 	-0.206097
beta_h	-0.0803268 	-0.066344 	-0.0527167 	-0.0671602 	-0.0686604
beta_10	0.279489	0.0968819	-0.0391797	-0.152221	0.0375334

M. Giersz

Data received 11 August.  Data are not completely checked.

Code: 		Monte Carlo Anisotropic Fokker-Planck, with apocentre escape criterion;
123720 shells;
energy generation using the standard Spitzer formula
with mass dependence as in eq 6-15 in his book with free factor which
gives correct result for equal masses
Hardware:	HyperSPARC 125 MHz
CPU Time:	166 hours

Times in Myr, masses in solar masses, radii in pc (gamma = 0.01).
time of core bounce determined from minimum of central potential.

time, total mass, r_10, r_h, r_tid, mbar_10, mbar_h, mbar_t, beta_10, beta_h, beta_t

t_cc:   1.3952E+04  4.6666E+04  1.3265E+00  7.0273E+00 2.7588E+01  7.5849E-01  3.4605E-01  2.6345E-01  5.1347E-03 -4.8672E-02 -8.5078E-02

t_h:    2.0928E+04  2.9948E+04  1.2472E+00  7.1098E+00 2.3797E+01  8.7243E-01  4.2504E-01  3.0465E-01  3.1935E-01  4.9457E-02 -1.0701E-01

t_end:  3.2989E+04  7.2822E+02  4.9012E-01  1.7426E+00  6.8942E+00  1.0098E+00  8.9202E-01  5.4183E-01 -1.0594E+00 -6.3682E-01 -3.1506E-01

Fig. 1 Anisotropy
Fig. 2 Average Masses
Fig. 3 Density
Fig. 5 Total Mass

Data received 1 October, corrected 13 October 1997

Code: 		As above, except with modified escape criterion (energy
and reduced overall time step;
Hardware:	2xHyperSparc 125MHz
CPU Time:	16 days

Times in Myr, masses in solar masses, radii in pc (gamma = 0.01).

time, total mass, r_10, r_h, r_tid, mbar_10, mbar_h, mbar_t, beta_10, beta_h, beta_t

t_cc: 1.3283E+04  4.5166E+04  1.2696E+00  6.9965E+00 2.7290E+01  7.7511E-01  3.5319E-01  2.6579E-01  4.1847E-02 -2.1978E-02 -8.5073E-02

t_h:  1.8970E+04  2.9990E+04  1.1973E+00  7.0305E+00 2.3809E+01  8.8929E-01  4.2210E-01  3.0267E-01  1.6082E-02 -6.9372E-02 -2.2761E-01

t_end:2.5998E+04  1.5184E+03  2.8041E-01  1.8342E+00 8.8077E+00  1.3001E+00  1.1456E+00  9.9762E-01 -2.5645E-01  9.4453E-02  2.9716E-02

Fig. 1 Anisotropy
Fig. 2 Average Masses
Fig. 3 Density
Fig. 5 Total Mass

Data received 13 October 1997, completed 17/10/97

Code: 		As above, except with corrected modified escape criterion
(energy above -M_tot/r_tide and radius above the tidal
radius"), and reduced overall time step;
Hardware:	Pentium II 266 MHz CPU
CPU Time:	180 hrs (incomplete)

Times in Myr, masses in solar masses, radii in pc (gamma = 0.01).

time, total mass, r_10, r_h, r_tid, mbar_10, mbar_h, mbar_t, beta_10, beta_h, beta_t

t_cc: 1.2842E+04  4.5128E+04  1.2678E+00  6.9820E+00 2.7283E+01  7.7400E-01  3.4764E-01  2.6445E-01  5.1925E-02 -1.5848E-02 -7.3359E-02

t_h: 1.8935E+04  2.9987E+04   1.2196E+00  7.0348E+00 2.3808E+01  8.6502E-01  4.1715E-01  3.0072E-01  4.7681E-02 -4.2881E-02 -1.9706E-01

t_end: 2.5902E+04  2.4619E+03  4.2470E-01  2.4169E+00 1.0347E+01  1.2046E+00  1.0625E+00  8.5422E-01 -3.0106E-01  2.5209E-02 -1.4839E-01

Fig. 1 Anisotropy
Fig. 2 Average Masses
Fig. 3 Density
Fig. 5 Total Mass

Code: 		As above, with corrected modified escape criterion
(energy above -M_tot/r_tide and radius above the tidal
Hardware:	Pentium II 266 MHz CPU
CPU Time:	591 hrs (incomplete)

Times in Myr, masses in solar masses, radii in pc (gamma = 0.01).

Collapse time

time        mass       density    potential     r_10         rh
r_t        _10      _h       _t        beta_10    beta_h
beta_t

1.2813E+04  3.9388E+04  4.2076E+06 -4.9784E+00  1.1224E+00  6.2009E+00
2.6073E+01  8.1120E-01  3.7336E-01  2.7938E-01  3.0240E-02 -3.6693E-02
-6.0971E-02

Half-mass time

time        mass       density    potential     r_10         rh
r_t        _10      _h       _t        beta_10    beta_h
beta_t

1.5774E+04  3.0004E+04  1.7497E+05 -2.7731E+00  1.1705E+00  5.9380E+00
2.3811E+01  8.4444E-01  4.2368E-01  3.0696E-01  3.6963E-02 -2.8838E-02
-7.7859E-02

End time (mass not 0)

time        mass       density    potential     r_10         rh
r_t        _10      _h       _t        beta_10    beta_h
beta_t

2.2209E+04  3.0847E+03  7.3392E+04 -7.9325E-01  2.7739E-01  2.3873E+00
1.1156E+01  1.2337E+00  1.0034E+00  8.1978E-01 -2.0359E-02  3.0083E-03
-5.3372E-02

S.L.W. McMillan et al

Data received 11 August.  Data provisional.

Authors:	S. McMillan, J. Makino,	P. Hut,	S. Portegies
Zwart, K. Engle
Software: 	a development version of kira, the Starlab N-body integrator,
which includes rudimentary handling of multiple systems.
Hardware:	500 MHz Alpha system with a 47-pipeline
GRAPE-4 board.
Timings:	9.5 hr for 4k, 43 hr for 10k, 154 hr for 20k.

Scaling assumes Coulomb logarithm = log(0.025N)
------------------------------------------------------------------------

N = 4096, self-consistent tidal field, stars removed at the Jacobi radius.

t_h              t_cc             t_end               M_cc

14.6 Gyr         10.0             26.7                4.22e4 M_sun

Conditions at t_h:                           Conditions at t_cc

r_h         r_10                              r_h         r_10

6.43 pc     1.58                              6.42        1.40

mbar        mbar_h      mbar_10               mbar        mbar_h       mbar_10

0.312 M_sun  0.426       0.603                0.270       0.354        0.724

beta        beta_h      beta_10               beta        beta_h       beta_10

-0.182      -0.150      -0.330                -0.142      -0.110       -0.119

-------------------------------------

N = 4096, no tidal field, stars removed at the Jacobi radius.

t_h              t_cc             t_end               M_cc

18.9 Gyr         11.0             35.7                4.76e4 M_sun

Conditions at t_h:                           Conditions at t_cc

r_h         r_10                              r_h         r_10

6.69 pc     1.57                              6.93        1.59

mbar        mbar_h      mbar_10               mbar        mbar_h       mbar_10

0.310 M_sun 0.443       0.803                 0.263       0.346        0.666

beta        beta_h      beta_10               beta        beta_h       beta_10

-0.0791     0.0627      0.256                 -0.0453     0.00562      -0.114

-------------------------------------

N = 4096, self-consistent tidal field, stars removed at twice the Jacobi radius.

t_h              t_cc             t_end               M_cc

15.5 Gyr         11.1             31.3                4.26e4 M_sun

Conditions at t_h:                           Conditions at t_cc

r_h         r_10                              r_h         r_10

7.09 pc     1.60                              6.43        1.48

mbar        mbar_h      mbar_10               mbar        mbar_h       mbar_10

0.303 M_sun 0.432       0.772                 0.267       0.375        0.625

beta        beta_h      beta_10               beta        beta_h       beta_10

-0.0356    -0.0267      -0.0349               -0.0654     -0.0346      -0.478

-------------------------------------

N = 10240, self-consistent tidal field, stars removed at twice the Jacobi radius.
(average of two runs)

t_h              t_cc             t_end               M_cc

12.5 Gyr         10.1             25.0                3.87e4 M_sun

Conditions at t_h:                           Conditions at t_cc

r_h         r_10                              r_h         r_10

6.81 pc     1.23                              6.57        1.18

mbar        mbar_h      mbar_10               mbar        mbar_h       mbar_10

0.301 M_sun  0.417      0.880                 0.276       0.366        0.795

beta        beta_h      beta_10               beta        beta_h       beta_10

-0.109      -0.0656     0.119                -0.153      -0.114       0.008

-------------------------------------

N = 10240, no tidal field, stars removed at the Jacobi radius.

t_h              t_cc             t_end               M_cc

17.0 Gyr         10.8             30.0                4.63e4 M_sun

Conditions at t_h:                           Conditions at t_cc

r_h         r_10                              r_h         r_10

06.90 pc    1.28                              6.95        1.29

mbar        mbar_h      mbar_10               mbar        mbar_h       mbar_10

0.305       0.449       0.867                 0.260       0.349        0.731

beta        beta_h      beta_10               beta        beta_h       beta_10

-0.225      -0.065       0.190               -0.0944     -0.0519      -0.0561

-------------------------------------

N = 20480, self-consistent tidal field, stars removed at the Jacobi radius.

t_h              t_cc             t_end               M_cc

10.8 Gyr         9.72             20.1                3.41e4 M_sun

Conditions at t_h:                           Conditions at t_cc

r_h         r_10                              r_h         r_10

5.67 pc     1.06                              5.86        0.950

mbar         mbar_h     mbar_10               mbar        mbar_h       mbar_10

0.304 M_sun  0.419      0.857                 0.288       0.392        0.862

beta        beta_h      beta_10               beta        beta_h       beta_10

-0.0170     -0.0528     -0.0315               -0.163      -0.0652      -0.0638

R. Spurzem

Software:	Anisotropic Gaseous Model Code, 8 mass groups
Hardware:	Ultra SPARC 166MHz
CPU time:	13 hours (On CRAY J90: 16 hours)

A model with 16 mass groups is in progress (75 hrs on UltraSparc, 50 on
J90).  Data are given only for the model with 8 mass groups.  All CPU
times are given until the mass is about 50% of the initial total mass.

See more comments at the end.

Quantities without subscript refer to 90 % of total mass,
(100% of mass) is given.
Gravothermal oscillations occurred and the run is still in progress.
Therefore, conditions at t_end cannot yet be given here.

Data are given for two runs, depending on the escape criterion.  The
first is for the standard energy criterion, the second is for an
"apocentric" criterion.  (See notes at end.)

Escape:	Energy	Apocentre

t_h      11.79	11.79
t_cc     9.41	9.77
t_end    13.4	15.43
M_cc     40560	44820

Conditions at t_cc

r_h      6.52	6.93
r_10     1.30	1.42
mbar     0.277	0.269
mbar_h   0.359	0.344
mbar_10  0.839	0.805
beta     0.186	0.148
beta_h   0.138	0.137
beta_10  0.115	0.121

Conditions at t_h

r_h      6.03	6.31
r_10     1.36	1.35
mbar     0.307	0.308
mbar_h   0.409	0.412
mbar_10  0.858	0.883
beta     0.226	0.179
beta_h   0.165	0.152
beta_10  0.128	0.124

Graphs of results for anisotropic gaseous model (ps files,
Energy escape criterion

Fig. 1 Total Mass
Fig. 2 Average Masses
Fig. 3 Anisotropies
Fig. 5 Central Densities

Abscissa is time in Gyr. Ordinates are in solar masses (Figs. 1,2) and pc
(Fig. 4). Fig. 3 shows beta in dimensionless values as required.

Apocentric escape criterion

Fig. 1 Total Mass
Fig. 2 Average Masses
Fig. 3 Anisotropies
Fig. 5 Central Densities

---------------------------------------------------------------------

Coulomb Logarithm ln(0.11N)

Radial mesh contained initially 139 points. The mesh is fixed,
number of mesh points reduces with time as mass decreases. Note that
some steps in the curves can be due to the discrete events of
reducing the number of meshes. 8 (16) mass bins were selected as
in Chernoff & Weinberg (1990) ranging in its stellar mass
from 0.118 to 1.27 solar masses (mmin=0.1, mmax=1.5 sol. m.) and
giving a total particle number of Ntot=245444
(for 16 comp: 0.109 to 1.38, mmin=0.1, mmax=1.5, Ntot=246939).
Max changes allowed per timestep: 2 per cent.

A Program used to generate initial
King models was kindly supplied by G. Quinlan.

Anisotropic gaseous model code using moment equations of
Boltzmann equation with Fokker-Planck collisional terms,
the latter self-sonsistently derived analytically from an
anisotropic multi-mass background by assuming a local Larson-type distribution
function (Legendre polynomial series up to order 2). Note that such treatment
standard Spitzer equipartition time.

Inclusion of escapers is very new feature and proved to be non-trivial.
Just to let mass stream across the outermost boundary is not the physically
correct treatment and leads to a much too small mass loss. Instead at
each radial shell the rate of escapers has to be computed by using some
approximate distribution function. The incorporation of this mechanism
was done only very recently and benefitted very much from comparisons with
K. Takahashi's equal mass data for W0=3,6,9. It is believed that
the treatment of escapers in the gaseous model
near the outer boundary is critical to the
anisotropy development and not yet a correct description of what happens
in an N-body system. Therefore the anisotropy values at the outer boundary
are still in not a good agreement with other models, whereas core and
core-collapse related quantities are reproduced rather well.

As an alternative to the standard energy escape criterion, results are
also given for an apocentre escape criterion (determining vesc,rad
and vesc,t locally, then adopting a King-like distribution function
in radial and tangential direction with different limiting velocities
to determined local escaper fraction.)

More details will be published soon, presumably in a common paper with
K. Takahashi.

The use of departmental computers at Dept. of Astronomy,
Kyoto Univ., Japan and at Dept. of Earth and Space Science,
Univ. of Osaka, Japan,
is gratefully acknowledged. Furthermore CRAY J90 of
Hoechstleistungs-Rechenzentrum (HLRZ) Juelich, Germany, was used.

Summary of results

Postscript files

Key to Symbols (170K) The order of symbols corresponds to the order in which the models appear under each author.

Time scales  (75K) All models have been reduced to a common value of the Coulomb logarithm

Time and mass at core collapse (90K)

Mean mass at core collapse (36K)

Anisotropy at core collapse (26K)

G.A. Drukier
Software: 	Cohn's isotropic Fokker-Planck code, as modified and described
by Drukier (1995), with stellar evolution turned off.
Hardware:	HP 9000 K260 (180 Mhz).
CPU time:	See below

The models had 25 mass components and the masses of the stars involved
ran from 0.1056 to 1.449 solar masses. I ran two runs, one with small
time steps and one with large time steps to suppress oscillations. The
results for the two runs are similar, except that in the latter case
core collapse is slightly later and shallower. The coulomb logarithm is
ln(0.11N).

Results
(times in Gyr, masses in solar masses, radii in pc)

Short time step
316,210 time steps, 104 CPU-hours to 1% of the initial mass.

core collapse:          Half mass lost:       End:
T       = 10.54            11.74               17.78 (18.0*)
M       = 3.54e+4          3.00e+4
r_10    = 1.07             1.23
r_h     = 5.84             5.71
mbar    = 0.291            0.314
mbar_h  = 0.391            0.425
mbar_10 = 0.814            0.805

Long time step
1375 time steps, 39 CPU-minutes to 2.2% of the initial mass.

core collapse:          Half mass lost:       End:
T       = 10.72            11.77               17.63 (17.9*)
M       = 3.46e+4          3.00e+4
r_10    = 1.08             1.22
r_h     = 5.82             5.69
mbar    = 0.293            0.315
mbar_h  = 0.395            0.425
mbar_10 = 0.816            0.806

*Extrapolated to zero mass.

The amplitude of the oscillations decrease as the mass decreases and
finally stops.  I don't think that the oscillations can really be
classed as gravothermal, as energy input from the binaries continues
throughout the oscillation cycle, and the amplitude is quite small. In
other models we've run, we see the energy input go to zero during the
expansion phase, and this, I think, is the sign of a true gravothermal
oscillation. Further, I don't see any sign of a velocity inversion in
the oscillations, but my set of curves is rather coarsly spaced and I
may have missed a small inversion.

```