GPU-Enabled Particle-Particle Particle-Tree Scheme for Simulating Dense Stellar Cluster System

We describe the implementation and performance of the ${\rm P^3T}$ (Particle-Particle Particle-Tree) scheme for simulating dense stellar systems. In ${\rm P^3T}$, the force experienced by a particle is split into short-range and long-range contributions. Short-range forces are evaluated by direct summation and integrated with the fourth order Hermite predictor-corrector method with the block timesteps. For long-range forces, we use a combination of the Barnes-Hut tree code and the leapfrog integrator. The tree part of our simulation environment is accelerated using graphical processing units (GPU), whereas the direct summation is carried out on the host CPU. Our code gives excellent performance and accuracy for star cluster simulations with a large number of particles even when the core size of the star cluster is small.


Background
Direct N -body simulation has been the most useful tool for the study of the evolution of collisional stellar systems such as star clusters and the center of the galaxy [1]. The force calculations, of which the cost is O(N 2 ), are the most compute-intensive part of direct N -body simulations. Barnes and Hut [2] developed a scheme which reduces the calculation cost to O(N logN ) by constructing the tree structure and evaluating the multipole expansions. Dehnen [3,4] developed a scheme to reduce the calculation cost to O(N ) by combining the fast multipole method [5] and the tree code. Recently, the graphical processing units (GPU), which is a device originally developed for rendering the graphical image, began to be used for scientific simulations. The tree code is also implemented on GPUs and it is much faster than that on CPUs [6,7]. Bédorf et al. [8] parallelized the tree code on GPUs and showed good scalability up to 18600 GPUs. They also simulated the Milky Way Galaxy with N of up to 242 billion and reported that the average calculation time per iteration on 18600 GPUs was 4.8 seconds.
The tree schemes are widely used for collisionless system simulations. However, for collisional system simulations, the use of the tree code has been very limited. One reason might be that a collisional stellar system spans a wide range in timescales. Thus it is essential that each particle has its own integration timestep. This scheme is called the individual timestep or the block timestep [9]. However, when we use the tree code and the block timestep together, the tree structure is reconstructed at every block timestep, because the positions of integrated particle are updated. The cost of the usual complete reconstruction of the tree is O(N logN ) and not negligible.
To reduce the cost of the reconstruction of the tree, McMillan and Aarseth [10] introduced local reconstruction of tree. They demonstrated a good performance, but there seems to be no obvious way to parallelize their scheme.
Recently, Oshino et al. [11] introduced another approach to combine the tree code and the block timesteps which they called the P 3 T scheme. This scheme is based on the idea of Hamiltonian splitting [12][13][14][15][16][17][18]. In the P 3 T scheme, the Hamiltonian of the system is split into short-range and long-range parts and they are integrated with different integrators. The long-range part is evaluated with the tree code and is integrated using the leapfrog scheme with a shared timestep. The short range part is evaluated with direct summation and integrated using the fourth-order Hermite scheme [19] with the block timesteps. They investigated the accuracy and the performance of the P 3 T scheme for planetary formation simulations and showed that the P 3 T scheme achieves high performance.
In this paper, we present the implementation of the P 3 T scheme on GPUs and report its accuracy and performance for star cluster simulations. We found that the P 3 T scheme demonstrates a very good performance for star cluster simulations, even when the core of the cluster becomes small.
The structure of this paper is as follows. In section 2, we briefly describe the P 3 T scheme. In section 3, we report the accuracy and performance of the P 3 T scheme. We summarize these results in section 4.

