How to Run and Interpret NBODY4 Simulations

(Last modified July 2006. Prepared by Sverre Aarseth and compiled by Vicki Johnson.)

Contents:

• Introduction to NBODY4
• Guidance on initial parameters
• Research projects that have used NBODY4
• Suggested simulations
• NBODY4 input parameters
• Suggested input values for experimentation
• Output values printed and plotted
• Quantities to review in the full NBODY4 output printout
• Sample NBODY4 output and interpretation
• N-body units and astrophysical units
• Integration methods
• Stellar evolution
• Relationship of NBODY6 to NBODY4
• For more information and contents of define.f

Introduction to NBODY4

This introduction gives an overview of NBODY4, the NBodyLab adaptation of NBODY4, examples of simulation results, timings, etc.

Introduction to Running Simulations with NBODY4 and NBODY6 contains most of the information below in a document format (40 pages).

Guidance on initial parameters

A convenient cluster density model can be generated with N=1000 members and KZ5 = 1 which gives an isotropic Plummer model. Approximate equilibrium is achieved by setting Q = 0.5 while smaller values lead to initial collapse and an early core-halo configuration. As for suitable output times, try DTADJ = 2 (or proportionally smaller for N>5000) for energy check interval and DELTAT = 10 for producing main results. The termination time, TCRIT, should be consistent with the maximum CPU time (in minutes) given by TCOMP.

The initial mass function (IMF) may be generated in several ways using option KZ20, the simplest being a power-law distribution in decreasing order with exponent ALPHA = 2.3 and KZ20 = 1. Upper and lower mass limits are specified by BODY1 and BODYN in solar mass units which are adjusted to yield a specified average, ZMBAR.

More realistic star cluster simulations can be made by including a population of primordial binaries. The optional procedure KZ8 = 1 or 3 requires an extra input line after the virial ratio. With a length unit scale factor RBAR = 1 pc, the maximum semi-major axis in AU becomes a = 2.0D+05*SEMI which gives a period of (a**3/MTOT)**0.5 yr for binary mass MTOT in Msun. The minimum size, AMIN = SEMI/RANGE, should not be too small - the solar radius is 0.005 AU. If ECC < 0, the eccentricities are randomized with a 'thermal' distribution and average value = 2/3. Since the pericentre distance is given by a*(1 - e) and stellar radii tend to increase with time, the overlapping collision condition may be satisfied. In this case a new composite star is created (still conserving total energy). A variety of mass functions are available for KZ20 > 1, with the binary components selected independently from several types of IMF according to option values (KZ20 = 4 is recommended). When binaries are included, the most massive single star has identity 2*NBIN0 + 1 and is usually more massive than the first binary component with identity 1. Note the meaning of the input value N here. Thus the number of single stars is NS = N - 2*KS for KS hard binaries. Hint on input parameters: take RMIN = 2*SEMI if possible in order for all hard binaries (but RMIN < 1/N) to be regularized initially. The choice KZ8 = 3 instead of 1 produces a binary output file OUT9 containing the following quantities: binding energy, eccentricity, central distance, masses (solar), period (in days) and stellar type, ranging from 0 to 14 according to initial mass and age..

The merging of two clusters is an interesting topic for study. The relevant initial conditions are given by KZ5 = 2 which distributes two Plummer spheres in bound orbit. Now the initial separation, APO, and orbital eccentricity, ECC, must be specified and read after the IMF parameters. This produces a cluster binary orbit with semi-major axis APO/(1+ECC). The membership of the second Plummer model, N2, and scale factor, SCALE, are also read here. Such results are best viewed using the 3D animation Java applet which includes rotation and zoom.

