EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 511 (February 2010) A&A, 511 (2010) A69 Full HTML
Free Access
Issue
A&A
Volume 511, February 2010
Article Number A69
Number of page(s) 9
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/200912870
Published online 12 March 2010
A&A 511, A69 (2010)

The validity of the super-particle approximation during planetesimal formation

H. Rein - G. Lesur - Z. M. Leinhardt

University of Cambridge, Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, UK

Received 11 July 2009 / Accepted 22 December 2009

Abstract
The formation mechanism of planetesimals in protoplanetary discs is hotly debated. Currently, the favoured model involves the accumulation of meter-sized objects within a turbulent disc, followed by a phase of gravitational instability. At best, one can simulate a few million particles numerically as opposed to the several trillion meter-sized particles expected in a real protoplanetary disc. Therefore, single particles are often used as super-particles to represent a distribution of many smaller particles. It is assumed that small-scale phenomena do not play a role and particle collisions are not modelled. The super-particle approximation is not always valid when applied to planetesimal formation because the system can be marginally collisional (of order one collision per particle per orbit). The super-particle approximation can only be valid in a collisionless or strongly collisional system, although, in many recent numerical simulations this is not the case.

In this work, we present new results from numerical simulations of planetesimal formation via gravitational instability. A scaled system is studied that does not require the use of super-particles. This system is simplified for computational practicality and proper identification of important processes: 1) the evolution of particles is studied in a local shearing box; 2) the particle-particle interactions such as gravity, physical collisions, and gas drag are solved directly assuming a constant background shear flow without any feedback from the particles. We find that the scaled particles can be used to model the initial phases of clumping if the properties of the scaled particles are chosen such that all important timescales in the system are equivalent to what is expected in a real protoplanetary disc. Constraints are given for the number of particles needed in order to achieve numerical convergence.

We compare this new method to the standard super-particle approach. We find that the super-particle approach produces unreliable results that depend on artifacts such as the gravitational softening in both the requirement for gravitational collapse and the resulting clump statistics. Our results show that short-range interactions (collisions) have to be modelled properly.

Key words: accretion, accretion disks - turbulence - methods: numerical - planets and satellites: formation

1 Overview

Extrasolar planets have been observed around a variety of parent stars from pulsars to solar-type stars to M-dwarfs (see e.g., Chauvin et al. 2004; Wolszczan & Frail 1992) indicating that planet formation is common and successful in a broad range of environments. However, the process of planet formation itself is not directly observable, leaving theory and numerical simulations to fill in the blanks between observations of hot circumstellar discs around young stars and planets orbiting mature stars. One of the most important unanswered questions in the theory of planet formation is what is the mechanism for planetesimal formation, i.e., the process by which the building blocks of planets are formed.

There are two main theories for planetesimal formation: mutual collisions (e.g., Hayashi et al. 1977) and gravitational instability in the dust layer (Goldreich & Ward 1973). In the first hypothesis, dust particles grow as the result of accretion-dominated collisions. Although the formation of planetesimals by mutual collisions is consistent with meteoritic evidence, the collision speed between dust particles (or aggregates) must be much slower than the typical velocity dispersion in a standard protostellar disc to avoid destructive collisions (Blum & Wurm 2008). In addition, the planetesimal formation process is so slow that meter-sized particles are in danger of spiraling into the star before growing large enough to decouple from the gas (Weidenschilling 1977b). Even if km-sized planetesimals were able to form, they would be in danger of being ground down again by mutual collisions (Stewart & Leinhardt 2009).

Gavitational instability is often considered to be a solution to most of these problems because the intermediate sizes are avoided all together. In this theory, the dust layer becomes dense enough for the Keplerian shear and velocity dispersion of the dust particles to be unable to support the dust against its own self gravity. The dust then collapses into clumps that eventually cool via drag forces and mutual collisions into planetesimals. Many different authors have worked on this subject. Until recently, the focus has been on quiet, non-turbulent, and low density regions of the accretion disc (see e.g., Tanga et al. 2004; Michikoshi et al. 2007,2009). However, the turbulent gas in the protoplanetary nebula stirs the dust, which increases the velocity dispersion of the dust particles. Several ideas have been proposed to overcome the turbulence-induced mixing of the dust particles and create localized clumps. For example, Cuzzi et al. (2008) suggest that the same turbulence that stirs the dust on larger scales may also collect the dust particles on small scales. A similar idea was proposed by Johansen et al. (2007), in which dust particles are localized into clumps using both turbulence and the streaming instability (Youdin & Goodman 2005). The clumps then become gravitationally unstable. A third hypothesis suggests that large structures, such as vortices, may be able to collect and protect dust particles from the turbulent background (Barge & Sommeria 1995). In this paper, we focus on the gravitational collapse in a very dense and turbulent region of the protoplanetary disc. The overdensity in the dust layer, which is approximately two orders of magnitude higher than the standard minimum-mass solar nebula, might have been created by any of the processes described above.

Numerical simulations must be used to test models of planetesimal formation. However, the separation of scales in the problem is huge. A meter-sized object is more than 12 orders of magnitude smaller than the protoplanetary disc thickness forcing numericists to use super-particles. These super-particles have a mass much higher than the mass of a single dust particle, but the forces (e.g., drag forces) acting on them are equivalent to those of a single dust particle. Tanga et al. (2004) were able to use this approximation to succesfully model the gravitational collapse in a regime in which collisions are unimportant. However, when the surface density is increased by two orders of magnitude, as mentioned above, the system does become marginally collisional and therefore the super-particle approach breaks down.

Michikoshi et al. (2009) perform a similar set of simulations, also in a low density, non-turbulent region of the disc. Although the authors describe their particles as super-particles, the way in which they accurately treat the collisions is the same as the method presented in this paper. Other approaches use Monte Carlo type methods to resolve statistically the outcome of overlapping particles (see e.g., Lithwick & Chiang 2007).

In the following, we look carefully at the numerical requirements of modelling gravitational instability accurately and test the validity of using super-particles in a high density region. For simplicity, we begin with a system without turbulence and assume that the surface density has already been increased by turbulence or some equivalent process. We find that the super-particle approach gives erratic results when collisions become important.

