On the number of particles in N-body simulations of planet–disk interaction

In this paper we present a new mixed-variables symplectic tree code for planetesimal dynamics, which allows to use a very large number of particles in numerical simulations of planet–disk interaction. Furthermore, we discuss how a wrong election of the number of particles to model a planetesimal disk, in order to carry out realistic numerical simulations, can cause spurious result. If the particle number is not enough, the close encounter between the planet and unrealistic massive particles will introduce a stochastic component of planet migration. This has an impact on the migration rate and it can even change its sign. Finally, we describe how to determine the minimum number of disk particles that should be used in numerical simulations, in order to have a correct description of the interaction, especially when the dynamical friction on the planet is the relevant mechanism of migration.


Introduction
The planet-planetesimal disk interaction has been extensively studied by means of numerical simulations. The attention of most papers was focused at the end of the runaway phase, with special emphasis in the structure of the disk in the neighborhood of the protoplanet Makino, 1992, 1993). In special Ida and Makino (1993) have shown that even for a small (non-migrating) protoplanet the ''planet-planetesimal-planetesimal'' interaction evolves in a ''planet dominated stage'' (planet-planetesimal interaction) in which the planet become a strong disk perturbator scatter the disk and slow-down the runaway growth. Cionco and Brunini (2002) have simulated the dynamical interaction between a massive planet and a non-selfgravitating Keplerian 3D disk with a relative large number of planetesimals. The planet can migrate freely within the disk and the disk can evolve without limits. The N-body simulations with large number of particles enable them to study the problem in a resonant perspective, together with other interactions as scattering, ejection, encounters and accretion. However, Cionco and Brunini (2002) were not able to extend their investigation to small mass planets, because of the possible artifacts that the number of particles could introduce in the numerical simulations.
Although in most papers containing N-body numerical simulations of planet disk interaction some attention is paid to the number of particles used, in none of them a detailed investigation of this effect is carried out. Therefore, we do not have a criterion helping us to decide the number of particles that a simulation of planet-disk interaction should include, in order to furnish reliable results.
Recently, we have developed a tree code for simulations of large collisional systems (Brunini and Viturro, 2003). A comparison between the behavior of the tree code and a direct hybrid integrator (the one used in Cionco and Brunini, 2002) revealed that the tree codes are useful in numerical simulations of planetary accretion, specially during intermediate stages, where dynamical friction dominates the evolution of planetary embryos embedded in disks, which is the case we are interested in.
In this paper, we first present an improved version of our tree code, where we take advantage of the fact that the bodies evolve mostly under a central potential. We therefore have implemented a hybrid symplectic algorithm, where the evolution under the central potential is computed with the analytical formulation of the two-body problem, but the perturbations are computed through an octal tree code and therefore is substantially faster than the previous version. After some numerical checks of the code, we perform a series of numerical simulations in order to determine the effect of the number of particles in N-body simulations of planet-disk interaction.

The numerical code
An efficient numerical integrator designed for the simulation of planetesimal dynamics in Keplerian disks must fulfill the following characteristics: The integrator should be symplectic in order to preserve the dynamical evolution of the system. It should take advantage of the fact that the bodies evolve mostly under a central potential.
Close encounters should be detected and integrated with enough accuracy. This a very important point for the comprehension of many process that affect the structure of the system.
The computation of the mutual gravitational interactions in direct form should be avoided. The reason of this is that the number of operations in its computation increases proportional to the square of the number of bodies.
At present, an integrator that incorporates all these desirable properties does not exist in the literature. Therefore, we have constructed a new integrator by combining features from different already known integrators. We used the symplectic integration scheme described by Duncan et al. (1998), to satisfies the first two points. In order to resolve the third point, we have incorporated the scheme described by Chambers (1999). Finally, to avoid the direct computation of the mutual interaction, we used a tree code, as it is described by Brunini and Viturro (2003). For the sake of completeness, we will describe briefly the proposed solutions in the next section.

