© 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 amino-acids 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 star-clusters, galaxies, and groups of galaxies.

N-body 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 time-step, a total of N(N − 1) ≈ N² force calculations should be made, where N is the number of interacting particles. From that, it is easy to conclude that N-body codes have time complexity O(N²) and require fast computers or even special purpose hardware, such as GRAPE systems, to simulate realistically large systems. This is why the evolution of N-body 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 N-body 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 N-body 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 block-time-step 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 N-body codes such as Starlab¹ (Portegies Zwart et al. 2001; Hut 2003), φGRAPE (Harfst et al. 2008), and Aarseth’s nbody² series (nbody4, nbody6, nbody6++ (Aarseth 1999a,b; Spurzem 1999)) use block time steps. Another milestone in the history of N-body 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 GRAPE-DR³ (Makino 2008), which is a more general purpose hardware specialized for N-body 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 host-computer. Today, GPUs are slowly replacing CPUs and the older GRAPEs in N-body 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 N-body codes with today’s software and hardware power. First of all, N-body 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 so-called 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 unequal-mass star clusters containing up to N = 131072 stars and an intermediate-mass 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⁸ − 10⁹), a simulation of a typical globular cluster (N ≈ 10⁵ − 10⁶) 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⁵ stars up to and beyond core collapse is possible, using a GRAPE-6 machine attached to a fast host-computer. 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⁶ − 10⁷).

During the past few years, N-body 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 N-body 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 Amaro-Seoane et al. (2009) used post-Newtonian effects on their N-body simulations to study the dynamics and the gravitational radiation from IMBH-binaries 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⁴, a new C++ N-body code for collisional and collisionless simulations. The code is loosely based on the instructive on-line books found in The Art of Computational Science website⁵. It has many of the features of existing codes, such as the use of block time steps, the Hermite 4th order integrator, GRAPE-6 support, and special treatment for binary and multiple sub-systems, but introduces a new structure and some new ideas especially in using GRAPE-6 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 sub-systems 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 so-called N-body 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

The core of Myriad⁶ uses C++ classes to store and manipulate data. The simplest or lower-level class is the class particle, which contains all necessary information about a single star⁷ (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 higher-level 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 highest-level class. The cluster-functions, 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 star-system, 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, half-mass radius etc.) at every diagnostic time step dt_diag.

2.2. Integration method

The integration method used in Myriad is the Hermite 4th order (H4) (Makino & Aarseth 1992) predictor-corrector scheme (PEC: Predict-Evaluate force-Correct) with block time steps. Accelerations and their derivatives (jerks) are computed using GRAPE-6. Binaries, close encounters, and multiple sub-systems that require small time steps are detected using GRAPE-6 and evolved using the time-symmetric 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)³).

For each block i, t_ic 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⁸ 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)³ integrator.

• ii.

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

• iii.

  We calculate the forces between binary-members on the CPU (O(k²)).

• 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 sub-systems(O(m)).

• (c)

  We evaluate the accelerations and their derivatives for all the particles of the block using the GRAPE-6 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 N-body code, with its former members and erase it.

• (i)

  We then send the corrected values to GRAPE-6 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 GRAPE-6 and which parts of the calculation are performed on the CPU and the GRAPE-6. Initially positions and velocities, and estimates of both the accelerations and the jerks of all particles are stored in GRAPE-6 memory (This procedure is omitted in the graph). Then, a prediction is made for a set of i-particles using the CPU. The predicted values and the time of each particle is sent to the GRAPE-6 buffer, which directly sends them to the internal predictor of the GRAPE-6. The internal predictor takes the data for the rest j-particles from the GRAPE-6 memory and predicts their positions and velocities at the time of the i-particles. The GRAPE-6 then computes the acceleration, its derivative (jerk), and the nearest neighbors for the i-particles and returns them to the corrector operating on the CPU. The corrector calculates the corrected values of position and velocity for the i-particles and sends the new data to the GRAPE-6 memory. A new set of i-particles 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 GRAPE-6 occurs three times per time step. The most time-consuming calculation, the calculation of the gravitational forces (O(Nlog (N))), performed on the GRAPE-6. The rest of the calculations are of order O(N). If the force-calculation were to be performed on the CPU, then it would be far more time-expensive as it would scale as O(N²).

2.2.1. Block time steps

The use of block time steps have proven to be an ideal method for N-body 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 GRAPE-6, because it lessens the communication time between the host and GRAPE-6. In individual time-step algorithms, a single particle is sent to GRAPE-6 everytime, while in block time step algorithms, groups of particles, that share the same time step are sent simultaneously to GRAPE-6 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 N-body 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ⁿ, 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¹⁴) stars.

