Issue 
A&A
Volume 567, July 2014



Article Number  A38  
Number of page(s)  8  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201423708  
Published online  07 July 2014 
Modeling dust growth in protoplanetary disks: The breakthrough case
Heidelberg University, Center for Astronomy, Institute of Theoretical
Astrophysics,
AlbertUeberleStr. 2,
69120
Heidelberg,
Germany
email:
drazkowska@uniheidelberg.de
Received: 25 February 2014
Accepted: 20 May 2014
Context. Dust coagulation in protoplanetary disks is one of the initial steps toward planet formation. Simple toy models are often not sufficient to cover the complexity of the coagulation process, and a number of numerical approaches are therefore used, among which integration of the Smoluchowski equation and various versions of the Monte Carlo algorithm are the most popular.
Aims. Recent progress in understanding the processes involved in dust coagulation have caused a need for benchmarking and comparison of various physical aspects of the coagulation process. In this paper, we directly compare the Smoluchowski and Monte Carlo approaches to show their advantages and disadvantages.
Methods. We focus on the mechanism of planetesimal formation via sweepup growth, which is a new and important aspect of the current planet formation theory. We use realistic test cases that implement a distribution in dust collision velocities. This allows a single collision between two grains to have a wide range of possible outcomes but also requires a very high numerical accuracy.
Results. For most coagulation problems, we find a general agreement between the two approaches. However, for the sweepup growth driven by the “lucky” breakthrough mechanism, the methods exhibit very different resolution dependencies. With too few mass bins, the Smoluchowski algorithm tends to overestimate the growth rate and the probability of breakthrough. The Monte Carlo method is less dependent on the number of particles in the growth timescale aspect but tends to underestimate the breakthrough chance due to its limited dynamic mass range.
Conclusions. We find that the Smoluchowski approach, which is generally better for the breakthrough studies, is sensitive to low mass resolutions in the highmass, lownumber tail that is important in this scenario. To study the low number density features, a new modulation function has to be introduced to the interaction probabilities. As the minimum resolution needed for breakthrough studies depends strongly on setup, verification has to be performed on a case by case basis.
Key words: accretion, accretion disks / circumstellar matter / protoplanetary disks / planets and satellites: formation / methods: numerical
© ESO, 2014
1. Introduction
Dust coagulation is an important astrophysical process, particularly in the case of planet formation. To form a planet, micronsized dust particles have to collide and stick to form ever larger dust aggregates. The main problem with this picture is what is known as the “metersize barrier”: once particles reach a size of decimeter to meter, they acquire relative velocities that are too large to allow hitandstick growth. Instead, collisions at these speeds, which are on the order of tens of meters per second, lead to the destruction of the colliding aggregates. This is often called the “fragmentation barrier”. Moreover, already at sizes of a millimeter, we encounter the “nonsticking problem”, which is often called the “bouncing barrier”: while the aggregates do not fragment upon collision at these sizes, they do not stick together either, meaning that the growth stalls. All these barriers are posing a major problem in our understanding of planet formation. In recent years, however, several new ideas have appeared in the literature to overcome these problems.
For example, the sweepup growth scenario has been proposed as a solution for the growth barriers issue (Windmark et al. 2012a). In this scenario, a few boulders can grow large by colliding with numerous smaller pebbles through the fragmentation with mass transfer process (Wurm et al. 2005; Teiser & Wurm 2009; Meisner et al. 2013). The growth of the pebbles is suppressed due to bouncing and/or fragmentation between themselves, while a few larger aggregates, or seeds, might form if the distribution of impact velocities due to stochastic turbulence is taken into account (Windmark et al. 2012b; Garaud et al. 2013; Meru et al. 2013a). Thanks to the low velocity tail of the distribution, some grains can be “lucky” enough to not experience any destructive collisions and undergo only lowvelocity sticking collisions, breaking through the growth barriers. This scenario enables the formation of planetary embryos while still keeping the disk dusty, which is consistent with observations.
Key to the trustworthiness of the conclusions derived from numerical models is the reliability of the codes and algorithms used. The problem of coagulation is extremely complex and nonlinear, and with the exception of some very simple coagulation kernels, no analytic solutions exist. So how do we know if the results of our codes are indeed correct? One way is to treat the problem with at least two distinct methods and compare the results.
Over the years, different approaches have been developed to study the dust coagulation problem. Besides numerous semianalytic models, two main numerical approaches are used nowadays: direct numerical integration of the Smoluchowski equation and various Monte Carlo codes. The former is traditional approach, which has been used in different versions by Weidenschilling (1980); Nakagawa et al. (1981); Wetherill & Stewart (1989); Ohtsuki et al. (1990); Dullemond & Dominik (2005); Tanaka et al. (2005); Nomura & Nakagawa (2006); Brauer et al. (2008); Okuzumi et al. (2009); Birnstiel et al. (2010); Charnoz & Taillifet (2012) and many others. This approach is often used when comparing dust coagulation models to observations, as it allows us to model the dust evolution in the global disk over very long timescales. The Monte Carlo approach is based on work by Gillespie (1975) and is used in one form by Ormel et al. (2007) and in another by Zsom & Dullemond (2008), subsequently used by Johansen et al. (2012); Ros & Johansen (2013) and Drążkowska et al. (2013). The Monte Carlo approach is useful to test different coagulation models and to include different properties of dust particles, such as the internal grain structure. It is also better suited to use along with hydrodynamic grid codes.
The two methods are usually benchmarked using analytical solutions of the coagulation equation that are available for three idealized growth kernels (see, e.g., Ohtsuki et al. 1990; Wetherill 1990; Lee 2000). However, these kernels do not necessarily represent any realistic growth scenario in the protoplanetary disk. In this work, we perform an explicit comparison between the two approaches for the first time. In this comparison, we focus on the sweepup growth scenario, which is challenging to model for both of the methods. In particular, it was already asserted by Windmark et al. (2012a,b) that an artificial breakthrough may occur when a low mass resolution is used in the Smoluchowski method. We study this issue in more detail, and we show that not only the high resolution but also a careful treatment of interactions in low particle number density bins is needed to avoid the nonphysical growth.
Until now, the sweepup growth triggered by the “lucky” growth was modeled using the Smoluchowski code only. Using a twodimensional Monte Carlo code, Drążkowska et al. (2013) showed that sweepup growth can occur at the inner edge of dead zone, but it was triggered by radial transport of big bodies grown in a dead zone in this case. In this paper, we implement the relative velocity distribution in the Monte Carlo code and directly compare results of the two approaches. We also present some of their major features and differences.
This paper is organized as follows: we describe both of our numerical models in Sect. 2. In Sect. 3, we compare results obtained with both codes. We discuss issues related to numerical convergence of both methods in Sect. 4. We summarize our findings in Sect. 5.
2. The numerical models
We study the Smoluchowski approach using the code developed by Brauer et al. (2008) and Birnstiel et al. (2010), along with an impact velocity distribution implemented as described in Windmark et al. (2012b). In this code, we let the dustgrain number density n(m,r,z) be a function of the grain mass m, the distance to the star r, and the height above the midplane z, and give it in number of particles per unit volume per unit mass. The dust evolution can then be solved by integration.
However, discretization of the problem is necessary in the integration process, which can lead to a significant numerical diffusion in massspace, because having a finite number of grid points means that particle collisions do not necessarily lead to particle masses m_{p} that directly correspond to one of the logarithmically spaced sampling points. The approach by Brauer et al. (2008) was to implement an algorithm that distributes the mass of the resulting particle into two adjacent mass bins corresponding to grid points i and i + 1, m_{i} < m_{p} < m_{i + 1}, according to (1)where ϵ·m_{p} is put into mass bin m_{i + 1}, and (1 − ϵ)·m_{p} is put into mass bin m_{i}. This algorithm is based on the work of Kovetz & Olund (1969) and is adopted by most of the modern dust coagulation codes. This approach, however, leads to some numerical diffusion, as m_{i+1} > m_{p} means that mass is inserted into a mass bin that corresponds to a larger mass than the mass of physical particle created in the collision. If the spacing between the mass bins is too coarse, this leads to a significant, artificial growth rate speedup (Ohtsuki et al. 1990). We show that the same effect strongly affects the number of seeds formed in the “lucky growth” scenario; however, the problem is more severe and requires a careful approach to low number density regions in this case.
In the Monte Carlo method, we use the representative particles approach described by Zsom & Dullemond (2008) and implemented by Drążkowska et al. (2013). We have slightly modified the method to match the vertical treatment of the Smoluchowski code (described by Brauer et al. 2008). Instead of performing an explicit vertical advection, we redistribute the particles according to a Gaussian distribution with a width: (2)where H_{gas} is the pressure scale height of the gas, St is the particle’s Stokes number and α is the turbulence strength parameter. In this way, we account for the reduction in the collision rate between small and large grains due to their different vertical settling. This occurs because small particles are more strongly affected by the turbulent diffusion, and, thus, their density in the midplane of the disk is lower than in the case of large particles.
One of the assumptions of the representative particle approach that we implement is that one representative particle, or swarm, represents a constant amount of mass, which is equal to (3)where M_{tot} is the total mass of dust present in the computational domain and N_{swarms} is the number of swarms used. In other words, the total mass of dust M_{tot} is divided into N_{swarms} equalmass units. Each of these units represents physical particles of mass m_{i} (m_{i} ≪ M_{swarm}), but the number of these particles N_{i} has to be such that N_{i}·m_{i} = M_{swarm}. The algorithm fails, if, for example, a physical coagulation kernel would lead to formation of only one massive particle with mass m_{p}, while keeping all the other particles small. If m_{p} < M_{swarm}, the single big particle cannot be resolved, because it involves less mass than the smallest available unit M_{swarm}. Increasing the number of swarms, N_{swarms}, lowers M_{swarm} and thus improves the mass resolution of the method; however, the computation time increases quadratically with the number of swarms used. The dynamic mass range of the representative particle approach is limited by the number of swarms used. This issue can be overcome by implementing a more advanced algorithm, as the “distribution method” proposed by Ormel & Spaans (2008) that involves continuously adjusting M_{swarm} by splitting and merging the swarms. This method was later used in the context of accretion among planetesimals (Ormel et al. 2010), allowing the Monte Carlo method to resolve a runaway coagulation kernel. However, the method was not yet tested with velocity distributions or complicated collision models that are needed to break through the growth barriers.
3. Comparison of results obtained with both methods
Because our Monte Carlo method is not capable of the same dynamic mass range as the Smoluchowski method, to directly compare the two methods, we choose a setup where the particle breakthrough (i.e. where the “lucky” particles can start to grow by mass transfer) occurs at relatively high particle number density, which is possible to resolve with our representative particle approach (see the previous section).
For the disk model, we use the minimummass extrasolar nebula (Chiang & Laughlin 2013) at 1 AU. The gas surface density is Σ_{gas} = 9900 g cm^{2}, and the temperature is T = 280 K. We assume a turbulence of α = 10^{2} and a standard dust to gas ratio of 10^{2}. We take a relative velocity between the dust grains driven by Brownian motion and turbulence into account, by calculating the rootmeansquare impact velocity v_{rms} from the formulas derived by Ormel & Cuzzi (2007), and we assume a Maxwellian distribution of the impact velocity.
We consider sticking, fragmentation, and mass transfer as possible outcomes of collision, which we refer to as the SF+MT model (Güttler et al. 2010; Windmark et al. 2012b). The collision outcome is determined by taking the sticking and fragmentation/mass transfer probabilities P_{s}(v_{rms}) and P_{f/mt}(v_{rms}) into account, which are calculated analogically as in Windmark et al. (2012b). We take the fragmentation threshold velocity to be v_{f} = 50 cm s^{1}. If a collision should lead to fragmentation but the mass ratio between the colliding particles mass is m_{1}/m_{2}> 20 (m_{1}>m_{2}), we assume projectile fragmentation with mass transfer with a 10% efficiency; that is the larger particle gains 10% of mass of the smaller one during the event. Realistic values of the collision parameters are poorly constrained, as discussed in Sect. 5.
Implementing the same setup in both of the codes, we perform a number of runs by varying the numerical resolution by around one order of magnitude. In the Smoluchowski code, we use from 3 to 40 mass bins per decade. In the Monte Carlo code, we use from 12 000 to 120 000 representative particles, and we repeat each run ten times with different random seeds.
The Monte Carlo method relies on random numbers used to determine which particles are participating in the subsequent collisions and to calculate collision time steps (Zsom & Dullemond 2008). Thus, outcomes of the Monte Carlo runs performed with different random seeds vary despite using the same setup. In the velocity distribution case, this effect is even stronger, meaning that multiple runs with different random seeds are necessary. This causes our high resolution Monte Carlo models to need a few days on an 8 core 3.1 GHz AMD machine. For comparison, the Smoluchowski models with 40 bins per mass decade take about one hour on a single core processor.
Fig. 1 Comparison of the mass distribution evolution obtained using our Smoluchowski and Monte Carlo codes. In a standard case when an insurmountable fragmentation barrier is present, the two methods perfectly agree (three upper panels); however, we encounter some differences between the results obtained with the two approaches (three bottom panels) with the possibility of breakthrough. For the Smoluchowski code, the presented results were obtained with a resolution of 30 mass bins per decade. For the Monte Carlo code, the results are averaged from the simulations using 120 000 particles. The “check points”, m_{1} and m_{2}, which are indicated with the dotted lines, are used to quantitatively compare dust growth timescales in Sect. 4.1. 