Symplectic integrators
The N-body problem is a Hamiltonian problem, and as we know, in order to integrate a Hamiltonian system numerically over a long time span, the symplectic integrators present some advantages over the traditional numerical integrators: they preserves the phase-space structure of the system, and the accumulation of the truncation errors does not display a secular component in the energy of the system.
The Hamiltonian for ðN þ 1Þ bodies in mutual gravitational interaction is where q i ¼ ðx i ; y i ; z i Þ and p i ¼ ðp x i ; p y i ; p z i Þ are the position and momentum of body i (of mass m i ) relative to a Cartesian inertial frame. The simplest division of the Hamiltonian is to separate this in the sum of the kinetic and potential energy terms (2) In this case a second-order symplectic integrator is the wellknown leapfrog integrator used in many codes However, in a Keplerian disk, the central mass m 0 is bigger than the mass of all the other bodies m i ðia0Þ. Thereby the gravitational force due to the central body is usually the dominant force, and as a result it is better to write the Hamiltonian as where H Kep is the part that describes the central motion of the bodies around the central body, and H int is the part that describes the mutual interaction among the other bodies. Symplectic integrators based on this division of the Hamiltonian are known as mixed-variable symplectic integrators (MVS methods). The way of how this separation is done depends on the canonical variables chosen. The standard MVS methods use the Jacobian coordinates (Wisdom and Holman, 1992). A widely known second-order MVS integrator is the SWIFT code of Levison and Duncan (1994). However, it is more suitable (and computational easier) to use the so-called democratic heliocentric variables (Duncan et al., 1998). These variables consist of coordinates Q i relative to the central body and barycentric momenta P i . They are canonical variables too, and using them the Hamiltonian is split as where Here, a third term H arises in the Hamiltonian, due to the barycentric momentum of the central body. For this division of the Hamiltonian a second-order symplectic integrator can be easily found. Its action over the temporal evolution of the phase-space canonical coordinates u, for a time t is whereĤ ½ ; H is a differential linear operator defined by the Poisson bracket with the relevant part of the Hamiltonian. Examining the way like each Hamiltonian part performs, we can describe the integration algorithm as follows: (i) the coordinates remain fixed, and each body receives an impulse to its momentum owing to the other bodies (but not from the central body) during an interval t=2, (ii) the momenta remain fixed, and each body takes a shift in position by an amount ðt=2m 0 Þ P N i¼1 P i , (iii) the bodies move on an unperturbed Keplerian orbit for t, (iv) as step (ii), (v) as step (i).
This integrator scheme has two remarkable advantages regarding the classical leapfrog integrator. First, the Keplerian motion given by the action of H Kep has a wellknown analytic solution. Thereby, the advance of each body in a time step is performed with a minimum computational effort. Moreover, contrary to standard MSV methods, all bodies evolve around the same central body, with the same central mass.
Second, if all the bodies remain in well-separated orbits, it will happen jH int j5jH Kep j, jH j5jH Kep j by the ratio of the body's mass to the central mass, ¼ P N i¼1 m i =m 0 . Therefore, one step of the integrator has a truncation error Oðt 3 Þ instead of Oðt 3 Þ as a leapfrog integrator. This allows us, for the same precision, to use a longer time step by a factor À1=3 compared with a leapfrog method. This reduces the computational time.
However, this advantage remains valid only if close encounters between the ''small'' bodies do not take place. If two bodies have a close encounter, the mutual distance will become small and then the corresponding term in H int becomes comparable in size to H Kep . This increases the error per step from Oðt 3 Þ to OðtÞ. In order to maintain the accuracy the size of the time-step might be reduced. However, in a symplectic integrator a change of the timestep destroys its symplectic nature.

Resolving close encounters
It is not a trivial task to develop a symplectic integrator that accurately predicts and integrate close encounters. The standard MSV methods, using Jacobian coordinates, cannot handle arbitrarily close encounters. One solution is the one given by Duncan et al. (1998) in their SYMBA integrator. This code is a multiple time-step symplectic integrator which employs democratic heliocentric variables, and it splits up the perturbation term, giving to each part a separated fixed step size where weaker perturbations have bigger step sizes. This integrator, although truly symplectic, is troublesome to implement and to modify for our purpose.
Another way to maintain in the split Hamiltonian, H int smaller than H Kep , is by transferring the conflictive term from the first to the second part, during each close encounter. However, this naı¨ve exchange at each close encounter eventually destroys the symplectic nature of the integrator. In order to make the integrator truly symplectic, it never has to transfer the terms between the different parts of the Hamiltonian. The correct procedure is described by Chambers (1999). He rewrites the interaction terms in H Kep and H int introducing a changeover function K which form is such as it tends to one when the mutual distance jQ i À Q j j among the bodies i and j is large, while tending to zero when this distance is small: This form for K ensures H int remains smaller than H Kep , even during a close encounter. In the absence of encounters, H Kep can be integrated analytically as before. On another hand, when two bodies have a close encounter we have a three body problem (the central body and the bodies in close encounters). This situation is no longer integrable analytically, but in the practice, the equations can be integrated numerically with a traditional integrator (typically, we employ the Bulirsch-Stoer method). This contribution non-symplectic in the above scheme gives to the resulting integrator the category of hybrid.