We therefore go on and present a method that does not use super-particles. We simulate a scaled system in which the number of particles is dramatically lower than in a real protoplanetary disc. To keep the surface density the same, the mass of individual particles has to be increased. The drag force produced by the surrounding gas is scaled accordingly, so that the stopping time remains constant. Furthermore, we increase the physical size of each particle. This allows us to keep the collision timescales in simulations of different particle numbers exactly the same. These requirements lead to results that are independent of the number of particles, which is the only free parameter in the simulation, and we call those simulations converged. The results agree with other simulations performed in the context of Saturn's rings (Goldreich & Tremaine 1982; Daisaka & Ida 1999). In this paper, we argue that this method should also be used in current simulations of planetesimal formation.

The outline of this paper is as follows. In Sect. 2, we define the important timescales in the problem and show that the super-particle approach breaks down in high density regions. In Sect. 3, we discuss the numerical method and the initial conditions of our simulations. Results of numerical simulations using both the super-particle approach and our scaling method are presented in Sect. 4. We conclude the paper with resolution constraints for numerical simulations and discuss the implications for the planetesimal formation process through gravitational instability in Sect. 5.

2 Orders of magnitude

2.1 Definition of timescales

The aim of this work is to simulate the dynamics of dust (or particles) interacting gravitationally inside a protoplanetary disc. Each particle is subject to five different physical processes, each operating on various typical timescales:

Stopping time $\tau_{\rm s}$

Each particle feels the effects of the surrounding gas through a linear drag force. This force can be written as $\vec{F}_{\rm drag}=\frac{m_{\rm p}}{\tau_{\rm s}} (\vec{v}_{\rm g}-\vec{v}_{\rm p})$, where $\tau_{\rm s}$ is the stopping time of a particle of mass $m_{\rm p}$, $\vec{v}_{\rm g}$ is the velocity of the gas, and $\vec{v}_{\rm p}$ is the velocity of the particle (Weidenschilling 1977a).

Physical collision timescale $\tau_{\rm c}$

Dust particles will suffer a physical collision with another dust particle on a timescale of $\tau_{\rm c}=(\sigma_{\rm c} {\bar v}_{\rm p} n)^{-1}$, where $\sigma_{\rm c}$ is the geometrical cross-section of the particles, ${\bar v}_{\rm p}$ is the particle velocity dispersion, and n is the number density of particles.

Orbital timescale $\tau_{\rm e}$

The timescale associated with the particles orbiting the central object in a local shearing patch is the epicyclic period $\tau_{\rm e}=2\pi\Omega^{-1}$, where $\Omega$ is the angular velocity at the semi-major axis of the shearing box.

In this paper, we consider two limits for the gravitational interactions between dust particles, long-range and short-range interactions. The long-range interaction can be seen as a collective process involving interactions between clouds of particles. On the other hand, the short-range interaction is important for resolving close approaches between pairs of particles (i.e., gravitational scattering). For the sake of clarity, we separate these two processes.

Gravitational collapse timescale $\tau_{\rm Gl}$

We introduce a length scale $\lambda \gg \delta_r$, where $\delta_r$ is the average distance between neighbouring dust particles. In this limit, which we call the long-range interaction, the dust particle distribution can be approximated by a continuous density $\rho$, and one can define a gravitational timescale, which is of the order of the free fall time of the system, defined by $\tau_{\rm Gl}=1/\sqrt{G \rho}$ where G is the gravitational constant.

Gravitational scattering timescale $\tau_{\rm Gs}$

Short-range gravitational interactions can be seen as an interaction between a pair of single particles that can result in a scattering event. As with physical collisions, the timescale for a gravitational interaction depends on the cross-section, $\tau_{\rm Gs}=(\sigma_{\rm G} {\bar v}_{\rm p} n)^{-1}$, where $\sigma_{\rm G}$ is the gravitational cross-section of the particle. Using Kepler's laws, one finds that $\sigma_{\rm G}\sim G^2 m_{\rm p}^2/{\bar v}_{\rm p}^4$. The gravitational cross-section is velocity dependent, which makes it qualitatively different from physical (billard ball) collisions between two particles. We note that gravitational focusing can also lower the physical collision timescale.

As mentioned earlier (Sect. 1), the particles might also be affected by excitation from a turbulent background state on a timescale defined by the turbelence itself. In this paper, we model this excitation in a simplified way, described in Sect. 3. In our case, the turbulence is scale independent and has no timescale associated with it.

2.2 The real physical system

To identify the dominant physical processes in a protoplanetary disc, one has to quantify and compare the relevant timescales as defined above. In the following discussion and in the numerical simulations, we assume R0=1 AU, a gas disc thickness of H/R0=0.01, and a minimum mass solar nebula (MMSN) with surface gas density $\Sigma=890~{\rm g~cm}^{-2}$. One usually assumes a solid to gas ratio of  $\rho_{\rm s}/\rho_{\rm g}=0.01$. We, however, assume a solid to gas ratio of unity. As mentioned above, this value is justified by numerical simulations (Johansen et al. 2009,2007) that demonstrate overdensities of the order of 101-102 can easily occur because of the interaction with a turbulent gas disc, the streaming instability, or vertical settling. We note that significantly different models of the solar nebula have been proposed (Desch 2007). However, all our results can easily be scaled to different scenarios.

We simplify the system by assuming that all solid components of the disc are in meter-sized particles. We assume that these boulders have a typical velocity dispersion caused by gas turbulence ${\bar v}_{\rm p}\sim0.05 c_{\rm s}\sim 30~{\rm m}/{\rm s}$, where $c_{\rm s}$ is the local sound speed, in accordance with numerical results (Johansen et al. 2007). Since ${\bar v}_{\rm p}\ll c_{\rm s}$, the particles tend to sediment toward the midplane, forming a finite thickness dust layer due to the non-zero velocity dispersion.

