Effect of turbulence on collisions of dust particles with planetesimals in protoplanetary disks
Key Words.:
Planets and satellites: formation, Planetdisk interactions, Turbulence, Accretion, Accretion disks, hydrodynamics, methods: numericalAbstract
Context:Planetesimals in gaseous protoplanetary disks may grow by collecting dust particles. Hydrodynamical studies show that small particles generally avoid collisions with the planetesimals because they are entrained by the flow around them. This occurs when St, the Stokes number, defined as the ratio of the dust stopping time to the planetesimal crossing time, becomes much smaller than unity. However, these studies have been limited to the laminar case, whereas these disks are believed to be turbulent.
Aims:We want to estimate the influence of gas turbulence on the dustplanetesimal collision rate and on the impact speeds.
Methods:We used threedimensional direct numerical simulations of a fixed sphere (planetesimal) facing a laminar and turbulent flow seeded with small inertial particles (dust) subject to a Stokes drag. A noslip boundary condition on the planetesimal surface is modeled via a penalty method.
Results:We find that turbulence can significantly increase the collision rate of dust particles with planetesimals. For a high turbulence case (when the amplitude of turbulent fluctuations is similar to the headwind velocity), we find that the collision probability remains equal to the geometrical rate or even higher for , i.e., for dust sizes an order of magnitude smaller than in the laminar case. We derive expressions to calculate impact probabilities as a function of dust and planetesimal size and turbulent intensity.
Conclusions:
1 Introduction
Conventional models for planet formation involve the hierarchical growth by accretion of “planetesimals”. Such building blocks undergo collisions and gravitational binding to eventually reach planetary sizes (see, e.g., Safronov 1972; Hayashi et al. 1985). Still several serious uncertainties remain in the processes leading from subm dust grains, which follow the subKeplerian gas, to kmsized planetesimals that are massive enough to move on Keplerian orbits. One of them is known as the “metersize barrier”. Metersized bodies experience a strong drag force, which causes rapid orbital decay due to angular momentum loss. The timescale of this decay is less than 100 years at 1AU (see, e.g., Weidenschilling 1984; Nakagawa et al. 1986), which is shorter by several orders of magnitude than the observationally inferred disk lifetime (several million years). Consequently, the planetforming material are lost from the disk.
One scenario for overcoming the metersize barrier is by the combined effect of streaming instabilities and pebble accretion. Because of growth to mm/cmsized particles and settling, dust concentrates in the midplane and may become unstable to streaming instabilities (Youdin & Goodman 2005; Johansen et al. 2007, 2011). This forms relatively large planetesimals. After that, these large planetesimals can sweep up surrounding grains and migrating pebbles (Ormel & Klahr 2010; Ormel & Kobayashi 2012; Lambrechts & Johansen 2012). Because the orbital decay of the pebbles is also fast, a large number of pebbles are supplied from outer disk regions on relatively short timescales, resulting in rapid planet growth. This scenario utilizes the too rapid migration problem, rather than avoids it. The pebble accretion model is actively discussed for the formation of the solar system (Morbidelli et al. 2015) and of closein superEarth systems in exoplanetary systems (Chatterjee & Tan 2014, 2015). However, the details of the the accretion of grains and pebbles by (large) planetesimals have not been fully clarified.
Guillot et al. (2014) proposed detaile expressions for the accretion rates of dust grains in a protoplanetary gaseous disk for a wide range of dust and planetesimal sizes. They point out that, when using the numerical results by Sekiya & Takeda (2003), the accretion probability drops off by several orders of magnitude at relatively small dust grains (within what they refer to as the “hydrodynamical regime”). This strong depletion can be easily explained by qualitative physical arguments. The motion of small grains is strongly coupled to the gas flow. Since they follow gas streamlines they may avoid (headon) collisions with the planetesimal. However, the numerical simulations by Sekiya & Takeda (2003) assume laminar flow, while it is known that the disk gas is most likely to be turbulent. Observationally inferred accretion rate in TTauri disks is far higher than that caused by molecular viscosity. Consequently, protoplanetary disks are expected to be turbulent, with a turbulent eddy viscosity much higher than the molecular kinematic viscosity. Turbulence should affect the reduction in the dust accretion probability in hydrodynamical regime. Recently, Mitra et al. (2013) studied the influence of turbulence on the dust impact velocity by means of twodimensional direct numerical simulations. They found a velocity distribution with exponential tails and argued that most of the dust collides with a speed comparable to that of the head wind.
The purpose of this paper is to propose expressions for dust accretion probabilities on planetesimals in turbulent gas, based on threedimensional direct hydrodynamical simulations. Gravity effects are omitted in this study and are left for a future work. We follow the approach introduced by Homann & Bec (2015) and consider an idealized situation in which the planetesimal is assumed to be spherical with a smooth surface. The gas is modeled as a purely hydrodynamical and incompressible fluid so that other possible modes originating from either the magnetorotational instability (MRI) or compressibility are not taken into account. Additionally, we neglect any shear, disregard gravitational effects, and do not allow for rotation of the planetesimal. The aim of this work is twofold: First, it shall provide a detailed study of collisions between small inertial grains and a large spherical object and, secondly, it shall determine its consequences on the accretion of dust by planetesimals in the context of planet formation.
It is worth mentioning that our results, involving the hydrodynamical interactions and the collisions between a large spherical inclusion and small particles with inertia, have important applications in contexts that go beyond planet formation. In atmospheric physics, this problem is known as “inertial impaction” and is relevant for estimating rates in rain formation or wet deposition of aerosols where a falling water drop scavenges smaller cloud droplets or solid pollutants. Studies of such problems often rely on collision efficiencies and, again, mainly formulas from laminar flow conditions are used (Berthet et al. 2010). A popular formula is that of Slinn (1974), who proposed a fit of the collision efficiency for the inertial impaction regime. Later, he included molecular diffusion and extended his formula to smaller projectiles (Slinn 1976, 1983). Inertial impaction is also important for the design and improvement of industrial filters. Haugen & Kragset (2010) studied the impact of particles on a cylinder in a twodimensional laminar inflow. Later, Hydle Rivedal et al. (2011) investigated the case of a turbulent inflow and found that turbulent fluctuations yield up to ten times more collisions than a corresponding laminar inflow.
The paper is organized as follows. In Section 2, we explain the disk model and notations. In Section 3, we describe the numerical method that we use in our approach. In Section 4, we present the results of our hydrodynamical simulations. In particular, we propose an expression for the accretion probability in laminar and turbulent flows as a function of both the Stokes number, which compares dust inertia to the perturbation of the gas motion due to the planetesimal, and the turbulent intensity of the surrounding gas flow. Section 5 is devoted to astrophysical discussions on the obtained results. Section 6 contains a summary and some concluding remarks.
2 Protoplanetary disk turbulence model & notations
2.1 Disk model
We consider a minimum mass solar nebula (MMSN) model (Weidenschilling 1977b; Hayashi 1981), where the surface density and temperature of the gas in the disk are given in terms of power laws:
(1)  
(2) 
and are the values at 1 a.u., for which we choose and K. (A list of symbols and their definitions used in this paper is summarized in the appendix A.) Assuming the disk is vertically isothermal, with the isothermal sound speed and a mean molecular weight of , the density reads
(3) 
where is the disk scale height and the orbital (Keplerian) frequency at distance from the star. It follows that , so that the disk is flared.
Due to the density and temperature gradients, the disk is slightly supported by pressure. The amount of pressure support compared to the solar gravity is often expressed in terms of a parameter defined as
(4)  
where we introduced the pressure logarithmic gradient , the Keplerian velocity , and the disk aspect ratio . The numbers are obtained from an MMSN disk with the above power law profiles for density and temperature. It follows that the motion of the disk is less than Keplerian (Weidenschilling 1977a); the gas flows at a speed equal to . Also, the headwind that is faced by a big body (a.k.a. planetesimal), which moves at a Keplerian speed through the subKeplerian rotating gas, then reads
(5) 
which, for an MMSN disk is independent of the distance to the star.
2.2 Turbulence model
The gas has a meanfree path of where (Chapman & Cowling 1970) is the molecular cross section. The kinematic molecular viscosity of the gas is where is the mean molecular thermal speed.
However, for these parameters the resulting diffusion timescales are too long to explain, for example, the measured accretion rate in TTauri disks. Consequently, protoplanetary disks are expected to be mildly turbulent, but still subsonic, with a turbulent eddy viscosity assumed to be much larger than the molecular kinematic viscosity. An often used parameterization is the alphaviscosity model of Shakura & Sunyaev (1973):
(6) 
The dimensionless constant is a diskdependent parameter. Its value may range from a minimum of , when the turbulence originates from the KelvinHelmholtz instability of a dust layer at the midplane, to perhaps for the most violent disks with a strong magnetorotational instability (the MRI; Balbus & Hawley 1991). Generally, the saturation level of the MRI turbulence depends on the magnitude of the magnetic field that threads the disks, the level of ionization, and the amount of dust (Turner et al. 2014). Consequently, a variety of values can be expected (Ormel & Okuzumi 2013). But even for a weak level of turbulence – meaning, – the Reynolds number of the flow, , is very large by virtue of the very low densities that characterize astrophysical environments. Indeed, for an MMSN density profile:
(7) 
For such expected large values of the Reynolds number, one can reasonably assume that the gas flow is in a developed threedimensional turbulent regime which can be described using Kolmogorov (1941) phenomenology. The kinetic energy is injected at a large (integral) scale by a hydro or magnetohydrodynamical instability. In both cases the injection mechanism is associated to the subKeplerian shear. The associated shear rate sets the typical turbulent largeeddy turnover time to (Cuzzi et al. 2001). The eddy viscosity then reads leading, together with Equation (6), to and . Note that the assumptions of incompressibility () and threedimensionality () of the turbulent flow in the disk are both fulfilled as long as . According to Kolmogorov (1941) phenomenology, the kinetic energy cascades downscale with a constant rate until it reaches the smallest active scales, the Kolmogorov scale below which it is dissipated by molecular viscosity. The typical velocity of eddies whose size lies in the inertial range reads . The Kolmogorov length is defined as the scale where the scaledependent Reynolds number is unity, so that and the associated timescale . The integral Reynolds number gives the extension of the inertial range: and .
2.3 Dimensionless quantities
Let us now turn back to our original problem, that is the
aerodynamic interactions between the disk gas flow and a solid
planetesimal of size , diameter . There are three dimensionless quantities characterizing
the system: (the numerical values below are those obtained for an MMSN disk)

