Issue 
A&A
Volume 522, November 2010



Article Number  A70  
Number of page(s)  22  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/200913890  
Published online  04 November 2010 
MYRIAD: a new Nbody code for simulations of star clusters
^{1}
Department of PhysicsAristotle University of Thessaloniki,
54124
Thessaloniki,
Greece
email: simos@astro.auth.gr
^{2}
Theoretical Astrophysics, EberhardKarls University of Tübingen,
72076
Tübingen,
Germany
Received:
16
December
2009
Accepted:
13
June
2010
Aims. We present a new C++ code for collisional Nbody simulations of star clusters.
Methods. The code uses the Hermite fourthorder scheme with block time steps, to advance the particles in time, while the forces and neighboring particles are computed using the GRAPE6 board. Special treatment is used for close encounters, and binary or multiple subsystems that form either dynamically or exist in the initial configuration. The structure of the code is modular and allows the appropriate treatment of more physical phenomena, such as stellar and binary evolution, stellar collisions, and evolution of close blackhole binaries. Moreover, it can be easily modified so that the part of the code that uses GRAPE6, could be replaced by another module that uses other acceleratinghardware such as the graphics processing units (GPUs).
Results. Appropriate choice of the free parameters give a good accuracy and speed for simulations of star clusters up to and beyond core collapse. Simulations of Plummer models consisting of equalmass stars reached core collapse at t ≃ 17 halfmass relaxation times, which compares very well with existing results, while the cumulative relative error in the energy remained below 10^{3}. Comparisons with published results of other codes for the time of core collapse for different initial conditions, show excellent agreement. Simulations of King models with an initial massfunction, similar to those found in the literature, reached core collapse at t ≃ 0.17 ± 0.05 halfmass relaxation times, which is slightly earlier than expected from previous work. Finally, the code accuracy becomes comparable to and even better than the accuracy of existing codes, when a number of close binary systems is dynamically created in a simulation. This is because of the high accuracy of the method that is used for close binary and multiple subsystems.
Conclusions. The code can be used for evolving star clusters containing equalmass stars or star clusters with an initial mass function (IMF) containing an intermediate mass black hole (IMBH) at the center and/or a fraction of primordial binaries, which are systems of particular astrophysical interest.
Key words: gravitation / globular clusters: general / galaxies: star clusters: general / methods: numerical
© ESO, 2010
1. Introduction
Systems of numerous individual members interacting with each other are found on both small and large scales. Gas molecules interact with each other by means of electrostatic van der Waals forces, the same forces are responsible for the interactions and collisions of aminoacids that lead to the construction of protein molecules. On the other hand, planets and planetesimals interact with both each other and the central star in a planetary system, with the gravitational force. The gravitational force is also the interaction force between members of greater systems such as open or globular starclusters, galaxies, and groups of galaxies.
Nbody codes provide one way of simulating systems like those mentioned above. Those codes require calculations of all the forces between all the individual particles at every integration time step. This means that at every timestep, a total of N(N − 1) ≈ N^{2} force calculations should be made, where N is the number of interacting particles. From that, it is easy to conclude that Nbody codes have time complexity O(N^{2}) and require fast computers or even special purpose hardware, such as GRAPE systems, to simulate realistically large systems. This is why the evolution of Nbody codes that could integrate these systems in a reasonable time, is closely connected to the hardware evolution.
The first attempt to perform a simulation, albeit without the use of a computer, but instead a system of 37 light bulbs represending two interacting model galaxies, was made by Holmberg (1941). The first computer Nbody simulations were developed by von Hoerner (1960, 1963) using initially 16 and later 25 particles. Later, as simulation methods and computers became faster, the number of particles increased to 100 (Aarseth 1963) and 250 (Aarseth 1968). Even in these early attempts some of the physical processes that occur in true star clusters, such as mass segregation and formation of hard binaries, were found in the simulations, but simulating systems closer to realistic star clusters required improvements in both software, by developing efficient algorithms, and hardware.
A significant step in the evolution of Nbody algorithms was the introduction of individual time steps: all particles do not share the same time step, but each one of them has its own, depending on its motion. In blocktimestep schemes (McMillan 1985), the time steps are quantized, usually as powers of 2, permitting blocks of particles to be updated simultaneously. This idea made the algorithms faster, because only the necessary calculations need to be made for a given accuracy. Today all sophisticated Nbody codes such as Starlab^{1} (Portegies Zwart et al. 2001; Hut 2003), φGRAPE (Harfst et al. 2008), and Aarseth’s nbody^{2} series (nbody4, nbody6, nbody6++ (Aarseth 1999a,b; Spurzem 1999)) use block time steps. Another milestone in the history of Nbody simulations was the introduction of special purpose hardware, responsible only for the calculation of the gravitational forces. The GRAPE series of computers (Makino et al. 1997, 2003a) were designed to perform only this type of calculation. The next step was GRAPEDR^{3} (Makino 2008), which is a more general purpose hardware specialized for Nbody codes which operates faster than a normal CPU. Simulations of realistically large star systems can be made using a single GRAPE machine attached to a fast hostcomputer. Today, GPUs are slowly replacing CPUs and the older GRAPEs in Nbody simulations (Portegies Zwart et al. 2007a; Belleman et al. 2008; Hamada & Iitaka 2007), because they are faster and considerably cheaper, than either CPUs or GRAPEs.
In this paper, we assess the limitations of Nbody codes with today’s software and hardware power. First of all, Nbody simulations are eithercollisionless or collisional. In collisionless simulations, the gravitational law is modified in close encounters to avoid collisions of particles and disencourage the dynamical formation of hard binary systems. The modified gravitational force acting on a single particle i due to all other particles is calculated to be: (1)where ϵ is the socalled softening parameter. The presence of close binary systems slows down the simulations, thus when handling large numbers of particles, one prefers to use collisionless simulations.
In simulations of collisional systems, the minimum number of stars that can be modeled is much lower and depends strongly on the fraction of stars that belong to binary systems in the initial configuration (primordial binaries). This is because the evolution of binary systems requires special treatment, including small time steps. Baumgardt et al. (2005) presented the results of simulations of unequalmass star clusters containing up to N = 131072 stars and an intermediatemass black hole (IMBH) at the center. Those simulations were oparated for 12 Gyr and no primordial binaries were included. In Portegies Zwart et al. (2007b), the results of simulations of star clusters with a total of N = 144179 stars including 10% primordial binaries were presented. The simulations were run for 100 Myrs. From the above numbers, it is easy to conclude that although it is not possible to simulate a typical galaxy with all its stars (N ≈ 10^{8} − 10^{9}), a simulation of a typical globular cluster (N ≈ 10^{5} − 10^{6}) for timescales of the order of some Gyr is possible, with today’s computing power. If primordial binaries are not included, the evolution of a star cluster with up to 10^{5} stars up to and beyond core collapse is possible, using a GRAPE6 machine attached to a fast hostcomputer. In the near future, using new algorithms and the coming generations of hardware, it should become possible to simulate the dense and massive centers of galaxies (N ≈ 10^{6} − 10^{7}).
During the past few years, Nbody simulations have provided a realistic tool for explaining the dynamics and structure of star clusters. Portegies Zwart et al. (1999) and Portegies Zwart & McMillan (2003) explained one way of forming IMBH in some globular clusters, by means of the process of runaway collisions of stars. Baumgardt et al. (2005) used detailed Nbody simulations to discuss the structure of a globular cluster containing an IMBH, while Gualandris & Merritt (2008) discussed the possibility of the ejection of supermassive black holes from galactic centers. Kupi et al. (2006) and AmaroSeoane et al. (2009) used postNewtonian effects on their Nbody simulations to study the dynamics and the gravitational radiation from IMBHbinaries in star clusters. Finally, Baumgardt et al. (2008) discussed the effects of primordial mass segregation on the evolution of globular clusters.
Here we introduce Myriad^{4}, a new C++ Nbody code for collisional and collisionless simulations. The code is loosely based on the instructive online books found in The Art of Computational Science website^{5}. It has many of the features of existing codes, such as the use of block time steps, the Hermite 4th order integrator, GRAPE6 support, and special treatment for binary and multiple subsystems, but introduces a new structure and some new ideas especially in using GRAPE6 to find neighboring stars and perturbers of binary systems.
These ideas and a description of the basic structure of the code are presented in Sect. 2. In Sect. 2.2, the detailed steps of the code are described. More attention is given in describing the method for dynamical evolution of binary or multiple subsystems that may appear in a simulation, in Sect. 2.2.3. In Sect. 3 we compare our preliminary results with the results of other codes. In the final two sections, some of the future directions of this code are presented with a discussion about the source of numerical errors and the CPU time needed for the simulations. The code uses the socalled Nbody units, which, together with some basic equations for calculating the physical parameters of a star cluster are given in Appendix A.
2. The new code
2.1. Code structure
Fig. 1 Class hierarchy in Myriad. A class cluster is a set of several block classes. A block contains a variable number of particle and binary classes. A binary is a set of two or more particle classes. A class particle contains information (mass, position, velocity etc.) for a single star and functions that act on them. The class cluster contains information about the system as a whole (total mass, energy, halfmass radius, core radius etc.) and the required functions for calculating them. 
Fig. 2 Simplified graphical representation of Myriad. The arrows show the flow of the data in the code. Boxes represent input or output datafiles or outside applications, while circles represent sets of codefunctions. Myriad applications are those inside the large dashed box. The other boxes and circles refer to satellite programs that create initial data or make animations from Myriadoutput. 
The core of Myriad^{6} uses C++ classes to store and manipulate data. The simplest or lowerlevel class is the class particle, which contains all necessary information about a single star^{7} (mass, radius, position, velocity, acceleration etc.). The class particle also contains a set of functions that act on a star following its evolution in time, finding its neighbors, its time step and all the required data for the star itself. A set of particle classes that share the same time step at a given time, forms a higherlevel class, the class block. As with the class particle, for the class block there is a set of functions that act on both the class itself and its particle members.
Finally, the complete set of block classes (it might only be one if all stars evolve with the same time step) forms a cluster, which is the highestlevel class. The clusterfunctions, act on all particle classes of the system, and mainly finding information about the cluster of stars itself (total mass, center of mass, escapers, binaries etc.). There is also another class that is a member of a class block, the class binary, which consists of two or more particle classes representing stars that lie close to each other in space. The binary class has the required functions for the initialization of a binary or multiple starsystem, its evolution in time, and its termination, when necessary, and also returns the required data back to the cluster.
Figure 1 shows the hierarchy of classes in the code described above. Figure 2 is a simplified graphical representation of Myriad. Initial data are provided to Myriad from a single data file. This initial data files can be constructed by the Starlab package. The main output of Myriad are snapshots of the evolving system at every output time step dt_{out} and information about the whole cluster (total mass, core radius, halfmass radius etc.) at every diagnostic time step dt_{diag}.
Fig. 3 Schematic data flow between the CPU and the GRAPE6 for the Hermite 4thorder integrator. See text for description. 
2.2. Integration method
The integration method used in Myriad is the Hermite 4th order (H4) (Makino & Aarseth 1992) predictorcorrector scheme (PEC: PredictEvaluate forceCorrect) with block time steps. Accelerations and their derivatives (jerks) are computed using GRAPE6. Binaries, close encounters, and multiple subsystems that require small time steps are detected using GRAPE6 and evolved using the timesymmetric Hermite 4th order integrator (Kokubo et al. 1998), where the prediction (P) step is the same as in the simple H4, but the force evaluation (E) and correction (C) steps are repeated 3 times (P(EC)^{3}).
For each block i, t_{i}_{c} is the current time of its particle classes, dt_{i} is the time step, and t_{if} = t_{ic} + dt_{i} is the time after a single update.
The integration procedure consists of the following steps (the time complexity of the slowest steps is displayed in parenthesis):

