EDP Sciences
Free Access
Issue
A&A
Volume 597, January 2017
Article Number A69
Number of page(s) 10
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201629561
Published online 03 January 2017

© ESO 2017

1. Introduction

One of the greatest problems in the theory of planet formation is to explain how millimeter- or centimeter-sized solid particles – in the following referred to as pebbles – grow to kilometer-sized planetesimals. Micron-sized dust grains can grow to pebble sizes by coagulation, but larger particles bounce or fragment under mutual collisions (Güttler et al. 2010; Zsom et al. 2010; Birnstiel et al. 2011). Growth might continue despite this so-called bouncing barrier for very porous ice particles (Wada et al. 2008, 2009) or by mass transfer in high-speed collisions (Wurm et al. 2005; Windmark et al. 2012).

A fundamental problem of planetesimal formation is the time constraint inflicted by the radial drift of solid particles, a problem that persists even under the assumption of perfect sticking. The orbital velocity of the gas in a protoplanetary disk is sub-Keplerian because the gas is supported against the gravity of the central star by a radial pressure gradient. The gas exerts a drag force on the particles in the disk, whose orbital speed would be equal to the Keplerian speed if the gas was not present, causing them to lose angular momentum and drift radially towards the star. The drift velocity depends on the size of the particles, but is in general highest for meter-sized particles in the inner regions of the disk. A particle with a size of 1 m, initially orbiting at a distance of 1 au from the star, drifts towards the star and sublimates in less than 100 years (Adachi et al. 1976; Weidenschilling 1977; Brauer et al. 2007). Hence, there needs to be some mechanism assisting the growth of pebbles into planetesimals, which are sufficiently large for the effect of the drag force exerted on them by the gas to be negligible, on a timescale shorter than the radial drift timescale.

The streaming instability provides a mechanism to concentrate solid materials and form planetesimals despite the poor sticking efficiency of the particles and their radial drift. It was discovered analytically by Youdin & Goodman (2005) and confirmed numerically by Youdin & Johansen (2007), Johansen & Youdin (2007), and Bai & Stone (2010a). The radial drift speed of solid particles decreases with increasing solid-to-gas density ratio because of the drag force exerted by the particles on the gas. A locally enhanced solid-to-gas ratio causes the local orbital velocity of the gas to be closer to Keplerian, and thus a reduction of the local drift speed of the particles. Hence, clusters of particles drift more slowly than isolated particles, and downstream clusters can accumulate upstream isolated particles, further reducing the drift speed of the clusters. Owing to this positive feedback loop, particles can be concentrated into filaments reaching maximum densities of up to several thousand times the local gas density (Bai & Stone 2010a; Johansen et al. 2012; Yang & Johansen 2014; Johansen et al. 2015), sufficient to undergo gravitational collapse and form planetesimals (Johansen et al. 2007; Simon et al. 2016). For this strong clustering of particles to occur, the solid-to-gas column density ratio needs to exceed a critical value (Johansen et al. 2009b; Bai & Stone 2010b), which depends on the radial pressure gradient supporting the gas (Bai & Stone 2010c) and the particle size (Carrera et al. 2015; Yang et al. 2016).

Although several studies have shown that the streaming instability can lead to the formation of planetesimals, their birth mass distribution has not been comprehensively investigated. However, the initial mass function of planetesimals is essential for the study of the formation of planetary systems because it determines the initial conditions for the evolution of the bodies that planetesimals evolve into, including planets, asteroids, and Kuiper belt objects. The asteroids in the asteroid belt provide a natural sample distribution that can be fitted with a broken power law. Bottke et al. (2005) argue that the current size distribution of asteroids larger than 120 km in diameter represents the birth size distribution of the planetesimals that formed in the asteroid belt (but have been strongly depleted by resonances with Jupiter, independent of their sizes), while smaller asteroids are largely fragments of collisions between the larger ones.

Both Johansen et al. (2015) and Simon et al. (2016) performed numerical simulations of planetesimal formation by the streaming instability and find that the differential distribution of the planetesimal birth masses is well-fitted with a power law with an exponent of about −1.6, albeit with the difference that, while the former observe an exponential tapering of the power-law distribution that constitutes the physical upper mass cutoff, the latter do not include such a tapering in their fits. In this paper, we compare power-law fits with and without exponential cutoff to evaluate how well the high-mass end of the initial mass function is described by an exponential cutoff.

Johansen et al. (2015) find the shape of the initial mass function to be relatively independent of the resolution of the simulations and the solid particle column density. They show that a higher resolution leads to the formation of planetesimals with a wider range of sizes, between 30 km and 120 km in radius in their simulation with the highest resolution because the size of the smallest planetesimal declines with increasing resolution. On the other hand, they observe the size of the largest planetesimal to mainly depend on the particle column density, with smaller column densities yielding smaller sizes. Simon et al. (2016) also studied the dependence of the shape of the birth mass distribution on the resolution of the simulations and obtain the same result. They further find the shape of the distribution to be largely independent of the strength of the self-gravity and the simulation time at which it is initiated, although the masses of the planetesimals are shifted to higher values with the increasing strength of self-gravity. The planetesimals that formed in their simulations typically range in radius from 50 km to a few hundred kilometers.