Formulation
In this section, we describe the P 3 T scheme. The Hamiltonian H of a gravitational N -body system is given by where p i , m i and q i are momentum, mass and position of the particle i, respectively. To avoid the singularity of the 1/r potential, we use the Plummer softening ǫ [1]. With the P 3 T scheme, H is split into H hard and H soft as follows [11]: Here W (s ij ) is a smooth transition function. A suitable form of W (s ij ) should be zero when a distance between two particles is smaller than the inner cutoff radius r in and should be unity if the distance is larger than the outer cutoff radius r cut . This splitting is introduced by Chambers [15] to avoid undesirable energy error from close encounters between particles. Similar splitting has been used with P 3 M (Particle-Particle Particle-Mesh) scheme, in which the long-range part of the interaction is evaluated by using FFT [20].
Forces derived from H hard and H soft are given by We call K(s ij ) the cutoff function. The tree algorithm is used for the evaluation of F soft,i to reduce the calculation cost.
The formal solution of the equation of motion for the phase space coordinate w = (q, p) at time t + δt for the given Hamiltonian H is Here the braces {, } stand for the Poisson bracket. In the P 3 T scheme, we use the second order approximation; Here, the formal solution for the H soft term is the simple velocity kick, since H soft contains the potential only. We numerically integrate the H hard term, since it cannot be solved analytically. We use the fourth-order Hermite scheme with the block timestep [19]. The fourth-order integrator requires K(s ij ) to be three-times differentiable with respect to position. We use the following formula: This K(x) is the lowest-order polynomial which satisfies the requirement that derivatives up to the third order is zero for x = 0 and 1 (i.e. The highest-order term of the lowest-order polynomial is the seventh, because there are eight boundary conditions at x = 0 and x = 1).
Here η is the accuracy parameter of the timestep and its typical value is 0.1. ∆t max is the maximum timestep which should be smaller than ∆t soft , a (n) i is the nth time derivative of the acceleration of particle i, a 0 is a constant introduced to prevent ∆t i from becoming too small when the distance to the nearest neighbor is close to r cut and α is a parameter to control a 0 . In this case, the acceleration from H hard becomes very small and there is no need to use very small ∆t i . According to [11], when we chose α ≤ 1, α hardly affects the energy error. Thus we set α = 0.1 for all simulations.
In our Hermite implementation, a  i , and as a consequence we cannot use equation (18) for the first step.
We use: This criterion dose not contain the 2nd and 3rd time derivatives of the acceleration.
To prevent the timestep derived by equation (20) from becoming too large, we set η s to be the one-tenth of η for all simulation in this paper. We summarize all accuracy parameters in table 1.

Implementation on GPUs
Even with the Barnes-Hut tree algorithm, obtaining F soft,i is still costly and dominates the total calculation time [11]. To accelerate this part, we use GPUs, by modifying the sequoia library (Bédorf, Gaburov and Portegies Zwart, submitted to ComAC), on which the high-performance tree code for parallel GPUs Bonsai [7] is based. Our library calculates the long range forces on all particles, F soft,i by the Barnes-Hut tree algorithm (up to the quadrupole moment). On the other hand, we calculate F hard,i on the host computer. The library also returns, for each particle, the list of particles within the distance h from it. We use this list of neighbors to calculate F hard,i . The value of h should be sufficiently larger than r cut to guarantee that the particles which are not on the list of the neighbors of particle i do not enter the sphere of the radius r cut around particle i during the time interval ∆t soft . We call the sphere with a radius of r cut the neighbor sphere and the shell between the sphere with a radius of h and the neighbor sphere the buffer shell. The particles of which the nearest neighbor is outside the sphere with radius h are considered isolated and the particles on the list of neighbors are considered neighbor particles. We denote the width of the buffer shell as ∆r buff (i.e. h = r cut + ∆r buff ).
The compute procedures of our implementation of the P 3 T scheme on GPU is as follows: 1 Evaluate long range forces on all particles F soft,i using GPU. 2 Particles are divided into two groups; isolated and non-isolated, by using the neighbour list made on GPU. 3 For non-isolated particles, F hard,i are calculated on the host computer. 4 All particles receive a velocity kick through F soft,i for ∆t soft /2. 5 Isolated particles are drifted by r i ← r i + ∆t soft v i . 6 Non-isolated particles are integrated with the fourth-order Hermite scheme for ∆t soft . 7 Evaluate F soft,i and make the neighbour list in the same way as in step 1-2. 8 All particles obtain the velocity kick again for ∆t soft /2. 9 go back to step 3.