1.
We define the universal current timet_{c} and the time after the step is completed t_{f} to be equal to the time and time forward of block 0, the block containing stars with the minimum time step, respectively (2)(3)

2.
We determine which blocks have time forward equal to t_{f} and label them as requiring update.

3.
We update each labeled block to its time forward. To update a single block with m members (N being the total number of stars):

(a)
If binary^{8} classes exist, update them to t_{f} To update a binary class with k members:

i.
We determine the time step dt_{b} for the P(EC)^{3} integrator.

ii.
We predict the position and velocity of each member (O(k)).

iii.
We calculate the forces between binarymembers on the CPU (O(k^{2})).

iv.
We calculate the perturbation forces from p neighboring stars on the CPU (O(kp)).

v.
We then correct the position and velocity of each member (O(k)). Processes (iii) to (v) are repeated n = 3 times.

vi.
We update the current time t_{c} = t_{c} + dt_{b} of the binary.

vii.
We continue from (i) until t_{c} equals t_{f}.

(b)
We predict the positions and velocities of all the particles of theblock, including normal stars and “virtual” stars that representthe centers of mass of binary or multiple subsystems(O(m)).

(c)
We evaluate the accelerations and their derivatives for all the particles of the block using the GRAPE6 hardware (O(mlog N)).

(d)
We find the new time step for each particle of the block. This will be used for the next update (O(m)).

(e)
We check for encounters between single stars or binaries.

(f)
We correct the positions and velocities of all the particles of the block (O(m)).

(g)
If an encounter is detected in step (e), we construct the binary, and replace its members with a “virtual” star represending their center of mass.

(h)
We check the binary classes. If a binary is terminated, replace its center of mass in the Nbody code, with its former members and erase it.

(i)
We then send the corrected values to GRAPE6 memory.

4.
We move particles between block classes according to theirnew time step.

5.
We update t_{c} and t_{f} for each updated block.

6.
We continue with step 1.
From the above, it is obvious that the speed of the code depends on:

1.
The total number of stars N.

2.
The number of block classes and the average number of stars per block.

3.
The average number p of stars responsible for significant perturbations on a binary system (perturbers).
Figure 3 shows how data flows between the CPU and the GRAPE6 and which parts of the calculation are performed on the CPU and the GRAPE6. Initially positions and velocities, and estimates of both the accelerations and the jerks of all particles are stored in GRAPE6 memory (This procedure is omitted in the graph). Then, a prediction is made for a set of iparticles using the CPU. The predicted values and the time of each particle is sent to the GRAPE6 buffer, which directly sends them to the internal predictor of the GRAPE6. The internal predictor takes the data for the rest jparticles from the GRAPE6 memory and predicts their positions and velocities at the time of the iparticles. The GRAPE6 then computes the acceleration, its derivative (jerk), and the nearest neighbors for the iparticles and returns them to the corrector operating on the CPU. The corrector calculates the corrected values of position and velocity for the iparticles and sends the new data to the GRAPE6 memory. A new set of iparticles then, is sent to the predictor and the procedure continues until all particles are evolved to the expected time. Data flow between the CPU and the GRAPE6 occurs three times per time step. The most timeconsuming calculation, the calculation of the gravitational forces (O(Nlog (N))), performed on the GRAPE6. The rest of the calculations are of order O(N). If the forcecalculation were to be performed on the CPU, then it would be far more timeexpensive as it would scale as O(N^{2}).
2.2.1. Block time steps
Fig. 4 Distribution of equal mass particles in different time steps (powers of 2). Aarseth’s criterion is used for the time step with η = 0.01. A small number of particles lies in time step 2^{14}. This is the smallest allowed time step for the Nbody code. Stars with this time step are candidates for close encounters or virtual particles representing the centers of mass of binary or multiple subsystems 
The use of block time steps have proven to be an ideal method for Nbody simulations. The advantage of block time steps over constant (shared) time steps for all particles is that a smaller number of operations is needed to achieve a given accuracy. The advantage over individual time steps is that particles are monitored with time in groups and not individually. The implementation of block steps is ideal with GRAPE6, because it lessens the communication time between the host and GRAPE6. In individual timestep algorithms, a single particle is sent to GRAPE6 everytime, while in block time step algorithms, groups of particles, that share the same time step are sent simultaneously to GRAPE6 for force evaluation.
At t = 0, all particles, unless there are binaries in the initial configuration, share the same time step, which is the smallest time step allowed by the Nbody code. This time step is given by (Aarseth 2003) (4)where is the mean mass of the system, η_{I} the initial accuracy parameter with typical value 0.01, and R_{cl} is the close encounter distance (5)where G is the gravitational constant, σ is the rms velocity dispersion given by (6)at equilibrium, N is the number of stars, and R_{V} is the virial radius of the cluster (7)where U is the total potential energy.
The time step Dt_{min} is rounded to the closest power of 2 and the resulting Δt_{min} is the time step of block 0 and the particles that belong to it.
Each time a particle i is updated, its next time step is calculated according to Aarseth’s criterion (Aarseth 2003)(8)where η is the accuracy parameter with a typical value 0.01, and a_{i,1}, ȧ_{i,1}, , and are the acceleration of particle i and its time derivatives, calculated in next section. Subscript “1” refers to the values of the parameters at the end of the time step (final values), while subscript “0” refers to the values at the beginning (initial values).
The second and third time derivatives of the acceleration at the end of the time step are given by (9)and (10)where Δt_{i,0} is the previous time step of particle i.
Time steps are quantized according to the rule Δt_{i,1} = 2^{n}, where n is a negative integer defined in such a way that (11)In this way, all particles have time steps that are powers of 2. All particles with the same time step are grouped in a block so that they are updated simultaneously. The time step of “virtual” particles replacing binaries or multiple systems (see below) is set to be Δt_{min}. In addition, Δt_{min} is the time step of particles that are relatively close to and approaching each other rapidly. To improve the efficiency of the code, every new time step cannot be longer than 2 times its previous value, jumps to higher time steps with step greater than 2 are not allowed.
The diagram of Fig. 4 shows an example of a particle distribution in different time steps in a cluster of 16 384 (2^{14}) stars.
2.2.2. The Hermite 4th order predictorcorrector scheme
The H4 scheme is the integrator used in many Nbody codes. The most important feature of this integrator is that it needs only one calculation of the acceleration and its derivatives (force calculation) in every time step. This makes H4 superior to other popular integrators of the same accuracy, which require more than one force calculations per time step. A RungeKutta 4th order algorithm evaluates the force three times every time step, which makes this integrator much slower than H4 and consequently inappropriate for Nbody simulations with large number of particles.
Another reason for choosing the H4 integrator is that GRAPE6 is specifically designed for use with this kind of integrator. This is because GRAPE6 returns the acceleration and its first derivative for each particle, which are the only data H4 needs. One may use a higher order integrator of the Hermite family (Nitadori & Makino 2008), but this would require extra calculations to be made of higher order derivatives of the acceleration on the CPU. Higher accuracy, also is not necessary for this kind of problem.
The integration steps for advancing a particle forward in time with the H4 scheme are given in the following lines (see also Makino & Aarseth 1992).
We assume that t_{i,0} is the current time of the particle, t_{i,1} = t_{i,0} + Δt_{i,0} is the time after a step is taken, and Δt_{i,0} is the current time step.

1.
We predict the position and velocity ofthe particle using already known data:(12)(13)This operation has a time complexity of O(1).

2.
We evaluate the acceleration a_{i} and its derivative (jerk) of the particle. To achieve this, we have first to predict the positions and velocities of all the other particles to time t_{f}. The acceleration and jerk are evaluated on either the GRAPE6 or the CPU using (14)(15)where (16)(17)If softening ϵ is included in the calculations (for collisionless simulations), then the denominators of Eqs. (14) and (15) are replaced with and , respectively.
We also, evaluate numerically higher order derivatives of the acceleration using (18)(19)The evaluation of a_{i,1} and has time complexity O(N), when it is performed on the CPU and O(log N), when it is performed on the GRAPE6. The other calculations have complexity O(1).Fig. 5 Representation of Hermite 4th order with block time steps. In the figure, a set of 7 particles distributed in 3 timeblocks is shown. See text for description.