It remains unclear if the planetesimal initial mass function depends on the dimensions of the simulation box. Both Johansen et al. (2015) and Simon et al. (2016) employed only one box size of 0.2 gas scale heights in the radial, azimuthal, and vertical directions. However, Yang & Johansen (2014) find that, while in the simulations with this box size, the solid particles are concentrated by the streaming instability into only one axisymmetric filament, multiple of these filaments form in simulations with larger box dimensions. This raises the question of whether the mass budget of planetesimal formation, and thus the shape of the initial mass function, is different when not only one filament is observed. In this paper, we study simulations with three different box sizes, the smallest of which is equal to that employed by Johansen et al. (2015) and Simon et al. (2016), while the others are two and four times larger, respectively, in the radial and azimuthal directions, which permits investigating planetesimal formation in several filaments. Furthermore, in simulations with larger box sizes, more planetesimals emerge, yielding better statistics in particular for the determination of the initial mass function.

The paper is structured as follows: in Sect. 2, the simulation setup, i.e. the initial conditions and the parameters that govern the evolution of the simulations are described. In Sect. 3, we present our results regarding the formation of planetesimals by the streaming instability and their radial migration. Further, we comment on the issue of permitting the mutual accretion of sink particles, which we use to model planetesimals. In Sect. 4, we discuss whether the planetesimal birth mass distribution is exponentially tapered and how its shape depends on the dimensions of the simulation box as well as the resolution. We conclude in Sect. 5.

2. Simulation setup

Table 1

Simulation specifications.

We conduct three-dimensional computer simulations with the Pencil Code1, a hybrid code for gas, for which the magnetohydrodynamic equations are solved on a fixed grid, with Lagrangian particles representing solid bodies. The code employs sixth-order finite differences in space and third-order Runge-Kutta steps in time.

We use the shearing box approximation (Goldreich & Lynden-Bell 1965), i.e. we assume that the size of the simulation box is small compared to the distance to the central star of the protoplanetary disk. Hence, the curvature of the disk is neglected and the stellar gravity is linearized. The rectangular simulation box is aligned such that the x-, y-, and z-directions correspond to the radial, azimuthal, and vertical directions, respectively, and co-rotates with the Keplerian velocity at its origin. For both gas and particles, sheared periodic boundary conditions are employed at the radial and azimuthal boundaries and periodic boundary conditions at the vertical boundaries (Hawley et al. 1995; Brandenburg et al. 1995; Youdin & Johansen 2007; Johansen et al. 2009a).

In total, we perform five simulations with three different simulation box dimensions and two different resolutions, as listed in Table 1. The two smaller boxes have a size of 0.2 and 0.4 gas scale heights, respectively, in the radial and azimuthal directions with a resolution of either 320 or 640 grid cells per scale height, while the largest box has a radial and azimuthal size of 0.8 scale heights with a resolution of 320 grid cells per scale height. All simulation boxes have a vertical size of 0.2 scale heights. The names of the simulations are composed of the radial and azimuthal dimension as the first number and the resolution as the second number.

2.1. Gas

The simulation box is filled with an isothermal, non-magnetized gas with an isothermal equation of state , where pg and ρg are the pressure and density, respectively, and cs is the (constant) sound speed. While the gas density is initially constant in the radial and azimuthal direction, it is stratified in the vertical direction because we take into consideration the vertical gravity of the central star, which causes both gas and solid particles to sediment to the mid-plane at z = 0. This background density stratification is determined by the equilibrium between vertical gravity and vertical pressure gradient, and is given by (1)where Hg = cs/ ΩK is the gas scale height, ΩK = 2π/PK the Keplerian orbital frequency, and PK the Keplerian orbital period. Here and in the following, the subscript zero refers to the mid-plane. As formulated in Yang & Johansen (2014), we subtract the background density stratification from the equations of the motion for the gas to numerically balance this equilibrium state down to machine precision.

Since the gas density is initially radially constant, there is no radial pressure gradient to support the gas and cause it to orbit with sub-Keplerian speed. Hence, a background pressure gradient set by the dimensionless parameter (2)where R is the orbital distance, is imposed. We refer to this parameter as Π adopting the notation by Bai & Stone (2010b). The resulting sub-Keplerian orbital velocity of the gas is given by vg = vK − Δv, where vK = ΩKR is the Keplerian orbital velocity and Δv = Πcs = 0.05cs, which is a representative value at orbital distances of the order of 1 au in a typical protoplanetary disk (Hayashi 1981; Bai & Stone 2010b; Bitsch et al. 2015).

2.2. Particles

Two types of Lagrangian particles are employed in the simulations: super-particles and sink particles to model pebbles and planetesimals, respectively. To achieve a good load balancing among the processors, we use the particle block domain decomposition algorithm implemented by Johansen et al. (2011). For the calculation of the mutual drag forces between the gas and the super-particles and the mutual gravitational forces between the super-particles and the sink particles, we apply the triangular shaped cloud scheme to map the particle masses and velocities onto the grid, and similarly interpolate the back-reaction drag forces and the self-gravitational forces onto the particles (Hockney & Eastwood 1981; Youdin & Johansen 2007; Johansen et al. 2007).

The gravitational potential of the super-particles as well as the sink particles is computed by solving Poisson’s equation using the fast Fourier transform algorithm (Gammie 2001), which entails gravitational softening. The softening length is of the order of the grid cell edge length. Even though the gravitational potential of the particles is vertically periodic, the particles are concentrated in a thin layer around the mid-plane such that their dynamics are not affected by the periodic potential away from the mid-plane.

We neglect the contribution of the gas to the gravitational potential under the assumption that the gas component of the protoplanetary disk is not in the gravitationally unstable regime. Indeed, in our simulations the gas density deviates by at most 2% from the mid-plane density ρg,0, the gas density perturbations are thus very small compared to the super-particle densities in the filaments when planetesimals form, which are of the order of the Roche density (see below).