A tree code
The direct evaluation of the inter-particle interaction is time-consuming due to the fact that it scales as N 2 . A way to speedup the computation is to employ a tree code because it scales as N logðNÞ.
An octal tree code specifically designed for planetesimal dynamics was implemented by Brunini and Viturro (2003). This code effectively predicts and integrates close encounters. However, the leapfrog algorithm is chosen to integrate the equations of motion and therefore it does not take advantage of the fact the bodies evolve mostly under a central potential. Our implementation is a modified version of the Brunini and Viturro code, which is based itself on the basic algorithm described by Barnes andHut (1986, 1989), Hernquist (1987), Barnes (1995) and Pfalzner and Gibbon (1996). A brief description of the tree code is as follows.
The first step in the tree code is the construction of the tree structure. This construction begins by creating a cube which contains the entire system of bodies and identifying this cube with the largest cell (called the root) of the tree. Then, this cell is subdivided into eight equal cubic subcells, and each one of these subcells are also recursively subdivided, until the ending cells (the leaves of tree) contain either one body or it is empty. With this spatial division of the system, the total masses, center-of-mass coordinates and quadrupoles moments are calculated for each cell to rove the tree using fast recursive formulae for the calculus. This provides a second-order approach to the gravitational potential of the system to be used in the computing of the interaction amongst the bodies.
The interaction system computing on a body starts traversing the tree from the root to the leaves. A cell will contribute as a whole to this interaction if it is far enough from the body. Otherwise, the process is repeated over the descendant subcells. The criterion to decide when a cell must be subdivide or not is related to the aperture angle defined as the ratio between the size of the cell and its distance from the body in consideration. In the simplest case a cell is subdivided when the aperture angle is greater than a certain threshold value (typically, in the range 0.5-0.7).

Detecting close encounters
As a bonus, the tree structure may be used to determine the particles in close encounters. As part of the process of integration, the code requires to know those particles in close encounters. A close encounter condition is defined as follows. For each particle i, we compute its Hill's radius, defined as where r i is the heliocentric distance of the particle and m i its mass. The particles in close encounter with i are the nearest-neighbors whose distance to i satisfies the condition where a is a proportionality factor. We adopted a ¼ 3 according to others authors (Levison and Duncan, 1994;Duncan et al., 1998;Chambers, 1999;Brunini and Melita, 2002). Hence, the problem requires an efficient method for the nearest-neighbor searching. A direct procedure, based on the computing of the mutual distance between pairs, is time-consuming. A better algorithm may be constructed using the tree structure. Our code implements the algorithm employed by Brunini and Viturro (2003), which is based on the same ideas described by Hernquist and Katz (1989). In order to search for neighbors at a distance rph from a body, it is put in a cube of size 2h (the ''search cube''). Starting with the root cell, the tree is traversed through in a similar manner as in the computation of the gravitational interactions. The question is if the ''search cube'' overlaps the cell that is being examined. If it does not overlap, then the branch of tree behind such cell will not contain nearest-neighbor and the search of it will be stopped. Otherwise, if the overlapping occurs, the cell will be divided into subcells and we will continue the search on the next hierarchy level of such branch. In the case that the cell that overlaps the ''search cube'' is a leave (this means, a body) it must be verified that such body really lie along the distance h. If this is the case, the body is added to the list of nearest-neighbors. This procedure is repeated over each tree branch.