Open with DEXTER 
Figure 1 shows the time evolution of the dust mass distribution obtained with both of the codes. For the Smoluchowski code, we display results obtained with the resolution of 30 bins per mass decade. For the Monte Carlo code, we display averaged results, along with their scatter, of our highest resolution simulations with 120 000 particles, which additionally showed a breakthrough before 100 yrs; otherwise, the figure gets indecipherable due to high noise.
The early evolution, as seen in the upper panels of Fig. 1, is very similar in both codes, and the point at which the dust distribution hits the fragmentation barrier, where the growth of the majority of particles is hindered due to fragmentation that occurs when the impact velocities exceed the fragmentation threshold velocity v_{f}, is identical (m ≅ 10^{2} g, t ≅ 20 yrs). If the fragmentation with mass transfer collisions are not included and there is no possibility of growth beyond the fragmentation barrier, the steady state is represented by the third panel of Fig. 1, and the two methods agree perfectly. However, as we include the breakthrough possibility, we encounter some differences between the two approaches, which are visible at later stages of the evolution. Both of the codes reveal that some of the particles can grow beyond the fragmentation barrier thanks to the lowvelocity sticking collisions. The growth is very quick in general and metersized bodies are formed within <1000 yrs. However, the breakthrough generally occurs later in the Monte Carlo code, and the population of big grains lacks the characteristic waves seen in the distribution obtained from the Smoluchowski code, which can be seen in the bottom panel of Fig. 1. The differences in the late stages of evolution are caused by the restricted dynamic mass range of our representative particle approach that was discussed in Sect. 2. We further discuss issues related to resolution dependencies of both methods in the following section.
4. Resolution dependence
The results obtained with both codes exhibit resolution dependence. In this section, we discuss specific issues connected to the numerical convergence of both Smoluchowski and Monte Carlo methods.
4.1. Growth timescales
To quantitatively compare dust growth obtained in both codes, we establish “check points” m_{1} and m_{2}, as marked with the dotted lines in Fig. 1. We arbitrarily choose m_{1} = 2 × 10^{5} g and m_{2} = 1 g, and mark the time at which the peak of the mass distribution reaches a mass corresponding to one of the check points.
Figures 2 and 3 show the results obtained for m_{1} and m_{2}, respectively. We show the times of crossing the “check points” as a function of the number of bins used in the Smoluchowski code and a number of representative particles used in the Monte Carlo method. We note that the scaling of xaxes of these figures is arbitrary, as there is no method to directly connect the number of bins in the Smoluchowski method to the number of representative particles in the Monte Carlo method.
Mass m_{1} is located before the fragmentation barrier. This part of the evolution corresponds to a standard growth scenario that does not include the “lucky” breakthrough possibility and is well resolved by both of the codes. In the Monte Carlo code case, the time of crossing m_{1} depends very weakly on the number of particles used. The difference between individual runs, as marked by the shaded region in Fig. 2, is also very low. The Smoluchowski code exhibits much stronger resolution dependence. In the case of our lowest resolution of 3 bins per mass decade, the growth is more than four times faster than with 40 bins. The results obtained with resolution of 40 bins per mass decade converge to a value that is consistent with the one given by the Monte Carlo code. The critical resolution agrees with findings of Lee (2000) and Okuzumi et al. (2009).
Fig. 2 Time at which the peak of mass distribution reaches the mass m_{1} = 2 × 10^{5} g for both the Smoluchowski and Monte Carlo code. The result changes with mass resolution, given by the number of bins per mass decade in the case of Smoluchowski code and the number of particles used in the Monte Carlo code. The scatter of results obtained in different runs with the same number of particles in the Monte Carlo code is marked by the shaded region around an averaged dependence. The Monte Carlo approach does not exhibit a strong resolution dependence. In contrast, the Smoluchowski algorithm overestimates the growth rate by a factor of few when we do not use enough mass bins. 