Another optional procedure (KZ5 = 3) includes the effect of an external massive perturber on a thin planetesimal disk in the XY-plane. The planetesimals are distributed within a ring from 0.5 to 1.0 in circular orbits around the Sun, with a small total mass 0.001*Msun. The initial separation, APO, and predicted minimum distance, DMIN, are in N-body units. To avoid possible numerical problems due to large mass ratios, it is recommended to choose DMIN well outside the disk radius of 1. Since the length unit may be taken as 1 AU (Earth-Sun distance) here, the scale factor RBAR is converted from pc to AU by 1.0/2.0D+05 in order to obtain the correct physical time. The choice of elliptic, parabolic or hyperbolic perturber orbit is governed by the eccentricity. Thus for hyperbolic two-body motion, ECC > 1 and the binding energy in all cases is proportional to -1/SEMI = (ECC - 1)/DMIN. The initial velocity is adjusted to the specified periastron distance DMIN ( < APO) for the given energy. Since the disk system is displaced from the centre of mass, an external tidal field is not appropriate and this option is turned off. Also note that the standard scaling by the virial ratio Q is bypassed. Some suggested input values: APO = 6, ECC = 1.1, DMIN = 3, SCALE = 0.5 and TCRIT=20. The total time for a parabolic passage can be estimated from Tp = 2*DMIN**1.5/(2*MTOT)**0.5, with MTOT = 1 + SCALE. The time in yr is then Tp/2*pi or TIME/2*pi. This calculation is best done with conservative regularization parameters: DTMIN = 1.0D-05, RMIN = 1.0D-04, KZ16 = 0. Plots of the total disk energy, EDISK, and average eccentricity, DISP, are shown as functions of time. From the theory of stellar dynamics, eccentricity is more sensitive to change than energy in fly-by systems.

An option is also included for studying the evolution of a dominant binary. As in the other special cases, KZ5 (= 4) is used to read a simple extra input line defined below. Note that the initial value of the semi-major axis is changed a bit by the final scaling of coordinates and velocities. A good choice of dominant masses (in units of the average mass) might be around 25. The binary is placed at the cluster centre with no overall velocity. The semi-major axis (a), binding energy (-M1*M2/(2*a)), cluster energy and eccentricity are written to a file every output time (DTADJ) provided at least one of the components is original. The evolution of the binding energy and effective cluster energy is plotted at the end, as well as the eccentricity. It is recommended to exclude stellar evolution mass loss (KZ19 = 0) for heavy binary components. Also note that an energetic binary may escape from the cluster after a strong interaction. The case of two free-floating bodies may be studied by taking ECC > 1, when the masses M1 and M2 are assigned to the two first cluster members. If significantly heavier than the average, a dominant binary may form following a stage of mass segregation. However, for some mass choices and values of TCRIT the binary plot may not appear.

It is useful to distinguish between the integration itself and the purpose. For a standard test run, the main parameters in the first category may be taken as ETA = 0.02 (dimensionless time-step factor), QE = 2.0D-05 (relative energy tolerance), ETAU = 0.2 (regularization time-step factor). It may be a good idea to reduce ETA for small systems (e.g. ETA = 0.005 for N<100). Other regularization parameters are assigned by the code with option KZ16 > 0. For a standard isolated cluster simulation, choose KZ14 = 0, while KZ14 = 1 specifies a circular orbit in a galactic tidal field with cluster half-mass radius near RBAR pc. Note that the value of RBAR should be consistent with the total cluster mass. A reasonable value would be 1 pc for N = 1000.

More experienced users may want to employ their own initial conditions. This can be achieved by generating a complete data file containing mass, position and velocity for each member and uploading the file. The data are scaled to the standard total energy (-0.25) and total mass (1.0), with velocity magnitudes adjusted by the input parameter Q for KZ22 = 2. There is no scaling with KZ22 = 3 which therefore requires the data set to be consistent (but without restriction on the masses).

Research projects that have used NBODY4

The Cambridge GRAPE-6 has been used with NBODY4 for a variety of N-body simulations. The study of realistic star clusters is concerned with modelling relevant aspects of stellar evolution combined with consistent dynamics. Typical models of rich open clusters studied have 30,000 single stars. Such calculations may require a month's dedicated effort in order to describe the behaviour until complete dispersal. Results of more complete models with up to 12,000 stars and 12,000 binaries have been compared with an actual cluster in order to understand the complex interplay between astrophysical processes and dynamics (Hurley et al., MNRAS, 323, 630, 2001). Most of these calculations were made at the American Museum of Natural History (Hurley et al. MNRAS 355, 1207, 2004).