Accuracy and Performance
We performed a number of test calculations using the P 3 T scheme on GPUs, to study its accuracy and performance. In this section, we describe the result of these tests. For most of them we adopted a Plummer model [21] with 128K (hereafter K=2 10 ) equal-mass particles as the initial condition. We use the so-called N -body unit or Heggie unit, in which total mass M=1, the gravitational constant G=1 and total energy E = −1/4 [22]. To avoid the singularity of the gravitational potential, we use the Plummer softening and set ǫ = 4/N . Since this value is a typical separation of a hard binary in the N -body unit, we can follow the evolution of the system up to the moment of the core collapse.
Note, in this paper, we use the energy errors as an indicator of the accuracy of the scheme. However, energy conservation dose not guarantee accuracy of simulations (though, it is necessary). Thus we will perform realistic simulations in section 3.2 and check the statistical character of stellar systems by comparing the results with the Hermite scheme, which is widely used in collisional stellar system simulations. As we will see later, for simulations of the core collapse of the star cluster, when the relative energy error is 10 −3 at the moment of the core collapse, the behavior of the core collapse with the P 3 T scheme agreed with that with the Hermite scheme very well.

Accuracy
With the P 3 T scheme, we have six accuracy parameters. In sections 3.1.1.1-3.1.1.3, we discuss how each parameter controls the accuracy of the P 3 T scheme. In section 3.1.1.4, we describe the accumulation of the energy error in a long-term integration. To measure energy errors accurately, we calculate potential energies by the direct summation instead of the tree code for all runs in this paper.
3.1.1.1 Effect of r cut , ∆t soft and θ In figure 2, we present the maximum relative energy error |∆E max /E 0 | over 10 N -body time units as a function of r cut and ∆t soft for several different values of the opening criterion of the tree, θ. Here ∆E max is the maximum energy error and E 0 is the initial energy. We chose η = 0.1, ∆t max = ∆t soft /4 and ∆r buff = 3σ∆t soft , where σ is the global three dimensional velocity dispersion and we adopt σ = 1/ √ 2. We can see that the error is smaller for smaller θ, smaller ∆t soft , or larger r cut . Roughly speaking, the error depends on two terms, ∆t soft /r cut σ and θ. If ∆t soft /r cut σ is large, it determines the error. In this regime, the error is dominated by the truncation error of the leapfrog integrator. If it is small enough, θ determines the error, in other words, the tree force error dominates the total error. Even for a very small value of θ like 0.2, the tree force error dominates if ∆t soft /r cut σ 0.05.
In figure 3, we plot the maximum energy error as a function of θ. We use the same η, ∆t max and ∆r buff as in figure 2. For the runs with r cut = 1/256 and ∆t soft = 1/512, the energy error does not drop below 10 −6 because the error of the leapfrog integrator is larger than the tree force error. In an chaotic system like the model used in our simulations such energy error is sufficient to warrant a scientifically reliable result [23]. On the other hand, for the run with r cut = 1/128 and ∆t soft = 1/1024, integration error is smaller than the tree force error.
3.1.1.2 Effect of ∆r buff In figure 4, we show the maximum relative energy error as a function of ∆r buff for the runs with ∆t max = ∆t soft /4, η = 0.1, θ = 0.2, for (∆t soft , r cut ) = (1/512, 1/128) and (1/1024, 1/256). The energy error is almost constant for ∆r buff 2∆t soft σ, which indicates that the energy error for ∆r buff < 2∆t soft σ is caused by particles that are initially outside the buffer shell (with radius r cut + ∆r buff ) and plunge into the neighbour sphere (with radius r cut ) during the timestep ∆t soft . We can prevent this by adopting ∆r buff 2∆t soft σ.

Long term integration
In figure 6, we show the time evolution of the relative energy error until T =500. We compare the accuracy of our P 3 T scheme with two other schemes, the direct fourth-order Hermite scheme and the leapfrog scheme with the Barnes-Hut tree code. The calculations with the direct Hermite scheme are performed by using the Sapporo library on GPU [24], and the calculations with the leapfrog scheme are performed by using the Bonsai library on GPU [7]. The energy error of the P 3 T scheme behaves like a random walk whereas that of the leapfrog and the Hermite schemes grow monotonically. In the right-hand panels of figure 6, we show the same evolution of the error as in the left panels, but time is plotted with a logarithmic scale. This allows us to realize that the error growth of Hermite and tree schemes are linear, whereas the error in the P 3 T scheme grows as ∝ T 1/2 . This latter proportionality is caused by the short-term error of the P 3 T scheme, which is dominated by the randomly changing tree-force error. For long-term integration the P 3 T scheme conserves energy better than the Hermite or leapfrog schemes.

