Initial mass function of planetesimals formed by the streaming instability
^{1} Hamburg Observatory, University of
Hamburg, Gojenbergsweg
112, 21029 Hamburg, Germany
email: urs.schaefer@hs.unihamburg.de
^{2} Lund Observatory, Department of
Astronomy and Theoretical Physics, Lund University, PO Box 43, 22100
Lund,
Sweden
Received:
20
August
2016
Accepted:
2
November
2016
The streaming instability is a mechanism to concentrate solid particles into overdense filaments that undergo gravitational collapse and form planetesimals. However, it remains unclear how the initial mass function of these planetesimals depends on the box dimensions of numerical simulations. To resolve this, we perform simulations of planetesimal formation with the largest box dimensions to date, allowing planetesimals to form simultaneously in multiple filaments that can only emerge within such large simulation boxes. In our simulations, planetesimals with sizes between 80 km and several hundred kilometers form. We find that a power law with a rather shallow exponential cutoff at the highmass end represents the cumulative birth mass function better than an integrated power law. The steepness of the exponential cutoff is largely independent of box dimensions and resolution, while the exponent of the power law is not constrained at the resolutions we employ. Moreover, we find that the characteristic mass scale of the exponential cutoff correlates with the mass budget in each filament. Together with previous studies of highresolution simulations with small box domains, our results therefore imply that the cumulative birth mass function of planetesimals is consistent with an exponentially tapered power law with a powerlaw exponent of approximately −1.6 and a steepness of the exponential cutoff in the range of 0.3–0.4.
Key words: hydrodynamics / instabilities / methods: numerical / planets and satellites: formation / protoplanetary disks
© ESO 2017
1. Introduction
One of the greatest problems in the theory of planet formation is to explain how millimeter or centimetersized solid particles – in the following referred to as pebbles – grow to kilometersized planetesimals. Micronsized 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 socalled bouncing barrier for very porous ice particles (Wada et al. 2008, 2009) or by mass transfer in highspeed 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 subKeplerian 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 metersized 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 solidtogas density ratio because of the drag force exerted by the particles on the gas. A locally enhanced solidtogas 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 solidtogas 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 wellfitted 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 powerlaw distribution that constitutes the physical upper mass cutoff, the latter do not include such a tapering in their fits. In this paper, we compare powerlaw fits with and without exponential cutoff to evaluate how well the highmass 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 selfgravity 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 selfgravity. 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
Simulation specifications.
We conduct threedimensional computer simulations with the Pencil Code^{1}, 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 sixthorder finite differences in space and thirdorder RungeKutta steps in time.
We use the shearing box approximation (Goldreich & LyndenBell 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 zdirections correspond to the radial, azimuthal, and vertical directions, respectively, and corotates 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, nonmagnetized gas with an isothermal equation of state , where p_{g} and ρ_{g} are the pressure and density, respectively, and c_{s} 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 midplane 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 H_{g} = c_{s}/ Ω_{K} is the gas scale height, Ω_{K} = 2π/P_{K} the Keplerian orbital frequency, and P_{K} the Keplerian orbital period. Here and in the following, the subscript zero refers to the midplane. 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 subKeplerian 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 subKeplerian orbital velocity of the gas is given by v_{g} = v_{K} − Δv, where v_{K} = Ω_{K}R is the Keplerian orbital velocity and Δv = Πc_{s} = 0.05c_{s}, 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: superparticles 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 superparticles and the mutual gravitational forces between the superparticles 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 backreaction drag forces and the selfgravitational forces onto the particles (Hockney & Eastwood 1981; Youdin & Johansen 2007; Johansen et al. 2007).
The gravitational potential of the superparticles 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 midplane such that their dynamics are not affected by the periodic potential away from the midplane.
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 midplane density ρ_{g,0}, the gas density perturbations are thus very small compared to the superparticle densities in the filaments when planetesimals form, which are of the order of the Roche density (see below).
The sedimentation of the superparticles to the midplane that is due to the vertical stellar gravity induces turbulence as a result of either the streaming instability or the KelvinHelmholtz instability, which are both caused by the mutual drag forces between the gas and the superparticles (Bai & Stone 2010b). This turbulence stirs up the superparticles, and hence counteracts the sedimentation. To give sedimentation and turbulence time to attain an equilibrium, in our simulations selfgravity is not introduced until t = 25 P_{K}. Its strength is then gradually increased from zero over 10 P_{K} until it reaches its final value at t = 35 P_{K} since initiating selfgravity instantaneously with full strength could cause significant impulses on the particles. While we initiate selfgravity at a simulation time at which the superparticles have already formed filaments, Johansen et al. (2015) introduced selfgravity 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 selfgravity by substituting (3)where the dimensionless selfgravity parameter (4)and G is the gravitational constant, into the righthand 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 selfgravity is introduced and the strength of the selfgravity. The Roche density depends on the selfgravity parameter γ and is given by (5)Each superparticle represents a large number of equally sized pebbles because it is computationally infeasible to simulate the pebbles individually. While the mass of a superparticle 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 superparticles is determined by the initial solidtogas column density ratio and their initial number. We set the solidtogas ratio (6)where Σ_{p,init} and (7)are the initial column densities of the superparticles and the gas, respectively, to Z = 0.02. This value corresponds to the critical solidtogas 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 superparticles is set equal to the total number of grid cells. The superparticles are randomly distributed among the entire simulation box to seed the streaming instability.
The Stokes number of the superparticles τ_{f} = Ω_{K}t_{f}, where t_{f} 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 particlegas dynamics are dominated by the most massive particles, and that the critical solidtogas ratio required for strong particle clustering owing to the streaming instability is determined by the total mass of all particles.
After selfgravity has attained its full strength at t = 35 P_{K}, every superparticle comprised in a cluster whose superparticle 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 superparticle 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 superparticles 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.
Superparticles 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 superparticle mass and momentum are added to the sink particle mass and momentum, respectively, and the superparticle 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 superparticles 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 superparticle accretion, and might as well be partially artificial. Nevertheless, it is required because in a superparticle cluster exceeding the sink particle creation threshold, all superparticles are replaced by sink particles. That is, the gravitational collapse of a superparticle cluster results in the creation of a cluster of sink particles, of which only one should persist to represent one newborn 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 P_{K}, the gas scale height H_{g}, and the midplane 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 selfgravity is initiated at t = 25 P_{K}. 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 M_{S} 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 P_{K} and H_{g}: 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 M_{S} = 1 M_{⊙}, and the chosen strength of the selfgravity γ = 1 /π to convert simulation units into physical units. For example, the midplane 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 selfgravity 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 P_{K}, when selfgravity is introduced (compare the evolution of the pebble density in the simulation time span 20 P_{K} ≤ t ≤ 25 P_{K} shown in Fig. 1 with their Fig. 3). We thus only report on the evolution of our simulations between t = 25 P_{K} and the end of the simulations, t = 40 P_{K}.
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 P_{K}, 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 selfgravity 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 selfgravity causes these filaments to fragment into pebble clusters that undergo gravitational collapse and form planetesimals. In the first 5 P_{K} after selfgravity 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 P_{K}, when the selfgravity 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 P_{K}, this sink particle merging process is for the most part completed (lower right panel). However, a few lowmass pebble clusters do not turn into sink particles. All sink particles emerge instantly at t = 35 P_{K}, as can be seen from Fig. 1, but a couple of clusters remain at t = 35.1 P_{K} (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.
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 P_{K} (top left panel), t = 30 P_{K} (top right panel), t = 35 P_{K} (bottom left panel), and t = 35.1 P_{K} (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 selfgravity has been initiated at t = 25 P_{K}, the pebbles aggregate into clusters and the axisymmetric filaments disperse (upper panels). When the selfgravity attains its full strength at t = 35 P_{K}, 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 P_{K} (lower right panel). 

Open with DEXTER 
3.2. Mutual sink particle accretion
For the sink particles to realistically represent newborn 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 superparticle 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 P_{K} after their formation, i.e. until t = 35.1 P_{K}, because in all simulations all sink particles emerge togather at t = 35 P_{K}. 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 m_{i} are the masses of the two sink particles, r_{i} their radii, which are calculated from m_{i} 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 ⟨b_{max}⟩, averaged over all mutual sink particle accretions in all five simulations where the lifetime of the accreted sink particle is greater than 0.1 P_{K} 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 P_{K}. 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. Bestfitting parameters
Fig. 3 Radial locations x of sink particles which are not accreted before t = 36 P_{K}, colorcoded 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, N_{tot} is their total number and M_{max} their maximum mass, and M_{pow} and α are fitting parameters. Since the formation of planetesimals by the streaming instability is a stochastic process and the actual M_{max} in a simulation might differ significantly from the M_{max} of the ensembleaveraged mass distribution of the sink particles for this model, we treat M_{max} as a fitting parameter.
The exponentially tapered power law is given by (17)where M_{exp} and β are fitting parameters. The condition N_{>}(M_{min}) /N_{tot} = 1, where M_{min} 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 powerlaw part, M_{pow}, 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 M_{max} in Eq. (16), M_{min} 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 leastsquares method to fit Eqs. (16) and (19) to the sink particle data at each Keplerian orbital period between t = 36 P_{K} and the end of the simulations, t = 40 P_{K}. In all five simulations, all sink particles emerge at once at the simulation time at which their formation is initiated, t = 35 P_{K}, but we begin the fitting at t = 36 P_{K} 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 P_{K}. In the legends for both fits, the standard deviation σ of the actual N_{>}(M) /N_{tot} for the sink particle masses M from the fitted N_{>}(M) /N_{tot} 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 powerlaw fits without exponential cutoff are larger, not only at t = 40 P_{K} as shown in the figure, but also at other simulation times. We further note that the shallower exponential cutoff represents the highmass 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 powerlaw fits with exponential cutoff.
4.2. Powerlaw fits with and without exponential tapering
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 P_{K} 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_{>}/N_{tot} (black crosses) from the fitted N_{>}/N_{tot} 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 bestfitting parameters M_{min}, α, M_{exp}, and β in Eq. (19) are shown for t = 36 P_{K} to t = 40 P_{K} 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 P_{K} does not permit to properly fit a power law and an exponential cutoff. The parameters α, M_{exp}, and β are approximately constant in time for all simulations apart from M_{exp} 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 P_{K} to t = 40 P_{K}.
Similar to M_{exp} for run_0.4_320, the fitted minimum mass M_{min} 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 M_{min} and M_{exp} 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, M_{min} is roughly constant. To provide a comparable value for each simulation despite these differences in the time dependence, in Table 2 we list the bestfitting values at t = 40 P_{K} for all simulations (except for run_0.2_320, for which only the value at t = 36 P_{K} is available). A comparison of these values shows that M_{min} declines with both increasing box dimensions and increasing resolution if the value for run_0.2_320 is disregarded.
The bestfitting values of the exponent of the powerlaw 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 lowmass sink particles that would constitute the powerlaw part to emerge. Both Johansen et al. (2015) and Simon et al. (2016) fit the differential mass distribution with a powerlaw 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 powerlaw distribution higher resolutions than seem to be required.
We find the mass scale of the exponential cutoff M_{exp} to correlate with the mass budget in every filament. The parameter M_{exp} 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 N_{f} 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 H_{g}, 0.4 H_{g}, and 0.8 H_{g}, 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 M_{exp} 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 bestfitting values for run_0.4_320 increase with time with a range enclosing the mean value for run_0.4_640. Hence, we find M_{exp} to be largely independent of the resolution.
Bestfitting 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 bestfitting 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 highmass 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 bestfitting 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 R_{min} and the radius scale of the exponential tapering R_{exp} can be calculated from M_{min} and M_{exp}, respectively, using a solid density of 3 g cm^{3}.
Sink particle statistics.
In Table 3, we list the number N_{tot}, minimum mass M_{min}, maximum mass M_{max}, and mean mass ⟨M⟩ of the sink particles at the end of our five simulations, t = 40 P_{K}. In the case of the two simulations with the smallest box dimensions, we observe less than ten sink particles, and the bestfitting parameters for these simulations are therefore afflicted with comparably large errors (see Fig. 5). We find the actual M_{min} and ⟨M⟩ to be of the order of the fitted M_{min} and M_{exp}, 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, N_{tot} increases and both M_{min} and ⟨M⟩ decline with increasing resolution. Like Johansen et al. (2015) and Simon et al. (2016), we find the maximum mass M_{max} to be relatively independent of the resolution.
Fig. 5 Bestfitting parameters M_{min} (upper left panel), α (upper right panel), M_{exp} (lower left panel), and β (lower right panel) of the exponentially tapered power law (Eq. (19)) at every Keplerian orbital period between t = 36 P_{K} and t = 40 P_{K} 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 P_{K} 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 M_{exp}, 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, M_{min} and M_{max} in general increase with the box dimensions, but we find M_{min} to considerably decrease if the box size is increased from 0.2 H_{g} to 0.4 H_{g} 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 ensembleaveraged values of M_{min} and M_{max} 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 selfgravity, the streaming instability concentrates pebbles into axisymmetric filaments. After selfgravity 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 selfgravity and thus on the solid particle column density, range from 80 km to 620 km. We have compared powerlaw 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 powerlaw part and to investigate its dependence on the box dimensions because the planetesimals that formed in our simulations are too large to constitute a powerlaw distribution at the lowmass 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 wellfitted 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 shortlived, 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 selfgravity, and the simulation at which selfgravity is initiated. We have complemented these parameter studies with an analysis of the boxsize dependence, but how, for instance, the solidtogas 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 278675PEBBLE2PLANET. A.J. is thankful for financial support from the Knut and Alice Wallenberg Foundation and from the Swedish Research Council (grant 20145775).
References
 Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Prog. Theor. Phys., 56, 1756 [NASA ADS] [CrossRef] [Google Scholar]
 Armitage, P. J. 2007, ArXiv Astrophysics eprints [arXiv:astroph/0701485] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010a, ApJS, 190, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010b, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010c, ApJ, 722, L220 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bottke, W. F., Durda, D. D., Nesvorný, D., et al. 2005, Icarus, 175, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Nordlund, Å., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Brauer, F., Dullemond, C. P., Johansen, A., et al. 2007, A&A, 469, 1169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chabrier, G. 2003, PASP, 115, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Gammie, C. F. 2001, ApJ, 553, 174 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & LyndenBell, D. 1965, MNRAS, 130, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1981, Computer Simulation Using Particles (New York: McGrawHill) [Google Scholar]
 Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Johansen, A., Youdin, A., & Klahr, H. 2009a, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A., & Mac Low, M.M. 2009b, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Klahr, H., & Henning, T. 2011, A&A, 529, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, A&A, 537, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Mac Low, M.M., Lacerda, P., & Bizzarro, M. 2015, Sci. Adv., 1, 1500109 [NASA ADS] [CrossRef] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Massey, P. 1998, in The Stellar Initial Mass Function, 38th Herstmonceux Conference, eds. G. Gilmore, & D. Howell, ASP Conf. Ser., 142, 17 [Google Scholar]
 Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salpeter, E. E. 1955, ApJ, 121, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2008, ApJ, 677, 1296 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Windmark, F., Birnstiel, T., Güttler, C., et al. 2012, A&A, 540, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wurm, G., Paraskov, G., & Krauss, O. 2005, Icarus, 178, 253 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Yang, C.C., & Johansen, A. 2014, ApJ, 792, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, C.C., Johansen, A., & Carrera, D. 2016, A&A, submitted [arXiv:1611.07014] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Zinnecker, H., & Yorke, H. W. 2007, ARA&A, 45, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
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 P_{K}, 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 
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 P_{K} (top left panel), t = 30 P_{K} (top right panel), t = 35 P_{K} (bottom left panel), and t = 35.1 P_{K} (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 selfgravity has been initiated at t = 25 P_{K}, the pebbles aggregate into clusters and the axisymmetric filaments disperse (upper panels). When the selfgravity attains its full strength at t = 35 P_{K}, 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 P_{K} (lower right panel). 

Open with DEXTER  
In the text 
Fig. 3 Radial locations x of sink particles which are not accreted before t = 36 P_{K}, colorcoded 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 
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 P_{K} 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_{>}/N_{tot} (black crosses) from the fitted N_{>}/N_{tot} 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 
Fig. 5 Bestfitting parameters M_{min} (upper left panel), α (upper right panel), M_{exp} (lower left panel), and β (lower right panel) of the exponentially tapered power law (Eq. (19)) at every Keplerian orbital period between t = 36 P_{K} and t = 40 P_{K} 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 P_{K} 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 