Testing the code
We have performed a series of numerical experiments in order to test the new code. Essentially they are the same we have already presented in Brunini and Viturro (2003) and therefore we consider redundant to include them here. We have obtained excellent agreement with the previous experiments, even using a step size considerably larger. This is because of the hybrid nature of our symplectic code, which, in contrast to the previous version, exploit the fact that most of the motion is performed in the central, Keplerian potential. The details of these tests are available in the web site indicated at the end of this paper.

On the number of particles in planet-disk numerical simulations
A remarkable point in planet-disk numerical simulations is the number of particles used to describe the disk of planetesimals. As the planet's mass M of a simulations decreases, more particles are required to obtain reliable results for a given disk mass. Thus, if the disk does not contain enough particles, the planet-disk interaction will be dominated by the strong interaction at close encounter between the planet and unrealistic massive ''particles''.
We have designed a series of numerical simulations in order to determine which is the right number of particles for the same scenario described in Cionco and Brunini (2002). Our model consists in one planet of mass M immersed in a disk of N equal mass planetesimals, each one of mass m, orbiting around the Sun. According to Wahde et al. (1996), the self-interaction between the planetesimals was not considered. This assumption is also supported by the results of Ida and Makino (1993). They have shown that the mutual interaction between planetesimals is weaker than the effect due to the scattering by the protoplanet in the case M=m410 2 . In all our simulations we have M=mb10 2 . Planet-planetesimal interaction was fully accounted for. Collisions with the planet were treated as perfect accretion. All those planetesimal reaching a hyperbolic orbit were eliminated from the system. In our simulations there is no particle that hits the Sun, indeed an insignificant number (p2) reaches a heliocentric distance less than 1 AU. The simulations were carried out for five different planet masses (10, 8, 6, 4 and 2M È ), placed initially at 5 AU on a circular orbit. In all the simulations we have used the same planetesimal disk with constant surface density of 10 g cm À2 consistent with a model of minimum mass solar nebula ð0:01M Þ. In all our simulations, the initial orbital eccentricities and inclinations of the particles were generated at random with uniform probability distribution in the range ð025Þ Â 10 À3 . The inner disk limit was setted at 2 AU, and the outer limit at 12.5 AU, thus covering all the relevant resonances. The time step of the integrations was t ¼ 0:1 yr ðp5 Â 10 À2 of the smallest orbital period in each simulation). As a result of some pre-tests, this step size allows us a reasonable quicker computational time and also a suitable close encounter detection.
One of the most interesting results in simulations of planet-disk interaction, is the variation of the planetary semi-major axis, which determines the migration speed and the torque of the disk on the planet (Cionco and Brunini, 2002). As it was shown by Cionco and Brunini (2002), for the simulated disk, the planet-disk interaction should result in an inward migration of the planet, with a migration rate compatible with the predictions of the linear theory of density waves, as representative of the dynamical friction in disks. With the aim of checking the sensitivity of the results with respect to the number N of planetesimals, we have run simulations with the same disk mass, but with different number of particles in the range of 1000-250 000, for each planetary mass M. In total we have run 16 simulations for five planet masses.
When the ''granularity'', defined by the relation M=m, is bigger than a critical value, we expect that the effect of resonant dynamical friction among the particle disk and the planet, is completely superseded by the effect of the close encounters, which drive us toward a wrong model. When the number of particles is increased, the noise of granularity is reduced and the action of dynamical friction becomes relevant, obtaining an appropriate migration rate. An example of the this situation is shown in Fig. 1. This figure shows the temporal evolution of the semi-major axis for the case of a planet with 2M È initially placed at 5 AU in a simulated disk with three different initial number N of planetesimals. In the first run, N ¼ 50 000, as a result the mass of each particle is m ¼ 0:00327M È (and M=m ' 610). In the second run, N ¼ 150 000 and m ¼ 0:00164M È (M=m ' 1830). In the last run, N ¼ 250 000, m ¼ 0:00065M È (M=m ' 3000). The above results show clearly that the disk models of the first and second run are not reliable while the third run shows a migration rate compatible with the scenario of our model.
In order to verify that the wrong results are a consequence of a poor modeling of the dynamical friction, we use a statistical description of planetesimal orbits (Stewart and Wetherill, 1988), the so-called ''particle-ina-box'' approach. This model conceives the planetesimals in nearly co-planar, nearly circular concentric orbits. The mutual gravitational interactions perturb this orbits, drive to small non-negligible values of their eccentricities and inclinations. The statistical method ignores the individual orbits and instead of uses a probability density to describe the distribution of orbital elements in the planetesimal population. Assuming the orbital perihelia and longitudes of the node are randomly oriented, the eccentricity and inclination distribution of the population can be described by a single variable, the mean square random velocity, v, with respect to a circular zero-inclination reference orbit where v K is the local circular Keplerian velocity, e is the eccentricity and i is the inclination (measured in radians). Several mechanisms can affect the random velocity evolution of the population. A single equation that describes the evolution of a test body of mass M interacting with a population of bodies of mass m was given by Stewart and Wetherill (1988). In our gas-free disk of particles, this equation splits into four terms. The rate of change in V , the rms velocity of body of mass M, is given by where the four terms on the right-hand side are:
2. B viscous stirring caused by inelastic collisions, 3. C velocity damping due to energy dissipated by inelastic collision, 4. D energy transfer from large bodies to small ones via dynamical friction, The proportional factors can be found, for example, in Lissauer and Stewart (1993). The dynamical friction term (and the second term in A) tends to drive the system toward an equipartion of random kinetic energy, where the random velocities of the smaller bodies increase while the random velocities of the larger bodies decrease. The particle-in-a-box approach described above only works well locally. Therefore, we have carried out the computation of the four terms in the rate of change of V taken only the planetesimals within an annulus of width 3R H centered at the distance of the protoplanet. The results show the C term does not strongly affect the velocity evolution. Thus the A þ B and D terms lead the random velocity evolution. If the dynamical friction must be significant, the D term must take control over the A þ B term and we expect the jA þ Bj=jDj ratio stays lower than one time along. Fig. 2 shows the relative contribution of A þ B and D terms on the three runs and we can confirm that just only in the third run the action of dynamical friction is relevant. In fact, for the first and second run, the condition jA þ Bj4jDj is satisfied by 12% of times, while only by 4% in the third run. Table 1 and Fig. 3 present the planetesimals number needed to obtain reliable results in our simulations. As we mention above, for our particular disk conditions, inward planetary migration is expected and, therefore, the dynamical friction term should be mv 2 À MV 2 o0.
In a planetesimal disk, we can assume that there is a mass dependence of random velocities of the form (Safronov, 1969) v where the exponent q is negative, and q ¼ À 1 2 if energy equipartition is completely achieved (Wetherill and Stewart, 1989). Combining these two last equations, and using the fact that the disks of all our simulations all have the same mass, we found that the number of planetesimals should be where M 0 is a proportionality constant.We observe in Fig. 3, that, in our simulations, the necessary planetesimal  number to obtain suitable results scales as M À1:3 , which is compatible with the prediction above. As it is evident, the exact value of q is not relevant to determine the scale relation between the number of planetesimals in a given disk and the mass of the planet to get a reasonable characterization of the dynamical friction on the planet.

Conclusions
In this paper, we have presented a new improved version of a tree code for planetesimal dynamics, which is completely symplectic and exploit the fact that in a planetesimal disk, most of the motion is performed in a central potential. The performance of the code has allowed us to perform simulations with very large number of planetesimals, and to explore the importance of using many particles in planet-disk simulations. We have shown that the use of less planetesimals than a certain number could get unreliable results. These results are potentially dangerous, because, as it is shown in Fig. 1, the semi-major axis of the planet could experience a smooth migration in the wrong direction, even if the dynamical interaction with the disk is poorly modeled. Although our results have been obtained with a particular disk structure (constant surface density), the scale relation would be valid in more general cases. Therefore, a possible strategy to use in a simulation with a given disk structure, is to perform a series of test simulations with a massive planet of mass M 0 , where a small number N 0 of planetesimals is required (and therefore representing a relative fast numerical experiment). Once the minimum N 0 furnishing reliable migration direction and rate is obtained, we could determine the number of particles N that should be used in the simulation of interest, with the given planetary mass M, through the scaling relation The DAEDALUS integrator: The symplectic integrator described in this paper, named DAEDALUS, is publicly available in the Grupo de Ciencias Planetarias site of the Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, at the address: http:// gcp.fcaglp.unlp.edu.ar. Copies of the source code, instructions for how to compile and run the program, and example integrations can be obtained here.