The planetesimal Reynolds number
(8) which corresponds to the strength of inertia with respect to molecular viscous forces for the gas flow surrounding the planetesimal and characterizes how turbulent is its wake.

The turbulent intensity
(9) which measures how strong the turbulent velocity fluctuations are compared to the speed of the planetesimal slip.

the ratio between the planetesimal size and the Kolmogorov scale
(10) which measures the range of turbulent eddies that might interfere with the planetesimal perturbation of the gas flow. It is indeed known that the flow perturbation occurs on scales of the order of , so that all turbulent eddies of sizes between and can potentially modify the gas flow around the spherical planetesimal.
Figure 1 shows the dependence of the various parameters with the orbital distance for two values of compared to the assumptions used for the numerical calculations. We can see that decreases with increasing orbital distance (kinematic viscosity goes up). Similarly, decreases because increases (turbulent scales are getting larger). The turbulent intensity is not a very strong function of disk radius, but it depends quite strongly on the turbulence parameter .
In the laminar case, corresponds to a 1 km planetesimal located at au in the MMSN disk, or e.g. to a 10 km planetesimal at au.
For the turbulent case, ideally we should find parameters for which all three lines match with parameters of the numerical experiments. However, because of the mismatch between the Reynolds number in the disk and the maximum one that can be attained by presentday numerical simulations, this is not yet possible. We come back to this issue when applying the results of numerical models to real disk conditions.
3 Numerical model
3.1 Gas flow
We focus on the collision rate of a stream of dust particles with one spherical planetesimal. To study this situation we perform 3D direct numerical simulations (DNS) of a hydrodynamic flow around a spherical object with a noslip boundary conditions at its surface. The overall method consists in a combination of a standard pseudoFourierspectral solver with a penalty method that is explained in this section.
We integrate the incompressible NavierStokes equations
(11) 
for the gas velocity , where is the gas density, its kinematic molecular viscosity and a force. The latter maintains a uniform inflow speed and an eventually ambient turbulent flow. Its form is described later. The pseudoFourierspectral approach consists in computing spatial derivatives in Fourierspace and convolutions arising from the nonlinear terms in real space. A FastFourier transform is used to switch between the two spaces. We use the P3DFFT library (see Pekurovsky 2012) that is very efficient on massive parallel computers.
The Navier–Stokes equation is integrated in the reference frame of the planetesimal. The spatial average of the velocity field is thus fixed and given by the planetesimal speed relative to the gas. Equation (11) is associated with a noslip boundary condition at the surface of the planetesimal, which is assimilated to a spherical object at rest whose diameter is denoted by and its center by . We thus have
(12) 
Numerically, this noslip condition is enforced by an immersed boundary technique (IBM). The latter technique was first used to simulate the blood in flow in the context of a human heart by Peskin (1972). Today, a variety of different approaches exists (Mittal & Iaccarino 2005). Generally speaking, IBM consists in solving fluid equations in domains with complex boundary conditions such as moving heart valves on Cartesian grids. As such boundaries are generally not grid conform, their effect on the flow has to be modeled. For our problem of a spherical obstacle at rest, the idea consists in defining the velocity field in the full domain enforcing a vanishing velocity in the entire object, that is for all such that . The NavierStokes equation (11) then has to be modified by introducing in its righthand side a penalty force , which acts as a Lagrange multiplier associated to the constraint defined by the boundary condition (12). The full problem (11)(12) can then be rewritten as
(13) 
where denotes the right hand side of (11). In order to compute we make use of a direct forcing method introduced by Fadlun et al. (2000) where we directly impose the planetesimal velocity to the grid using the technique of a pressure predictor and an improved modeling of the spherical inclusion. Benchmarks (see Homann et al. 2013) of this method for a fixed sphere show good agreement with existing data. This method has also been used for the study of moving neutrally buoyant particles in homogeneous isotropic turbulence in Homann & Bec (2010) and Cisse et al. (2013). A similar IBM has been used by Uhlmann (2005) together with second order finitedifference NavierStokes solver to simulate turbulent suspension involving many particles.
3.2 Dust particles
Dust particles are modeled by spherical inertial particles with a radius much smaller than the smallest scales of the flow, so that they can be approximated by point particles. Further, we assume that these particles move sufficiently slow with respect to the gas and that their mass density is much higher than the gas density . With these assumptions the dominant hydrodynamic force exerted by the gas is a drag force, which is proportional to the velocity difference between the particle and the gas flow (see Maxey & Riley 1983; Gatignol 1983)
(14) 
where the dots stand for time derivatives. is called the response (or stopping) time and is a measure of particle inertia. It is the relaxation time of the particle velocity to that of the gas (see Guillot et al. 2014, for a description of how the particle size relates to the stopping time). Finally, we also assume that the particles are sufficiently diluted to neglect any interaction among them and any backreaction on the flow.
Usually, particle inertia is quantified in terms of the Stokes number defined by nondimensionalizing their response time by a characteristic time scale of the carrier flow. The present problem involves different relevant time scales. , the time it takes a dust particle to pass the planetesimal, , the turbulent largeeddy turnover time and , the turbulent dissipation time scale. If not otherwise specified, we use as this time scale rules the collision efficiency in the laminar regime: Particles with small St are closely coupled to the flow and are swept around the obstacle. Large St particles preferentially collide with it.^{1}^{1}1Note that in the context of astrophysical disks, a different Stokes number is often defined from . The definition of St that is used here is denoted by in e.g., Sekiya & Takeda (2003) and Guillot et al. (2014).
In studies concerned with the smallscale dynamics of inertial particles one usually uses (Bec et al. 2006) and we denote the associated Stokes number .
3.3 Simulation setup
We performed simulations of a planetesimal moving through a uniform and different turbulent flows. The physical situation is illustrated in Fig. 2 where the planetesimal, together with dust particles are shown in a small slice. These are snapshots taken from our threedimensional simulations. The planetesimal experiences a headwind speed from the right that is advecting the dust particles. The upper panel is taken from a simulation with a uniform inflow while the two others include turbulence with two different intensities. Dust particles that have collided with the planetesimal are missing downstream, resulting in an empty region behind the planetesimal that is clearly visible for the laminar flow, but much less in the turbulent cases.
These figures show important differences between laminar and turbulent gas flows. The spatial distribution of the dust as well as its local velocity strongly depend on the carrier flow type. Note that while for all cases shown in Fig. 2, for and for so that the clustering properties of the particles change from one value of to the other.
The runs involving a turbulent flow are set up in the following way (see table 2 for parameters and definitions). An initially smooth large scale flow is integrated according to (11) without any forcing . Once a turbulent flow has developed (after approx. ) the velocity field is forced by keeping constant the energy content of the two lowest wave number shells ( in Fourier space. This leads to a statistically stationary turbulent flow to which the planetesimal is added. For this, the rootmeansquare value of the velocity fluctuations is normalized to the values given in table 2 for each simulation. A mean velocity of in one direction is imposed by keeping constant the zero Fourier mode of the corresponding component of the velocity. The planetesimal is modeled via the mentioned immersed boundary method. The integration is continued until a statistically stationary state is reached again. At this point, inertial particles are seeded at a constant rate into the flow in a plane sufficiently far from the planetesimal so that they have enough time () to relax to the flow. During the simulation we remove and record all the particles that are touching the spherical planetesimal or reaching the end of the computational domain. On average the domain is filled with approximately ten million particles.
The laminar flow simulations (see table 1 for parameters and definitions) start with a uniform flow () in which the planetesimal is placed. Disturbances produced by its wake are removed at the end of the computational domain via another application of the penalty method, so that they are not reinjected upstream the planetesimal by the periodic boundary conditions. The dust particles are introduced into the flow once a (statistically) stationary state is reached.
The time integration of (11) uses a RungeKutta scheme of third order. The grid resolution is chosen in order to resolve all small scales of the problem: those of the spherical planetesimal boundary layer, all turbulent scales in its wake and the smallest scales of the possibly turbulent ambient flow.
All physical flow parameters are determined by three dimensionless parameters: the Reynolds number of the planetesimal , the Reynolds number of the gas , being the integral scale of the ambient turbulent flow, and the turbulent largescale intensity . The latter measuring the strength of the largescale turbulent fluctuations compared to the mean flow velocity. In the laminar case (, ), we varied the planetesimal Reynolds number from 100 to 1000. We analyze the effect of turbulent fluctuations for one specific choice of the planetesimal Reynolds number, namely (turbulent wake) and vary from 0.14 to 1.18. The corresponding flow Reynolds numbers (listed in Table 2) are a consequence of our particular choice of the external force. Freezing the energy content of the lowest shells in spectral space does not allow for changing but only that enters in the definitions of both and .
The particle dynamics is characterized by the Stokes number St. In all simulations we consider streams of heavy dust particles with response times in between 0.04 and 81.92, corresponding to Stokes numbers St in the range . The main parameters of all simulations are summarized in table 2 and table 1.
d  

400  200  1  0.14  0.002  0.8  2.92  21.0  
400  450  1  0.29  0.002  0.8  3.04  10.6  
400  600  1  0.39  0.002  0.8  3.09  7.9  
400  900  1  0.60  0.002  0.8  3.10  5.2  
400  3000  1  1.18  0.002  0.8  3.40  2.9 
100  0.8  

400  0.8  
1000  0.65 
4 Results
4.1 Laminar settings
In a uniform gas flow all dust particles have the same velocity far from the planetesimal. This physical situation is fully determined by the Reynolds number of the planetesimal and the Stokes number St of the dust particles. A typical flow pattern is illustrated in the upper panel of Fig. 2 where a planetesimal flies through a stream of dust particles while creating a moderately turbulent wake. If an approaching dust particle collides with the planetesimal or not is only determined by its Stokes number and impact parameter. Without the hydrodynamic flow that deflects dust around the obstacle it would just be the impact parameter ruling the collision rate so that the hydrodynamic forces reduce the collision rate below that of the geometric cross section.
The colliding dust particles preferentially hit the planetesimal on the axis of symmetry with a decreasing probability to its edge. Small inertial particles only touch the planetesimal in a central region leaving eventually an outer noncollisional ring (see Fig. 3). But this region already disappears for Stokes numbers around unity, so that dust collisions fill the complete front of the spherical planetesimal. In the limit of infinite the distribution become uniform. Rear collisions virtually never happen (see Fig. 4) in a laminar flow that is because inertia prevents those particles that moved around the planetesimal from getting into the recirculation region of the wake. The few records observed might be numerical noise.
A central quantity of the accretion problem is the collision efficiency (or collision rate) that is the ratio between the actual collision rate and that expected for freestreaming particles that do not feel any hydrodynamic force from the gas flow. The latter rate is obtained from the geometrical crosssection of the planetesimal and reads , where designates the number density of dust particles. We use the following notation: denotes the uniform gas flow case that we are going to analyze in this section and denotes the turbulent case parameterized by the turbulent intensity that is discussed in the following section. The efficiency is a monotonously increasing function of the dust particle Stokes number St that asymptotically reaches unity for particles with very large inertia (see Fig. 5). The higher is the more probable are collisions, especially at small St. In this low inertia limit, sharply drops for all . The reason for this is the well established existence of a critical Stokes number below which no collisions occur at all (Taylor 1940).
Indeed, in the case of an inviscid Euler flow (zero viscosity, ) it can be shown that drops to zero for small but finite (Ingham et al. 1990). The viscous boundary layer of a finite planetesimal that shrinks as reduces the collision probability and increases . For the limiting case of a Stokes flow Michael & Norey (1970) computed numerically a critical Stokes number of .
Slinn (1974, 1976, 1983) proposed a fitting function of the collision rate of the form
(15) 
is thus a critical number below which the (extremely low) collision probability is determined by different physical mechanisms. These expressions appear to have been fitted empirically to both experimental and numerical data with uncertainties of order on the determination of (see Slinn 1974).
It is useful to look in more detail at the small and large St efficiencies separately. In Fig. 6 the collision efficiency is shown as a function of and not simply as a function of St. We hence focus on the behavior of the collision probability when approaching the critical Stokes number. is chosen in such a way that all curves fall on top of each other so that differences for different Reynolds numbers disappear. Our estimated are close to that of Slinn (15) (see inset of Fig. 6). The case of large but finite values of has also been addressed by Phillips & Kaye (1999), using matched asymptotics together with numerical simulations of the particle dynamics inside the obstacle boundary layer. They found critical Stokes numbers slightly larger than those of Slinn (15), so that our value are in between these two predictions.
increases linearly at small (see Fig. 6), that is in contradiction to Slinn’s formula (15) that predicts . To incorporate this smallStokes behavior we propose the fitting function
(16) 
which is in excellent agreement with our data. Here we mention (as already remarked by Slinn) that an exponential fit of the form also works quite well for not too small St but this form does evidently not conform to a critical Stokes number.
We also note that Slinn (1974) had also proposed a similar fit (linear instead of with a exponent), but opted for eq. (15) which showed a marginal improvement to the experimental and numerical data available at the time. The uncertainties that we obtain here are much smaller and allow to discriminate against the fit from eq. (15).
To study the largeStokes number behavior in details it is useful to analyze that measures how the freestreaming particle limit is recovered as a function of St. All curves for different (see Fig. 7) fall on the top of each other and reveal a behavior at large St. Again, we observe slight deviations to Slinn’s fitting formula, while our proposed expression (16) fairly matches the data.
4.2 Turbulent settings
One can expect that turbulent velocity fluctuations of the gas, resulting in dust velocity fluctuations, alter the collision statistics of dust and planetesimals. An analysis of these changes is the subject of this section.
The turbulent accretion problem involves more parameters than the laminar problem presented in the previous section. It is determined by four dimensionless parameters. The planetesimal and dust properties are of course still described by the Reynolds number and the dust Stokes number St, but turbulence adds two additional dimensionless parameters specifying the ambient turbulent flow. They are the gas Reynolds number and the turbulent intensity that compares the amplitude of turbulent fluctuations with the mean velocity of the gas.
In this work, we explicitly vary St and the turbulent intensity and fix the planetesimal Reynolds number to (weakly turbulent wake) in order to limit computational costs. The gas Reynolds number is implicitly varied as it is coupled to due to our specific forcing scheme that prescribes to approximately half of the domain size.
The physical situation for two different turbulent intensities is illustrated in the mid and bottom panel of Fig. 2. Dust particles are in these cases advected by a chaotic and irregular flow possessing coherent structures, i.e. structures eventually persisting for a long time. This has two important implications: First, the velocity of dust is fluctuating spatially and temporally (see Fig. 2 b) and c)) and in turn (as we study in details below) modifies the collision statistics. Second, from studies of hydrodynamic turbulence it is known that inertial particles tend to escape from coherent rotating regions of the flow and tend to cluster in straining regions (Squires & Eaton 1991). These agglomerations are called preferential concentrationsShaw (2003); Balachandar & Eaton (2010). However, we expect that these concentrations play a negligible role for the present study in which we are concerned with averaged accretion quantities such as the collision efficiency. The reason is that the temporally averaged dust concentration in the ambient flow is approximately the same as in the laminar case (not shown) so that particle density fluctuations average out.
4.2.1 Collision efficiency
Turbulent fluctuations significantly increase the collision probability. Figure 8 shows the collision efficiency for various . One observes that higher is the turbulent intensity , more collisions happen. This increase is the strongest at small St and disappears asymptotically at large St. As is discussed later, the collision efficiency around remarkably exceeds unity.
For the smallest Stokes numbers and the largest turbulent intensity we observe more than one hundred times more collisions than in the reference laminar flow. This relative increase of the turbulent efficiency compared with the laminar one is shown in Fig. 9. It is represented as a function of that is the distance from the critical Stokes number as diverges when . The large values of at small St decreases as a power law with an exponent close to . The dependence of on is nearly quadratic. Indeed, the different curves almost fall on the top of each other when normalized by as can be seen from the inset of Fig. 9.
Turbulent collision efficiencies (compare Fig. 8) are well fitted by the function
(17) 
containing two parameters and . The dependence of these parameters on the turbulent intensity is shown in Fig. 10 together with simple fitting formulas. The parameter is decreasing roughly exponentially as a function of . is close to zero for small intensities and strongly increasing starting from .
In the small Stokes number limit, the turbulent collision efficiency displays a clear exponential falloff (see Fig. 11) indicating that the critical Stokes number disappears when turbulence influences the dust motion. Additionally, curves for different fall approximately on top of each other once shown as a function of . It is thus a turbulent time scale, namely the largeeddy turnover time, that replaces the advection time relevant for the laminar problem.
There are different ways in which turbulence enhances the collision probability. Fluctuations of the headwind speed mix collision efficiencies for different and St and allow for collisions of dust particles beyond the critical Stokes number . In turbulent gas flows, the dust heads onto the planetesimal with a fluctuating velocity. The probability distribution of the headwind is given by the onepoint turbulent velocity distribution that is known to be close to a Gaussian. For the present problem, its standard deviation is given by the turbulent intensity . These headwind variations lead to fluctuating planetesimal Reynolds number and dust Stokes numbers. Especially close to the laminar value this results in higher collision probabilities and an absence of a critical Stokes number. Another consequence of a variable headwind are fluctuations of the streamwise velocity gradient in the upstream boundary layer of the planetesimal. They modify the local Stokes number that can be defined by . The probability for a collision of a dust particle with stopping time is then , and thus relates to the probability of observing a large velocity gradient at the particle surface. The probability density function (PDF) of shown in the inset of Fig. 11. The tails are close to exponential, although we cannot rule out stretched exponential tails as observed in homogeneous isotropic turbulence. Exponential tails mean which is consistent with the formally observed exponential fall off of at small values of St.
Another mechanism to increase the collision probability relates to turbulent diffusion. This enhances the mobility of dust and increases its flux onto the planetesimal. As a simplified gedankenexperiment one can think of a sphere moving in a sinusoidal velocity field representing the largescale turbulent fluctuations. Let us assume (in the reference frame of the sphere) a velocity of the form . The volume swept by the sphere is evidently larger in the turbulent case () than in the laminar (). According to inertia, dust particles are more or less coupled to the gas which makes the swept volume additionally depending on St. Small St particles stick to the gas trajectories so their swept volumes equal. This effect can explain the observation of collision efficiencies exceeding unity in turbulent flows. For large St, particles move ballistically so that they only sweep through the laminar () volume.
We conclude this section with a study of the spatial distribution of dust impacts on the surface of a planetesimal in a turbulent disk. In Fig. 4 we saw that dust collisions preferentially happen close to the stagnation point of the flow in a quiescent disk. This is still true when turbulent fluctuation agitate the dust (see Fig. 12). But the added randomness leads to a homogenization of the impact position. The larger is the more the collisions fill the entire planetesimal surface. And especially for small St particles, backward collision become frequent.
4.2.2 Impact velocity
The relative velocity of dust and a planetesimal at impact is crucial for the dust accretion problem as it determines, together with the angle of impact, the outcome of a collision. Low collision speeds lead to sticking of dust on the target surface, while high speeds lead to bouncing, fragmentation with mass transfer or erosion (e.g., Blum & Wurm 2008; Windmark, F. et al. 2012).
In laminar disks the mean impact speed of smallsize dust (with a stopping time smaller than the orbital period) is a monotonously increasing function of inertia (Weidenschilling 1977a). Dust particles with inertia close to the critical Stokes number only slightly touch the planetesimal surface while large St particles collide with the full headwind speed. Dust particles in turbulent disks experience gas velocity fluctuations and in turn drag variations. They follow preferentially turbulent structures with characteristic timescales that equal their response time. This coupling of inertial particles is known to create nontrivial phenomena such as the mentioned smallscale preferential concentrations that are the most effective for particles (Reade & Collins 2000; Bec et al. 2006).
We observe such a eddydust coupling also for the collision velocity (the norm of the dust velocity vector at impact) of dust particles with a planetesimal (see Fig. 13). Asymptotically, smallSt dust (small in the sense of particles with a small collision efficiency) still only mildly touches the surface, while largeSt dust collides with the speed of the headwind. However, at intermediate values of St, turbulent velocity fluctuations lead to an increase of the collision speed that even exceeds the headwind speed. For the highest turbulent intensity that is studied here (), the average impact speed is approximately 75% higher than the mean headwind speed. We remark that once particle inertia is measured in terms of the characteristic time scale of turbulent structures of size (planetesimal diameter) ( for ) all maxima of the curves (located at ) in Fig. 13 align at . Turbulent eddies of the size of the planetesimal are thus responsible for this increase of the impact velocity.
To estimate the broadness of the velocity distributions we measured in Fig. 14 their standard deviation. Here, the peak is even more important than for the mean impact speed. For a dust particle in a headwind, the velocity distribution is up to six times broader than in a nonturbulent disk.
The impact speed has, besides its average value and standard deviation, a probability distribution that varies with both and St. It reveals that highimpact speeds of the order of several times the headwind speed are quite probable. For a fixed value of St the distribution becomes monotonously broader with increasing (see Fig. 15). The numerical data is there compared to the noncentral chisquared distribution that would be obtained if were the norm of a threedimensional random Gaussian vector with prescribed mean and variance. Up to statistical accuracy, it seems from Fig. 15 that such an approach gives a rather good description of actual fluctuations of the impact velocity. This approximation is however valid only if the Stokes number St is sufficiently large. Figure 16 indeed represents the same distributions for a fixed value of at varying St. One observes deviations from the chisquared prediction in both tails at the smallest value of the Stokes number. It seems nevertheless that the distribution still belongs to the same family and can be approximated by a chisquared distribution with a smaller number of degrees of freedom. Everything happens as at small St, the dust velocity fluctuations with respect to the planetesimal were constrained in a space with dimension less than three.
The outcome of a collision also depends on the angle of impact. Figure 17 shows the average collision angle with respect to the surface normal direction . Small St dust preferentially touches the planetesimal with an mean impact angle close to . High inertia dust heads straight onto the planetesimal and experiences an average impact angle of . Turbulent fluctuations randomize the impact angle and favor this way to the same mean angle of .
5 Astrophysical application
5.1 Linear cross section
We apply these results to realistic disk conditions. In order to study filtering of dust by planetesimals, Guillot et al. (2014) had applied the results of numerical simulations by Sekiya & Takeda (2003) for a laminar disk and a fixed planetesimal Reynolds number . Our simulations in the laminar case cover a range of from to . As seen in Fig. 1 and Eq. (8), this corresponds to planetesimals with diameters between 4 m and 40 m at 1 au and between 1 km and 12 km at 10 au. Since the outcome of the simulations is only weakly dependent on (the critical Stokes number is a function of –see eq. (16)), we expect to be able to extrapolate the results outside this range.
In most cases however, turbulence is expected to be important. Our simulations in the turbulent case have been calculated for various intensities of the turbulence , but a fixed planetesimal Reynolds number . However, when turbulence becomes important, we expect the results to become very weakly dependent on . This is for two reasons: First as seen in Fig. 2, for values of approaching unity, the flow around planetesimals is perturbed very significantly and becomes controlled by the turbulence of the disk instead of by the planetesimal properties. Second, turbulence perturbs the boundary layer around the planetesimal independently of its properties to offer new possibilities for dust particles to impact.
But for small planetesimals and/or low turbulent intensity, the planetesimal size can become smaller than the Kolmogorov scale, i.e., the minimum scale for turbulent eddies. In that case, the planetesimals experience a headwind of variable intensity and direction. It is expected that, in the limit of , the situation becomes similar to the laminar case, but with a headwind that is increased by . We write
(18) 
where is the collision efficiency as defined by Guillot et al. (2014), and and the fitting formula (16) and (17) . With this definition, is the planetesimal linear collisional cross section and its surface collisional cross section. Thus, for a planetesimal smaller than the smallest turbulent eddy, the flow is considered laminar, but we account with the use of for a flow velocity that is slightly higher on average. On the other hand, when the planetesimal is larger than the Kolmogorov scale, we use the results of the simulations in the turbulent case directly.
Figure 18 illustrates how the factor varies as a function of particle size in an MMSN disk for a 1 kmradius planetesimal either at 7.1 au or at 1 au. The first case corresponds to a planetesimal Reynolds number equal to the one used in the hydrodynamical simulations with turbulence. The second case corresponds to a much higher Reynolds number outside the range of our simulations.
In all cases, particles which are larger than a critical value (i.e., with a Stokes number higher than unity) are accreted with a crosssection approximately equal to the geometric one (i.e., ). In laminar disks or when turbulence is small, the cross section drops for small particles such that . If turbulence is large enough, this drop occurs at even smaller sizes, with an offset that corresponds to one to two orders of magnitude for the case with .
The comparison of the laminar cases shows a relatively good agreement between our work and the previous results of Guillot et al. (2014) who used results from Sekiya & Takeda (2003) for a fixed planetesimal Reynolds number of . When turbulence is added, it is worth noticing that while the case with stands out and allows much smaller particles to collide, the cases with and are almost indistinguishable from the laminar case. This is a direct consequence from the fact that for these: the smallest turbulent cell is expected to be larger than the planetesimal size which implies that we switch to the laminar case in eq. (18). Our approach is thus discontinuous in , but resolving this issue would require dedicated simulations beyond the scope of the present work
At high Reynolds number (i.e., the 1 au case in Fig. 18), the Kolmogorov parameter is generally high which implies a nearly continuous behavior from high to low values of . A small issue seen for low values of the viscosity and dust sizes corresponding to Stokes number close to unity is that the value of for a disk with low turbulence (e.g., ) can become smaller than the laminar value, which according to our simulations is unlikely. Clearly, this is a consequence of the fact that our expressions have been derived for a relatively low planetesimal Reynolds number and are applied very far from that value. Again, dedicated simulations would be needed, but may be out of reach with presentday computing power.
5.2 Collision probabilities in disks
We now examine the consequences for collisions of dust grains with planetesimals with the same approach as Guillot et al. (2014). In protoplanetary disks, collisions between drifting dust particles and planetesimals occur with a probability that is a function of the planetesimal cross section, the scale height of the dust disk and of the drift velocity of the dust particles. The latter depends on orbital distance , orbital (keplerian) frequency and stopping time . In the limit that gravitational effects and gas drift may be neglected and for circular orbits, this probability can be shown to write (Guillot et al. 2014):
(19) 
where accounts for hydrodynamical effects discussed previously (the purely geometrical limit is recovered for ).
In reality, gravity becomes important both for median to large planetesimals (kilometer size and more) and for large grains (above meter size) and seriously complicates the picture. Several interaction regimes may be defined as follows (see Ormel & Klahr 2010; Guillot et al. 2014):