Calculation cost
In this section, we discuss the calculation cost of the P 3 T scheme and its dependence on the number of particles N , required accuracy, and other parameters.
We first construct a simple theoretical model of the dependence of the calculation cost on parameters of the integration scheme such as N , ∆t soft , θ and r cut in section 3.1.2.1. In section 3.1.2.2 we derive the optimal set of parameters from the model and compare this model with the result of the numerical tests. We found that the calculation cost per unit time is proportional to N 4/3 .

Theoretical model
The calculation cost for the force evaluations in P 3 T is split into the tree part and the Hermite part. For the tree part, the calculation cost of evaluating forces for all particles per tree step is proportional to O(θ −3 N logN ). Since we use constant timestep for the tree part, the calculation costs of the integration of particles per unit time for the tree part is proportional For the Hermite part, since each particle has its own neighbour particles and timesteps, the number of interactions for all particles per unit timstep is given by Here N ngh,i is the number of the neighbour particles around particle i, N step,i is the number of timesteps required to integrate particle i for one unit time, n i is the local density around particle i, ∆t i is the average timestep of particle i over one unit time and ∆t is the average of ∆t i over all particles. Here we assume n i is constant within the radius of r cut + ∆r buff around particle i.
Next we express the ∆t as a function of N and r cut . To simplify the discussion, we define the timestep of the particle through the relative position and velocity from its nearest neighbour particle; ∆t ∝ r NN /v NN , where r NN and v NN are the relative position and the velocity of the nearest neighbour particle. We can replace v NN to the velocity dispersion σ. Thus average timestep is given by To further simplify the derivation we assume that the number density of particles in the system is uniform. If r cut is larger than the mean inter-particle distance r (i.e. if most particles have neighbour particles), the average timestep is roughly given by where R is the typical size of the system. In this case, the average timestep depend only on N (dose not depend on r cut ). If r cut is small compared to r , most particles are isolated and most of the nonisolated particles have only one neighbour particle. In this case, ∆t is given by In figure 7 we show the number of steps per particle per unit time N step for a plummer sphere as a function of N (top panel) and as a function of r cut (bottom panel). In the top panel, we can see that N step is roughly proportional to N 1/3 for large N (i.e. r is small). On the other hand when N is small N step is almost constant because r is large [see equation (26)].
The bottom panel of figure 7 shows that all curves eventually approach to constant values for both of large and small r cut . For large r cut , the timesteps of the nonisolated particles are determined by N , not by r cut [see equation (25)], whereas for small values of r cut the non-isolated particles have a timesteps ∆t max . This is because most neighbouring particles are in the buffer shell and not in the neighbour sphere. For runs with ∆t soft =1/2048, 1/1024 and 1/512, we can see bumps of N step at r cut ∼ 1/512 due to the dependence on r cut shown in equation (26).
Using above discussions, the number of interactions for all particles per unit time of the Hermite part N int,hard and the tree part N int,soft are given by 3.1.2.2 Optimal set of accuracy parameters In this section, we derive the optimal values of r cut and ∆t soft from the point of view of the balance of the calculation costs between the tree and the Hermite parts, in other words we express r cut and ∆t soft as functions of N such that N int,hard /N int,soft is independent of N . Following the discussion in section 3.1.1.1 and because the energy errors can be controlled through ∆t soft /r cut , r cut should be proportional to ∆t soft . From section 3.1.1.2, ∆r buff should be also proportional ∆t soft . The requirements are met for N int,hard ∝ N 7/3 (r cut + ∆r buff ) 3 (or N 2 (r cut + ∆r buff ) 3 ), ∆t soft ∝ N −1/3 and r cut ∝ N −1/4 and both N int,hard and N int,soft are proportional to N 4/3 (or N 5/4 ). Here we have neglected the log N dependence in the tree part. This is illustrated in figure 8, where we plot N int,hard for a plummer sphere as a function of N . Following above discussions, we use the N -dependent tree timestep: ∆t soft = (1/256)(N/16K) −1/3 and N int,hard as well as N int,soft are proportional to N 4/3 .
In figures 9 and 10, we plot the wall-clock time of execution T cal and the maximum relative energy errors |∆E max /E 0 | for the time integration for 10 N -body units against N . Top (bottom) panel in figure 9 shows the results of the runs with r cut /∆t soft =2 (top panel) and 4 (bottom panel). All runs in these figures are carried out on NVIDIA GeForce GTX680 [1] GPU and Intel Core i7-3770K CPU. For each run, we use one CPU core and one GPU card.
We also perform the simulations using the direct Hermite integrator with the same η and the standard tree code with the same θ and ∆t soft . These calculations are performed with the Sapporo GPU library [24] and a standard tree code with the same θ and ∆t soft using the Bonsai GPU library [7]. The calculation time for our P 3 T implementation is also proportional to N 4/3 , as we presented in section 3.1.2.1, while for the Hermite integrator it is proportional to N 7/3 . The P 3 T scheme is faster than the direct Hermite integrator for N > 16K and when N =1M (M=2 20 ), the P 3 T scheme is about 50 times faster than the direct Hermite scheme. The pure tree [1] GTX680 does not have ECC (Error Check and Correct) memories. However, as we will see later, we do not observe any large energy error in all of our runs, which means the hardware error dose not affect our result. Betz, DeBardeleben and Walker [25] performed Molecular Dynamics simulations, in order to investigate the rate of bit-flip error events. They observed a single bit-flip error event in about 4700 GPU*hours without ECC and conclude that the bit-flip error is exceedingly rare. code is slightly faster than the P 3 T scheme, but the integration errors are worse by several orders of magnitude (see figure 6 and 10).