An application to binary black hole evolution in stellar systems containing some 240,000 particles proved interesting when relativistic effects were included (Aarseth, Astrophys. Space. Sci., 285, 367, 2003). Other simulations have been concerned with studying the evolution of core radii in small globular clusters (Wilkinson et al., MNRAS, 343, 1025, 2003) and the formation of stable hierarchical systems (Aarseth, IAU Colloq. 191, 2003).

Suggested simulations

1. Two unequal Plummer models with higher mean density in the second.
2. Comparison of escape rate for equal masses versus general IMF.
3. Study the remnant bound core for positive total energy (Q = 1.1).
4. Time of significant binary formation as function of N (EB/E > 1/N).
5. Compare energy errors by varying random seed (DTADJ = 2, TCRIT = 10).
6. Plot radii of mass fractions for initial collapse (Q = 0, KZ7 = 3).
7. Mass segregation (two specific mass groups; upload by KZ22 = 2 or 3). (Hint: look at mean mass in the core as function of time: MC/NC.)
8. Study the effect of one dominant binary (KZ5 = 4, extra input).
9. Binary formation of two free-floating bodies (KZ5 = 4, ECC > 1).

NBODY4 Input parameters

Parameters (detailed descriptions below)	NBodyLab Default Values
KSTART TCOMP GPID	1 10.0 0
N NFIX NCRIT NRAND NRUN	1000 1 5 50000 1
ETA DTADJ DELTAT TCRIT QE RBAR ZMBAR	0.02 2.0 10.0 100.0 2.0D-05 1.0 0.5
KZ1 KZ2 KZ3 KZ4 KZ5 KZ6 KZ7 KZ8 KZ9 KZ10	0 0 1 0 1 0 5 0 0 0
KZ11 KZ12 KZ13 KZ14 KZ15 KZ16 KZ17 KZ18 KZ19 KZ20	0 0 0 0 1 0 0 0 3 0
KZ21 KZ22 KZ23 KZ24 KZ25 KZ26 KZ27 KZ28 KZ29 KZ30	1 0 2 0 1 2 0 0 0 2
KZ31 KZ32 KZ33 KZ34 KZ35 KZ36 KZ37 KZ38 KZ39 KZ40	0 0 0 0 0 0 1 0 0 0
DTMIN RMIN ETAU ECLOSE GMIN GMAX	1.0E-05 1.0D-04 0.2 1.0 1.0E-06 0.001
ALPHA BODY1 BODYN NBIN0 ZMET EPOCH0 DTPLOT	2.3 10.0 0.2 0 0.02 0 10.0
APO ECC N2 SCALE (Line inserted only when KZ5=2)	6.0 0.5 500 .5
APO ECC DMIN SCALE (Line inserted only when KZ5=3)	6.0 1.1 3 .5
SEMI ECC M1 M2 (Line inserted only when KZ5=4)	1.0D-03 0.5 25.0 25.0
Q   0 0 0	0.5   0 0 0
SEMI ECC RATIO RANGE   0 0 0       (Above line inserted only when KZ8=1 or 3)	0.0005 -1.0 0.0 100.   0 0 0

Definition of input parameters for experimentation

If KZ8=1: For primordial binary information,
an additional input line   SEMI ECC RATIO RANGE 0 0 0
is inserted as shown in red in the above table of input parameters. Also set NBIN0 (described above).

If KZ5=2: For generation of two Plummer spheres in orbit,
an additional input line   APO ECC N2 SCALE
is inserted as shown in red in the above table of input parameters.

If KZ5=3: To model the effect of a passing perturber on a planetesimal disk,
an additional input line   APO ECC DMIN SCALE
is inserted as shown in red in the above table of input parameters.

If KZ5=4: To study a massive binary or two free-floating bodies,
an additional input line   SEMI ECC M1 M2
is inserted as shown in red in the above table of input parameters.

Output values printed and plotted

Quantities printed and plotted on the main simulation web page include:

