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.
The boundary conditions I propose are as follows:
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 r_h half-mass radius (i.e. radius of shell containing 0.5*M) r_10 10% radius 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.
Comments on the specification
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.
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.eduInterested replies have been received from those with a cross (x)
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.0Graphs of results for N = 32768 Abscissa is time in Myr. Clockwise from top left these show
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 - - - - - -
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 ------------------------------------------------- Comments on simulations 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 pcH.-M. Lee
These data were received 22/7/97----------------------------------------------------------------------- 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
Energy escape criterion Fig. 1 Total Mass Fig. 2 Average Masses Fig. 3 Anisotropies Fig. 4 Radii 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. 4 Radii Fig. 5 Central Densities --------------------------------------------------------------------- Comments on simulations: 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 leads to additional anisotropy dependent terms as compared to the 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.
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)
Radii at core collapse (90K)
Mean mass at core collapse (36K)
Anisotropy at core collapse (26K)
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). Further comments at end. 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.