For meter-sized boulders, the physical cross-section is $\sigma_{\rm c}\sim 1~{\rm m}^2$, whereas the gravitational scattering cross-section is $\sigma_{\rm G}\sim10^{-21}~{\rm m}^2$ showing clearly that gravitational scattering is irrelevant to the dynamics of dust particles embedded in a disc. However, all other timescales are roughly equivalent. Meter-sized boulders are weakly coupled to the gas with a stopping time $\tau_{\rm s}\sim \tau_{\rm e}\sim \Omega^{-1}$ (Weidenschilling 1977a). The physical collision timescale is $\tau_{\rm c}=7.3 ~ \Omega^{-1}$ and the long-range gravitational interaction timescales is $\tau_{\rm Gl} \sim \Omega^{-1}$.

This physical system is expected to become gravitationally unstable, according to the Toomre criterion (Toomre 1964). The instability occurs when the gravitational collapse timescale $\tau_{\rm Gl}$ is shorter than the transit time due to random particle motions $\lambda /{\bar v}_{\rm p} $ and the orbital timescale. Assuming that the particles can be modelled by an isotropic gas[*] with a sound speed ${\bar v}_{\rm p}$, the system becomes unstable when

\begin{displaymath}Q\equiv \frac{{\bar v}_{\rm p} \Omega}{\pi G \Sigma}<Q_{\rm crit}\simeq 1. \end{displaymath} (1)

In that case, the most unstable wavelength is given by

\begin{displaymath}\lambda_{\rm T}=\frac{2\pi^2G\Sigma}{\Omega^2}\cdot \end{displaymath} (2)

In the following, we compare the physical parameters from this section to their counterpart in numerical simulations. We first summarize the super-particle approach before presenting the scaling method.

2.3 Super-particle approximation

A super-particle represents many smaller particles. The dynamics of the small particles are not calculated exactly. It is assumed that all particles behave similarly and in a collective manner. To simulate gravitational interaction between super-particles, one has to use a softened potential. Without this, the gravitational scattering cross-section of super-particles becomes too large and super-particles undergo gravitational scattering events, which is unphysical since these events never occur in a real system. Individual particle-particle collisions are not modelled.

There are various examples where this approach is used successfully. For example, smoothed particle hydrodynamics (SPH) uses the super-particle approximation to simulate many gas molecules (Lucy 1977). These systems are often assumed to be strongly collisional to ensure thermodynamical equilibrium inside each super-particle. One can therefore assign collective properties to clouds of particles such as pressure and temperature. Another example is the evolution of galaxies. When two galaxies collide, individual particles (stars) will usually not undergo gravitational scattering events or physical collisions. In that case, the super-particle approach models a collisionless system in which collective dynamics is the only important physical process. In both cases, the super-particle approximation is a valid approach for simulating the system numerically, but will break down as soon as the system is marginally collisional.

As an example, we consider two clouds of particles undergoing a ``collision''. In the strongly collisional case, the clouds will slowly merge and the thermodynamic variables (e.g., temperature) will diffuse between the clouds, the particles inside each cloud following a random walk trajectory due to numerous collisions. On the other hand, in the collionless regime, the two clouds simply do not see each other because there is no short-range interaction present between the particles. In a marginally collisional regime, some particles will collide with particles of the other cloud, leading to a partial thermalisation of the velocity distribution, but some other particles will not have any collision at all and will follow approximately a straight line. Evidently, the outcome of this event cannot be described using a super-particle approach.

2.4 Scaling method

The idea of our scaling method is that one should keep all important timescales in a numerical simulation as close to those of the real physical system as possible and model all particle collisions explicitly.

The numerical system consists of $N_{\rm num}$ particles in a box with length H, simulating a patch of the disc at a fixed radius R0. The particle mass is $m_{\rm num}=\Sigma H^2/N_{\rm num}$ (H is the box size and one scale height). Thus, the density in the box and the long-range gravitational interactions are unchanged compared to the real system.

The gravitational scattering cross-section of the numerical system is then given by

\begin{displaymath}\sigma_{\rm G,~{\rm num}}=\frac{G^2 \Sigma^2 H^4}{N_{\rm num}^2 {\bar v}_{\rm p}^4}\cdot \end{displaymath} (3)

For the initial surface density and velocity dispersion used in the simulation (Sect. 2.2), one finds

\begin{displaymath}\sigma_{\rm G,~{\rm num}}\simeq \frac{ H^2}{N_{\rm num}^2} \simeq \frac{0.0001{\rm AU}^2}{N_{{\rm num}}^2}\cdot \end{displaymath} (4)

The physical collision cross-section in the simulation is $\sigma_{\rm c,~{\rm num}}~=~\pi a_{\rm num}^2$ where $a_{\rm num}$, is the radius of the particles in the simulation. We derive two constraints from the physical and gravitational collision cross-sections:
1.
Same mean free path in simulation and real physical systems

That means $N_{{\rm num}} \sigma_{\rm c,~{\rm num}} =N \sigma_{\rm c} $, where $\sigma_{\rm c}$ and N are the geometrical cross-section and the number of particles in a region of size H3 in a real disc, respectively. In other words, this condition ensures that the physical collision timescale in the simulation is exactly the same as in the real system.

2.
Negligible gravitational scattering

When two particles approach each other, the outcome should be a physical collision, which means that $\sigma_{\rm c,~{\rm num}}~\gg~\sigma_{\rm G,~{\rm num}}$. The gravitational scattering timescale remains long compared to the physical collision timescale.

The first condition places a constraint on the particle size in the simulation. With $N=4.7\times 10^{18}$ and a particle radius of $a=1~{\rm m}$ as found in a real disc within a box of volume H3, using the parameters from above finds that

\begin{displaymath}a_{\rm num} = \sqrt{\frac{\sigma_{\rm c,~{\rm num}}}\pi} = \s... ...rm num}}\pi}} \simeq \frac{0.014}{\sqrt{N_{\rm num}}}{\rm AU}. \end{displaymath} (5)

We note that once this condition is satisfied in the initial conditions it will be automatically satisfied at all times. In simulations, we typically find that the collision timescale is reduced by more than one order of magnitude during the gravitational collapse.