Open with DEXTER 
Fig. 3 Same as Fig. 2 but for the “check point” m_{2} = 1 g. The resolution dependence of the Smoluchowski method is the same as for m_{1}, but it changes significantly in the case of the Monte Carlo code. The Monte Carlo code tends to underestimate the breakthrough possibility when used with a low number of particles due to the limited dynamic mass range. 

Open with DEXTER 
The reason for the high diffusion of the Smoluchowski method is the linearity of the algorithm described in Eq. (1), which is a necessity for solvers that employ implicit integration schemes to overcome the numerical stiffness of the equations (see Brauer et al. 2008, for further details). Explicit solvers, on the other hand, are capable of implementing higher order mass distribution schemes, which lowers the numerical diffusion. We find that steady states that arise when a bouncing or fragmentation barrier is met and that have no possibility of breaking through are significantly less dependent on mass resolution. Thus, the conclusions of most of the papers that did not include breakthrough are not affected, even if a lower resolution was used.
As can be seen in the Fig. 3, the situation changes significantly in the case of the second check point m_{2}. This point is located beyond the fragmentation barrier and breakthrough point. Due to limited dynamic mass range, the resolution dependence of the Monte Carlo code is much stronger. The dispersion of results obtained in individual runs is much higher and can reach even an order of magnitude. Generally, the less particles we use, the longer we have to wait for the breakthrough. In the case of the Smoluchowski code, the resolution dependence is the same as for m_{1}, meaning that the lower resolution we use, the faster the breakthrough occurs. This is consistent with the findings of Garaud et al. (2013). The results obtained with both of the codes roughly converge for our highest accuracy. Both of the methods become computationally inefficient when used with even higher resolutions.
The different resolution dependence seen in the Fig. 3 is a result of a fundamental difference between the two approaches. In a real protoplanetary disk, the number of physical particles is so high that even when the breakthrough probability is low, some particles will be able to be “lucky” and overcome the fragmentation barrier quickly. However, the number of representative particles is restricted in the Monte Carlo code and the breakthrough probability is additionally reduced. In contrast, the breakthrough is resolved much easier because the Smoluchowski code deals with number densities instead of discrete particles. However, this introduces another problem, which is discussed in the following section.
4.2. Breakthrough probability
The “lucky growth” scenario has introduced a new issue in how the numerical diffusion affects the global dust evolution in the Smoluchowski method. In this section, we discuss this issue and introduce a way to limit its effect by including a modulation function to the coagulation algorithm that suppresses the interactions with mass bins containing unrealistically low particle numbers.
As discussed above, when velocity distributions are introduced, the collision barriers are naturally smeared out. Because a particle that is more massive than the grains at the mass distribution peak must be “lucky” and has to grow by only interacting with other particles in lowvelocity collisions, it becomes necessary to accurately resolve the highmass tail of the distribution. Otherwise, if not all sticking events are resolved properly, the slope of the tail becomes incorrect, creating artificially large mass ratios between the luckiest grains and those in the peak.
As an example of this, we can consider the extreme case where the entire highmass tail is represented by a single mass bin m_{i}. If two grains in the peak undergo a single sticking event, forming a particle of mass m ≪ m_{i + 1}, some mass will still be put into the mass bin i + 1, even though the particles would need to undergo several consecutive sticking events to reach a mass m_{i + 1} in reality. Such a badly resolved largeparticle tail could cause an artificial breakthrough of growth barriers, as unrealistically large particles can form that continue to grow by the sweepup process where no such particles would form in a better resolved case.
To show this issue clearly, we perform additional simulations with the Smoluchowski code using a critical mass transfer ratio, that is the mass ratio above which a fragmenting collision leads to mass transfer (see Sect. 3), of 500. The point of breakthrough then occurs at a very low number density, which is impossible to resolve with our implementation of the Monte Carlo method, as the particles that should break through would involve less mass than the mass of a single swarm in our simulation (see also Sect. 2). We have therefore performed a resolution study with the Smoluchowski code only. Simulations were run with resolutions between 3 and 40 bins per decade of mass with the results shown in Fig. 4. As can be seen, the SF+MT case is extremely sensitive to the mass resolution. Between the highest and the lowest resolutions, the point of breakthrough differs by more than 25 orders of magnitude in surface density. For resolutions above 30 bins per decade, the point of breakthrough would have occurred at a density lower than the density corresponding to 1 particle in an annulus of 0.1 AU width.
Fig. 4 Effect of mass resolution in the case of SF+MT setup with a critical mass ratio of 500 for the Smoluchowski code. The resolutions included a range from 3 to 40, and the mass distribution is shown after 100 yrs of evolution. The point of breakthrough is very sensitive to the resolution used, and the breakthrough happens only in the low resolution cases, which is clearly nonphysical. The dashed lines mark the density corresponding to 1 and 10^{6} particles of a given mass in an annulus of 0.1 AU width. 