2.2.2. The Hermite 4th order predictor-corrector scheme

The H4 scheme is the integrator used in many N-body 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 Runge-Kutta 4th order algorithm evaluates the force three times every time step, which makes this integrator much slower than H4 and consequently inappropriate for N-body simulations with large number of particles.

Another reason for choosing the H4 integrator is that GRAPE-6 is specifically designed for use with this kind of integrator. This is because GRAPE-6 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 GRAPE-6 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.

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

  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 GRAPE-6. The other calculations have complexity O(1).

• 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 time-blocks are updated using the H4 scheme. We note that in a typical star cluster, the number of time-blocks 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 predict-evaluate-and-correct every dt. Particles of block 1 predict every dt, but subjet to evaluate-and-correct every 2dt. Finally, particles of block 2 predict every dt and evaluate-and-correct 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 three-body 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 N-body 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 time-symmetric 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. Time-symmetry is achieved by applying the force evaluation-correction (EC) steps n times, after the prediction (P(EC)ⁿ) and by appropriate choice of the time step. The value n = 3 is used in the code because this is sufficient to achieve time-symmetry. Time-symmetry 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, non-symplectic 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 sub-systems is given at the beginning of Sect. 2.2, between items (i) and (vii).

The choice of the time step for the (P(EC)³) integrator is important to preserve time-symmetry. 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₀ and the end dt₁ 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 N-body 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 N-body 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 N-body system is a slow procedure with time complexity O(N²), unless all particle classes are members of a linked list, as in a binary search tree, or GRAPE-6 is used for acceleration. Here we use the neighbor-module of GRAPE-6 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 neighbor-module of GRAPE-6, 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 equal-mass 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 time-step 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 non-equal-mass system. Finally, the distance R_hi that is sent to GRAPE-6 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 N-body 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 rest-frame of the N-body 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)

In Eq. (30), , where M is given by Eq. (26), is the mean mass of the system, R_b is the size of the binary defined to be either its semimajor axis for hard binaries (i.e., binaries of binding energy greater than the average energy of a particle) or the distance R between its members in the case of soft binaries (i.e., binaries with binding energy lower than the average energy of a particle).

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 equal-mass 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 equal-mass 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 sub-systems 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 sub-system occurs when the following three criteria are met:

• 1.

  When the dimensionless perturbation γ_j acting on the particle-member of the sub-system 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 so-called “inner binary” is the hardest double sub-system inside the multiple (the one with the greatest binding energy).

In some rare cases, two particles-members 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 sub-system is controlled by the internal forces acting between its member stars. External perturbations from stars passing close to the center of mass of the sub-system 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 star-member 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 sub-system is done by using the neighbor-module of GRAPE-6. During the force calculation for an “virtual” star i, a distance (40)is sent to the neighbor-module of GRAPE-6. 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 GRAPE-6, each “virtual” star is sent alone to GRAPE-6, after cleaning the previous neighbor lists from GRAPE-6 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₆, where r₆ 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₆ 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₆. 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 GRAPE-6 memory does not overflow and GRAPE-6 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.

Figure 6 shows schematically the distribution of perturbers and non-perturbers 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 half-mass 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 equal-mass 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 N-body time unit. The speed of the code is measured by the CPU time required for the simulations. We also performed simulations of systems with equal-mass 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 mass-loss, and binary evolution are also not included. If close black-hole binaries are formed, they evolve using purely Newtonian dynamics, without any post-Newtonian 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 32-chip GRAPE-6 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 N-body 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⁴). In codes that use block time steps, such as Myriad, the parameter that controls the error is the accuracy parameter η is used in the time-step 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 GRAPE-6. Myriad uses double precision for all calculations on the CPU, while the accuracy of the calculation of the accelerations and derivatives is controlled by GRAPE-6 (Makino et al. 2003b).