The second condition is then satisfied by changing the number of particles. An interesting result is that if $N_{\rm num}>10$, the second condition is easily satisfied if the first one is satisfied. We note however that $\sigma_{\rm G,{\rm num}}$depends strongly on the velocity dispersion. In particular, a velocity dispersion 5 times smaller than the initial value (as found in some simulations) leads to an increase by a factor of 600 in the gravitational scattering cross-section. Therefore, we suggest using a large safety factor ( $N_{\rm num}>10^5$) to ensure that the second condition is always satisfied, even for significantly smaller  ${\bar v}_{\rm p}$. This condition also allows us to estimate when our approach breaks down, namely when the number of clumps in the system is so low ( $N<10\sim100$) that gravitational scattering becomes important.

Although it would be helpful as a further simplification, it is not possible to perfom the above mentioned simulation in two dimensions. In that case, the filling factor, which is defined to be the ratio of the volume (or area) of all particles to the total volume (or area), is of the order of one for the above parameters. We note that in 2D an increase in particle number (and decrease in particle size as required by Eq. (5)) does not decrease the filling factor if the collisional lifetime remains constant.

3 Methods

Two different kinds of simulation are considered in this paper:

1.
Particles are assumed to be point masses and have no physical size, the gravitational field of the particles being approximated with a smoothing length to avoid numerical divergences (see Sect. 3.1).

2.
Particles have a physical size and, therefore, no gravitational smoothing is required but physical collisions must be included.
We refer to the particles as super-particles and scaled particles, respectively.

We perform our simulations in a cubic box with shear periodic boundary conditions in the radial (x) and perdiodic boundary conditions in the azimuthal (y) and vertical (z) directions, as illustrated in Fig. 1 (Wisdom & Tremaine 1988). In the local approximation, the force per mass on each particle is a sum of the contributions from Hill's equations and the interaction terms. Hill's equation can be written as

\begin{displaymath}\bm{F}_{\rm hill} = -2\Omega \;\vec{e}_z\times \vec{v}_{\rm p}+3 \Omega^2 \;x\;\vec{e}_x-\Omega^2z\;\vec{e}_z. \end{displaymath} (6)

The interaction term $\bm{F}_{\rm int}$ is divided into components related to self gravity, physical collisions between particles, or drag and excitation forces:

\begin{displaymath}\bm{F}_{\rm int} = \bm{F}_{\rm grav} + \bm{F}_{\rm col} + \bm{F}_{\rm drag} + \bm{F}_{\rm turb}. \end{displaymath} (7)

We solve the resulting equations of motion with a leap frog (kick drift kick) time stepping scheme. In the following subsections, we describe the numerical methods used to compute each of these terms and their physical relevance.

The box size of $0.01~{\rm AU}$was chosen such that there are always several unstable modes in the box, when the system is pushed into the regime of gravitational collapse, as estimated by Eq. (2).

\begin{figure} \par\includegraphics[width=8.9cm,clip]{12870f1.eps} \end{figure} Figure 1:

Shearing box, simulating a small patch of the protoplanetary disc. The grey box represents the shearing box of interest, the dashed boxes surrounding the central box represent one ghost ring.

Open with DEXTER

3.1 Self gravity


In the N-body problem that we consider, we have to solve Newton's equations of universal gravitation for a large number of particles N. The gravitational force on the i-th particle is given by

\begin{displaymath}\bm{F}_{\rm grav} = \sum_{j\neq i} G \frac{m_im_j}{\left( r_{ij} +b \right)^2}\vec{e}_{ij}, \end{displaymath} (8)

where b is the smoothing length used to avoid divergences in numerical simulations (in simulations that include physical collisions b = 0) and $\vec{e}_{ij}$ is the unit vector in the direction of the gravitational force between the i-th and j-th particle. Calculating the gravity for each particle from each other particle results in $O\left(N^2\right)$ operations. To reduce the number of operations, one can use different approximations to Eq. (8). In this paper, we use a Barnes-Hut (BH) tree code (Barnes & Hut 1986).

The BH tree in three dimensions divides the original box into eight smaller cells with half the length of the original box. This process is continued recursively until there is only one particle per cell left. The depth of the tree in a homogeneous medium is approximately $\log_8 N$. One then calculates the total mass and the centre of mass for every cell at every level of the tree.

To calculate the force acting on a particle, one starts at the top of the tree and descends into the tree as far as necessary to achieve a given accuracy. If the current cell is far away from the particle for which the force is calculated, then the detailed density structure within this cell is not important. All that matters is the box's monopole moment (total mass and centre of mass). One therefore does not have to descend into this tree branch any further. The BH tree reduces the number of calculations to $O\left(N\log N\right)$.

We use 8 rings of ghost boxes in the radial and azimuthal direction (Fig. 1 shows one ring). A ghost box is simply a shifted copy of the main box. The gravity on each particle is then calculated by summing over contributions from each (ghost) box. This setup approximates a medium of infinite horizontal extent and avoids large force discontinuities at the boundaries. We do not need any ghost boxes in the vertical direction because the disc is stratified.

A finite number of ghost rings can only act as an approximation of an infinite medium and the gravitational force will tend to concentrate particles horizontally in the centre of the box. Due to the 8 ghost rings, the asymmetry of the gravitational force between the centre and the faces of the box is reduced by a factor of 20 compared to using no ghost rings at all. We note, however, that a small asymmetry remains in our simulations and leads to higher concentrations in the box centre. In other simulations such as Tanga et al. (2004), this effect is not seen because the system is initially gravitationally unstable and integrated only for approximately one orbit. The system that we are interested in is marginal gravitationally unstable and this slight asymmetry will become important after several orbits. In the beginning of our simulations, we integrate for many orbits in order to reach a stable equilibrium. We then push the system into a gravitationally unstable regime (see Sect. 3.5 for details). During the stable phase, horizontal over-concentrations are to be expected and are indeed observed because of this asymmetry. However, since in a real protoplanetary disc the gravitational instability is considered to appear locally (e.g., inside a vortex or a similar structure), this effect of having a preferred location for gravitational instability is not unphysical and does not affect any conclusions made in this paper. We tested this by applying a linear cut-off to the gravitational force at a distance of one box length (A. Toomre, private communication). In this case, there was no preferred location in the box and the simulation evolved in exactly the same way.

3.2 Physical collisions