The sedimentation of the super-particles to the mid-plane that is due to the vertical stellar gravity induces turbulence as a result of either the streaming instability or the Kelvin-Helmholtz instability, which are both caused by the mutual drag forces between the gas and the super-particles (Bai & Stone 2010b). This turbulence stirs up the super-particles, and hence counteracts the sedimentation. To give sedimentation and turbulence time to attain an equilibrium, in our simulations self-gravity is not introduced until t = 25 PK. Its strength is then gradually increased from zero over 10 PK until it reaches its final value at t = 35 PK since initiating self-gravity instantaneously with full strength could cause significant impulses on the particles. While we initiate self-gravity at a simulation time at which the super-particles have already formed filaments, Johansen et al. (2015) introduced self-gravity with full strength at the start of their simulations and observe qualitatively the same planetesimal birth mass distribution as we do.

We achieve this gradual initiation of the self-gravity by substituting (3)where the dimensionless self-gravity parameter (4)and G is the gravitational constant, into the right-hand side of Poisson’s equation. We choose , and thus γ = 1 /π = 0.318. We note that Simon et al. (2016) find the shape of the initial mass function of the planetesimals formed in their simulations of the streaming instability to be relatively independent of both the simulation time at which self-gravity is introduced and the strength of the self-gravity. The Roche density depends on the self-gravity parameter γ and is given by (5)Each super-particle represents a large number of equally sized pebbles because it is computationally infeasible to simulate the pebbles individually. While the mass of a super-particle is equal to the total mass of the pebbles it models, its friction time is the same as that of an individual constituent pebble.

The mass of the super-particles is determined by the initial solid-to-gas column density ratio and their initial number. We set the solid-to-gas ratio (6)where Σp,init and (7)are the initial column densities of the super-particles and the gas, respectively, to Z = 0.02. This value corresponds to the critical solid-to-gas ratio necessary for strong clustering of pebbles due to the streaming instability to occur (Johansen et al. 2009b; Bai & Stone 2010b,c; Carrera et al. 2015), and is slightly higher than the solar metallicity. The initial number of super-particles is set equal to the total number of grid cells. The super-particles are randomly distributed among the entire simulation box to seed the streaming instability.

The Stokes number of the super-particles τf = ΩKtf, where tf is the friction time, is set to τf = π/ 10 = 0.314, which, at an orbital distance of 2.5 au in the Minimum Mass Solar Nebula (MMSN), corresponds to a size of approximately 25 cm (Bai & Stone 2010b; Johansen et al. 2015). While we employ only one fixed particle size, Bai & Stone (2010b) performed simulations with particles of a range of sizes. They find that the particle-gas dynamics are dominated by the most massive particles, and that the critical solid-to-gas ratio required for strong particle clustering owing to the streaming instability is determined by the total mass of all particles.

After self-gravity has attained its full strength at t = 35 PK, every super-particle comprised in a cluster whose super-particle density ρp exceeds a threshold value ρp,thres is replaced by a sink particle. This sink particle creation threshold is set to ρp,thres = 200 ρg,0, i.e. about seven times the Roche density (see Eq. (5)). We note that Johansen et al. (2015) compared simulations similar to ours with three different threshold values and find the masses of the sink particles emerging in their simulations to be largely independent of the threshold above a value of five times the Roche density.

The simulation time at which the formation of sink particles is introduced is arbitrary since both the gravitationally bound super-particle clusters that exist beforehand and the sink particles that emerge from them afterwards represent planetesimals. Owing to the limited resolution of the gravitational forces the behavior of many super-particles inside one grid cell is comparable to that of a few sink particles. Nevertheless, the computational expense of the simulations is lowered substantially by the introduction of sink particles.

Super-particles within the accretion radius of a sink particle, which is set equal to one grid cell edge length, are accreted by it, i.e. the super-particle mass and momentum are added to the sink particle mass and momentum, respectively, and the super-particle is removed. This accretion, however, might in parts be artificial because the physical accretion radius, i.e. the Bondi radius for pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012), could be smaller than the simulated accretion radius, especially in the case of less massive sink particles. On the other hand, since the mutual gravitational forces between the super-particles and the sink particles within one grid cell can only be computed inaccurately, the chosen accretion radius corresponds to the highest accuracy our numerical simulations can offer.

We further permit sink particles to accrete one another. This accretion is handled analogously to the super-particle accretion, and might as well be partially artificial. Nevertheless, it is required because in a super-particle cluster exceeding the sink particle creation threshold, all super-particles are replaced by sink particles. That is, the gravitational collapse of a super-particle cluster results in the creation of a cluster of sink particles, of which only one should persist to represent one new-born planetesimal. See Sect. 3.2 for further discussion of this topic.

2.3. Units and scaling relations

We report our results using the Keplerian orbital period at the origin of the simulation box PK, the gas scale height Hg, and the mid-plane gas density ρg,0 as the units of time, length, and density, respectively. We note that a shearing box freely scales with these units until self-gravity is initiated at t = 25 PK. Afterwards, (8)(see Eq. (4)), and hence the unit of mass .