Open with DEXTER 
To correctly simulate the lowest dust densities, it becomes necessary to include a modulation function, f_{mod}, that limits the interactions between mass bins that have particle numbers that are too low. The reasoning behind this is the following. For breakthrough to occur, we need at least 1 real particle within the simulation domain that is large enough to trigger a sweepup. Because the Smoluchowski code deals with number densities in a single point in space, the limiting particle density becomes somewhat arbitrary and relies on choosing a reasonable physical domain size.
The modulation function works in a similar way to the active bin method introduced by Lee (2000). Instead of speeding up the code by deactivating lowpopulated mass bins; however, it is a continuous method that prevents growth of bins with unrealistically low particle numbers. The collision frequency between particle species i and j can then be written as (4)where n_{i} and n_{j} are the number densities, and K_{ij} is the coagulation kernel. We choose f_{mod} = exp(−1/N_{i} − 1 /N_{j}), where N_{i} and N_{j} are numbers of particles in a given surface of the disk, which we take to be an 0.1 AU wide annulus at 1 AU. The result of this is that mass is allowed to be put into a mass bin with a nonunity amount of particles by coagulation, for example, but the mass inside is unable to coagulate further until the bin contains sufficient mass.
For our setup, the modulation function f_{mod} suppresses the breakthrough in the case of resolution larger than 30 bins per mass decade. The breakthrough that occurs in the runs with lower resolution is a result of the numerical diffusion introduced by the mass distribution algorithm (Eq. (1)), which is discussed in Sect. 2. This is seen by the sharp cutoff at masses between 10^{2}−10^{1} g for the highest resolutions. We stress that this particular resolution dependency varies very strongly between setups, and it is necessary to confirm numerical convergence individually.
We want to note that even a high resolution run can lead to a nonphysical breakthrough if f_{mod} is not included. In such a case, breakthrough would initially occur at densities corresponding to less than one particle in an annulus of 0.1 AU width, but the density of the highmass tail could increase over the threshold density during further evolution. Thus, the combination of sufficient numerical resolution and the modulation function is necessary to ultimately confirm the possibility of forming planetesimals in the breakthrough scenario under given conditions.
In this section, we clarified the issue of too crude mass resolutions found by Windmark et al. (2012a,b). Although the lower resolutions commonly used in previous papers might work well for less extreme cases, which did not include the possibility of breakthrough, an artificial growth might occur for problems of the kind discussed here. Thus, we stress that careful convergence tests are necessary to confirm breakthrough possibility. We find that the convergence depends strongly on the strength of the growth barrier, and therefore, a separate convergence study is needed for each setup.
As was shown in the previous section, the Monte Carlo code underestimates the breakthrough chance with low particle numbers in contrast to the Smoluchowski code, and, thus, the time at which breakthrough happens increases. If the Monte Carlo code is used with not sufficient number of particles, the breakthrough might be completely suppressed. This happens when the mass that should be involved in the breakthrough is lower than M_{swarm}, the mass represented by a single swarm, that limits the dynamic range of the Monte Carlo method (see Sect. 2). The minimum number of representative particles required is dependent on the breakthrough mechanism. Here, the breakthrough is driven by distribution in the impact velocities, which smear out the fragmentation barrier by changing the slope of the mass distribution highmass end. The breakthrough is only possible if the largest particles are more massive than the particles in the mass distribution peak by a factor defined by the critical mass ratio. Thus, the number of breakthrough particles is defined by the mass distribution function slope and the value of the critical mass transfer ratio. The higher the slope and the critical mass ratio, the lower number of particles can break through. However, it was shown by Drążkowska et al. (2013) that the breakthrough can also be driven by radial mixing of dust aggregates between regions of different grain sizes. In such case, the number of breakthrough particles follows different constraints than here.
5. Discussion and conclusions
The issue of numerical resolution of the Smoluchowski code has been discussed for a long time in the context of dust coagulation in protoplanetary disks (Ohtsuki et al. 1990; Wetherill 1990; Lee 2000). These codes used different algorithms that do not necessarily result in similar resolution dependencies as the implicit integration scheme with the linear mass distribution algorithm we use.
In this paper, we extended the prior studies by implementing a possibility of breaking through growth barriers with impact velocity distributions and including a direct comparison to the Monte Carlo algorithm with representative particle approach. The two methods have never before been explicitly compared. Our work showed that the two methods give consistent results when applied to usual coagulation problems. However, we find that modeling of the recently discovered planetesimal formation via “lucky growth” is much more challenging. Although the results obtained with the two methods converge for sufficient resolution, the approaches are fundamentally different and their limitations have to be realized when performing scientific models.
In agreement with previous studies, we find that simulations with the Smoluchowski code require a sufficiently high mass resolution to avoid an artificial speedup of the growth rate. This problem arises from numerical diffusion and our implementation of the algorithm, which is required to describe how the resulting mass is distributed after a collision event. Additionally, we show that the introduction of the modulation function that prevents interactions between mass bins containing less than one physical particle is necessary to study the dust coagulation at low number densities. In Sect. 4.2, we have shown that the numerical issues can change both the quantitative and the qualitative result in the case of the breakthrough scenario.
The Monte Carlo approach used to study the breakthrough scenario, in which only a few “lucky” particles break through the growth barriers, results in a high noise. In the presented tests, the individual runs show times of breakthrough that are orders of magnitude different and their averaged value depends on mass resolution. However, contrary to the Smoluchowski approach, the Monte Carlo approach does not suffer a strong resolution dependence when the dust aggregates grow with a relatively narrow size distribution, which usually is the case for small aggregate evolution when no breakthrough or runaway growth are possible. The Monte Carlo methods are generally computationally expensive, and they require numerous runs with different random seeds to reduce the noise.
The convergence of each of the methods can be very different for every setup. We cannot present a general recipe for the minimum resolution required to study dust growth, because this can vary enormously from case to case. Thus, it is important to run resolution tests for every new physical model until convergence of results is obtained.
Since we find that the numerical convergence is sensitive to the collision model parameters, it is important to use realistic values; however, these are poorly constrained, as the amount of data we have from laboratory experiments is restricted, and it is still not necessarily reproduced by direct numerical simulations. Laboratory experiments show that bouncing collisions start to occur at velocities between 0.1 and 10 cm s^{1} (Kothe et al. 2013), while numerical simulations claim that such collisions only rarely occur, if ever (Wada et al. 2011; Seizinger & Kley 2013). Fragmentation is also hotly discussed, as laboratory experiments find fragmentation to occur at velocities as low as a few cm s^{1} between 5 cm grains (Schräpler et al. 2012) and up to a few m s^{1} for mmsized grains (Lammel 2008). Numerical simulations, on the other hand, predict significantly higher threshold velocities for fragmentation, ranging between 1 and 12 m s^{1} for 6 to 10 cmsized grains (Wada et al. 2011; Meru et al. 2013b). The value of the fragmentation threshold velocity determines the maximum size of particles that we are able to obtain before the breakthrough happens. Mass transfer experiments are even more uncertain with the mass transfer efficiency that ranges between 0 and 60% (Teiser & Wurm 2009; Kothe et al. 2010; Beitz et al. 2011; Meisner et al. 2013), and the critical mass ratio, which is in principle unexplored in the laboratory and only for small ratios numerically (Meru et al. 2013b; Wada et al. 2013). In all of these cases, additional material and collisional properties, such as porosity, composition, structure and impact angle also greatly influence the outcome, which adds to the uncertainty. As discussed in this work, the critical mass transfer ratio might need to be significantly lower than estimated in prior studies, due to the need for a relatively high numerical resolution of Smoluchowski solvers to accurately represent the highmass tail.
The collision models used in sweepup modeling attempt to simplify the very complex physics of collisions between dust agglomerates. There have been a few attempts at more rigorous models (Zsom et al. 2010; Windmark et al. 2012a) to consolidate the recent progress in laboratory and numerical collision experiments (Wada et al. 2007; Blum & Wurm 2008; Güttler et al. 2010), but these are still not necessarily more correct than the simple models that narrow the modeling down to only a few key parameters (Windmark et al. 2012b; Garaud et al. 2013; Meru et al. 2013a; Drążkowska et al. 2013). As the critical mass ratio plays crucial role in determining if the breakthrough is possible, restricting its realistic values should be a priority for the future laboratory studies.
We find the breakthrough scenario to be more sensitive to resolution issues than other problems we have tested, when the dust growth is utterly stopped by bouncing or fragmentation. At the same time, this scenario is of particular importance for the current planetesimal formation theory. Regardless of the method used, modeling of the “lucky growth” requires extreme computational force, and restricting the resolution to make the models more efficient can lead to serious numerical artifacts: the nonphysical breakthrough in the Smoluchowski approach case and lack of breakthrough in the case of the Monte Carlo approach.
To conclude, we list features that characterize the two main coagulation methods. One should keep these in mind when deciding which method to use for a particular scientific application.
The Smoluchowski equation solver with implicit integration scheme