Physical collisions between particles are treated in the following way. After each timestep, we check whether any two particles are overlapping. Using the already existing tree structure from the gravity calculation, one can again reduce the computational costs from $O\left(N^2\right)$ to $O\left(N\log N\right)$. In these simulations, a collision between particles is defined as an overlap between two particles that are approaching each other. When a collision is detected, it is resolved assuming energy and momentum conservation (perfectly elastic collisions). We note that the timestep has to be small enough such that no collision is missed and the overlap is always small compared to the particle size.

The particle radius is given by Eq. (5). In order to keep the collision timescale $\tau_{\rm c}$ close to unity, which is numerically the worst case scenario because all timescales are equivalent, we increase the radius by a factor of 4 in all simulations.

3.3 Drag force

Each particle feels a drag force. The background velocity of the gas $\vec{v}_{\rm g}$ is assumed to be a steady Keplerian profile

\begin{displaymath}\vec{v}_{\rm g}=-\frac32 \Omega x\;\vec{e}_y. \end{displaymath} (9)

Although accretion flows are turbulent, we choose to use a simple velocity profile so that the behaviour of self-gravitating particles can be understood in conditions that can be easily controlled.

3.4 Random excitation

The system of particles described above is dynamically unstable even without self-gravity, as the coefficient of restitution and the stopping time do not depend on the particle velocities[*]. The particles can either settle down in the midplane and create a razor-thin disc if the stopping time is too short, or they can expand vertically forever if the stopping time is above some critical value and the excitation mechanism is provided only by collisions. Dust particles in accretion discs are found in the settling regime. However, a complete settling never occurs in accretion disks because the background flow is always turbulent due to the Kelvin-Helmoltz instability (Johansen et al. 2006), the MRI, or other hydrodynamic instabilities (see e.g., Lesur & Papaloizou 2010) which diffuse particles vertically (Fromang & Papaloizou 2006). This stirring process of turbulence is not present in the gas velocity field of the simulations presented here (see Eq. (9)). To approximate the turbulent mixing, we added a random excitation (white noise in space and time) to the particles in the simulation. This allows us to have a well defined equilibrium in which the system is stable rather than starting from unstable initial conditions that might influence the final state.

We perturb the velocity components of each particle after each timestep on a scale $\Delta {v}_i =\sqrt{\delta t}~ \xi$, where $\delta t$ is the current timestep and $\xi$ is a random variable with a normal distribution around 0 and variance s. This excitation mechanism heats up the particles and allows us to have a well defined three dimensional equilibrium as shown in Sect. 4. In our simulations, we use a value of $s=1.3\times 10^{15}\;{\rm m}^2~{\rm s}^{-3}$ and $s=8.9\times 10^{14}\;{\rm m}^2~{\rm s}^{-3}$for the simulations without and with collisions, respectively. These values were chosen such that the equilibrium state is approximately equivalent for both types of simulations.

3.5 Initial conditions

All particles, which have equal mass, smoothing length, and physical size, are placed randomly inside the box in the x-y plane. We note that particles have either a physical size or a smoothing length associated with them, depending on the type of simulation we perform (Sect. 3). In the z-direction, the particles are placed in a layer with an initially Gaussian distribution about the midplane and a standard deviation of 0.05 H. We allow the system to reach equilibrium by integrating it until $t=30\Omega^{-1}$ (4.8 yrs).

We tested various ways of pushing the system into the unstable regime, either by reducing the turbulence stirring or by shortening the stopping time, but did not see any qualitative difference as long as we start from a well defined equilibrium state. In the simulations presented in this paper, we switch off the turbulence stirring. Following this modification, the system becomes gravitationally unstable and bound clumps form within a few orbits.

4 Results

4.1 Super particles

\begin{figure} \par\includegraphics[angle=270,width=8.8cm,clip]{12870f2.eps} \vspace*{-2mm} \end{figure} Figure 2:

Super-particles: velocity dispersion as a function of time in four different runs. The simulation labeled LS has a ten times larger smoothing length. The curves do not overlap because the simulations are not converged.

Open with DEXTER

We first present simulations without physical collisions (super-particles) which rely on the smoothing length b to avoid any divergencies. Although we use a tree code and therefore a smoothed potential for the force calculation is assumed, the results are equivalent to an FFT-based method where the grid length acts as an effective smoothing length. In general, to check the numerical resolution, the particle number $N_{\rm num}$ is increased and the smoothing length b is reduced independently. A simulation is resolved when the result is independent of both $N_{\rm num}$ and b. This turns out to be impossible in the present situation. The main reason is that the smallest scale in a gravitational collapse will ultimately depend on the smoothing length. By varying both parameters at the same time, an empirical scaling of  $b\sim1/\sqrt{N_{\rm num}}$ works fairly well if $N_{\rm num}$is large enough. However, this procedure is not justified and is unphysical. The smoothing length was introduced to avoid divergent terms and not to model any small-scale physical process. Therefore, it should not have any impact on the physical outcome of the simulation. Incidentally, the existence of this dependency indicates that short-range interactions are important for the result of these simulations and should be modelled with care.

\begin{figure} \par\includegraphics[angle=270,width=8.8cm,clip]{12870f3.eps} \vspace*{-2mm} \end{figure} Figure 3:

Scaled-particles: velocity dispersion as a function of time in three different runs that include physical collisions. All curves overlap because the simulations are converged.

Open with DEXTER

In Fig. 2, we plot the velocity dispersion as a function of time. The simulations begin from the stable equilibrium described in Sect. 3. After $t=30~\Omega^{-1}$, we switch off the excitation mechanism. Because the system continues to cool, it becomes gravitationally unstable within one orbit, as estimated by Eq. (1). Once clumping occurs, the velocity dispersion begins to rise again. All simulations use the particle radius given by Eq. (5) as a smoothing length except the ones labeled LS, which use a ten times larger value. The larger softening length is approximately a 128-th of the box length and illustrates the kind of evolution expected from an FFT method using a 1283 grid.

\begin{figure} \par\includegraphics[width=16.3cm,clip]{12870f4.ps} \end{figure} Figure 4:

Super-particles: snapshots of the particle distribution in the xy plane. Both simulations use 160 000 particles with different smoothing lengths. The simulation on the bottom uses a ten times larger smoothing length than the one on the top. The snapshots were taken ( from left to right) at $t=0, 30, 37, 40~\Omega^{-1}$. With a large smoothing length, the outcome looks very different, the system is more stable, more stripy structure can be seen and clumps form later, if at all.

Open with DEXTER
\begin{figure} \includegraphics[width=16.3cm,clip]{12870f5.ps} \end{figure} Figure 5:

Scaled particles: snapshots of the particle distribution in the xy plane. The simulations ( from top to bottom) use 40 000, 320 000, 640 000 particles with their physical size given by Eq. (5). In all simulations, the collision timescale was kept constant. The snapshots were taken ( from left to right) at $t=0, 30, 37, 40~\Omega^{-1}$. The simulation with 40 000 particles has a large filling factor by the last frame, which prevents clumps from forming. The intermediate resolution simulation ( middle row) and the highest resolution simulation ( bottom row) have more and smaller particles and thus a smaller filling factor. The results (i.e., number of clumps in the last frame) are very similar in the intermediate and high resolution simulations.

Open with DEXTER

Snapshots of the particle distribution of two simulations are shown in Fig. 4. Both simulations use 160 000 particles and all parameters, except the smoothing length, are the same. The top row is the simulation with a smoothing length given by Eq. (5), whereas the bottom row uses a smoothing length that is ten times larger. The simulations correspond to the blue (medium dashed) and red (solid) curve in Fig. 2. We show the snapshots to illustrate the importance of the smoothing length in simulations without physical collisions. One can see that the simulations differ already before clumps form. Stripy structures appear on a scale given by Eq. (2). As soon as we enter the unstable regime at  $t\sim32\Omega^{-1}$ (i.e., the velocity dispersion begins to rise, see also Fig. 2) the simulations evolve very differently. The simulation on the top row of Fig. 4 forms many clumps at an early time, whereas the bottom row simulation forms only a few, more massive clumps at later times.

This can be confirmed by looking at spectra of the same two simulations as shown in Fig. 6. These spectra were generated by mapping the particles onto a $128\times 128$ grid in the x-y plane and computing the fast Fourier transform of the mapping in the x direction. The resulting spectra were finally averaged in the y direction (see Tanga et al. 2004, for a complete description of the procedure). As suggested by the snapshots, the nonlinear dynamics of the gravitational instability strongly depends on the smoothing length used. In particular, one observes that smaller scales are amplified more slowly for larger smoothing lengths. The resulting spectra at $t=40\Omega^{-1}$ also differ significantly. With a small enough smoothing length, the spectrum looks almost flat, whereas a large smoothing length introduces a cutoff at $k/2\pi\simeq 10$. We also note that the smoothing length in the latter case is of the order of the grid size ( $k/2\pi\sim100$). The cutoff observed clearly demonstrates that the smoothing length modifies the dynamics on scales up to 10 times larger than the smoothing length itself.

4.2 Scaled particles

The velocity dispersion evolution of simulations with N=160 000, 320 000, 640 000 particles including physical collisions is presented in Fig. 3. In these runs, no gravitational smoothing length is needed. Again, the simulations start from a stable equilibrium and after $t=30~\Omega^{-1}$, we change the same parameters as in the non-collisional case to push the system into the gravitationally unstable regime. Snapshots of the particle distributions for three runs are also plotted in Fig. 5. The middle and bottom rows of snapshots correspond to the purple (medium dashed) and dark blue (long dashed) lines of Fig. 3.

One can see in Fig. 3 as well as Fig. 5 that the intermediate and high resolution runs produce very similar results. We call those simulations converged, as a change in particle number does not change the outcome. An additional, lower resolution run (N=40 000) is not converged because the filling factor is too large in dense regions, preventing clumps from being gravitationally bound. The problem does not occur in the non-collisional cases because particles are allowed to overlap. Since the filling factor scales with particle number as  $1/\sqrt{N_{\rm num}}$, this issue is resolved for runs with more than a few hundred thousand particles. We note that the velocity dispersion might differ slightly at later stages for the converged runs because the final clump radius is still determined by the physical particle radius and therefore by the particle number.

\begin{figure} \par\includegraphics[angle=270,width=9cm,clip]{12870f6} \end{figure} Figure 6:

Super-particles: density distribution spectrum with 160 000 particles and two different smoothing lengths that differ by a factor 10 at times t=30 (red, solid curve), t=36 (green, long dashed curve), and t=38 (blue, short dashed curved). The spectra at similar times are not on top of each other because the simulations are not converged.

Open with DEXTER

\begin{figure} \par\includegraphics[angle=270,width=9cm,clip]{12870f7} \end{figure} Figure 7:

Scaled particles: density distribution spectrum with 160 000, 320 000, and 640 000 particles (colors are equivalent to Fig. 6). At each time, the spectra are converged, i.e., the spectra are on top of each other and independent of the particle number, besides the noise level on very small scales. At late times ( $t\gtrsim 38$, blue curve), one begins to see more structure on small scales in simulations with larger number of particles (see text).

Open with DEXTER

We show the density distribution spectrum in Fig. 7 as a function of the number of particles. The resulting spectra do not depend on the number of particles in the system, as expected. In comparison with Fig. 6, the simulations with collisions are converged since the dynamic is not controlled by the numerical parameter N, which is the only free parameter in the system.

At very late times ( $t\gtrsim 38\Omega^{-1}$), the spectra begin to show a systematic trend towards more structure on smaller scales for runs with larger N. This is expected and due to the filling factor, already mentioned above. As soon as one clump in the simulation is only determined by the size of the individual particles it consists out of, our approach breaks down. One can also see these clumps forming by looking at the velocity dispersion in Fig. 3. The particles inside a clump begin to dominate the velocity dispersion over the background after $t\sim38\Omega^{-1}$. A clear indication of this is the spiky structure with a typical correlation time of ${\sim}\Omega^{-1}$ corresponding to different clumps interacting and merging with each other. One way around this issue, which will be considered in future work, is to allow particles to merge (Michikoshi et al. 2009). Using this accretion model, the mass of the clump can be used as an upper limit.

5 Discussion and implications