In the following, relations for the scaling of relevant quantities and units with the orbital distance R, the temperature T, and the mass of the central star MS are given. A mean atomic weight of μ = 2.33 is used and the scaling relations for ρg,0 and [M] are calculated applying those for PK and Hg: Hereafter, we use the above scaling relations, the properties of the asteroid belt, i.e. an orbital distance of R = 2.5 au, a temperature of T = 180 K, and a stellar mass of MS = 1 M, and the chosen strength of the self-gravity γ = 1 /π to convert simulation units into physical units. For example, the mid-plane gas density ρg,0 = 9.4 × 10-10 g cm-3, which is almost one order of magnitude greater than the corresponding value in the MMSN, ρg,0 = 1.1 × 10-10 g cm-3 (Hayashi 1981; Bai & Stone 2010b). The streaming instability has been shown to form planetesimals for pebble column densities similar to that in the MMSN (Johansen et al. 2015). Nevertheless, we choose this comparably high gas density and thus high solid density to promote planetesimal formation, enabling us to better constrain the initial mass function of planetesimals.

3. Evolution of the simulations

Apart from the use of self-gravity and sink particles to model planetesimals, our simulation setup is identical with that applied by Yang & Johansen (2014). Hence, we find our simulations to be consistent with their simulations until t = 25 PK, when self-gravity is introduced (compare the evolution of the pebble density in the simulation time span 20 PKt ≤ 25 PK shown in Fig. 1 with their Fig. 3). We thus only report on the evolution of our simulations between t = 25 PK and the end of the simulations, t = 40 PK.

thumbnail Fig. 1

Pebble column density Σp, integrated over the vertical dimension and averaged over the azimuthal dimension of the simulation box, as a function of radial location x and simulation time t for the simulation with the largest box size, run_0.8_320. Locations at which sink particles emerge are indicated with white dots. Though every pebble cluster forms in one of the filaments, they migrate up to the distance to one of the adjacent filaments. All sink particles emerge as soon as their formation is initiated at t = 35 PK, and they emerge nearly evenly distributed among the entire radial extent of the box owing to the radial migration of the pebble clusters.

Open with DEXTER

3.1. Planetesimal formation and migration

If self-gravity is not taken into account the streaming instability concentrates pebbles into axisymmetric filaments that are elongated in the azimuthal direction (Johansen & Youdin 2007; Johansen et al. 2009b; Bai & Stone 2010b; Yang & Johansen 2014). In Figs. 1 and 2, we show how self-gravity causes these filaments to fragment into pebble clusters that undergo gravitational collapse and form planetesimals. In the first 5 PK after self-gravity has been initiated, the filaments begin to disperse because the pebbles accumulate into clusters owing to their mutual gravitational attraction (upper panels of Fig. 2). At t = 35 PK, when the self-gravity reaches its full strength, these clusters contain most pebbles available in the simulation, and the filaments are no longer observable (lower left panel).

At this point, we commence the formation of sink particles. Almost all pebble clusters exceed the sink particle creation threshold, consequently the pebbles in each of these clusters are replaced by sink particles which merge into one massive sink particle that represents the gravitationally collapsed cluster. At t = 35.1 PK, this sink particle merging process is for the most part completed (lower right panel). However, a few low-mass pebble clusters do not turn into sink particles. All sink particles emerge instantly at t = 35 PK, as can be seen from Fig. 1, but a couple of clusters remain at t = 35.1 PK (three such clusters can be spotted in the lower right panel of Fig. 2). Although they are not sufficiently dense to exceed the sink particle creation threshold, these clusters probably represent gravitationally bound planetesimals with low masses.

We observe that the planetesimals on average move through more than half of the radial extent of the simulation boxes. Figure 1 shows that each pebble cluster forms in one of the filaments, but that they migrate in the radial direction, some of them only marginally, others the entire distance to one of the adjacent filaments. As a result, the sink particles emerge almost evenly distributed among the whole radial dimension of the box. From Fig. 3, it can be seen that the sink particles continue this radial migration, they on average pass through over half of, a few of them even through the whole radial extent of the box. The mean standard deviation of the radial displacement of the sink particles from the locations at which they emerge, averaged over and weighted by the lifetime of every sink particle and then averaged over all sink particles in a simulation, amounts to between 26% and 36% of the radial box size in each of the five simulations. We note that the extent of the migration in the radial direction increases with the box size without converging for the box sizes we consider. The radial motions result from the mutual gravitational scattering of sink particles that closely pass by each other. It remains to be investigated whether planetesimals are composed of not only pebbles from the filament they form in, but also of an appreciable amount of pebbles from filaments they migrate to.

thumbnail Fig. 2

Pebble column density Σp, integrated over the vertical dimension of the simulation box, as a function of radial location x and azimuthal location y at four different simulation times t = 25 PK (top left panel), t = 30 PK (top right panel), t = 35 PK (bottom left panel), and t = 35.1 PK (bottom right panel) for the simulation with the largest box size, run_0.8_320. In the lower right panel, sink particles are plotted as white dots and three pebble clusters are indicated using white circles. After self-gravity has been initiated at t = 25 PK, the pebbles aggregate into clusters and the axisymmetric filaments disperse (upper panels). When the self-gravity attains its full strength at t = 35 PK, most pebbles are concentrated in clusters and the filaments are no longer visible (lower left panel). At this point, we introduce the formation of sink particles, and the pebbles comprised in all but the three encircled clusters are replaced by sink particles. The sink particles emerging from the same pebble cluster undergo a merging process until only one of them remains. This process is largely completed at t = 35.1 PK (lower right panel).

Open with DEXTER

3.2. Mutual sink particle accretion