Examples of practical applications
In sections 3.1.1 and 3.1.2, we presented a detailed discussion on the accuracy and performance of our P 3 T scheme. However, we performed simple simulations, where the stellar systems are in the dynamical equilibrium. In this section, we study the performance of our P 3 T scheme when applied to more realistic, or more difficult, simulations by comparing the results of the Hermite scheme. In section 3.2.1, we discuss the case of the simulation of star clusters up to core collapse. In section 3.2.2, we discuss the case of a galaxy model with massive central black hole binary.

Star cluster down to core collapse
In this section, we discuss the performance of our P 3 T scheme for the simulation of the core collapse of a star cluster. In section 3.2.1.1, we describe the initial condition and parameters of the integration scheme. In section 3.2.1.2 we compare the calculation results obtained by the P 3 T and Hermite schemes, and in section 3.2.1.3 the calculation speed.

Initial conditions
We apply the P 3 T scheme to the evolution of a star cluster consisting of 16K stars to the moment of the core collapse [26]. We use an equal-mass plummer model as an initial density profile and we adopt η = 0.1. We apply the Plummer softening ǫ = 4/N = 1/4096. The simulations are terminated when the core number-density exceeds 10 6 , at which point the mean interparticle distance in the core is comparable to ǫ. Next, we set θ. We must chose θ so that the tree force error is smaller than the force due to the two-body relaxation. Hernquist et al. [27] pointed out that, for θ = 0.5 with monopole and quadrupole, the treeforce error is much smaller than the force due to the two-body relaxation. Thus we chose θ = 0.4 with quadrupole as a standard model. For comparison, we also perfrome a run with θ = 0.8.
To resolve the motions of the particles in the core, we impose ∆t soft to be smaller than 1/128 of the dynamical time of the core (∼ 3π/16ρ core , where ρ core is the core density). To reduce the calculation cost for the Hermite part we require r cut ∝ ρ −1/3 core and set the initial value of r cut = 1/64. We also change ∆r buff = 3σ core ∆t soft , where σ core is the velocity dispersion in the core, and ∆t max = ∆t soft /4, as ∆t soft and σ core are changing. Here, to calculate ρ core and σ core , we use the formula proposed by Casertano and Hut [28]. The same simulation is repeated using the fourth-order Hermite scheme with the block timesteps with the same value of η = 0.1.