In this paper, we have shown that convergence in N-body simulations of planetesimal formation via gravitational instability can be achieved when taking into account all relevant physical processes. It is absolutely vital to simulate gravity, damping, excitation, and physical collisions simultaneously, as the related timescales are of the same order and therefore all effects are strongly coupled.

A set of simulations is defined to be converged when the results do not depend on the particle number or any other artificial numerical parameter such as a smoothing length. As a test case, a box with shear periodic boundary conditions was used, containing hundreds of thousands of self-gravitating (super-)particles in a stratified equilibrium state. The particles were then pushed into a self gravitating regime that eventually led to gravitational collapse. Simulations with and without physical collisions were studied for a range of particle numbers to test convergence.

In cases without physical collisions, convergence can not be achieved. However, it was possible to change multiple free simulation parameters at the same time, namely the particle number N and the smoothing length b, such that in special cases the results did not depend on the particle number. We note however that there is a free parameter in the simulation (i.e., the smoothing length) that effects the outcome and makes it impossible to find the real physical solution.

In the protoplanetary disc, physical collisions dominate over gravitational scattering. In the case of planetesimal formation, the system is marginally collisional. Physical collisions will partially randomise the particle distribution, whereas a smoothing length does not because the softening length is very large compared to the gravitational cross-section. Even if the gravitational particle-particle scattering becomes important, it can be seen from Eq. (3) that the velocity dependence of the cross-section $\sigma_{\rm G}$ differs fundamentally from the velocity independent cross-section $\sigma_{\rm C}$. We therefore argue that the super-particle approach should not be used when collisions become important.

With collisions, the simulation outcome is both quantitatively and qualitatively very different from simulations with no physical collisions. The initial size of the clumps is larger and the number of clumps is smaller. As there is no smoothing length, there is effectively one fewer free parameter. Changing the particle number while keeping the collision time constant does not change the outcome. We therefore call these simulations converged.

Additional tests including inelastic collisions with a normal coefficient of restitution of 0.25 have been performed to confirm that the results do not depend on the way physical collisions are modelled. As expected, the qualitative outcome in terms of number of clumps and clump size is different because the physical properties of the system have been changed. However, once individual collisions are resolved the results converge in exactly the same manner as in those simulations presented above (which have a coefficient of restitution of 1).

We have also explored other gravity solvers by comparing a fast Fourier (FFT) gravity solver with the BH tree code (using a smoothing length) and confirmed our previous results. The grid in the FFT code acts as a smoothing length and a direct comparison between the two gravity solvers gave an approximate effective smoothing length of a quarter of the grid length. Because one has to solve individual particle-particle interactions, the grid size has to be smaller than the physical size of the particles. As a result, the FFT method is always slower than a tree code. In addition, the tree structure can be reused to search efficiently for collisions.

This paper focused on the numerical requirements to study gravitational instability and the formation of planetesimals in protoplanetary discs. The initial conditions were chosen such that the equilibrium is a well defined starting point for the convergence study. The collision, damping, collapse and orbital timescales are all of the same order, which is effectively the worst case scenario. To provide any constraints on planetesimal formation itself, one would have to properly simulate the turbulent background gas dynamics which is beyond the scope of this paper.

However, we were able to determine the numerical resolution needed to resolve dust particles in a protoplanetary disc. Assuming that one wishes to simulate dust particles realistically, including collisions, and that the particles are uniformly distributed in a box of base L2 and height H, the size of each particle is determined by the collision rate in the real disc (see Eq. (5))

\begin{displaymath}a_{{\rm num}} = \sqrt{\frac{N}{N_{{\rm num}}}} a, \end{displaymath} (10)

where a and N are the size and number of particles in the volume L2H of the real disc. The number of particles in the simulation $N_{\rm num}$ has to be large enough to ensure that the filling factor is low until clumps form. The requirement that the filling factor is less than unity at t=0 is given by

\begin{displaymath}\frac43 \pi \left(a_{{\rm num}}\right)^{3} \cdot N_{{\rm num}} \le C_{\rm f}\; L^2 H, \end{displaymath} (11)

where $C_{\rm f}$ is a safety factor. Although $C_{\rm f}\la 1$ would be enough to resolve collisions initially, it is insufficient to simulate the collapse phase. During the collapse, the filling factor rises rapidly. Assuming that one wishes to resolve collisions correctly when the particles have contracted by one order of magnitude, one has to include a safety factor of $C_{\rm f}=10^{-3}$. If the collapse occurs mainly in two dimensions, as in the simulations presented here, a factor of $C_{\rm f}=10^{-2}$ is sufficient. Furthermore, one has to ensure that the box contains at least a few unstable modes (see Eq. (2)). All this together places tight numerical constraints on numerical simulations of planetesimal formation via gravitational instability.

In future work, we will expand the discussion to include more physics, namely a proper treatment and feedback of the background gas turbulence. This will allow us to further clarify the behaviour of planetesimals in the early stages of planet formation and eventually test the different formation scenarios.

Acknowledgements
We thank J. Papaloizou, G. Ogilvie, S.-I. Inutsuka, A. Johansen, Y. Lithwick, and A. Youdin for helpful discussions and comments. Hanno Rein is supported by an Isaac Newton Studentship and St John's College Cambridge. Geoffroy Lesur acknowledges support from STFC. Zoë M. Leinhardt is an STFC postdoctoral fellow. Simulations were performed on the Astrophysical Fluids computational facilities and on Darwin, the Cambridge University HPC facility. The authors are grateful to the Isaac Newton Institute for Mathematical Sciences in Cambridge where the final stages of this work were carried out during the Dynamics of Discs and Planets research programme. A simplified version of the simulation used in this paper can be downloaded free of charge as an application for the iPhone at http://itunes.com/apps/gravitytree.