For the sink particles to realistically represent new-born planetesimals, only one of them should be allowed to form from every pebble cluster that undergoes gravitational collapse. However, this requires a precise determination of the extent of each cluster, which entails several issues, for instance the treatment of overlapping clusters. Therefore, we replace every super-particle that is part of a cluster that exceeds the sink particle creation threshold by a sink particle and allow the sink particles to accrete one another until only one of them remains. We observe that in our five simulations, on average 81% of all accreted sink particles are accreted within 0.1 PK after their formation, i.e. until t = 35.1 PK, because in all simulations all sink particles emerge togather at t = 35 PK. Therefore, the merging process of sink particles that emerged from the same pebble cluster is probably largely completed at this point (see the lower right panel of Fig. 2), and most of the mutual sink particle accretions are a part of this process.

On the other hand, we note that the merging of sink particles afterwards may be artificial. The accuracy of the calculation of the gravitational forces between sink particles is limited by the resolution, hence we cannot determine whether two sink particles that encounter one another inside a grid cell collide or pass by each other. Nevertheless, we find the latter to be more probable. Taking into account gravitational focusing, the maximum impact parameter leading to a collision of two sink particles (15)where mi are the masses of the two sink particles, ri their radii, which are calculated from mi using a solid body density of 3 g cm-3, and Δv is their relative velocity (see, e.g., Armitage 2007). The mean maximum impact parameter bmax, averaged over all mutual sink particle accretions in all five simulations where the lifetime of the accreted sink particle is greater than 0.1 PK and weighted by the lifetimes of the accreted sink particles, amounts to only 3.5% of the grid cell edge length. Mutual accretions of three or more sink particles at the same simulation time are not included in this statistic because we cannot infer their outcome. We note, though, that the sink particle data are not written out after each simulation time step, but every 0.01 PK. Thus, we can only imprecisely determine the simulation time at which a sink particle merging occurs and the maximum impact parameter for this encounter.

4. Initial mass function

4.1. Best-fitting parameters

thumbnail Fig. 3

Radial locations x of sink particles which are not accreted before t = 36 PK, color-coded according to the locations at which they emerge, as functions of the simulation time t for the simulation with the largest box size, run_0.8_320. On average, the motions of the sink particles span more than half of the radial extent of the simulation box, a couple of them even migrate through the entire radial dimension of the box.

Open with DEXTER

We fit the cumulative mass distribution of the sink particles that emerge in our simulations using two functional forms: an integrated power law and a power law with an exponential cutoff. We choose to fit the cumulative mass distributions because, in particular for small numbers of sink particles, they are less affected by noise than the differential mass distributions and can thus be fitted more accurately.

The integrated power law can be expressed as (16)where N> is the number of sink particles with masses greater than M, Ntot is their total number and Mmax their maximum mass, and Mpow and α are fitting parameters. Since the formation of planetesimals by the streaming instability is a stochastic process and the actual Mmax in a simulation might differ significantly from the Mmax of the ensemble-averaged mass distribution of the sink particles for this model, we treat Mmax as a fitting parameter.

The exponentially tapered power law is given by (17)where Mexp and β are fitting parameters. The condition N>(Mmin) /Ntot = 1, where Mmin is the minimum sink particle mass, can be used to eliminate one of the fitting parameters in Eq. (17). We choose to eliminate the characteristic mass of the power-law part, Mpow, which is then given by (18)Substituting Eq. (18) into (17) yields (19)We use Eq. (19) to fit the sink particle data because we find the resulting fits to be better than the ones we obtain applying Eq. (17). For the same reason as Mmax in Eq. (16), Mmin is treated as a fitting parameter.

To investigate the dependence of the shape of the initial mass function on the resolution and particularly the dimensions of the simulation box, we determine an individual initial mass function for every simulation. At first, we employ the least-squares method to fit Eqs. (16) and (19) to the sink particle data at each Keplerian orbital period between t = 36 PK and the end of the simulations, t = 40 PK. In all five simulations, all sink particles emerge at once at the simulation time at which their formation is initiated, t = 35 PK, but we begin the fitting at t = 36 PK to give the sink particles that emerge from the same pebble cluster time to merge into one. Averaged over this period, we then calculate mean values of the fitting parameters that do not vary significantly with time and are thus probably relatively unaffected by artificial sink particle merging.

In the left panel and the right panel of Fig. 4 we show the sink particle mass distributions as well as the fitted integrated power laws and power laws with exponential tapering for the two simulations with the largest numbers of sink particles, the one with the largest box size, run_0.8_320, and the one with the middle box size and the higher resolution, run_0.4_640, respectively, at t = 40 PK. In the legends for both fits, the standard deviation σ of the actual N>(M) /Ntot for the sink particle masses M from the fitted N>(M) /Ntot is given.

We find that the exponential tapered power laws fit the mass distributions better than the integrated power laws. The standard deviations for the power-law fits without exponential cutoff are larger, not only at t = 40 PK as shown in the figure, but also at other simulation times. We further note that the shallower exponential cutoff represents the high-mass end of the distributions better than the steeper cutoff of the integrated power law, and that the exponential tapering better reproduces the smooth change of the slope of the distributions. Hence, we limit our further analysis to the power-law fits with exponential cutoff.

4.2. Power-law fits with and without exponential tapering

thumbnail Fig. 4

Cumulative mass distributions of sink particles (black crosses) and fitted integrated power laws (Eq. (16), green lines) and power laws with exponential cutoff (Eq. (19), red lines) at t = 40 PK for the two simulations with the largest numbers of sink particles: the one with the largest box dimensions, run_0.8_320, (left panel) and the one with the middle box dimensions and the higher resolution, run_0.4_640 (right panel). Standard deviations of the actual N>/Ntot (black crosses) from the fitted N>/Ntot are given in the legends. The standard deviations for the integrated power laws are up to twice as large as those for the exponentially tapered power laws. In addition, it can be seen that the shallower exponential cutoffs fit the actual cutoffs more accurately than the steeper cutoffs of the integrated power laws and better replicate the smooth change of the slope of the mass distributions. We thus find the power laws with exponential tapering to represent the mass distributions better than the integrated power laws.