Results
In figure 11 we present the evolution of the core densities ρ core (top panel) and the core radii r core (bottom panel) for P 3 T and Hermite schemes. For each scheme, we perform three runs, changing the initial random seed for generating the initial conditions of the Plummer model. The behaviors of the cores for all runs are similar. The differences between two schemes are smaller than run-to-run variations. Figure 12 shows the relative energy errors of the runs with the same initial seed as functions of the core density (top panel) and the time (bottom panel). The energy errors of the runs with P 3 T scheme change randomly, whereas those of the Hermite code grow monotonically. As a result, the P 3 T scheme with θ = 0.4 conserves energy better than the Hermite scheme in the long run. The errors for the P 3 T scheme with θ = 0.8 is slightly worse than that of the Hermite scheme, but the behavior of the core are similar with other runs. Thus the choice of θ = 0.4 is enough to follow the core collapse simulations. Figure 13 shows the calculation time of the P 3 T scheme (θ = 0.4) and Hermite scheme on GPU. As shown in this figure, the calculation time of the P 3 T scheme is dominated by the tree (soft) part calculation.

Calculation speed
Initially the P 3 T scheme is much faster than the Hermite scheme, but after the time when ρ core ∼ 10 4 , the P 3 T scheme is slightly slower than the Hermite scheme because in the P 3 T scheme, ∆t soft is proportional to ρ −1/2 core . However, even for the P 3 T scheme, the CPU time spent after ρ core reaches 10 4 is small. As a result, the calculation time to the moment of the core collapse of the P 3 T scheme is smaller than that of the Hermite scheme by a factor of two.

Orbital evolution of SMBH binary
In this section, we also discuss the performance of the P 3 T scheme applied to simulations of a galaxy with a supermassive black hole (SMBH) binary. In section 3.2.2.1, we describe the initial conditions and parameters of the integration scheme. In section 3.2.2.2 we compare the calculation results obtained by the P 3 T and Hermite schemes, and in section 3.2.2.3 the calculation speed.

Initial conditions and methods
We use the Plummer model with N =16K, 128K and 256K as the initial galaxy model. Two SMBH particles with a mass of 1 % of that of the galaxy are placed at the positions (± 0.5, 0.0, 0.0) with the velocities (0.0, ± 0.5, 0.0). We use three values for the cut off radius with respect to three different kinds of interactions. For the interaction between field stars (FSs), we set r cut,FS−FS = 1/256. For the interaction between SMBHs, the force is not split and F soft = 0. In other words, the force between SMBHs is integrated with the pure Hermite scheme. We set the cut off radius between SMBH and FS r cut,BH−FS = 1/32 which is large enough that ∆t soft is smaller than the Kepler time of a particle in orbit around the SMBH binary at a distance of r cut,BH−FS . We use the Plummer softening ǫ = 10 −4 for the interactions between FS-FS and FS-SMBH. For the SMBH-SMBH interaction, we do not use the softening. The accuracy parameter of timestep criterion for FS η FS is 0.1, and for SMBH η BH is 0.03. We adopt ∆r buff = 3σ∆t soft , ∆t max = ∆t soft /4 and θ = 0.4.
We use ∆t soft = 1/1024 at T = 0, and as the binary becomes harder, we decrease ∆t soft to suppress the aliasing error of the binary. As a standard model, we set ∆t soft to be less than half of the Kepler time of the SMBH binary t kep . Only for N =128K, we also perform two other runs, where ∆t soft < t kep /4 and t kep .
We also perform the same simulations by the Hermite scheme with the same η FS and η BH . Figure 14 shows the evolution of the semi-major axis (top panel) and eccentricity (middle panel) of the SMBH binary and the relative energy error (bottom panel) as functions of time for our standard models (∆t soft < t kep /2). The behaviors of the semi-major axis of the SMBH binary for the runs with the same N agree very well. The hardening rate of the binary depends on N because of the loss-cone refilling through the two-body relaxation [29][30][31]. The evolution of the eccentricity has large variation, because this evolution is sensitive to small N fluctuation [32]. In the cases of N =16K with the Hermite scheme, the relative energy error increases dramatically after T = 150 because the binding energy and the eccentricity of the binary are very high. Figure 15 is the same as figure 14 but for several different values of ∆t soft . Thick solid, dashed and dotted curves indicate the results for ∆t soft < t kep /4, t kep /2 and t kep , respectively. The orbital parameters show similar behaviors for all runs. The absolute value of the energy errors of P 3 T runs (∼ 10 −5 ) are small compared with the binding energy of SMBH binary, which is roughly 0.05. Figure 16 shows the calculation time for runs for several different values of N with ∆t soft < t kep /2. Initially, the P 3 T scheme is much faster than the Hermite scheme. As the SMBH binary becomes harder, the P 3 T scheme slows down more significantly than the direct Hermite scheme does. We can see that T cal of the Hermite scheme is roughly proportional to a −1 for a −1 > 300, whereas that of the P 3 T scheme is roughly proportional to a −5/2 , because ∆t soft is proportional to the Kepler time of the binary (∝ a 3/2 ). However, the calculation time for all runs with the P 3 T scheme is shorter than that with the Hermite scheme by a = 1/800. We can also confirm that as we use more N , the ratio of the calculation time of the P 3 T scheme to the Hermite scheme become larger. The reason why the P 3 T scheme becomes slower for large a −1 is simply that we force the timestep of all particles to be smaller than the orbital period of the SMBH binary. For the Hermite scheme, we do not put such constraint. Thus, in the Hermite scheme, particles far away from the SMBH have the timestep much larger than the orbital period of the SMBH binary. This large timestep can cause accuracy problem [33]. With P 3 T, it is possible to apply the perturbation approximation to F soft between the SMBH binary and other particles. Such a treatment should improve the accuracy and speed of the P 3 T scheme when SMBH binary becomes very hard.