3.
We correct the position and velocity of the particle using higher order terms (20)The complexity of this operation is O(1).
Figure 5 shows how particles distributed in 3 timeblocks are updated using the H4 scheme. We note that in a typical star cluster, the number of timeblocks usually used is 10–15. All particles lie initially at the bottom, where t = 0. Initially block 0 contains 3 particles, while blocks 1 and 2 each contain 2 particles. A black dot represents a particle after the correction step (a straight blue arrow represents the correction step), while a red dot represents a particle after a prediction step (a curved red arrow represents the prediction step). Force evaluation is necessary for every particle, only before a correction step, i.e., force evaluation happens for a particle only when it is black. For the force evaluation of a particle, a prediction of all other particles to the current time of the particle is needed. It is obvious that particles in block 0 are considered for predictevaluateandcorrect every dt. Particles of block 1 predict every dt, but subjet to evaluateandcorrect every 2dt. Finally, particles of block 2 predict every dt and evaluateandcorrect every 4dt. The curved brown arrows show how particles are allowed to move between blocks depending on the time step calculated using Eq. (8). A jump to a higher block is allowed only if it doubles the time step of a particle (particle 3 jumps from block 0 to block 1). In contrast, a jump to any lower block is allowed if necessary (particle 6 jumps immediately from block 2 to block 0).
2.2.3. Binaries and multiples
During the evolution of the system, close encounters between stars may occur. Close binaries also may be formed dynamically, after threebody encounters. Those subsystems require a small time step to be integrated accurately in time and this time step may be much shorter than the shorterst time step Δt_{min} used in the Nbody code. This means that this subsystem has to be specially treated to ensure that it evolves in time with the required accuracy and speed. The special treatment used for close encounters, binaries, and small multiple systems is described in the following lines.
The timesymmetric Hermite 4th order integrator (Kokubo et al. 1998) is the numerical integrator for binaries, multiples, and close encounters. This integrator is more accurate than the simple H4 integrator with the penalty of being a few times slower. Timesymmetry is achieved by applying the force evaluationcorrection (EC) steps n times, after the prediction (P(EC)^{n}) and by appropriate choice of the time step. The value n = 3 is used in the code because this is sufficient to achieve timesymmetry. Timesymmetry is responsible for transforming the H4 scheme into a symplectic integrator, which is more appropriate for the time evolution of small systems such as binaries. Symplectic integrators (Kinoshita et al. 1990) conserve the characteristics of the integrated system (energy, semimajor axis, eccentricity, period) with high accuracy, paying of course the penalty of time. In reality, when using symplectic integrators the energy of the integrated system is bounded and changes periodically with time. In contrast, nonsymplectic integrators are faster, but the energy of the integrated system increases monotonically with time. A brief description of the integration method for binaries and multiple subsystems is given at the beginning of Sect. 2.2, between items (i) and (vii).
The choice of the time step for the (P(EC)^{3}) integrator is important to preserve timesymmetry. The time step at each time t must be such that if the time were to be reversed and the particles move backwards, the time step choice would remain the same. We followed the logic of the algorithm presented in Kaplan et al. (2005) for the choice of the time step. According to this algorithm, the time step is only allowed to either increase or decrease by a factor of 2 or remain the same with respect to its last value and an increase in the time step is allowed only at even times. In this way, the time step remains symmetric and quantized as a power of 2. The current time is characterized as either even or odd with respect to the last time step taken by the system. We apply this algorithm using as a time step criterion the shortest among the time steps dt_{bi} of the members of the binary (21)where η_{b} is a parameter to control the accuracy, and a_{i} and ȧ_{i} are the acceleration and its derivative for particle i, respectively.
At each time t, we begin by finding whether t is even or odd with respect to the last time step, Δt_{old}, taken. If t is odd, we maintain the same time step and by using Δt_{old} we determine the time step at both the beginning dt_{0} and the end dt_{1} of that time step and find their average (22)We constrain the time step used Δt_{old} to be shorter than dt and longer than dt/2. If it satisfies this condition, then the chosen time step is the one to use in the integration. If not, then the time step must be halved. In the case of even time, we follow the same procedure, testing whether the doubling of the previous time step satisfies the criterion of being smaller than dt and greater than dt/2, where dt is given by Eq. (22). If the criterion is satisfied, then 2Δt_{old} is the time step for the next step. If not, then we follow the same logic using Δt_{old} as a test time step and follow exactly the same steps as if t was odd.
For this algorithm to work properly, an appropriate choice for η_{b} must be made. This choice needs to be small enough for dt not to change by more than a factor of 2 with respect to its previous value. In addition, η_{b} is responsible for providing a smooth transition to the time step from the Nbody code to the binary module of the code, to ensure that a particle that has either just joined a binary system or just left it does not change its time step by more than a factor of 2. This is required because major changes in the time step introduce significant errors in the integration. To guarantee this behavior, η_{b} is calculated for every binary at the time of its formation. To do this, we determine the particle i with the smallest value of a_{i}/ȧ_{i} and using Eq. (21) we determine η_{b} by requiring that dt = Dt_{min}/2, which is given by Eq. (4) and is the smallest time step of the Nbody code; it is also the last time step taken by all the future members of the binary or the next time step for all the particles that leave a binary. In the case of a binary with more than 2 members, we use a shared time step for all its members. To avoid major changes in the time step criterion, we use an even smaller value for η_{b}. The problem here is that if a single particle escapes for example from a triple system, leaving behind a close binary system, the time step of the escaping particle is forced to jump to Dt_{min} even if it was much smaller when the particle was still a member of the triple. This necessary change in the time step may introduce errors in the integration. To solve this problem, an individual time step must be used for multiple subsystems. In this way, all members of the binary would have their own time steps and an escaping particle would increase its time step by just a factor of 2.
Finding the binary systems in an Nbody system is a slow procedure with time complexity O(N^{2}), unless all particle classes are members of a linked list, as in a binary search tree, or GRAPE6 is used for acceleration. Here we use the neighbormodule of GRAPE6 to search more rapidly for binaries and close encounters. Each star i with mass m_{i} is associated with a distance (23)where m_{g} is the mass of the most massive star of the system, N the total number of stars, and R_{cl} is given by Eq. (5). During the calculation of the force for the star i, this distance is sent to the neighbormodule of GRAPE6, which returns the identities of the stars that lie in distances r ≤ R_{hi} from the star i. Those stars are considered as “neighbors” of the star i. Then, a search between the neighbors is performed to find out whether one or more of them is about to have a close encounter or to form a close binary system with star i. The criterion for a close encounter between the star i and its neighbor j is (24)(25)where R_{ij} is the separation between stars i and j, m_{i}, and m_{j} their masses, N is the total number of stars, and V_{ij} is the relative velocity of the two stars. The choice of the distance R_{critij} comes from our finding that in an equalmass system, if the distance between two stars is more than R_{cl}, then they can be integrated accurately using the H4 integrator with a time step greater than or equal to Dt_{min}, given by Eq. (4), which is the shortest time step allowed by the block timestep scheme. If the separation is smaller than R_{cl}, then a smaller time step is needed, so integrating the two particles in time with acceptable accuracy, needs some special treatment. The factor in Eq. (24) is used to generalize the criterion for a nonequalmass system. Finally, the distance R_{hi} that is sent to GRAPE6 during the calculation of the force for star i, is simply 5 times the critical distance between this star and the most massive star of the system, to ensure that every close encounter between two stars is detected.
If the above criteria are fulfilled, a binary with the two particle classes as members is created. The two particle classes are removed from the Nbody system and replaced by a new particle that corresponds to a “virtual” star represending the center of mass of the removed stars. This “virtual” star has mass (26)and is placed at position (27)with velocity (28)and behaves as if it were a normal star in the restframe of the Nbody system.
A third star k at distance R_{k} from the center of mass of the binary, becomes a member of the binary transforming it to a triple system, if one of the criteria below is satisfied:

1.
In all cases, if the distance of the star k from at least one of the members of the binary is smaller than 2/3R_{cl}

2.
If the relative motion of the two stars is a simple close encounter (e > 1, where e the eccentricity of the system), when (29)where M is the total mass of the binary and m_{k} the mass of the third star.

3.
If the relative motion of the two stars is a close binary (e < 1), then the dimensionless perturbation γ_{k}, acting on the binary from star k(30)becomes greater than a critical value, i.e., (31)
The same criteria are used when a fourth particle comes close enough to a triple system. The critical value of γ for soft binaries is γ_{crit} = 0.125. For hard binaries, we use the value γ_{crit} = 0.015625. In an equalmass system, according to Eq. (30), the first value corresponds to a distance R_{k} = 2 R_{b} for the third star, while the second to a distance R_{k} = 4 R_{b}. This means that in an equalmass system, a third star would become a member of a soft binary when it lies at a distance twice the size of the binary. If the binary is hard, the critical distance is 4 times the size of the binary.
Another way of forming multiple subsystems with number of stars n ≥ 4 is when two binaries or multiples approach each other. When the distance between their centers of mass becomes smaller than three times the sum of their sizes, then the two binaries merge forming a larger system where all particles interact strongly with each other. The size of a binary is defined to be its semimajor axis, for hard binaries or the distance between its members, for soft binaries.
Termination of a triple, quadruple, or higher order subsystem occurs when the following three criteria are met:

1.
When the dimensionless perturbation γ_{j} acting on the particlemember of the subsystem j on the inner binary, becomes less than γ_{crit}(32)

2.
The relative velocity of particle j, v_{j}, with respect to the inner binary is positive, i.e., particle j moves outwards.