Open with DEXTER

In Fig. 5 the best-fitting parameters Mmin, α, Mexp, and β in Eq. (19) are shown for t = 36 PK to t = 40 PK for all of our simulations. In the case of the simulation with the smallest box dimensions and the lower resolution, run_0.2_320, the small number of sink particles persisting after t = 36 PK does not permit to properly fit a power law and an exponential cutoff. The parameters α, Mexp, and β are approximately constant in time for all simulations apart from Mexp for the simulation with the middle box dimensions and the lower resolution, run_0.4_320, which seems to converge to an upper limit. In Table 2, we list mean values of these three parameters, averaged over the simulation time span from t = 36 PK to t = 40 PK.

Similar to Mexp for run_0.4_320, the fitted minimum mass Mmin increases with time for the simulation with the largest box size, run_0.8_320, and the one with the middle box size and the higher resolution, run_0.4_640, but appears to saturate at an upper limit. The increase of Mmin and Mexp is likely due to both pebble accretion by the sink particles – as less and less pebbles remain towards the end of the simulations (see Fig. 1), both masses converge to an upper limit – and possibly artificial merging of sink particles. In contrast to this, for run_0.4_320 and the simulation with the smallest box dimensions and the higher resolution, run_0.2_640, Mmin is roughly constant. To provide a comparable value for each simulation despite these differences in the time dependence, in Table 2 we list the best-fitting values at t = 40 PK for all simulations (except for run_0.2_320, for which only the value at t = 36 PK is available). A comparison of these values shows that Mmin declines with both increasing box dimensions and increasing resolution if the value for run_0.2_320 is disregarded.

The best-fitting values of the exponent of the power-law component α for all simulations vanish. This is because a lower limit for the sizes of the pebble clusters, and thus the sink particle masses, is set by the resolution, which even in the case of the higher resolution is too large for low-mass sink particles that would constitute the power-law part to emerge. Both Johansen et al. (2015) and Simon et al. (2016) fit the differential mass distribution with a power-law exponent of about −1.6, corresponding to α = 0.6. Since the higher resolution we employ, , is the lowest one that is considered in these papers, to properly study the power-law distribution higher resolutions than seem to be required.

We find the mass scale of the exponential cutoff Mexp to correlate with the mass budget in every filament. The parameter Mexp should increase with the distance between the filaments because, if the distance is larger, more pebbles can be accreted by the planetesimals forming in each filament. In Col. 6 of Table 2, the numbers of filaments Nf we observe in our simulations are listed. One, three, and four filaments form in the simulation boxes with radial and azimuthal dimensions of 0.2 Hg, 0.4 Hg, and 0.8 Hg, respectively. Therefore, the mass reservoir of pebbles in each filament is similar for the smallest and the largest box sizes, but smaller for the middle box sizes. We indeed see that the mean values of Mexp for run_0.2_320 and run_0.8_320 are similar, but larger than the one for run_0.4_320. Likewise, the mean value for run_0.2_640 is greater than that for run_0.4_640. Even though the mean values for the two simulations with the middle box size differ by more than one standard deviation, we note that the best-fitting values for run_0.4_320 increase with time with a range enclosing the mean value for run_0.4_640. Hence, we find Mexp to be largely independent of the resolution.

Table 2

Best-fitting parameters.

With the mean values of the exponent of the exponential cutoff β ranging from 0.28 to 0.38, the exponential cutoff is rather smooth. Johansen et al. (2015) fit their data by eye using β = 4 / 3, which is a significantly steeper cutoff, but this steep cutoff might be an artefact of the small box size they employed, in which only a few massive planetesimals formed. The best-fitting values for the two simulations with the largest number of sink particles, run_0.8_320 and run_0.4_640, are nearly equal, but somewhat greater than the values for run_0.4_320 and run_0.2_640, which are also roughly equal. This indicates that only in the former two simulations enough sink particles emerge to completely capture the high-mass end of the mass distribution. However, the mean values for all simulations lie in a rather small range of 0.1, hence we find β to be relatively independent of the box size and the resolution.

Substituting the best-fitting parameters listed in Table 2 into the cumulative (Eq. (19)) or the differential mass distribution, (20)yields an initial mass function for each simulation. The cumulative mass distribution can be converted to a cumulative size distribution, (21)where R is the radius of every sink particle and the minimum radius Rmin and the radius scale of the exponential tapering Rexp can be calculated from Mmin and Mexp, respectively, using a solid density of 3 g cm-3.

Table 3

Sink particle statistics.

In Table 3, we list the number Ntot, minimum mass Mmin, maximum mass Mmax, and mean mass M of the sink particles at the end of our five simulations, t = 40 PK. In the case of the two simulations with the smallest box dimensions, we observe less than ten sink particles, and the best-fitting parameters for these simulations are therefore afflicted with comparably large errors (see Fig. 5). We find the actual Mmin and M to be of the order of the fitted Mmin and Mexp, respectively (compare with Cols. 2 and 4 of Table 2). For run_0.8_320 and run_0.4_640, the 10% of the sink particles which are most massive contain 66% and 70%, respectively, of the total sink particle mass. That is, in our simulations the most massive sink particles dominate the total mass.