NSTEPS		Number of integration steps (direct Hermite and regularized pairs)
DE		Relative energy error (every energy check, summed with sign main output)
<R>		Half-mass radius (N-body units or pc as appropriate)
DMIN		Closest two-body separation (first is all-time and 2nd is since last main ouput time)
N		Cluster membership (in case of escape removal)
DETOT		Final sum of all energy changes (relative and total)
E		Total energy (should be constant except for escape removal)
EB/E		Energy in binaries (relative to total energy)
MEMBERS		Identity of most energetic binary components
AMIN		Semi-major axis of dominant binary (in AU: 2.0D+05*RBAR*AMIN)
RC		Core radius (size of high-density region, membership NC)
RTIDE		Tidal radius (N-body units; typically around 10)
DISP		Eccentricity dispersion (rms value) for KZ5 = 3
EDISK		Total energy of disk particles for KZ5 = 3

Quantities to review in the full NBODY4 output printout

KS		Number of regularized binaries (i.e. energetically important or hard)
MC		Core mass (N-body units; average is MC/NC)
NS		Number of single stars (N = NS + 2*KS except for triples)
Q		Ratio of kinetic and potential energy (=0.5 for equilibrium)
T6		Time in Myr (RBAR assumed = 1 nominally for isolated cluster)

N-body units and astrophysical units