The geometric regime, corresponds to the most simple case in which drag, hydrodynamical and gravity effects may be neglected.

We define the hydrodynamical regime as an extension of this regime at small dust sizes when we must account for the deflection of dust grains around planetesimals.

The Safronov regime corresponds to the case when large “dust” (effectively, boulders) which are very weakly affected by gas drag migrate inward so slowly that they feel a gravitational focusing by the planetesimal which increases the collision probability.

In the threebody regime, the gravity fields of the planetesimal and that of the central star must be taken into account.

The settling regime corresponds to the case when gravitational acceleration from the planetesimal and gas drag on the dust particles lead to an enhanced capture probability.
Accounting for the complexity of the problem we thus calculate the collision probability between dust and planetesimals, , from eq. (43) of Guillot et al. (2014), assuming monodisperse size distributions^{2}^{2}2In doing so, we correct for the fact that in Guillot et al. (2014), the 3D collision probability in the hydro mode was overestimated because it neglected the reduction in the vertical cross section, i.e., had been assumed instead of . Because this affected the hydro mode with an already very low collision probability this had negligible effect on the qualitative results., but including from eq. (18). In doing so, we also adopt an important modification stemming from the work of Johansen et al. (2015) and Visser & Ormel (2015): Instead of limiting the extent of the settling regime to when the capture radius is larger than the physical size of the planetesimal as in Guillot et al. (2014), we instead look for solutions of the settling regime equations outside of this range and adopt for the collision probability the maximum of the probabilities obtained in the settling and geometric+hydro regimes. This is important in regions where is extremely low but gravity and gas drag can still affect the trajectories of the dust particles.
Figure 19 shows how varies with dust and planetesimal size for a fixed orbital distance of 1 au. We focus on planetesimals smaller than 100 km and down to 10 m with the caveat that for planetesimals smaller than about 1 km, gas drag should be included. The top panel shows the previous results from Guillot et al. (2014), which correspond to the case of a laminar flow and no extension of the settling regime. The bottom panel shows the results for a turbulent flow with full account for gravity effects even for lowplanetesimal sizes.
The comparison between the top and bottom panels of fig. 19 shows that even a weak planetesimal gravity effectively limits the decrease of the collision probabilities in the extended settling regime for dust smaller than m and planetesimals between one and 100 km. The inclusion of turbulence effects also shifts the hydrodynamic regime to smaller dust sizes. The shift is about one order of magnitude for all planetesimal sizes considered when comparing the results for to those for a laminar disk. Particles of 0.1 mm can hence be accreted relatively efficiently by planetesimals for all the sizes considered. However, smaller particles still end in the hydrodynamical regime with a strongly reduced collision efficiency. For example micronsized dust particles are very inefficiently captured by planetesimals larger than a few kilometers in size.
5.3 The inefficient capture of small dust grains
We now turn to the examination of how (in)efficiently individual small dust grains may have been collected by planetesimals as a function of their sizes and orbital distance. A full model would require studying the size distribution of dust and planetesimals and is beyond the scope of the present paper. However, we can identify the parameter space for which this collection of dust is inefficient by identifying when the collision cross section becomes smaller than of the geometrical one (i.e., corresponding to in the limit when gravity effects are not important). Because the drop in collision probability in the hydro region of fig. 19 is quite abrupt, we expect that if dust is indeed collected individually by planetesimals, this process should leave its imprint on the size distribution of individual grains in meteorites.
Figure 20 identifies these regions as a function of orbital distance and planetesimal size, either in the case of a high turbulence level (top panel) or a low turbulence level (bottom panel). In both cases, the collection of very small particles (nanometer sizes) is found to be very inefficient, at least inside of 10 au. Dust particles of progressively larger sizes can be collected up to shorter orbital distances, but the efficiency then strongly depends on the turbulence level.
For 1 mm particles (corresponding to a typical size of chondrules), we do not see in fig. 20 a region of strongly inefficient collection when turbulence is high (), but for , these particles avoid collisions with planetesimals between about 0.3 and 30 km within a fraction of an au. Onemicron particles are collected inefficiently inside a region extending from about 0.1 to 3 au depending on planetesimal size for , but this region extends from 0.8 to 6 au for . Smaller particles can collide with planetesimals only at larger orbital distances, when the gas density has decreased and the stopping time increased for a given particle size.
A larger turbulence level can therefore allow collisions of smallsize particles which would otherwise be avoided due to the hydrodynamical flow around the planetesimals. However, this is also balanced by the fact that higher turbulence means a thicker dust subdisk which lowers the collision probability (see Guillot et al. 2014). Due to the form of eq. 18, the effect of turbulence becomes weaker at large orbital distances, when the size of the smallest turbulent eddies becomes larger than the planetesimals. This thus explains why the contour lines for very small dust particles are identical for the and cases.
The change in behavior of the contour plots for , dust sizes between 1 m and mm and orbital distances from 0.05 to 0.3 au is due to a change in drag behavior for these particles: At short orbital distances, the gas density is so high that they are in the Stokes regimes and they switch to an Epstein drag beyond about 1 au.
For the planetesimal sizes considered, particles smaller than about 10 nm have a collision probability that is independent of alpha. This is because collisions can occur only in the outer disk where , i.e., the smallest Kolmogorov scale is still larger than the planetesimals considered.
Small particles such as the presolar grains present in meteorites, which can have sizes of only a few nanometers (e.g., Clayton & Nittler 2004) must have either collided with planetesimals far out in the disk or be incorporated into larger grains which would have themselves collided with planetesimals (e.g., Ormel et al. 2008). For some of the presolar grains, given their very low abundance, it remains possible that they were incorporated directly in planetesimals, although this would have to be quantified. In any case, this should have occurred without leading to any melting or dissociation of these grains which kept their identify throughout.
6 Conclusions
We have derived the accretion probability of small particles by a planetesimal in a turbulent gas. In order to do so, we performed highresolution hydrodynamical simulations of the flow around a spherical planetesimal of diameter moving with a velocity , assuming incompressibility. We studied both the case of a laminar flow and that of a turbulent one, the intensity of the turbulence being related to the turbulent viscosity of the disk. Dust particles of variable size were implemented in the flow to determine collision rates.
For laminar flows, we confirm that small particles with a Stokes number (corresponding to stopping times shorter than the time to cross the planetesimal) see a significant drop in their collision rate with the planetesimal. For turbulent flows however, this drop occurs for sizes that can be significantly smaller, i.e., turbulence helps accreting dust particles with sizes up to one to two order of magnitudes smaller than for laminar disks.
We thus derived collision probabilities both in the laminar case [eq. (16)] and in the turbulent case [eq. (17)]. These expressions, even if limited to limited to , can be used for a wide range of situations. We propose an approximate recipe to use either the laminar case if the planetesimal size is smaller than the Kolmogorov scale and the turbulent case otherwise [eq. (18)].
When applied to real disks, our new expressions shift the boundary with the hydro regime – where accretion rates are greatly suppressed – to smaller sizes. For example, for , the upper limit dust size in the hydrodynamical regime is decreased by a factor 100 and even subm size particles collide efficiently with onekilometer planetesimals. They also show that the accretion of extremely small particles is difficult and generally requires to be done by small planetesimals (less than km size) at large orbital distances (beyond 1 au) and/or late in time, when the disk has become less massive. We believe that these results are important to interpret, among other things, the presence and characteristics of presolar grains in meteorites since these vary in size from several microns down to only a few nanometers.
In order to apply our results to protoplanetary disks, we had to approximate the effect of gravity, often by extrapolations far from the regime in which numerical experiments were conducted. Future efforts will be directed towards including the gravitational force directly in our hydrodynamical simulations.
Appendix A Definition of symbols
We summarize here the symbols used in this paper and their definitions.
Symbol  Description  Equation 