is capable of simulating dust evolution over long timescales (even at high resolutions, 0D simulations are finished within minutes);

resolves equilibrium states well, as the implicit integration scheme allows for very large timesteps once the solution approaches steadystate (Birnstiel et al. 2010);

has a very high dynamical range that allows phenomena involving a single physical particle of tiny mass (compared to the total dust mass) to be resolved, which makes it ideal for breakthrough studies (Windmark et al. 2012a), and to produce synthetic observations, as the opacity is dominated by small particles, which may contain only a low fraction of the dust mass;

is slowed by a factor of , where n is the number of mass bins, for each additional dust property beyond mass (Windmark et al. 2012a), although numerical tricks exist to circumvent this (see, e.g., Okuzumi et al. 2009);

suffers from high numerical diffusion that affects both growth timescale and breakthrough likelihood. The growth timescale can be benchmarked against the analytical kernels (e.g. Ohtsuki et al. 1990), but it depends strongly on the strength of the barrier in the breakthrough case.
The Monte Carlo coagulation algorithm with equalmass representative particles

makes it easy to implement additional particle properties, such as porosity (Zsom & Dullemond 2008; Zsom et al. 2010), because the computation time does not depend significantly on the number of properties that are evolved, but only on the number of collisions performed;

is straightforward to develop to further spatial dimensions (Zsom et al. 2011; Drążkowska et al. 2013), as the representative particles can be treated as Lagrangian tracer bodies;