3.
The distance of particle j from any of the other particles/members of the multiple system, is greater than R_{cl}.
The socalled “inner binary” is the hardest double subsystem inside the multiple (the one with the greatest binding energy).
In some rare cases, two particlesmembers of a multiple system with number of stars n ≥ 4 are removed together from the system and form a close binary.
Termination of a simple binary occurs when the distance between its members becomes greater than R_{cl}. Hard binaries are not terminated, unless another particle becomes a member and interacts strongly with their members. The outcome of these interactions may be a softer or harder binary, a collision, or even a dissolved system.
Collisions:
Another rare case for termination of a binary is a collision between its members. A collision between two stars occurs when their distance r_{ij}(33)where R_{i}, R_{j} the radii of the two stars.
The stellar radius R_{i} for a main sequence star is given by (34)where m_{i} is the mass of the star, R_{⊙} is the radius, and M_{⊙} is the mass of the sun. The parameter a depends on the mass m_{i} (Bowers & Deeming 1984) such that (35)If the star is a black hole, then R_{i} corresponds to its Schwarchild radius (36)where c is the speed of light.
Perturbers:
The internal motion of a binary or multiple subsystem is controlled by the internal forces acting between its member stars. External perturbations from stars passing close to the center of mass of the subsystem are also included. A star k is considered as a perturber of a binary system if the dimensionless perturbation γ_{k}, acting from the star to the binary, is greater than a critical value (37)A typical value of 10^{7} is used for γ_{pert}, but this value is usually modified to limit the total number of perturbers for a binary system according to the total number of stars. During the force calculation, the external perturbations for every starmember of a binary or multiple system are calculated and added to their internal values. The perturbation in the acceleration a_{i} and first derivative acting on a star i, a member of a binary or multiple system, caused by an external star j are calculated using (38)(39)where a_{cm,j} is the contribution of particle j to the acceleration of the center of mass, is the contribution of particle j to the first derivative of the acceleration of the center of mass, and a_{i,j} and are the direct contributions of particle j to the total acceleration and jerk of particle i, respectively.
Identifying the perturbers of a binary or multiple subsystem is done by using the neighbormodule of GRAPE6. During the force calculation for an “virtual” star i, a distance (40)is sent to the neighbormodule of GRAPE6. Where m_{g} is the mass of the most massive star in the system, m_{i} the mass of the “virtual” star, i.e., the total mass of the binary or multiple subsystem, and d the size of the subsystem i.e. the maximum distance between its members. To avoid the case of overflow of the neighbor lists of GRAPE6, each “virtual” star is sent alone to GRAPE6, after cleaning the previous neighbor lists from GRAPE6 memory. To avoid overflow and limit the number of perturbers per binary system in systems with tens of thousands of stars, R_{pi} is also modified to ensure that ~100 neighbors are returned. This is achieved by preventing R_{pi} becoming greater than 2.5 × r_{6}, where r_{6} is the distance of the sixth nearest neighbor of star i, used to find the core radius of the cluster, as we see below, and 2.5 × r_{6} is an estimate of the distance of the 100th nearest neighbor of star i, assuming that the density around star i is constant and unchanged from the time of the last calculation of r_{6}. This means that the calculation of the core radius should be repeated relatively often, to estimate the true number density around each star at any time. This calculation is usually iterated a few times every crossing time, which is defined to be (41)as described in Appendix A. With the above restrictions, the GRAPE6 memory does not overflow and GRAPE6 returns the identities of the “neighbors” of star i, which are the perturbers of the corresponding binary system. Nevertheless, if an overflow occurs, then the perturbers are found after searching among all the stars of the system and the result is used to avoid future overflows.
Fig. 6 Perturbers. Stars with γ ≥ γ_{pert} are considered as perturbers of a binary system. Those are the stars represented as empty (red) circles in the above image. The filled (black) circles represent the members of the binary or multiple subsystem. If the dimensionless perturbation γ, given by Eq. (30), becomes greater than γ_{crit}, then the star is included into the binary or multiple subsystem. More stars could be added forming a multiple system with up to ~10 stars. Those systems are usually dispersed quickly. On the other hand, a star with γ becoming smaller than γ_{pert} is no longer a perturber. Those are the stars represented with filled (blue) squares. For an equal mass system, we typically set γ_{pert} ~ 10^{7} and γ_{crit} = 0.125 if the inner binary is considered as “soft”, while γ_{crit} = 0.015625 is assumed if it is considered as “hard”. 
Figure 6 shows schematically the distribution of perturbers and nonperturbers for a close binary system.
We note here a couple of definitions used regularly in the remainter of the paper. In all simulations we refer to as escapers, the star clusters that are assumed to be isolated i.e., for which no external galactic field exists. Stars that are moving outwards and reach a distance two times the initial size of the cluster, are considered as escapers from the system. When a star escapes, the energy of the system is reinitialized. In addition we refer to as hard binary a close binary system whose energy per unit mass h < − 1 or whose semimajor axis a < R_{h}(m_{i} + m_{j})/2, where m_{i} and m_{j} the masses of its members and R_{h} its halfmass radius. In an equal mass system, the limit for a becomes a < R_{h}/N, where N is the number of stars.
3. Results and performance
We measure the accuracy and performance of Myriad by using it to evolve different systems in time. The most simple tests are those where an equalmass Plummer sphere (Plummer 1911) is used as an initial configuration. We evolve systems with different numbers of stars and measure the accuracy and the speed of Myriad. As a measure of the accuracy of the code, we use the relative error in the total energy, after evolving the systems for one Nbody time unit. The speed of the code is measured by the CPU time required for the simulations. We also performed simulations of systems with equalmass stars distributed in a Plummer sphere and systems with different initial mass functions distributed in either a Plummer sphere or according to a King density profile (King 1966), up to core collapse. In these simulations, we measure the time of core collapse and compare the results with other codes and with results found in the literature.
In all simulations, the stars are considered as interacting point particles. The real dimensions of a star are taken into account only when it collides with another star. Stellar evolution, stellar massloss, and binary evolution are also not included. If close blackhole binaries are formed, they evolve using purely Newtonian dynamics, without any postNewtonian terms in the equations. Finally, the effects of the Galactic tidal forces on the star clusters are ignored i.e., star clusters are assumed to be isolated. Initial data in all simulations were provided by the Starlab software environment.
The host computer used for these tests is an AMD Athlon(tm) 64 Processor 3500+ operating at 2.2 GHz with 2GB of RAM, connected to a 32chip GRAPE6 board and running Fedore Core 1 GNU/Linux with the 2.4.22 kernel.
3.1. Accuracy
There are many different sources of errors in an Nbody code. Two of them, which can be easily controlled, are the choice of integration scheme and its order. Here, as mentioned earlier, the integrator is 4th order. If the time step were constant and common for all stars, then the errors would scale as O(dt^{4}). In codes that use block time steps, such as Myriad, the parameter that controls the error is the accuracy parameter η is used in the timestep criterion of Eq. (8). Another source of errors, which is clear when small values of η are chosen, is the accuracy of the hardware, both CPU and GRAPE6. Myriad uses double precision for all calculations on the CPU, while the accuracy of the calculation of the accelerations and derivatives is controlled by GRAPE6 (Makino et al. 2003b).
Fig. 7 Cumulative relative energy error as a function of η for a Plummer model with N = 8192 (top left), N = 16384 (top right), N = 32768 (bottom left), and N = 65536 (bottom right). The integrations ended at t = 1 Nbody units. 
We performed a series of experiments using Myriad to evolve equalmass Plummer models, with different particle numbers N from t = 0 to t = 1 Nbody time units. For all the integrations, we recorded the cumulative relative energy error (ΔE/E = (E_{t = 1} − E_{t = 0})/E_{0}) and studied its dependence on the accuracy parameter η and particle number N. For this study, we varied the accuracy parameter η from 10^{4} to 0.2 and the particle number from 8192( ≡ 8K) to 65536( ≡ 64K), where 1K ≡ 1024. The results are shown in Fig. 7.
For small values of η (η ≤ 5 × 10^{3}), the relative energy error is almost constant. This is because for those values of η the integration error is small and the total error is dominated by the hardware precision. For values of the accuracy parameter η ≥ 5 × 10^{3}, the error increases with η as expected. When η ≥ 0.1, the errors are too large, so these choices of accuracy parameter are inappropriate for simulations. The typical choices for the accuracy parameter are (0.001 ≤ η ≤ 0.01) and the relative energy error between these limits is ΔE/E_{0} ≤ 10^{9} for N ≤ 16 K and ΔE/E_{0} ≤ 10^{7} for greater values of the particle number N.
The dependence of the error on the particle number N is weak and evident only at the lower values of η. For greater values of N, the relative energy error saturates to slightly smaller limits.
3.2. Performance
The elapsed time T for a simulation depends on the accuracy parameter η and particle number N. The elapsed time (wallclock time) as a function of η and N for the simulations of Fig. 7 is shown in Fig. 8. As expected, T grows with smaller η and greater N. We note that smaller η indicates that shorter time steps are used on average.
To study the speed of Myriad, we performed a set of integrations of equalmass Plummer models with different particle numbers N from t = 0 to t = 1 Nbody time units. For all simulations, we recorded the wallclock time as a function of N and the results are shown in Fig. 9, where N is chosen to be between 512( = 2^{9}) and 262144( = 2^{18}). The accuracy parameter used for all simulations was η = 0.01. In the same figure, the slopes of N^{2}, N^{3}, and Nlog _{10}N are shown for comparison.
The complete set of calculations inside Myriad has a time complexity of O(Nlog N). There is only one calculation that has a time complexity O(mN), where 1 ≤ m ≤ N is the average number of neighbors per particle returned by GRAPE6 or, if GRAPE6 returns an overflow, calculated by the CPU. This is the calculation of the core radius of the star cluster. This calculation is repeated every diagnostic time step dt_{diag}. For the experiments of Fig. 9, this calculation took place only once, before t = 1. As the number of stars increases, the calculation of the core radius takes a comparable amount of time to the integration time. This is why initially the slope in Fig. 9 is similar to the slope of Nlog N, but as N increases, it tends to more closely resemble the slope of N^{2}. We note that all stars at t = 0 begin with the shortest time step Dt_{min} and then, as time passes by, they are slowly distributed in all the available time blocks. Because of this, the simulation initially operates more slowly. We note that primordial binary systems were not included in all runs, while no close binary systems were dynamically formed until t = 1 Nbody time units. At later times, dynamically formed binary systems would force the simulation to run slower.
Fig. 8 Wallclock elapsed time as a function of N and η for the simulations of Fig. 7. All simulations ended at t = 1 in Nbody units. 
Fig. 9 Wallclock elapsed time as a function of the number of stars N. All simulations ended at t = 1 in Nbody units. The accuracy parameter η had the same value 0.01 in all simulations. The slopes of O(N^{3}), O(N^{2}), and O(Nlog N) are also shown for comparison. 
3.3. Binaries and multiples
Fig. 10 Cumulative energy error (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as a function of time for the first 10 periods for the simulations of equalmass binary systems with 4 different eccentricities. The eccentricities of the binaries are e = 0.19 (red line), e = 0.51 (green dashed line), e = 0.75 (blue dotted line), and e = 0.91 (cyan dasheddotted line). 
Fig. 11 Variation in the time step used for the evolution of the binary systems of Fig. 10 during 2 periods of each binary. The time step parameter η_{b} is chosen so that the error in all the binaries is of the same order. 
Fig. 13 Cumulative energy error (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as a function of time for 10^{5} periods for the simulation of a binary system with eccentricity e = 0.89 and mass ratio m_{heavy}/m_{light} = 15. 
Fig. 14 The trajectory of the two stars in the simulation of Fig. 13. The two stars are orbiting around their common center of mass, which is located at O(0,0). 
Fig. 15 Cumulative energy error (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as a function of time for the Pythagorean threebody system. 
Fig. 16 Trajectories of the particles in the simulation of the Pythagorean threebody system. The initial positions of the particles are marked with black circles. 
Before testing Myriad in simulations of clusters up to core collapse, we investigated the accuracy of the the code in handling close binaries and multiples. We recall that the algorithm used in the binary module of Myriad is the timesymmetric Hermite 4th order scheme described in Sect. 2.2.3 and note that the force calculation is done on the CPU.
Figure 10 shows the cumulative relative error in the energy for the first 10 periods of binaries consisting of equal mass particles and having different eccentricities. The choice of the accuracy parameter η_{b} in the time step calculation and the maximum allowed time step are such that for all simulations the initial energy errors are of the same order. The number of time steps per period for the binary with eccentricity e = 0.19 is 720, for the binary with e = 0.51 is 2200, for the binary with e = 0.75 is 2200, and for the binary with e = 0.91 is 4300. Figure 11 shows how the time step changes according to our time step criterion, during the first 2 periods of each of the binaries. As shown in the figure, the time step is quantized as a power of 2. Finally, the cumulative energy error for evolving the binaries for 10^{5} periods is shown in Fig. 12. For a clearer visual result, we have printed a fraction of the points in the figure. It is obvious that the algorithm used for the dynamical evolution of close binary systems is timesymmetric and preserves the energy since there is no linear growth in the energy errors. The energy error oscillates between a lower and higher value. The lower error is observed during the apastron passage, while the higher value is observed during the periastron passage. We would have the same behavior for the energy error even if we decrease the number of steps per period for the simulations. The error would be greater but still have an upper limit and no linear growth with time.
Figure 13 shows the cumulative relative energy error in the case of a binary with unequal mass members (mass ratio = 15) and eccentricity e = 0.89. We followed the simulation for 10^{5} periods. The energy error does not grow linearly with time, but instead follows some kind of random walk between some limits. This result shows that the binary module of Myriad is capable of evolving binaries with high mass ratios and high eccentricities for a long period without significant energy errors. The trajectories of the two particles around their common center of mass for the full simulation of 10^{5} periods are shown in Fig. 14.
We finally tested the performance of Myriad in handling triple systems by evolving the well known Pythagorean threebody system. This system consists of three different masses initially located at rest at the corners of a right triangle whose sides have lengths 3, 4, and 5. Each of the three masses is located in such a way that its mass is equal to the length of the opposite side. Mass m_{1} = 3 is opposite the side of length 3, m_{2} = 4 is opposite the side of length 4 and m_{3} = 5 is opposite the side of length 5. The three masses interact strongly with each other, and their trajectories after some time are shown in Fig. 16. The final result is a hard binary consisting of particles with masses m_{3} = 5 and m_{2} = 4 and an escaping particle. The evolution of the cumulative relative energy error for the simulation is shown in Fig. 15. We note the two “jumps” in the energy error. These “jumps” occur when two of the particles come very close to each other and as a result more time steps are needed to perform an accurate integration. We used a constant accuracy parameter η_{b} = 0.0001. The cumulative relative energy error at t = 100 Nbody units is 2.23776 × 10^{8}.
3.3.1. Choice of the parameters
In Myriad, assuming that a good choice for the accuracy parameter η is 0.01 and adjusting the accuracy parameter for binary evolution η_{b} so that there is a smooth transition between the binary module and the Hermite integrator, there is another free parameter that balances between speed and accuracy in collisional simulations up to core collapse. This is the number of perturbers for every binary or multiple subsystem controlled by the critical dimmensionless perturbation γ_{pert}. This parameter is discussed in Sect. 2.2.3 (Eq. (37)). For the rest of this work, we choose γ_{pert} = 10^{7}, which is small enough to ensure that the errors introduced into the simulations by replacing two or more particles that lie close to each other, with their center of mass, are not significant. Smaller values of γ_{pert} would slightly increase accuracy, but the speed of the code would be lowered significantly, since the forces of the perturbers on the members of a binary are computed on the CPU and not on the GRAPE6. Finally, if the value of γ_{pert} is responsible for large numbers of perturbers per binary system, this value is increased so that no more than 100 perturbers per binary are used.
3.4. Core collapse
The evolution of a star cluster up to and beyond core collapse is one of the challenges for any Nbody code. During core collapse, the density at the center reaches its maximum value, while the core radius, defined by the equation (42)which is described in Appendix A, reaches its minimum value. Figures 23 and 17 clearly show this behavior as a star cluster of 1024 equalmass stars evolves beyond core collapse. As the cluster evolves, close binary or multiple systems are formed and interact with single stars or other binary or multiple systems. Those interactions become more frequent and violent as the cluster approaches core collapse, because the density of stars around the center of the cluster reaches its maximum. All of these encounters have to be resolved with an acceptable accuracy. In addition, some stars approach the outer bounds of the system and consequently escape the system.
We tested the behavior of Myriad in evolving star clusters up to core collapse and compared the results with those produced by Starlab or found in the literature. For these tests, we used 3 different initial configurations. We initially used an equalmass Plummer model and studied the time of core collapse and the evolution of the error in the total energy. The same simulation with similar parameters was performed using Starlab for comparison. The next test was the evolution of a star cluster with an initial mass function up to core collapse. The stars were again distributed initially according to a Plummer density profile. Finally, we performed several simulations of star clusters with an initial mass function of a broad range of masses, distributed initially according to a King density profile. The result of the latter simulations was a measure of the core collapse time (t_{cc}) as a function of the halfmass relaxation time (t_{rlx}). These results were compared with existing results found in the literature.
3.4.1. Equalmass plummer models
We performed a simulation of an equalmass Plummer model up to core collapse to test the stability, accuracy, and speed of Myriad. Figure 17 shows the evolution of the Lagrangian radii that corresponds to fractions of the total mass of the system and the evolution of the core radius of a system, which was initially an equalmass Plummer model with N = 1024 stars. It is obvious that the core radius and the smallest Lagrangian radii reach a minimum at t_{cc} ≃ 340 Nbody time units, when core collapse occurs. From equation
The core collapse times for different initial models.
The comparison between Myriad and Starlab is also shown in Fig. 20, where the respective relative energy error ratio Myriad/Starlab is presented. The error of Myriad at any time step is smaller by 300 times than the error of Starlab, while most of the time the error of Starlab is greater. After t ~ 250 and until core collapse, the error of Myriad becomes several times smaller than that of Starlab. After core collapse, both Starlab and Myriad have large instantaneous errors and their ratio scatter from 10^{5} to 300. Figure 21 shows the cumulative relative energy error for the simulation as it evolves with time. The same error of Starlab again is presented for comparison. It is obvious that before core collapse the two codes show identical behavior and after core collapse the error of Myriad remains smaller than that of Starlab.
The simulation ended at t ≃ 402 Nbody time units. In Fig. 22, the CPU times of Myriad and Starlab are compared. For Myriad, the total CPU time and the CPU time spent in force calculation and binary evolution are shown. Before core collapse, the CPU time of Myriad is constantly greater than that of Starlab. This is because of the small differences between the two codes in the time step criterion, which on average causes Myriad to store more particles in small timestep blocks. A search for neighbors also occurs every time these particles are updated, which slows down the simulation significantly. The CPU time of the Myriad simulation until core collapse is about T ≃ 1.7 h, while the same time for Starlab is about T ≃ 1 h. After core collapse, the CPU time of Starlab becomes longer because of the formation of hard binaries. Hard binaries lead to an increase in the CPU time of Myriad cause it to finish in about T ≃ 8 h. The simulation done with Starlab, with similar parameters as those used in Myriad, took about 13 h, spending most of the time in the postcollapse part. Finally, in Fig. 23 we show the evolution of the mass density and the number of stars at the core. As expected, the core density increases with time and reaches a maximum at the time of core collapse. On the other hand, the number of stars at the core drops with time and reaches a minimum value at core collapse.
Fig. 17 Core radius (heavy black line) and Langrangian radii containing 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, and 3% (blue dashed lines from top to bottom) of the total mass. The half mass radius is indicated by a heavier dashed line. Core collapse is reached at t_{cc} ≃ 340 Nbody time units or t_{cc} ≃ 17t_{rlx}. 
Fig. 18 Lagrangian radii containing from top to bottom 90%, 50%, 10%, and 1% of the total mass. The continuous lines are the results produced by Myriad, while the blue dashed lines are produced using Starlab. 
Fig. 19 Error in the energy per quarter of the time unit (δE = E_{(t)} − E_{(t − 0.25)}) as it evolves in time for the simulation of an equalmass Plummer model of 1024 stars. The simulation ended at t ≃ 402 Nbody time units. The dashed line is the energy error of Starlab when it was used for evolving the same system. Before core collapse, the energy errors are identical. 
3.4.2. Plummer models with a Salpeter initial mass function
A simulation of a Plummer model with N = 16384 stars was performed assuming a (Salpeter 1955) initial mass function (IMF) and limiting masses of m_{lower} = 0.5 M_{⊙} and m_{upper} = 5 M_{⊙}. In Fig. 24 we show the evolution of the Lagrangian radii and core radius of the system. In Fig. 25, we illustrate the evolution of the error in the total energy of the system for this simulation. Core collapse is reached at t ≃ 2.2t_{rlx}, which is much earlier than expected for an equalmass system, given empirically by Eq. (A.15). This is expected because heavier stars tend to sink to the core faster and because of that core collapse occurs earlier. The time taken for a star of mass m to sink to the center from a circular orbit at r ≫ r_{c} is given by (Portegies Zwart & McMillan 2003) (44)According to Eq. (44), the most massive stars of the above system would sink to the center in a time of about t_{s} ≃ 0.73t_{rlx}. Until t = 2.2t_{rlx}, most of the massive stars are located close to the center and this leads to the observed core collapse of the cluster.
3.4.3. King models with an initial mass function
We performed simulations of star clusters with initial density distribution of a W_{0} = 6 King model (King 1966) up to core collapse. We used a (Scalo 1986) initial mass function (IMF) with lower and upper limits of 0.1 M_{⊙} and 100 M_{⊙}, respectively. We varied the number of stars from N = 6122 to N = 12288 and found the core collapse time t_{cc} as a fraction of the half mass relaxation time t_{rlx} of the systems. The systems chosen are similar to those used in Portegies Zwart & McMillan (2003), where the simulations were performed using the Starlab package including stellar mass loss from stellar evolution. The results are presented in Table 1. As the core collapse time we consider the time of the creation of the first dynamically formed hard binary system in agreement with Portegies Zwart & McMillan (2003). We performed 18 simulations and our result for the core collapse time is t_{cc} ≃ (0.17 ± 0.05)t_{rlx} . The result is slightly smaller than Portegies Zwart & McMillan (2003), which is t_{cc} ≃ (0.19 ± 0.08)t_{rlx} . This difference in results may be due to the lack of stellar mass loss and stellar evolution in Myriad which, according to Portegies Zwart & McMillan (2003), tends to delay core collapse.
The evolution of Lagrange radii and core radius in one of the simulations containing N = 6122 stars is shown in Fig. 26. For a clearer visual result, we have applied a cubic spline interpolation to the core radius. The cumulative error in the energy for the same simulation is presented in Fig. 27. In Figs. 28 and 29, we present the same results for one of the simulations containing N = 12288 stars. The cumulative energy error in all simulations remained below 10^{3}.
Fig. 20 Respective energy error ratio Myriad/Starlab as a function of time. 
4. Discussion
Sources of energy errors:
There are unavoidable numerical errors in all Nbody simulations. Here we discuss the source of numerical errors in simulations with Myriad and also how the total CPU time is divided between the different parts of the code.
The precollapse error does not depend on the binary evolution since only some shortlived close encounters and soft binaries need the special treatment of the binary module of the code. This error depends mainly on the accuracy parameter η as described in Sect. 3.1. The postcollapse error is more complex to analyze, because longliving hard binaries form and interact with individual particles and other binaries. In some cases, multiple subsystems form. All of these uncertainties need to be treated very carefully. Since the error in the evolution of isolated hard binaries and multiples can be very small as described in Sect. 3.3, the relatively large errors in the energy after core collapse (see Figs. 19, 21, 25, 27, and 29) come mainly from the interactions of binaries and multiples with the rest of the system. The external perturbations in a binary are controlled by the dimensionless parameter γ_{pert} discussed in Sect. 2.2.3. This parameter plays role in the total energy error, since if it is large enough, some of the perturbations are ignored and systematic errors are introduced into the binary evolution and the evolution of stars that lie close to the binary. Our choice of γ_{pert} is 10^{7}, which probably needs adjustment to systems containing stars with large mass ratios.
Fig. 21 Cumulative energy error (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as a function of time for the simulation described in Fig. 19. The dashed line is the energy error of Starlab, while the heavy black line is the energy error of Myriad. 
When a binary or multiple subsystem forms according to the rules presented in Sect. 2.2.3, the transformation of coordinates from the center of mass of the cluster to the local center of mass of the binary or multiple, introduce another source of energy error. If that binary is a soft binary, it is usually required to form close to periastron passage and is deformed close to apastron passage once every period, until it transforms into a hard binary or its eccentricity becomes more than 1. This continuous transformation of coordinates introduces systematic errors into the simulation. Systematic errors are also introduced when a third star orbits with high eccentricity around the center of mass of a hard binary. According to the rules of Myriad, the triple system would have to form and deform often depending on the hardness of the hard binary and the eccentricity of the third star. In this case, another source of error may introduce large errors into the simulation: the third star would have to jump from small time steps, when it is a member of the triple, to large time steps, when it is considered to be an escaper from the triple and evolves as a single star in the whole system.
Fig. 22 CPU time for the simulation as a function of simulation time. The continuous line is the total CPU time of Myriad, while the heavy dashed line is the total CPU time of Starlab. The dasheddotted line shows the CPU time taken for the force calculation and neighbor list creation done on GRAPE6. The dotted line is the CPU time required for the evolution of binary systems. 
Finally, it is not clear what the optimal parameters are for the construction and destruction of a binary or a multiple system, because there are many different cases. A multiple subsystem for example may have formed around a soft binary or hard binary or even a close encounter between two stars may occur, and a particle can escape the subsystem under different conditions. The parameters used in Myriad are those described in Sect. 2.2.3, but they need further examination to reduce systematic energy errors. For example, it is unclear when an interaction between 2 hard binaries can lead to the formation of a multiple subsystem with 4 stars and if such a multiple is formed, it is not easy to determine under the conditions under which it would be destroyed forming smaller subsystems.
Fig. 23 Core massdensity (in Nbody units) (top) and number of stars in the core (bottom). At the time of core collapse the core density reaches a maximum value. 
From the above, it is clear that the main source of error in the postcollapse part of a simulation, where hard binaries form in a dense environment, depends on the treatment of the interactions between the binaries and their surrounding environment and on the choices we make for the rules of creation and dissipation of multiples. The sudden jumps in the time step of particles that interact with hard binaries are another important source of error. Finally, the dense environment itself introduces large errors, because there are relatively many stars sharing small time steps and more time steps per time unit, introduce large errors in the energy.
CPU time:
The total CPU time T_{CPU} taken for a simulation with Myriad can be divided into the CPU time needed for the different modules of the code (45)where T_{herm} is the CPU time for the predictor, corrector, and finding the time steps of the Hermite 4 scheme. In the simulation of an 1K system up to and beyond core collapse presented in Sect. 3.4.1, this time is less than 1% of the total CPU time.
The parameter T_{force} is the CPU time required by GRAPE6 to calculate the accelerations and identify the neighbors of particles. We reduced this CPU time by asking GRAPE6 to search for and return neighbors only for particles that have small time steps. We restricted the search for neighbors to the particles that belong to the first 3 time blocks, without any change in the accuracy of the code. For the simulation of Sect. 3.4.1, T_{force} ≃ 51%T_{cpu}. As the number of particles in a simulation increases, T_{HERM} is expected to increase according to the rule Nlog N.
The parameter T_{bin} is the CPU time spent in the evolution of binary and multiple systems that form in a cluster. It depends on the number of binaries formed, on the time step used for their evolution, and on the number of stars that are considered as perturbers of a binary system. If primordial binaries exist, this CPU time is expected to decelerate the simulation considerably. In the simulation of Sect. 3.4.1 where no primordial binaries exist, and hard binaries formed just before core collapse, T_{bin} ≃ 42%T_{cpu}. In the same simulation, T_{bin} before core collapse is insignificant.
T_{check} is the CPU time needed to check for encounters and binaries in the system. We limited this check to particles that have small time steps, since the time step given by Eq. (8) is an indicator of the dynamics around a particle. If at a given time the time step is small, that means that there is a probability that the particle is about to undergo a close encounter with some other particle, and a search for neighbors is required. If the timestep is large, then the particle probably lies in a relatively sparce environment and a close encounter with other particles is improbable, at least at the current time. T_{check} depends on the density of stars at the cluster center and for this reason grows for clusters that are close to core collapse. In the simulation described in Sect. 3.4.1, we found T_{check} ≃ 4%T_{cpu}.
Finally, T_{diag} is the CPU time required to calculate the global parameters of the cluster, such as the core radius, remove escaping stars, and write output snapshots to a file. The most “expensive” in terms of CPU time is the calculation of the core radius of the system, because it evaluates the density around each particle. GRAPE6 is used for this calculation, but the process is still relatively slow. For the simulation of Sect. 3.4.1, it is T_{diag} ≃ 2%Tcpu.
Fig. 24 Core radius (heavy black line) and Lagrangian radii containing 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, 3%, 2%, 1%, and 0.1% (blue dashed lines from top to bottom) of the total mass. The halfmass radius is indicated by a heavier blue dashed line. The initial condition was a Plummer model of 16 384 stars. A Salpeter mass function with mass limits m_{lower} = 0.5 M_{⊙} and m_{upper} = 5 M_{⊙} was chosen as the initial mass function. Core collapse is reached at t_{cc} ≃ 634 Nbody time units or t_{cc} ≃ 2.2t_{rlx}. The simulation ended at t ≃ 660 Nbody time units and took ~5 days. 
Fig. 25 Error in the energy per quarter of the time unit (δE = E_{(t)} − E_{(t − 0.25)}) as it evolves in time for the simulation of Fig. 24. 
Fig. 26 Core radius (heavy black line) and Lagrangian radii containing 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, and 3% (dashed lines from top to bottom) of the total mass. The half mass radius is indicated by a heavier dashed line. The initial condition was a King model with W_{0} = 6 containing 6144 stars. A Scalo mass function with mass limits m_{lower} = 0.1 M_{⊙} and m_{upper} = 100 M_{⊙} was chosen as the initial mass function. Core collapse is reached at t_{cc} ≃ 20.6 Nbody time units or t_{cc} ≃ 0.21t_{rlx}, which is the time of formation of the first hard binary. 
Fig. 27 Cumulative error in the energy (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as it evolves in time for the simulation of Fig. 26. 
Fig. 28 Core radius (heavy black line) and Lagrangian radii containing 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, and 3% (dashed lines from top to bottom) of the total mass. The half mass radius is indicated by a heavier dashed line. The initial condition was a King model with W_{0} = 6 containing 12 288 stars. A Scalo mass function with mass limits m_{lower} = 0.1 M_{⊙} and m_{upper} = 100 M_{⊙} was chosen as the initial mass function. Core collapse is reached at t_{cc} ≃ 37.2 Nbody time units or t_{cc} ≃ 0.21t_{rlx}, which is the time of formation of the first hard binary. 
Fig. 29 Cumulative error in the energy (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as it evolves in time for the simulation of Fig. 28. 
5. Conclusions
We have developed a new C++ Nbody code called Myriad. The code can simulate the evolution of star clusters with excellent accuracy, while its computational speed is satisfactory. The accumulated relative error in the total energy in the simulation of N = 16384 stars with a Salpeter IMF (m_{upper} = 5 M_{⊙}, m_{lower} = 0.5 M_{⊙}), evolved up to and beyond core collapse, is smaller than 2.5 × 10^{5}, while the error in the energy per quarter of Nbody time unit is smaller than 6 × 10^{6}. The computational time needed for this simulation was less than 3 days, which is comparable to the computational time of most Nbody codes,
but could be improved. The results for the core collapse time of star clusters with an IMF or with equalmass stars are also comparable with those found in the literature or produced by other codes such as Starlab. In Anders et al. (2009), Starlab and NBODY4 are compared in detail, and Myriad successfully reproduces the result for the core collapse time of an equalmass Plummer model consisting of 1K particles.
We now summarize the most important innovations introduced by Myriad:

1.
The accuracy parameter η_{b} for the evolution of binary and multiple subsystems is adjusted for every one of them by the time of its formation, while in other codes it is fixed to some value. Variable η_{b} makes the transition of time steps between the H4 scheme and the binary module smooth. In the future, η_{b} may change during the binary evolution to retain the accuracy in desirable limits.

2.
The algorithms of Kaplan et al. (2005) and Kokubo et al. (1998) have been successfully applied to the P(EC)^{3} scheme for the evolution of binary and multiple subsystems. The error in the energy of these systems does not increase linearly with time, but remains bounded within some limits in most of the cases.

3.
We have introduced new rules for creating triple systems when single stars are approaching binaries. We do not use any distance criterion, but the dimensionless perturbation parameter γ defined in Eq. (30), which depends on the hardness of the binary, the distance of the third star and their masses. We use the same parameter to remove a star from a multiple subsystem.
There are 4 important parameters of Myriad that balance between speed and accuracy. These are:

1.
The accuracy parameter η, used for Aarseth’s time step criterion given in Eq. (8). This is set by the user before beginning a simulation with Myriad. A typical value of η is 0.01.

2.
The accuracy parameter η_{b} for the time step criterion in Eq. (21), used for the evolution of binary and multiple subsystems. This parameter is adjusted during runtime until a smooth transition between the time step of the binary and the Hermite integrator is achieved.

3.
The critical value of the perturbation parameter γ_{crit}, discussed in Sect. 2.2.3, above which a particle that perturbs a binary system, becomes a member of the same binary, forming a triple system. The value of this parameter is set in Myriad to be 0.125 for soft binaries and 0.015625 for hard binaries.

4.
The critical value of the perturbation parameter γ_{pert} above which a particle is considered as a perturber of a binary or multiple subsystem (see Sect. 2.2.3: Perturbers). A typical value of this parameter that is used throughout the present work is 10^{7}, but it is automatically increased if it results in a large number of perturbers per binary. The maximum number of perturbers per binary system is set to 100.
From the results of the present work, we conclude that Myriad has the basic structure of a code that could be used to model the evolution of star clusters with IMBHs at their centers and study the effects of the presence of an IMBH on the system and the possibility of its ejection due to encounters and collisions with other stellarmass black holes. To achieve this, additional improvements to the code will be required so that more physical phenomena are included in the simulations. One of the improvements is to include postNewtonian dynamics into the equations for the evolution of close binary black holes. Because of the modular structure of the code, this improvement could be made by simply adding a function under the binary class. This function would calculate the appropriate postNewtonian terms (up to 2.5 PN) for the acceleration and its first derivative, following the formalisms found in Blanchet (2002), Berentzen et al. (2008), and Kupi et al. (2006). Another extension that would make the binary class more complete in treating binary black holes is the inclusion of the relativistic effect of the asymmetric emission of gravitational radiation and recoil velocity, during the merging of two orbiting black holes. According to numerical relativity results, when two spinning black holes collide, gravitational radiation could be emitted asymmetrically. This would lead to a recoil velocity in the resulting black hole that might be as high as 4000 km s^{1} (Gonzalez et al. 2007; Campanelli et al. 2007b,a; Healy et al. 2009; Herrmann et al. 2007), depending on the mass ratio of the initial black holes and the directions of their spins, but this velocity might be significantly suppressed by the relativistic alignment of the spins (Kesden et al. 2010). It is not easy, of course to include numerical relativity in an Nbody code, but semianalytical formulae, coming from fitting between numerical relativity results and postNewtonian theory (Lousto & Zlochower 2008; Baker et al. 2008; Lousto et al. 2010) to determine the direction and magnitude of the recoil velocity. These semianalytical formulas could be easily included in the binary module of Myriad, making the code a tool capable of reproducing collisions of black holes realistically. Those two extensions include effects beyond the Newtonian theory of gravitation, and with extensive tests and comparisons to theoretical and already published results, will be presented in a future paper. After some development, Myriad should be able to study the behavior of BHs that come from mergers in a cluster of stars, the possible changes in the structure of the cluster, and the possibility of the escape of IMBHs from clusters (HolleyBockelmann et al. 2007), after collisions with stellarmass BHs.
Additional improvements to the code could be made by adding stellar and binary evolution and the influence of the Galactic tidal field. This would make the simulation more realistic, because mass loss from stellar evolution and the tidal force of the Galaxy play an important role in the evolution of the whole cluster. Because of the code structure, the addition of stellar evolution may be achieved by simply including the appropriate function as part of the particle class, while binary evolution could be included as a simple function under the binary class.
Finally, improvements could be made to the speed of the code, by connecting the particle classes in a binary search tree. The search for close encounters and perturbers would then become faster and more accurate, since GRAPE6 memory overflow would be avoided. Another change in the structure of the code would be to divide the calculations between the CPU, GRAPE6, and the GPU, such that all three parts of the hardware cooperate in making the code even faster. A version of the code that uses the GPUs instead of GRAPE6 to calculate the forces, is under development. Until now, this version has usee the Sapporo library (Gaburov et al. 2009) to emulate the GRAPE6 functions on the GPUs and a shared time step for all particles. Results of preliminary tests show that, in agreement with the introductory paper of Sapporo, the speed of the code using 4 Nvidia Tesla C1060 GPUs is comparable with the speed of the code when GRAPE6 is used. Parallelization of the code, so that it could run on clusters of computers attached to accelerating hardware of the GRAPE6 family or GPUs, is an improvement that needs to be made to be able to simulate realistically large starsystems.
http://www.ids.ias.edu/~starlab/ also see http://muse.li for a set of Nbody simulation tools.
Myriad is an ancient Greek unit for 10 000, which is the order of magnitude of the size of a star cluster that can be evolved using the code, in reasonable time. In modern English, the word refers to an unspecified large number http://en.wikipedia.org/wiki/Myriad
Myriad will be soon available for download on http://www.astro.auth.gr/~simos/mediawiki1.6.7/index.php/Myriad
Acknowledgments
We would like to acknowledge the authors of Starlab for making their code freely available, and the creators of the Art of Computational Science website for their freely available codes and instructive books. We also would like to thank Kleomenis Tsiganis, for his very helpful suggestions and the anonymous referee for his/her instructions that helped us improve and correct the manuscript. This work was supported by the German Science Foundation (DFG) via SFB/TR7 on “Gravitational Waves” and by the Greek GSRT Pythagoras I program. The use of the GRAPE6 system was supported by a grant by the Empirikion Foundation. S.K. would like to thank the EberhardKarls University of Tübingen, where part of the code was written, for hospitality and support.
References
 Aarseth, S. J. 1963, MNRAS, 126, 223 [NASA ADS] [Google Scholar]
 Aarseth, S. J. 1968, Bull. Astron, 3, 105 [Google Scholar]
 Aarseth, S. J. 1999a, PASP, 111, 1333 [NASA ADS] [CrossRef] [Google Scholar]
 Aarseth, S. J. 1999b, Celestial Mechanics and Dynamical Astronomy, 73, 127 [Google Scholar]
 Aarseth, S. J. 2003, Gravitational NBody Simulations (Cambridge, UK: Cambridge University Press) [Google Scholar]
 AmaroSeoane, P., Miller, M. C., & Freitag, M. 2009, ApJ, 692, L50 [NASA ADS] [CrossRef] [Google Scholar]
 Anders, P., Baumgardt, H., Bissantz, N., & Portegies Zwart, S. 2009, MNRAS, 395, 2304 [NASA ADS] [CrossRef] [Google Scholar]
 Baker, J. G., Boggs, W. D., Centrella, J., et al. 2008, ApJ, 682, L29 [Google Scholar]
 Baumgardt, H., Makino, J., & Hut, P. 2005, ApJ, 620, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Baumgardt, H., De Marchi, G., & Kroupa, P. 2008, ApJ, 865, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Belleman, R. G., Bédorf, J., & Portegies Zwart, S. F. 2008, New Astron., 13, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Berentzen, I., Preto, M., Berczik, P., Merritt, D., & Spurzem, R. 2008, Astron. Nachr., 329, 904 [NASA ADS] [CrossRef] [Google Scholar]
 Blanchet, L. 2002, Living Rev. Relativity, 5 [Google Scholar]
 Bowers, R. L., & Deeming, T. 1984, Astrophys., 1 – Stars (Boston, MA: Jones and Bartlett) [Google Scholar]
 Campanelli, M., Lousto, C. O., Zlochower, Y., & David, M. 2007a, ApJ, 659, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Campanelli, M., Lousto, C. O., Zlochower, Y., & David, M. 2007b, Phys.Rev.Lett., 98, (231102) [Google Scholar]
 Freitag, M., & Benz, W. 2001, A&A, 375, 711 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaburov, E., Harfst, S., & Portegies Zwart, S. 2009, New Astron., 14, 630 [Google Scholar]
 Giersz, M., & Heggie, D. C. 1991, MNRAS, 268, 257 [Google Scholar]
 Gonzalez, J. A., Hannam, M., Sperhake, U., Brügmann, B., & Husa, S. 2007, Phys. Rev. Lett., 98 [Google Scholar]
 Gualandris, A., & Merritt, D. 2008, ApJ, 678, 780 [NASA ADS] [CrossRef] [Google Scholar]
 Hamada, T., & Iitaka, T. 2007, submitted to New Astron., 5 Mar., http://arxiv.org/abs/astroph/0703100 [Google Scholar]
 Harfst, S., Gualandris, A., Merritt, D., & Mikkola, S. 2008, MNRAS, 389, 2 [Google Scholar]
 Healy, J., Herrmann, F., Hinder, I., et al. 2009, Phys. Rev. Lett., 102, (0411101) [CrossRef] [Google Scholar]
 Heggie, D. C. & Mathieu, R. D. 1986, in The Use of Supercomputers in Stellar Dynamics, ed. P. Hut, & S. McMillan, Proceedings of a Workshop Held at the Institute for Advanced Study, 267, 233 [Google Scholar]
 Heggie, D. C., Trenti, M., & Hut, P. 2006, MNRAS, 368, 677 [Google Scholar]
 Herrmann, F., Hinder, I., Shoemaker, D. M., Laguna, P., & Matzner, R. A. 2007, Phys. Rev. D, 76, (084032) [CrossRef] [Google Scholar]
 HolleyBockelmann, K., Gultekin, K., Shoemaker, D., & Yunes, N. 2007, ApJ, 686, 829 [Google Scholar]
 Holmberg, E. 1941, ApJ, 94, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 2003, in Astrophysical Supercomputing using Particle Simulations, IAU Symp. 208, held 10–13 July 2001 in Tokyo, Japan, 331 [Google Scholar]
 Kaplan, M., Saygin, H., Hut, P., & Makino, J. 2005, unpublished, [arXiv:astroph/0511304] [Google Scholar]
 Kesden, M., Sperhake, U., & Berti, E. 2010, ApJ, 715, 1006 [NASA ADS] [CrossRef] [Google Scholar]
 King, I. R. 1966, Astron. J., 71, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Kinoshita, H., Yoshida, H., & Nakai, H. 1990, Celest. Mech. Dyn. Astron., 50, 59 [Google Scholar]
 Kokubo, E., Yoshinaga, K., & Makino, J. 1998, MNRAS, 297, 1067 [NASA ADS] [CrossRef] [Google Scholar]
 Kupi, G., AmaroSeoane, P., & Spurzem, R. 2006, MNRAS, 371, L45 [NASA ADS] [Google Scholar]
 Lousto, C. O., & Zlochower, Y. 2008, Phys. Rev. D, 79 [Google Scholar]
 Lousto, C. O., Campanelli, M., & Zlochower, Y. 2010, Class. Quant. Grav., 27, 114006 [Google Scholar]
 Makino, J. 2008, in Dynamical Evolution of Dense Stellar Systems, Proceedings of the International Astronomical Union, IAU Symp., 246, 457 [NASA ADS] [Google Scholar]
 Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141 [NASA ADS] [Google Scholar]
 Makino, J., Fukushige, T., Koga, M., & Namura, K. 2003a, PASJ, 55, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Makino, J., Fukushige, T., Koga, M., & Namura, K. 2003b, PASJ, 55, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Makino, J., Taiji, M., Ebisuzaki, T., & Sugimoto, D. 1997, ApJ, 480, 432 [NASA ADS] [CrossRef] [Google Scholar]
 McMillan, S. L. W. 1985, in The Use of Supercomputers in Stellar Dynamics, ed. P. Hut, & S. McMillan, Proceedings of a Workshop Held at the Institute for Advanced Study, 156 [Google Scholar]
 Nitadori, K., & Makino, J. 2008, New Astron., 13, 498 [NASA ADS] [CrossRef] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [CrossRef] [Google Scholar]
 Portegies Zwart, S., & McMillan, S. 2003, New Horizons in Globular Cluster Astronomy, ASP Conf. Proc., 296, 85 [NASA ADS] [Google Scholar]
 Portegies Zwart, S. F., Belleman, R. G., & Geldof, P. M. 2007a, New Astron., 12, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., Makino, J., McMillan, S. L. W., & Hut, P. 1999, A&A, 348, 117 [NASA ADS] [Google Scholar]
 Portegies Zwart, S. F., McMillan, S. L. W., & Makino, J. 2007b, MNRAS, 374, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Portegies Zwart, S. F., McMillan Stephen, L. W., Hut, P., & Makino, J. 2001, MNRAS, 321, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
 Scalo, J. M. 1986, Fund. Cosmic Phys., 11, 1 [NASA ADS] [EDP Sciences] [Google Scholar]
 Spitzer, L. 1987, Dynamical evolution of globular clusters (Princeton, NJ: Princeton University Press) [Google Scholar]
 Spurzem, R. 1999, J. Comput. Appl. Math., 109, 407 [NASA ADS] [CrossRef] [Google Scholar]
 von Hoerner, S. 1960, Z. Astrophys., 50, 184 [Google Scholar]
 von Hoerner, S. 1963, Z. Astrophys., 57, 47 [Google Scholar]
Appendix A: Units and basic equations
The units adopted in the simulations are the usual Nbody units^{9} (Heggie & Mathieu 1986), where G = M_{t} = R_{V} = 1, and G is the gravitational constant, M_{t} is the total mass of the system, and R_{V} is its virial radius. In these units, the total energy is E_{t} = −1/4. Transformation to physical units can be made if the total mass m_{t} and the virial radius r_{v} are known in physical units. If the mass m_{t} is given in solar masses [M_{⊙}] and r_{v} in parsecs [pc], the mass m_{i} of a star i in solar masses (A.1)Any distance R in Nbody units is transformed to [pc] using (A.2)In addition, there are functions for transforming time and velocity from Nbody (T,V) to physical units (t,v) (A.3)and (A.4)where (A.5)and (A.6)The numerical factors 6.557 × 10^{2} and 14.94 come from expressing 1 × 10^{5}(GM_{⊙}/L^{ * }) and (L^{ * }^{3}/GM_{⊙})^{1/2} in cgs units, respectively, where L^{ * } is taken to equal 1 pc and G and M_{⊙} are expressed in cgs units.
In the following, we follow the definitions of Aarseth (2003).
Halfmass radius:
The half mass radius r_{h} is defined as the radius of the sphere that has its center at the center of mass of the star cluster and contains half of the mass of the total system.
Halfmass relaxation time:
The equation for the halfmass relaxation time is adopted from Aarseth (2003) and Spitzer (1987)(A.7)where is the mean mass and lnΛ is the Coulomb logarithm (Spitzer 1987) given by (A.8)where N is the number of stars and the factor γ is usually chosen to be 0.4 (Aarseth 2003), but in some cases values of around 0.1 may be more appropriate (Giersz & Heggie 1991). A typical value of lnΛ used in Nbody simulations is 10.
Crossing time:
The time to cross the cluster is given by (A.9)where is the mean mass of the system. In Nbody units, it is . Equation (A.9) is derived from (A.10)where σ is the rms velocity dispersion given at equilibrium by (A.11)
Core radius:
The core radius is defined by the equation (A.12)where ρ_{i} is the local density around star i given by (A.13)and r_{6} is the distance to the sixth nearest neighbor of star i, while r_{d} is the position of the density center whose definition follows.
Density center:
The density center of the star cluster is defined to be (A.14)
Core collapse time:
For an equalmass system, the time for core collapse is (Spitzer 1987) (A.15)For a system containing stars with different masses distributed according to a King density profile (King 1966), the core collapse time becomes (Portegies Zwart & McMillan 2003) (A.16)In an Nbody simulation, a clearer definition of the initial core collapse time is sometimes required since in some cases it is unclear when exactly the core shrinks to a minimum size and then starts its expansion. In many publications and for more accurate comparison between different codes (Anders et al. 2009), the core collapse time is considered to be the time of the formation of the first longliving hard binary with binding energy greater than 100 kT, where kT is the equivalent to the thermal energy of a gas given for a cluster of N stars by (A.17)and E_{kin} is the total kinetic energy of the cluster. In other publications (Heggie et al. 2006), linear fitting is applied to the diagram of r_{c} versus time and the core collapse time is determined by the intersection of two lines, one being the fitting line of data before the core reaches its first minimum and the other being the fitting line after the core reaches its minimum.
In some other publications, the core collapse time is simply found by observing the diagram r_{c} versus time where there is a clear and sharp minimum.
In this work, the time of core collapse is determined both by simply observing the diagram of the evolution of core radius with time and by using the formation of the first hard binary with binding energy greater than the 100 kT criterion.
All Tables
All Figures
Fig. 1 Class hierarchy in Myriad. A class cluster is a set of several block classes. A block contains a variable number of particle and binary classes. A binary is a set of two or more particle classes. A class particle contains information (mass, position, velocity etc.) for a single star and functions that act on them. The class cluster contains information about the system as a whole (total mass, energy, halfmass radius, core radius etc.) and the required functions for calculating them. 

In the text 
Fig. 2 Simplified graphical representation of Myriad. The arrows show the flow of the data in the code. Boxes represent input or output datafiles or outside applications, while circles represent sets of codefunctions. Myriad applications are those inside the large dashed box. The other boxes and circles refer to satellite programs that create initial data or make animations from Myriadoutput. 

In the text 
Fig. 3 Schematic data flow between the CPU and the GRAPE6 for the Hermite 4thorder integrator. See text for description. 

In the text 
Fig. 4 Distribution of equal mass particles in different time steps (powers of 2). Aarseth’s criterion is used for the time step with η = 0.01. A small number of particles lies in time step 2^{14}. This is the smallest allowed time step for the Nbody code. Stars with this time step are candidates for close encounters or virtual particles representing the centers of mass of binary or multiple subsystems 

In the text 
Fig. 5 Representation of Hermite 4th order with block time steps. In the figure, a set of 7 particles distributed in 3 timeblocks is shown. See text for description. 

In the text 
Fig. 6 Perturbers. Stars with γ ≥ γ_{pert} are considered as perturbers of a binary system. Those are the stars represented as empty (red) circles in the above image. The filled (black) circles represent the members of the binary or multiple subsystem. If the dimensionless perturbation γ, given by Eq. (30), becomes greater than γ_{crit}, then the star is included into the binary or multiple subsystem. More stars could be added forming a multiple system with up to ~10 stars. Those systems are usually dispersed quickly. On the other hand, a star with γ becoming smaller than γ_{pert} is no longer a perturber. Those are the stars represented with filled (blue) squares. For an equal mass system, we typically set γ_{pert} ~ 10^{7} and γ_{crit} = 0.125 if the inner binary is considered as “soft”, while γ_{crit} = 0.015625 is assumed if it is considered as “hard”. 

In the text 
Fig. 7 Cumulative relative energy error as a function of η for a Plummer model with N = 8192 (top left), N = 16384 (top right), N = 32768 (bottom left), and N = 65536 (bottom right). The integrations ended at t = 1 Nbody units. 

In the text 
Fig. 8 Wallclock elapsed time as a function of N and η for the simulations of Fig. 7. All simulations ended at t = 1 in Nbody units. 

In the text 
Fig. 9 Wallclock elapsed time as a function of the number of stars N. All simulations ended at t = 1 in Nbody units. The accuracy parameter η had the same value 0.01 in all simulations. The slopes of O(N^{3}), O(N^{2}), and O(Nlog N) are also shown for comparison. 

In the text 
Fig. 10 Cumulative energy error (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as a function of time for the first 10 periods for the simulations of equalmass binary systems with 4 different eccentricities. The eccentricities of the binaries are e = 0.19 (red line), e = 0.51 (green dashed line), e = 0.75 (blue dotted line), and e = 0.91 (cyan dasheddotted line). 

In the text 
Fig. 11 Variation in the time step used for the evolution of the binary systems of Fig. 10 during 2 periods of each binary. The time step parameter η_{b} is chosen so that the error in all the binaries is of the same order. 

In the text 
Fig. 12 Same as Fig. 10, but for the duration of 10^{5} periods for each binary. 

In the text 
Fig. 13 Cumulative energy error (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as a function of time for 10^{5} periods for the simulation of a binary system with eccentricity e = 0.89 and mass ratio m_{heavy}/m_{light} = 15. 

In the text 
Fig. 14 The trajectory of the two stars in the simulation of Fig. 13. The two stars are orbiting around their common center of mass, which is located at O(0,0). 

In the text 
Fig. 15 Cumulative energy error (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as a function of time for the Pythagorean threebody system. 

In the text 
Fig. 16 Trajectories of the particles in the simulation of the Pythagorean threebody system. The initial positions of the particles are marked with black circles. 

In the text 
Fig. 17 Core radius (heavy black line) and Langrangian radii containing 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, and 3% (blue dashed lines from top to bottom) of the total mass. The half mass radius is indicated by a heavier dashed line. Core collapse is reached at t_{cc} ≃ 340 Nbody time units or t_{cc} ≃ 17t_{rlx}. 

In the text 
Fig. 18 Lagrangian radii containing from top to bottom 90%, 50%, 10%, and 1% of the total mass. The continuous lines are the results produced by Myriad, while the blue dashed lines are produced using Starlab. 

In the text 
Fig. 19 Error in the energy per quarter of the time unit (δE = E_{(t)} − E_{(t − 0.25)}) as it evolves in time for the simulation of an equalmass Plummer model of 1024 stars. The simulation ended at t ≃ 402 Nbody time units. The dashed line is the energy error of Starlab when it was used for evolving the same system. Before core collapse, the energy errors are identical. 

In the text 
Fig. 20 Respective energy error ratio Myriad/Starlab as a function of time. 

In the text 
Fig. 21 Cumulative energy error (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as a function of time for the simulation described in Fig. 19. The dashed line is the energy error of Starlab, while the heavy black line is the energy error of Myriad. 

In the text 
Fig. 22 CPU time for the simulation as a function of simulation time. The continuous line is the total CPU time of Myriad, while the heavy dashed line is the total CPU time of Starlab. The dasheddotted line shows the CPU time taken for the force calculation and neighbor list creation done on GRAPE6. The dotted line is the CPU time required for the evolution of binary systems. 

In the text 
Fig. 23 Core massdensity (in Nbody units) (top) and number of stars in the core (bottom). At the time of core collapse the core density reaches a maximum value. 

In the text 
Fig. 24 Core radius (heavy black line) and Lagrangian radii containing 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, 3%, 2%, 1%, and 0.1% (blue dashed lines from top to bottom) of the total mass. The halfmass radius is indicated by a heavier blue dashed line. The initial condition was a Plummer model of 16 384 stars. A Salpeter mass function with mass limits m_{lower} = 0.5 M_{⊙} and m_{upper} = 5 M_{⊙} was chosen as the initial mass function. Core collapse is reached at t_{cc} ≃ 634 Nbody time units or t_{cc} ≃ 2.2t_{rlx}. The simulation ended at t ≃ 660 Nbody time units and took ~5 days. 

In the text 
Fig. 25 Error in the energy per quarter of the time unit (δE = E_{(t)} − E_{(t − 0.25)}) as it evolves in time for the simulation of Fig. 24. 

In the text 
Fig. 26 Core radius (heavy black line) and Lagrangian radii containing 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, and 3% (dashed lines from top to bottom) of the total mass. The half mass radius is indicated by a heavier dashed line. The initial condition was a King model with W_{0} = 6 containing 6144 stars. A Scalo mass function with mass limits m_{lower} = 0.1 M_{⊙} and m_{upper} = 100 M_{⊙} was chosen as the initial mass function. Core collapse is reached at t_{cc} ≃ 20.6 Nbody time units or t_{cc} ≃ 0.21t_{rlx}, which is the time of formation of the first hard binary. 

In the text 
Fig. 27 Cumulative error in the energy (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as it evolves in time for the simulation of Fig. 26. 

In the text 
Fig. 28 Core radius (heavy black line) and Lagrangian radii containing 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, and 3% (dashed lines from top to bottom) of the total mass. The half mass radius is indicated by a heavier dashed line. The initial condition was a King model with W_{0} = 6 containing 12 288 stars. A Scalo mass function with mass limits m_{lower} = 0.1 M_{⊙} and m_{upper} = 100 M_{⊙} was chosen as the initial mass function. Core collapse is reached at t_{cc} ≃ 37.2 Nbody time units or t_{cc} ≃ 0.21t_{rlx}, which is the time of formation of the first hard binary. 

In the text 
Fig. 29 Cumulative error in the energy (ΔE/E_{0} = (E_{t} − E_{0})/E_{0}) as it evolves in time for the simulation of Fig. 28. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.