References

  1. Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] [Google Scholar]
  2. Barnes, J., & Hut, P. 1986, Nature, 324, 446 [NASA ADS] [CrossRef] [Google Scholar]
  3. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Chauvin, G., Lagrange, A.-M., Dumas, C., et al. 2004, A&A, 425, L29 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
  5. Cuzzi, J. N., Hogan, R. C., & Shariff, K. 2008, ApJ, 687, 1432 [NASA ADS] [CrossRef] [Google Scholar]
  6. Daisaka, H., & Ida, S. 1999, Earth, Planets, and Space, 51, 1195 [Google Scholar]
  7. Desch, S. J. 2007, ApJ, 671, 878 [NASA ADS] [CrossRef] [Google Scholar]
  8. Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Goldreich, P., & Tremaine, S. 1982, ARA&A, 20, 249 [NASA ADS] [CrossRef] [Google Scholar]
  10. Goldreich, P., & Tremaine, S. D. 1978, Icarus, 34, 227 [NASA ADS] [CrossRef] [Google Scholar]
  11. Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [NASA ADS] [CrossRef] [Google Scholar]
  12. Hayashi, C., Nakazawa, K., & Adachi, I. 1977, PASJ, 29, 163 [NASA ADS] [Google Scholar]
  13. Johansen, A., Henning, T., & Klahr, H. 2006, ApJ, 643, 1219 [NASA ADS] [CrossRef] [Google Scholar]
  14. Johansen, A., Oishi, J. S., Low, M.-M. M., et al. 2007, Nature, 448, 1022 [Google Scholar]
  15. Johansen, A., Youdin, A., & Mac Low, M. 2009, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
  16. Lesur, G., & Papaloizou, J. C. B. 2010, A&A, in press, [arXiv:0911.0663] [Google Scholar]
  17. Lithwick, Y., & Chiang, E. 2007, ApJ, 656, 524 [NASA ADS] [CrossRef] [Google Scholar]
  18. Lucy, L. B. 1977, AJ, 82, 1013 [NASA ADS] [CrossRef] [Google Scholar]
  19. Michikoshi, S., Inutsuka, S.-i., Kokubo, E., & Furuya, I. 2007, ApJ, 657, 521 [NASA ADS] [CrossRef] [Google Scholar]
  20. Michikoshi, S., Kokubo, E., & Inutsuka, S.-i. 2009, ApJ, 703, 1363 [NASA ADS] [CrossRef] [Google Scholar]
  21. Stewart, S. T., & Leinhardt, Z. M. 2009, ApJ, 691, L133 [NASA ADS] [CrossRef] [Google Scholar]
  22. Tanga, P., Weidenschilling, S. J., Michel, P., & Richardson, D. C. 2004, A&A, 427, 1105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Toomre, A. 1964, ApJ, 139, 1217 [NASA ADS] [CrossRef] [Google Scholar]
  24. Weidenschilling, S. J. 1977a, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
  25. Weidenschilling, S. J. 1977b, Ap&SS, 51, 153 [Google Scholar]
  26. Wisdom, J., & Tremaine, S. 1988, AJ, 95, 925 [NASA ADS] [CrossRef] [Google Scholar]
  27. Wolszczan, A., & Frail, D. A. 1992, Nature, 355, 145 [NASA ADS] [CrossRef] [Google Scholar]
  28. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]

Footnotes

... gas[*]
This would be true for a strongly collisional system, but is notformally valid in the studied regime.
... velocities[*]
This instability is qualitatively similar to the instability described by Goldreich & Tremaine (1978) for Saturn rings, although in the latter case the drag forces are absent.

All Figures

  \begin{figure} \par\includegraphics[width=8.9cm,clip]{12870f1.eps} \end{figure} Figure 1:

Shearing box, simulating a small patch of the protoplanetary disc. The grey box represents the shearing box of interest, the dashed boxes surrounding the central box represent one ghost ring.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[angle=270,width=8.8cm,clip]{12870f2.eps} \vspace*{-2mm} \end{figure} Figure 2:

Super-particles: velocity dispersion as a function of time in four different runs. The simulation labeled LS has a ten times larger smoothing length. The curves do not overlap because the simulations are not converged.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[angle=270,width=8.8cm,clip]{12870f3.eps} \vspace*{-2mm} \end{figure} Figure 3:

Scaled-particles: velocity dispersion as a function of time in three different runs that include physical collisions. All curves overlap because the simulations are converged.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=16.3cm,clip]{12870f4.ps} \end{figure} Figure 4:

Super-particles: snapshots of the particle distribution in the xy plane. Both simulations use 160 000 particles with different smoothing lengths. The simulation on the bottom uses a ten times larger smoothing length than the one on the top. The snapshots were taken ( from left to right) at $t=0, 30, 37, 40~\Omega^{-1}$. With a large smoothing length, the outcome looks very different, the system is more stable, more stripy structure can be seen and clumps form later, if at all.

Open with DEXTER
In the text

  \begin{figure} \includegraphics[width=16.3cm,clip]{12870f5.ps} \end{figure} Figure 5:

Scaled particles: snapshots of the particle distribution in the xy plane. The simulations ( from top to bottom) use 40 000, 320 000, 640 000 particles with their physical size given by Eq. (5). In all simulations, the collision timescale was kept constant. The snapshots were taken ( from left to right) at $t=0, 30, 37, 40~\Omega^{-1}$. The simulation with 40 000 particles has a large filling factor by the last frame, which prevents clumps from forming. The intermediate resolution simulation ( middle row) and the highest resolution simulation ( bottom row) have more and smaller particles and thus a smaller filling factor. The results (i.e., number of clumps in the last frame) are very similar in the intermediate and high resolution simulations.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[angle=270,width=9cm,clip]{12870f6} \end{figure} Figure 6:

Super-particles: density distribution spectrum with 160 000 particles and two different smoothing lengths that differ by a factor 10 at times t=30 (red, solid curve), t=36 (green, long dashed curve), and t=38 (blue, short dashed curved). The spectra at similar times are not on top of each other because the simulations are not converged.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[angle=270,width=9cm,clip]{12870f7} \end{figure} Figure 7:

Scaled particles: density distribution spectrum with 160 000, 320 000, and 640 000 particles (colors are equivalent to Fig. 6). At each time, the spectra are converged, i.e., the spectra are on top of each other and independent of the particle number, besides the noise level on very small scales. At late times ( $t\gtrsim 38$, blue curve), one begins to see more structure on small scales in simulations with larger number of particles (see text).

Open with DEXTER
In the text

Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.