As stated above, a higher resolution enables us to observe the formation of smaller pebble clusters, and thus less massive sink particles. Hence, Ntot increases and both Mmin and M decline with increasing resolution. Like Johansen et al. (2015) and Simon et al. (2016), we find the maximum mass Mmax to be relatively independent of the resolution.

thumbnail Fig. 5

Best-fitting parameters Mmin (upper left panel), α (upper right panel), Mexp (lower left panel), and β (lower right panel) of the exponentially tapered power law (Eq. (19)) at every Keplerian orbital period between t = 36 PK and t = 40 PK for all five simulations. Standard errors are plotted as error bars. In the case of the simulation with the smallest box size and the lower resolution, run_0.2_320, only parameter values for t = 36 PK are plotted because afterwards only four sink particles persist, which prevents us from properly fitting the mass distributions with a power law and an exponential tapering.

Open with DEXTER

We expect the number of planetesimals to increase with the number of filaments, and thus with the radial box dimension, and with the length of the filaments, i.e. the azimuthal box dimension. Analogously to the mass scale of the exponential cutoff Mexp, as discussed above, we further expect M to increase with the distance between the filaments. Our findings are consistent with these expectations, with the exception of M for run_0.2_320. One and three filaments form in run_0.2_320 and run_0.4_320, respectively (see Col. 6 of Table 2), therefore the value for the former simulation should be greater than the one for the latter simulation, yet it is slightly smaller. This shows that, at least at the lower resolution, the smallest boxes, in which only one filament forms, might be too small to accurately capture the mass budget of each filament.

Furthermore, Mmin and Mmax in general increase with the box dimensions, but we find Mmin to considerably decrease if the box size is increased from 0.2 Hg to 0.4 Hg in the radial and azimuthal directions. This also indicates that the smallest box dimensions might not capture the sink particle mass distribution as well as the larger boxes. However, it may also be a stochastic effect because, especially for a small number of sink particles, the ensemble-averaged values of Mmin and Mmax might differ significantly from the actual values.

5. Summary and discussion

We have investigated the formation of planetesimals by the streaming instability in numerical simulations with three different box sizes and two different resolutions. In particular, we have studied the initial mass function of these planetesimals, employing the largest box dimensions to date with radial and azimuthal sizes of up to 0.8 gas scale heights. These large box sizes have enabled us to study planetesimal formation in multiple axisymmetric filaments formed by the streaming instability and to yield better statistics because more planetesimals emerge in simulations with larger box sizes.

In the absence of self-gravity, the streaming instability concentrates pebbles into axisymmetric filaments. After self-gravity has been introduced, these filaments disperse within about ten Keplerian orbital periods because the pebbles accumulate into clusters that undergo gravitational collapse and form planetesimals. We have observed that, after their formation, the planetesimals on average migrate through more than half of the radial dimension of the simulation box owing to mutual gravitational scattering. The extent of the radial migration does not converge for the box sizes we have taken into consideration. Further studies could provide insights regarding the implications of the migration through multiple filaments for the dependence of the composition of planetesimals on the orbital distance.

The radii of the planetesimals formed in our simulations, which depend on the strength of the self-gravity and thus on the solid particle column density, range from 80 km to 620 km. We have compared power-law fits to their cumulative mass distribution with and without exponential tapering and have found that a rather shallow exponential cutoff fits the distribution better than the steeper cutoff of an integrated power law. Johansen et al. (2015) also find the initial mass function to be represented best by an exponentially tapered power law, although they studied planetesimals that are smaller than the ones formed in our simulations. In their simulation with the highest resolution, the planetesimal radii amount to between 30 km and 120 km. In contrast to this, Simon et al. (2016) find that a power law without exponential tapering is suitable to fit the birth mass distribution of planetesimals with radii between 50 km and a few hundred kilometers.

We have found a value of the exponent of the exponential cutoff of about 0.3 to 0.4, which is largely invariant under changes in the box size and the resolution, but considerably smaller than the value of 4/3 that Johansen et al. (2015) determine, which is based on a much smaller number of massive planetesimals that formed in a small simulation domain. However, the resolutions we have considered are insufficient to constrain the shape of the power-law part and to investigate its dependence on the box dimensions because the planetesimals that formed in our simulations are too large to constitute a power-law distribution at the low-mass end of the initial mass function.

Both the characteristic mass of the exponential cutoff and the mean mass of the planetesimals correlate with the pebble mass budget in every filament. In this regard, we have found indications that a simulation box with a size of 0.2 gas scale heights in the radial, azimuthal, and vertical directions, in which only one filament emerges, may be too small to properly capture the mass reservoir. This is consistent with the observation by Yang & Johansen (2014) that box dimensions of 0.2 scale heights in the radial and azimuthal directions are too small to capture all scales relevant for the streaming instability, and is further supported by Li et al. (in prep.) finding that the density distribution function of solid particles is consistent only for boxes with radial and azimuthal sizes of at least 0.4 scale heights.

The current size distribution of the asteroids in the asteroid belt with diameters between 120 km and several hundred kilometers, corresponding to the sizes of the planetesimals that emerge in our simulations, is well-fitted with a power law. This power law, which Bottke et al. (2005) argue represents the primordial asteroid size distribution, is in contrast to the exponential cutoff we (and Johansen et al. 2015) find. Subsequent pebble accretion therefore appears to be necessary to convert the exponential tapering of the birth size distribution into the power law observed in the asteroid belt (Johansen et al. 2015).