average headwind velocity  
turbulent intensity  
planetesimal Reynolds number  
outer gas flow Reynolds number  
rootmeansquare velocity of the gas flow  
gas velocity  
kinetic energy of the gas flow  
mean kinetic energy dissipation rate  
kinematic viscosity  
planetesimal diameter  
planetesimal radius  
Kolmogorov length scale  
Kolmogorov time scale  
integral scale  
integral time scale  
response or stopping time of dust  
Stokes number  
critical Stokes number  
dust collision speed  
disk turbulence parameter  
disk turbulent viscosity  
disk surface density  
gas temperature  
speed of sound  
orbital distance in the disk  
vertical height in the disk  
disk scale height  
disk scale height for the dust  
orbital (Keplerian) frequency  
stellar mass 
Acknowledgements.
We thank Satoshi Okuzumi, Zoe Leinhart and Rico Visser for useful discussions. Most of the simulations were done using HPC resources from GENCIIDRIS (Grant i2011026174). Part of them were performed on the “Mésocentre de calcul SIGAMM” at the Observatoire de la Côte d’Azur. The research leading to these results has received funding from the Agence Nationale de la Recherche (Programme Blanc ANR12BS0901104).References
 Balachandar & Eaton (2010) Balachandar, S. & Eaton, J. K. 2010, Annual Review of Fluid Mechanics, 42, 111
 Balbus & Hawley (1991) Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214
 Bec et al. (2006) Bec, J., Biferale, L., Boffetta, G., et al. 2006, Journal of Fluid Mechanics, 550, 349
 Berthet et al. (2010) Berthet, S., Leriche, M., Pinty, J.P., Cuesta, J., & Pigeon, G. 2010, Atmospheric Research, 96, 325
 Blum & Wurm (2008) Blum, J. & Wurm, G. 2008, ARA&A, 46, 21
 Chapman & Cowling (1970) Chapman, S. & Cowling, T. G. 1970, The mathematical theory of nonuniform gases. an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases (Cambridge University Press)
 Chatterjee & Tan (2014) Chatterjee, S. & Tan, J. C. 2014, The Astrophysical Journal, 780, 53
 Chatterjee & Tan (2015) Chatterjee, S. & Tan, J. C. 2015, The Astrophysical Journal Letters, 798, L32
 Cisse et al. (2013) Cisse, M., Homann, H., & Bec, J. 2013, Journal of Fluid Mechanics, 735, R1
 Clayton & Nittler (2004) Clayton, D. D. & Nittler, L. R. 2004, ARA&A, 42, 39
 Cuzzi et al. (2001) Cuzzi, J. N., Hogan, R. C., Paque, J. M., & Dobrovolskis, A. R. 2001, ApJ, 546, 496
 Fadlun et al. (2000) Fadlun, E. A., Verzicco, R., Orlandi, P., & Mohd, J. 2000, Journal of Computational Physics, 60, 35
 Gatignol (1983) Gatignol, R. 1983, J. Méc. Théor. Appl., 1, 143
 Guillot et al. (2014) Guillot, T., Ida, S., & Ormel, C. W. 2014, A&A, 572, A72
 Haugen & Kragset (2010) Haugen, N. E. L. & Kragset, S. 2010, Journ. Fluid Mech., 661, 239
 Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
 Hayashi et al. (1985) Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, ed. D. Black & M. Matthews (University of Arizona Press, Tucson, AZ), 1100–1153
 Homann & Bec (2010) Homann, H. & Bec, J. 2010, J. Fluid Mech., 651, 81
 Homann & Bec (2015) Homann, H. & Bec, J. 2015, Physics of Fluids, 27, 053301
 Homann et al. (2013) Homann, H., J.Bec, & Grauer, R. 2013, J. Fluid Mech, 721, 155
 Hydle Rivedal et al. (2011) Hydle Rivedal, N., Granskogen Bjørnstad, A., & Haugen, N. E. L. 2011, ArXiv eprints
 Ingham et al. (1990) Ingham, D., Hildyard, L., & Hildyard, M. 1990, Journal of Aerosol Science, 21, 935
 Johansen et al. (2011) Johansen, A., Klahr, H., & Henning, T. 2011, Astronomy & Astrophysics, 529, A62
 Johansen et al. (2015) Johansen, A., Mac Low, M.M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 15109
 Johansen et al. (2007) Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022
 Kolmogorov (1941) Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301
 Lambrechts & Johansen (2012) Lambrechts, M. & Johansen, A. 2012, Astronomy & Astrophysics, 544, A32
 Maxey & Riley (1983) Maxey, M. R. & Riley, J. 1983, Phys. Fluids, 26, 883
 Michael & Norey (1970) Michael, D. H. & Norey, P. W. 1970, Canadian Journal of Physics, 48, 1607
 Mitra et al. (2013) Mitra, D., Wettlaufer, J. S., & Brandenburg, A. 2013, The Astrophysical Journal, 773, 120
 Mittal & Iaccarino (2005) Mittal, R. & Iaccarino, G. 2005, Annual Review of Fluid Mechanics, 37, 239
 Morbidelli et al. (2015) Morbidelli, A., Lambrechts, M., Jacobson, S., & Bitsch, B. 2015, Icarus, 258, 418
 Nakagawa et al. (1986) Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
 Ormel & Klahr (2010) Ormel, C. & Klahr, H. 2010, Astronomy & Astrophysics, 520, A43
 Ormel & Kobayashi (2012) Ormel, C. & Kobayashi, H. 2012, The Astrophysical Journal, 747, 115
 Ormel et al. (2008) Ormel, C. W., Cuzzi, J. N., & Tielens, A. G. G. M. 2008, ApJ, 679, 1588
 Ormel & Okuzumi (2013) Ormel, C. W. & Okuzumi, S. 2013, ApJ, 771, 44
 Pekurovsky (2012) Pekurovsky, D. 2012, SIAM Journal on Scientific Computing, 34, C192
 Peskin (1972) Peskin, C. S. 1972, Journal of Computational Physics, 10, 252
 Phillips & Kaye (1999) Phillips, C. & Kaye, S. 1999, Journal of Aerosol Science, 30, 709
 Reade & Collins (2000) Reade, W. C. & Collins, L. R. 2000, Phys. Fluids, 12, 2530
 Safronov (1972) Safronov, V. S. 1972, Evolution of the protoplanetary cloud and formation of the earth and planets. (Israel Program for Scientific Translations, Keter Publishing House)
 Sekiya & Takeda (2003) Sekiya, M. & Takeda, H. 2003, Earth, Planets, and Space, 55, 263
 Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
 Shaw (2003) Shaw, R. A. 2003, Annual Review of Fluid Mechanics, 35, 183
 Slinn (1983) Slinn, W. 1983, Atmospheric Sciences and Power Production, Chap.11
 Slinn (1974) Slinn, W. G. N. 1974, Proceedings of the USAEC Symposium
 Slinn (1976) Slinn, W. G. N. 1976, Geophys. Res. Lett., 3, 21
 Squires & Eaton (1991) Squires, K. & Eaton, J. 1991, Phys. Fluids A, 3, 1169
 Taylor (1940) Taylor, G. I. 1940, Collected Works G.I. Taylor, 3, 236
 Turner et al. (2014) Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI, 411
 Uhlmann (2005) Uhlmann, M. 2005, Journal of Computational Physics, 209, 448
 Visser & Ormel (2015) Visser, R. G. & Ormel, C. W. 2015, ArXiv eprints
 Weidenschilling (1977a) Weidenschilling, S. J. 1977a, MNRAS, 180, 57
 Weidenschilling (1977b) Weidenschilling, S. J. 1977b, Ap&SS, 51, 153
 Weidenschilling (1984) Weidenschilling, S. J. 1984, Icarus, 60, 553
 Windmark, F. et al. (2012) Windmark, F., Birnstiel, T., Güttler, C., et al. 2012, A&A, 540, A73
 Youdin & Goodman (2005) Youdin, A. N. & Goodman, J. 2005, ApJ, 620, 459