can be used along with hydrodynamic grid codes (Johansen et al. 2012);

experiences no numerical diffusion of the mass function in general, so there is no danger of encountering an artificial speedup of the growth;

has difficulty resolving features that include low fraction of total mass, which makes it less useful in the case of breaking through the growth barriers or runaway growth modeling, although the algorithm can be developed to overcome this issue (Ormel & Spaans 2008),

makes it hard to model evolution over long timescales in general, because it is impossible to use extremely long timesteps, as every collision needs to be resolved.
Acknowledgments
We thank the anonymous referee as well as Chris Ormel, Andras Zsom, Satoshi Okuzumi, Til Birnstiel, Sebastian Stammler and Sebastian Lorek for useful comments that helped to improve the manuscript. J.D. was partially supported by the Innovation Fund FRONTIER of the Heidelberg University. F.W. was funded by the Deutsche Forschungsgemeinschaft within the Forschergruppe 759 “The Formation of Planets: The Critical First Growth Phase”. J.D. would also like to acknowledge the use of the computing resources provided by bwGRiD (http://www.bwgrid.de), member of the German DGrid initiative, funded by the Ministry for Education and Research (Bundesministerium für Bildung und Forschung) and the Ministry for Science, Research and Arts BadenWuerttemberg (Ministerium für Wissenschaft, Forschung und Kunst BadenWürttemberg)
References
 Beitz, E., Güttler, C., Blum, J., et al. 2011, ApJ, 736, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Charnoz, S., & Taillifet, E. 2012, ApJ, 753, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444 [NASA ADS] [CrossRef] [Google Scholar]
 Drążkowska, J., Windmark, F., & Dullemond, C. P. 2013, A&A, 556, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Garaud, P., Meru, F., Galvagni, M., & Olczak, C. 2013, ApJ, 764, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Gillespie, D. T. 1975, J. Atm. Sci., 32, 1977 [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]
 Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, A&A, 537, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kothe, S., Güttler, C., & Blum, J. 2010, ApJ, 725, 1242 [NASA ADS] [CrossRef] [Google Scholar]
 Kothe, S., Blum, J., Weidling, R., & Güttler, C. 2013, Icarus, 225, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Kovetz, A., & Olund, B. 1969, J. Atm. Sci., 26, 1060 [NASA ADS] [CrossRef] [Google Scholar]
 Lammel, C. 2008, Bachelor’s Thesis, Technische Universität Carolo Wilhelmina Braunschweig [Google Scholar]
 Lee, M. H. 2000, Icarus, 143, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Meisner, T., Wurm, G., Teiser, J., & Schywek, M. 2013, A&A, 559, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meru, F., Galvagni, M., & Olczak, C. 2013a, ApJ, 774, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Meru, F., Geretshauser, R. J., Schäfer, C., Speith, R., & Kley, W. 2013b, MNRAS, 435, 2371 [NASA ADS] [CrossRef] [Google Scholar]
 Nakagawa, Y., Nakazawa, K., & Hayashi, C. 1981, Icarus, 45, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Nomura, H., & Nakagawa, Y. 2006, ApJ, 640, 1099 [NASA ADS] [CrossRef] [Google Scholar]
 Ohtsuki, K., Nakagawa, Y., & Nakazawa, K. 1990, Icarus, 83, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., & Sakagami, M.A. 2009, ApJ, 707, 1247 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., & Spaans, M. 2008, ApJ, 684, 1291 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., Dullemond, C. P., & Spaans, M. 2010, Icarus, 210, 507 [NASA ADS] [CrossRef] [Google Scholar]
 Ros, K., & Johansen, A. 2013, A&A, 552, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schräpler, R., Blum, J., Seizinger, A., & Kley, W. 2012, ApJ, 758, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Seizinger, A., & Kley, W. 2013, A&A, 551, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tanaka, H., Himeno, Y., & Ida, S. 2005, ApJ, 625, 414 [NASA ADS] [CrossRef] [Google Scholar]
 Teiser, J., & Wurm, G. 2009, MNRAS, 393, 1584 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2007, ApJ, 661, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weidenschilling, S. J. 1980, Icarus, 44, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Wetherill, G. W. 1990, Icarus, 88, 336 [NASA ADS] [CrossRef] [Google Scholar]
 Wetherill, G. W., & Stewart, G. R. 1989, Icarus, 77, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Windmark, F., Birnstiel, T., Güttler, C., et al. 2012a, A&A, 540, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012b, A&A, 544, L16; erratum, 548, C1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wurm, G., Paraskov, G., & Krauss, O. 2005, Icarus, 178, 253 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Zsom, A., & Dullemond, C. P. 2008, A&A, 489, 931 [NASA ADS] [CrossRef] [EDP Sciences] [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]
 Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Figures
Fig. 1 Comparison of the mass distribution evolution obtained using our Smoluchowski and Monte Carlo codes. In a standard case when an insurmountable fragmentation barrier is present, the two methods perfectly agree (three upper panels); however, we encounter some differences between the results obtained with the two approaches (three bottom panels) with the possibility of breakthrough. For the Smoluchowski code, the presented results were obtained with a resolution of 30 mass bins per decade. For the Monte Carlo code, the results are averaged from the simulations using 120 000 particles. The “check points”, m_{1} and m_{2}, which are indicated with the dotted lines, are used to quantitatively compare dust growth timescales in Sect. 4.1. 

Open with DEXTER  
In the text 
Fig. 2 Time at which the peak of mass distribution reaches the mass m_{1} = 2 × 10^{5} g for both the Smoluchowski and Monte Carlo code. The result changes with mass resolution, given by the number of bins per mass decade in the case of Smoluchowski code and the number of particles used in the Monte Carlo code. The scatter of results obtained in different runs with the same number of particles in the Monte Carlo code is marked by the shaded region around an averaged dependence. The Monte Carlo approach does not exhibit a strong resolution dependence. In contrast, the Smoluchowski algorithm overestimates the growth rate by a factor of few when we do not use enough mass bins. 

Open with DEXTER  
In the text 
Fig. 3 Same as Fig. 2 but for the “check point” m_{2} = 1 g. The resolution dependence of the Smoluchowski method is the same as for m_{1}, but it changes significantly in the case of the Monte Carlo code. The Monte Carlo code tends to underestimate the breakthrough possibility when used with a low number of particles due to the limited dynamic mass range. 

Open with DEXTER  
In the text 
Fig. 4 Effect of mass resolution in the case of SF+MT setup with a critical mass ratio of 500 for the Smoluchowski code. The resolutions included a range from 3 to 40, and the mass distribution is shown after 100 yrs of evolution. The point of breakthrough is very sensitive to the resolution used, and the breakthrough happens only in the low resolution cases, which is clearly nonphysical. The dashed lines mark the density corresponding to 1 and 10^{6} particles of a given mass in an annulus of 0.1 AU width. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.