The code employs so-called N-body units which are defined by the sum of masses =1, total energy of bound systems -0.25 and the gravitational constant of unity. This implies a mean square velocity 0.5 at equilibrium and a harmonic mean separation of 1. Hence the average mass is 1/N and the half-mass radius for a Plummer density distribution is 0.8 model units. In order to obtain results for astronomical length units in pc (1 pc = 3 light years), scale the coordinates by the input parameter RBAR and velocities by VSTAR. From dimensional analysis (Aarseth's book p.112), the physical time in Myr is given by Tphys = 15*(RBAR**3/MTOT)**0.5*TIME, where MTOT is the total mass in solar units and TIME is the N-body time. For the velocity scaling, VSTAR = 0.066*(MTOT/RBAR)**0.5 km/sec. In the standard case, MTOT = ZMBAR*N with ZMBAR specified as the mass unit (in Msun) at input. After the scaling, ZMBAR continues to have the meaning of solar mass unit and is therefore rescaled (divided) by 1/N such that individual masses in Msun become BODY*ZMBAR. Finally, the mean crossing time is simply 2*sqrt(2) which suggests a convenient output time DTADJ = 2 (reducing to 1, 0.5, 0.25 for N >= 5,000, 10,000, 12,000 (to a maximum allowable of 15,000) to maintain the browser connection).

Integration methods

The basic integration employs the Hermite scheme, as summarized in the excerpt from the GRAPE book by Makino & Taiji . The special-purpose GRAPE hardware evaluates the force and first derivative for up to 48 particles at each cycle from the predicted coordinates and velocities. These values are used to construct the two next force derivatives whose contributions to the predicted quantities are added as a corrector on the host. We take advantage of the parallel architecture by introducing hierarchical time-steps (divisible by 2) such that many particles can be advanced as a group. In general, there are few members at the smallest time-steps and about 12-15 different levels in the hierarchy, depending on N and the range in density.

NBODY4 does not rely on softening of the force and several powerful procedures are included on the host to deal with strong point-mass interactions. In the first instance, close two-body encounters are handled by the Kustaanheimo--Stiefel regularization method (as discussed in Aarseth's book). Two particles are selected for treatment when their separation becomes smaller than RMIN and the time-steps fall below DTMIN. With standard N-body units, typical values of the regularization parameters may be taken as RMIN = 1/N and DTMIN = 0.01/N for N > 1000 and a bit more conservative for smaller N. If option KZ16 = 1, these parameters are updated at every energy check according to the core density. An arbitrary number of binaries and hyperbolic encounters can be treated simultaneously. Likewise, stable triples are converted temporarily into two-body systems where the inner binary is replaced by the combined mass until such time as the stability condition is violated. The stability of a hierarchical system is defined in terms of the orbital elements of the inner and outer binary (Mardling & Aarseth 1999). Strong interactions involving three or four members (of the type B--S or B--B) are selected for chain regularization (KZ15 > 0, KZ30 > 0), provided similar time-step and distance criteria are satisfied as for the two-body case. Note that all close encounter procedures may be bypassed by setting KZ16 = 0 and DTMIN = 1D-20, although this is not recommended in general.

Stellar evolution

Stars more massive than the Sun undergo significant changes over the life-time of a typical rich open cluster. Synthetic stellar evolution procedures are used to model the effect of continuous wind mass loss during the giant stage. Full descriptions of the relevant algorithms are given in Tout et al. MNRAS 291, 732, 1997 and Hurley et al. MNRAS 315, 543, 2000. The most massive stars (M > 8 Msun) are given a velocity kick in supernovae explosions. Corrections are made for the mass loss from the cluster such that total energy is still conserved. Strong interactions involving binaries often lead to their escape by recoil which also modifies the total energy of the remaining system. Plots of the stellar luminosity and effective temperature are shown in the Hertzsprung-Russel (HR) diagram for single stars and binaries.

Relationship of NBODY6 to NBODY4

The code NBODY6 is intended for different types of star cluster simulations and runs on standard laptops and workstations. A closely related code, NBODY6++, is available for conventional parallel supercomputers (Spurzem 1999). The main difference with NBODY4 is that a neighbour scheme is used to speed up the integration while the GRAPE versions rely on superior computing power to obtain the particle accelerations. A description of a few additional parameters is given in AA2003. All these codes employ the full range of close encounter procedures as well as the Hermite integration scheme. Although NBODY6 advances the particles sequentially, it still relies on the hierarchical time-step algorithm.

NBODY6 requires a few extra input parameters to deal with combining two force polynomials. This includes a second time-step tolerance, ETAR, for the so-called regular (or distant) force which is only updated on a somewhat longer (factor of 10) time-scale. A neighbour list of maximum size NNBMAX for each member is employed to obtain the force and first derivative due to particles inside a given radius (denoted RS0). Each list is updated during the total force evaluation, using a neighbour radius criterion based on the local density contrast. Particles at larger central distances are therefore assigned a decreasing number of neighbours because of the weaker irregular force. Because of compensating factors, the actual performance of the scheme is not a sensitive function of the neighbour membership. Replacing summation of the background force by prediction at each irregular time-step leads to a considerable speed-up while the self-consistent (or collisional) nature of the total contribution is preserved. A convenient feature is that derivative corrections for both the force polynomials due to change of neighbours may be avoided for the single particles, provided that the output times are commensurate with the maximum time-step of one time unit. However, this still entails introducing three explicit force derivatives to be used for neighbour changes relating to perturbed c.m. motion.

The treatment of regularization is formally the same as in NBODY4, hence all the corresponding routines are employed. However, NBODY6 exploits the advantage of selecting perturbers for KS or chain regularization from the existing c.m. neighbour list. The frequent check on extending unperturbed KS motion in NBODY4 makes use of the fast GRAPE function in determining the closest particle only. Some care is needed to ensure that sufficient neighbours are available for selection by using a more generous value of NNBMAX but the situation improves with increasing N.

Output from NBODY6 is fairly similar but contains a few additional data points. The ratio of irregular to regular time-steps is a good indicator for the efficiency of the scheme and this tends to increase slowly with N. These step counters are printed as the first and third entries at main output, while the fourth gives the number of regularized steps. Provided the average neighbour number (denoted ) is well behaved (say around sqrt(N)), profiling shows that an increasing fraction of the CPU time is spent evaluating the total force in larger systems. Even so, the number of irregular time-steps are very similar to the actual time-steps in NBODY4. As for a direct performance comparison, the neighbour scheme already becomes favourable for quite small cluster memberships (e.g. N = 50).

NBODY4 and NBODY6 have been designed to use similar data structures and options. A few options relate to incompatible features. Further information can be found in the documentation accompanying the download page.

See NBODY6 output examples and binaries for Windows and Mac OS.

For more information

For more information on all the NBODY4 input parameters see define.f and the documentation and code at the download page.

To sign up for news, get help or request permission for longer simulations, write to nbodylab@interconnect.com.