It is interesting to compare the initial mass function of planetesimals to the classical concept of an initial mass function of stars. The formation of stars is comparable to the formation of planetesimals by the streaming instability insofar as both stars and planetesimals form by gravitational collapse, the former from molecular cloud cores, and the latter from pebble clusters. The differential mass distribution of stars with masses greater than about 1 M is given by a power law with an exponent of approximately −2.3 (Salpeter 1955; Massey 1998; Chabrier 2003). That is, the total stellar mass is dominated by the least massive stars, in contrast to the total mass of the planetesimals we have observed in our simulations, which is dominated by the most massive ones. It has been argued that there is a physical upper mass cutoff of the stellar initial mass function, but massive stars are rare and short-lived, and their mass distribution is therefore difficult to observe (Zinnecker & Yorke 2007).

Johansen et al. (2015) and Simon et al. (2016) investigated the dependence of the shape of the planetesimal initial mass function on the resolution, the pebble column density, the strength of the self-gravity, and the simulation at which self-gravity is initiated. We have complemented these parameter studies with an analysis of the box-size dependence, but how, for instance, the solid-to-gas ratio, the friction time of the pebbles, the radial gas pressure gradient, and the vertical box size affect the shape of the birth mass distribution remains to be investigated. Finally, the streaming instability has been shown to operate in protoplanetary disks with turbulence driven by the magnetorotational instability (Johansen et al. 2007). It remains to be seen how turbulence and magnetic fields influence the shape of the initial mass function.


Acknowledgments

We thank Piero Ranalli for his advice on how to calculate the standard errors of the fitting parameters. We further thank the anonymous referee for providing comments and questions that helped to improve the paper. The simulations presented in this paper were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC), the Alarik system at Lunarc at Lund University and the Beskow system at the PDC Center for High Performance Computing at the KTH Royal Institute of Technology. This research was supported by the European Research Council under ERC Starting Grant agreement 278675-PEBBLE2PLANET. A.J. is thankful for financial support from the Knut and Alice Wallenberg Foundation and from the Swedish Research Council (grant 2014-5775).

References

All Tables

Table 1

Simulation specifications.

Table 2

Best-fitting parameters.

Table 3

Sink particle statistics.

All Figures

thumbnail Fig. 1

Pebble column density Σp, integrated over the vertical dimension and averaged over the azimuthal dimension of the simulation box, as a function of radial location x and simulation time t for the simulation with the largest box size, run_0.8_320. Locations at which sink particles emerge are indicated with white dots. Though every pebble cluster forms in one of the filaments, they migrate up to the distance to one of the adjacent filaments. All sink particles emerge as soon as their formation is initiated at t = 35 PK, and they emerge nearly evenly distributed among the entire radial extent of the box owing to the radial migration of the pebble clusters.

Open with DEXTER
In the text
thumbnail Fig. 2

Pebble column density Σp, integrated over the vertical dimension of the simulation box, as a function of radial location x and azimuthal location y at four different simulation times t = 25 PK (top left panel), t = 30 PK (top right panel), t = 35 PK (bottom left panel), and t = 35.1 PK (bottom right panel) for the simulation with the largest box size, run_0.8_320. In the lower right panel, sink particles are plotted as white dots and three pebble clusters are indicated using white circles. After self-gravity has been initiated at t = 25 PK, the pebbles aggregate into clusters and the axisymmetric filaments disperse (upper panels). When the self-gravity attains its full strength at t = 35 PK, most pebbles are concentrated in clusters and the filaments are no longer visible (lower left panel). At this point, we introduce the formation of sink particles, and the pebbles comprised in all but the three encircled clusters are replaced by sink particles. The sink particles emerging from the same pebble cluster undergo a merging process until only one of them remains. This process is largely completed at t = 35.1 PK (lower right panel).

Open with DEXTER
In the text
thumbnail Fig. 3

Radial locations x of sink particles which are not accreted before t = 36 PK, color-coded according to the locations at which they emerge, as functions of the simulation time t for the simulation with the largest box size, run_0.8_320. On average, the motions of the sink particles span more than half of the radial extent of the simulation box, a couple of them even migrate through the entire radial dimension of the box.

Open with DEXTER
In the text
thumbnail Fig. 4

Cumulative mass distributions of sink particles (black crosses) and fitted integrated power laws (Eq. (16), green lines) and power laws with exponential cutoff (Eq. (19), red lines) at t = 40 PK for the two simulations with the largest numbers of sink particles: the one with the largest box dimensions, run_0.8_320, (left panel) and the one with the middle box dimensions and the higher resolution, run_0.4_640 (right panel). Standard deviations of the actual N>/Ntot (black crosses) from the fitted N>/Ntot are given in the legends. The standard deviations for the integrated power laws are up to twice as large as those for the exponentially tapered power laws. In addition, it can be seen that the shallower exponential cutoffs fit the actual cutoffs more accurately than the steeper cutoffs of the integrated power laws and better replicate the smooth change of the slope of the mass distributions. We thus find the power laws with exponential tapering to represent the mass distributions better than the integrated power laws.

Open with DEXTER
In the text
thumbnail Fig. 5

Best-fitting parameters Mmin (upper left panel), α (upper right panel), Mexp (lower left panel), and β (lower right panel) of the exponentially tapered power law (Eq. (19)) at every Keplerian orbital period between t = 36 PK and t = 40 PK for all five simulations. Standard errors are plotted as error bars. In the case of the simulation with the smallest box size and the lower resolution, run_0.2_320, only parameter values for t = 36 PK are plotted because afterwards only four sink particles persist, which prevents us from properly fitting the mass distributions with a power law and an exponential tapering.

Open with DEXTER
In the text

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.