Calculation speed
In figure 17, we plot the calculation time of the hard and soft parts for the standard model with N =128k. We can see that the soft parts dominate the calculation time.
In figure 18, we compare the calculation time for the runs with various ∆t soft ( < t kep , t kep /2, t kep /4). Since the most of the calculation time is spent after the binary becomes hard, the calculation time strongly depends on the criterion of the ∆t soft . From figure 15, the evolution of the orbital parameters for all runs with the P 3 T scheme are similar for various ∆t soft criterion. Thus we could chose larger ∆t soft t kep after the binary formation.

Conclusions
We described the implementation and performance of the P 3 T scheme for simulating dense stellar systems. In our implementation, the tree part is accelerated using GPU.
The accuracy and performance of the P 3 T scheme can be controlled through six parameters: ∆r cut , ∆r buff , ∆t soft , ∆t max , η and θ. We find that ∆r buff 2σ∆t soft is good choice to prevent non-neighbour particles from entering the neighbour sphere. The integration errors can be controlled through ∆t soft /∆r cut σ. For θ = 0.2, if we set ∆t soft to be less than 0.05∆r cut /σ, the integration error is smaller than the tree force error. For the Hermite part, if we chose η 0.2, the errors hardly depend on ∆t max .
From the point of view of the balance of the calculation costs between the tree and Hermite parts, we derive the optimal set of accuracy parameters, and found that the calculation cost is proportional to N 4/3 .
The P 3 T scheme is suitable for simulating large N stellar clusters with a high density contrast, such as star clusters or galactic nuclei. We demonstrate the efficiency of the code and show that it is able to integrate N -body systems to the moment of the core collapse. We also performed the simulations of the galaxy with the SMBH binary and found that the P 3 T scheme can be applied to these simulations.
Finally, we discuss the possibilities of implementation of two important effects on star cluster evolution to P 3 T. The first is an effect of a tidal field which dramatically change the collapse time and the evaporation time of a star cluster. The tidal field effect can be included in the soft part.
The other is an effect of the stellar-mass binary. A stellar-mass binary plays an important role in halting the core collapse. In this paper, we introduce the Plummer softening and neglect these binary effect. However, we could treat these effects by integrating stellar-mass binaries in the hard part.
Our P 3 T code is incorporated in the AMUSE frameworks and free for use [34,35].        In left and right panels, the x-axes are linear and logarithmic scales, respectively. Thin curves in right panels are proportional to T (solid) and T 1/2 (dashed).