We performed a series of experiments using Myriad to evolve equal-mass Plummer models, with different particle numbers N from t = 0 to t = 1 N-body time units. For all the integrations, we recorded the cumulative relative energy error (ΔE/E = (E_t = 1 − E_t = 0)/E₀) 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₀ ≤ 10^-9 for N ≤ 16 K and ΔE/E₀ ≤ 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 (wall-clock 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 equal-mass Plummer models with different particle numbers N from t = 0 to t = 1 N-body time units. For all simulations, we recorded the wall-clock time as a function of N and the results are shown in Fig. 9, where N is chosen to be between 512( = 2⁹) and 262144( = 2¹⁸). The accuracy parameter used for all simulations was η = 0.01. In the same figure, the slopes of N², N³, and Nlog ₁₀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 GRAPE-6 or, if GRAPE-6 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². 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 N-body time units. At later times, dynamically formed binary systems would force the simulation to run slower.

3.3. Binaries and multiples

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 time-symmetric 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⁵ 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 time-symmetric 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⁵ 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⁵ periods are shown in Fig. 14.

We finally tested the performance of Myriad in handling triple systems by evolving the well known Pythagorean three-body 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₁ = 3 is opposite the side of length 3, m₂ = 4 is opposite the side of length 4 and m₃ = 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₃ = 5 and m₂ = 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 N-body 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 GRAPE-6. 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 N-body 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 equal-mass 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 equal-mass 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 half-mass relaxation time (t_rlx). These results were compared with existing results found in the literature.

3.4.1. Equal-mass plummer models

We performed a simulation of an equal-mass 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 equal-mass 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 N-body time units, when core collapse occurs. From equation

(43)which is described in Appendix A, we find that for this system the initial half-mass relaxation time is t_rlx = 19.92 N-body time units, so t_cc ≃ 17t_rlx, which is close to what is expected from Eq. (A.15) and consistent with the results of other codes presented in Table 1 of Freitag & Benz (2001) and Anders et al. (2009), where a detailed comparison between the two major N-body codes Starlab and NBODY4 is presented. In Fig. 18, we compare the results of Myriad and Starlab for the evolution of selected Lagrangian radii of the same system. The small differences in the Lagrangian radii computed by the two codes, may be caused by the difference in the escaper-removal criterion and small differences in the time-step allocation criteria. Because of these differences, the simulation performed with Myriad ended with 983 stars remaining in the system, while that completed with Starlab contained 973. Figure 19 shows the evolution of the instantaneous energy error δE for this simulation. In the same figure, the simulation error of Starlab is presented for comparison. The error is initially small and unaffected by close encounters and binary systems. It grows slowly with time, as expected for the H4 integrator. As the system approaches core collapse, the interactions at the center become more frequent and violent and the energy error is contolled by them. Hard binaries form and their interactions with single stars or other binaries introduce another source of error. After core collapse, the average error becomes smaller, but some peaks, due to strong encounters between hard binaries, continue to appear.

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 N-body 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 time-step 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 post-collapse 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.

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 equal-mass 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₀ = 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.

4. Discussion

Sources of energy errors:

There are unavoidable numerical errors in all N-body 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 pre-collapse error does not depend on the binary evolution since only some short-lived 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 post-collapse error is more complex to analyze, because long-living 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.

When a binary or multiple sub-system 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.

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 sub-system 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 sub-system 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 sub-system 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 sub-systems.

From the above, it is clear that the main source of error in the post-collapse 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 GRAPE-6 to calculate the accelerations and identify the neighbors of particles. We reduced this CPU time by asking GRAPE-6 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 time-step 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. GRAPE-6 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.

5. Conclusions

We have developed a new C++ N-body 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 N-body 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 N-body codes,

but could be improved. The results for the core collapse time of star clusters with an IMF or with equal-mass 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 equal-mass 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 sub-systems 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)³ scheme for the evolution of binary and multiple sub-systems. 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 sub-system.

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 stellar-mass 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 post-Newtonian 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 post-Newtonian 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 N-body code, but semi-analytical formulae, coming from fitting between numerical relativity results and post-Newtonian theory (Lousto & Zlochower 2008; Baker et al. 2008; Lousto et al. 2010) to determine the direction and magnitude of the recoil velocity. These semi-analytical 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 (Holley-Bockelmann et al. 2007), after collisions with stellar-mass 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 GRAPE-6 memory overflow would be avoided. Another change in the structure of the code would be to divide the calculations between the CPU, GRAPE-6, 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 GRAPE-6 to calculate the forces, is under development. Until now, this version has usee the Sapporo library (Gaburov et al. 2009) to emulate the GRAPE-6 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 GRAPE-6 is used. Parallelization of the code, so that it could run on clusters of computers attached to accelerating hardware of the GRAPE-6 family or GPUs, is an improvement that needs to be made to be able to simulate realistically large star-systems.

---

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 GRAPE-6 system was supported by a grant by the Empirikion Foundation. S.K. would like to thank the Eberhard-Karls University of Tübingen, where part of the code was written, for hospitality and support.

References

Appendix A: Units and basic equations

The units adopted in the simulations are the usual N-body units⁹ (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 N-body units is transformed to [pc] using (A.2)In addition, there are functions for transforming time and velocity from N-body (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).

Half-mass 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.

Half-mass relaxation time:

The equation for the half-mass 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 N-body 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 N-body 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₆ 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 equal-mass 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 N-body 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 long-living 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