© ESO, 2016

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 sub-μm dust grains, which follow the sub-Keplerian gas, to km-sized planetesimals that are massive enough to move on Keplerian orbits. One of them is known as the “meter-size barrier”. Meter-sized 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 yr at 1 AU (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 planet-forming material are lost from the disk.

One scenario for overcoming the meter-size barrier is by the combined effect of streaming instabilities and pebble accretion. Because of growth to mm/cm-sized 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 close-in super-Earth systems in exoplanetary systems (Chatterjee & Tan 2014, 2015). However, the details of 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 (head-on) 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 T-Tauri disks is far higher than that caused by molecular viscosity. Consequently, protoplanetary disks are expected to be turbulent, with a turbulent eddy viscosity ν_t 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 two-dimensional 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 three-dimensional 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 magneto-rotational 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 two-dimensional 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 Sect. 2, we explain the disk model and notations. In Sect. 3, we describe the numerical method that we use in our approach. In Sect. 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 and 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: Σ₁ and T₁ are the values at 1 au, for which we choose Σ₁ = 1700 g cm^-2 and T₁ = 270 K. (A list of symbols and their definitions used in this paper is summarized in Appendix 6.) Assuming the disk is vertically isothermal, with the isothermal sound speed and a mean molecular weight of μ = 2.34m_H, the density reads (3)where H = c_s/ Ω_K is the disk scale height and the orbital (Keplerian) frequency at distance r from the star. It follows that H/r ∝ r^1/4, 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 Δg = (dP/ dr) /ρ compared to the solar gravity is often expressed in terms of a parameter η defined as (4)where we introduced the pressure logarithmic gradient ∇_logP ≡ −d(log P)/d(log r), the Keplerian velocity v_K ≡ Ω_Kr, and the disk aspect ratio h ≡ c_s/v_K = H/r. 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 (1−η) v_K. Also, the headwind that is faced by a big body (a.k.a. planetesimal), which moves at a Keplerian speed through the sub-Keplerian rotating gas, then reads (5)which, for an MMSN disk is independent of the distance r to the star.

2.2. Turbulence model

The gas has a mean-free path of where σ_mol ≈ 2 × 10^-15 cm² (Chapman & Cowling 1970) is the molecular cross section. The kinematic molecular viscosity of the gas is ν_mol = (1/2) v_thℓ_mfp where is the mean molecular thermal speed.

However, for these parameters the resulting diffusion timescales ~r²/ν_mol are too long to explain, for example, the measured accretion rate in T-Tauri disks. Consequently, protoplanetary disks are expected to be mildly turbulent, but still subsonic, with a turbulent eddy viscosity ν_t assumed to be much larger than the molecular kinematic viscosity. An often used parameterization is the alpha-viscosity model of Shakura & Sunyaev (1973): (6)The dimensionless constant α is a disk-dependent parameter. Its value may range from a minimum of α ≈ 10^-5, when the turbulence originates from the Kelvin-Helmholtz instability of a dust layer at the mid-plane, to perhaps α ≈ 0.1 for the most violent disks with a strong magneto-rotational 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, α ≪ 1 – the Reynolds number of the flow, Re = ν_t/ν_mol, is very large by virtue of the very low densities that characterize astrophysical environments. Indeed, for a 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 three-dimensional turbulent regime which can be described using Kolmogorov (1941) phenomenology. The kinetic energy is injected at a large (integral) scale L by a hydro- or magnetohydrodynamical instability. In both cases the injection mechanism is associated to the sub-Keplerian shear. The associated shear rate sets the typical turbulent large-eddy turnover time to (Cuzzi et al. 2001). The eddy viscosity then reads ν_t = L²/t_L = L² Ω_K leading, together with Eq. (6), to L = α^1/2H and v_L ≡ L/t_L = α^1/2c_s. Note that the assumptions of incompressibility (v_L ≪ c_s) and three-dimensionality (L ≪ H) of the turbulent flow in the disk are both fulfilled as long as α ≪ 1. According to Kolmogorov (1941) phenomenology, the kinetic energy cascades downscale with a constant rate until it reaches the smallest active scales, the Kolmogorov scale ℓ_Kol below which it is dissipated by molecular viscosity. The typical velocity v_ℓ of eddies whose size ℓ lies in the inertial range ℓ_Kol ≪ ℓ ≪ L reads v_ℓ ≃ (ε_Kolℓ)^1/3. The Kolmogorov length is defined as the scale where the scale-dependent Reynolds number Re(ℓ) ≡ v_ℓℓ/ν_mol is unity, so that and the associated timescale . The integral Reynolds number Re = ν_t/ν_mol gives the extension of the inertial range: L/ℓ_Kol = Re^3/4 and t_L/t_Kol = Re^1/2.

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 R_p, diameter d = 2R_p. 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 d, so that all turbulent eddies of sizes between ℓ_Kol and d 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 Re_p decreases with increasing orbital distance (kinematic viscosity goes up). Similarly, d/ℓ_Kol decreases because ℓ_Kol 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, Re_p = 400 corresponds to a 1 km planetesimal located at ~5 au in the MMSN disk, or e.g. to a 10 km planetesimal at ~12 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 ~10⁶ mismatch between the Reynolds number in the disk and the maximum one that can be attained by present-day 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 no-slip boundary conditions at its surface. The overall method consists in a combination of a standard pseudo-Fourier-spectral solver with a penalty method that is explained in this section.

We integrate the incompressible Navier-Stokes equations (11)for the gas velocity u, where ρ_g is the gas density, ν_mol its kinematic molecular viscosity and f a force. The latter maintains a uniform inflow speed and an eventually ambient turbulent flow. Its form is described later. The pseudo-Fourier-spectral approach consists in computing spatial derivatives in Fourier-space and convolutions arising from the non-linear terms in real space. A Fast-Fourier 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 U_c relative to the gas. Equation (11) is associated with a no-slip boundary condition at the surface of the planetesimal, which is assimilated to a spherical object at rest whose diameter is denoted by d and its center by X_S. We thus have (12)Numerically, this no-slip 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 x such that | x − X_S | ≤ d/ 2. The Navier-Stokes Eq. (11) then has to be modified by introducing in its right-hand side a penalty force f_b(x,t), 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 L(u) denotes the right hand side of (11). In order to compute f_b 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 a second order finite-difference Navier-Stokes solver to simulate turbulent suspensions involving many particles.

3.2. Dust particles

Dust particles are modeled by spherical inertial particles with a radius a 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 ρ_p is much higher than the gas density ρ_g. 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. t_s 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 back-reaction on the flow.

Usually, particle inertia is quantified in terms of the Stokes number St = t_s/t_c defined by non-dimensionalizing their response time by a characteristic time scale t_c of the carrier flow. The present problem involves different relevant time scales. T_d = d/U_c, the time it takes a dust particle to pass the planetesimal, t_L, the turbulent large-eddy turn-over time and t_Kol, the turbulent dissipation time scale. If not otherwise specified, we use t_c = T_d 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¹.

In studies concerned with the small-scale dynamics of inertial particles one usually uses t_c = t_Kol (Bec et al. 2006) and we denote the associated Stokes number St_Kol.

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 three-dimensional simulations. The planetesimal experiences a headwind speed U_c 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 for 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 St = 0.8 for all cases shown in Fig. 2, St_Kol = 1.25 for I = 0.29 and St_Kol = 9.8 for I = 1.18 so that the clustering properties of the particles change from one value of I 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 f. Once a turbulent flow has developed (after approx. 1−2 t_L) the velocity field is forced by keeping constant the energy content of the two lowest wave number shells (1 ≤ k ≤ 2) in Fourier space. This leads to a statistically stationary turbulent flow to which the planetesimal is added. For this, the root-mean-square value v_L of the velocity fluctuations is normalized to the values given in Table 2 for each simulation. A mean velocity of U_c 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 (>t_s) 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 (U_c = 1) 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 re-injected 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 Runge-Kutta 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 Re_p = U_cd/ν_mol, the Reynolds number of the gas Re = v_LL/ν_mol, L being the integral scale of the ambient turbulent flow, and the turbulent large-scale intensity I = v_L/U_c. The latter measuring the strength of the large-scale turbulent fluctuations compared to the mean flow velocity. In the laminar case (Re = 0, I = 0), we varied the planetesimal Reynolds number Re_p from 100 to 1000. We analyze the effect of turbulent fluctuations for one specific choice of the planetesimal Reynolds number, namely Re_p = 400 (turbulent wake) and vary I from 0.14 to 1.18. The corresponding flow Reynolds numbers Re (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 L but only v_L that enters in the definitions of both I and Re.

The particle dynamics is characterized by the Stokes number St. In all simulations we consider streams of heavy dust particles with response times t_s in between 0.04 and 81.92, corresponding to Stokes numbers St in the range 0.05 ≤ St ≤ 63. The main parameters of all simulations are summarized in Tables 2 and 1.

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 Re_p 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 non-collisional 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 St the distribution becomes 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) E(St) that is the ratio between the actual collision rate and that expected for free-streaming particles that do not feel any hydrodynamic force from the gas flow. The latter rate is obtained from the geometrical cross-section of the planetesimal and reads n_dU_cπd²/ 4, where n_d designates the number density of dust particles. We use the following notation: E₀ denotes the uniform gas flow case that we are going to analyze in this section and E_I denotes the turbulent case parameterized by the turbulent intensity I that is discussed in the following section. The efficiency E₀ 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 Re_p the more probable are collisions, especially at small St. In this low inertia limit, E₀ sharply drops for all Re_p. The reason for this is the well established existence of a critical Stokes number St_c below which no collisions occur at all (Taylor 1940).

Indeed, in the case of an inviscid Euler flow (zero viscosity, Re_p = ∞) it can be shown that E₀ drops to zero for small but finite (Ingham et al. 1990). The viscous boundary layer of a finite Re_p planetesimal that shrinks as ~ reduces the collision probability and increases St_c. For the limiting case of a Stokes flow Michael & Norey (1970) computed numerically a critical Stokes number of 0.605.

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 0.1 on the determination of E₀(St) (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 St − St_c and not simply as a function of St. We hence focus on the behavior of the collision probability when approaching the critical Stokes number. St_c 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 St_c are close to that of Slinn (15) (see inset of Fig. 6). The case of large but finite values of Re_p 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.

E₀ increases linearly at small St − St_c (see Fig. 6), that is in contradiction to Slinn’s formula (15) that predicts . To incorporate this small-Stokes 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 exp(−a/St) 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 3/2 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 large-Stokes number behavior in details it is useful to analyze 1 − E₀(St) that measures how the free-streaming particle limit is recovered as a function of St. All curves for different Re_p (see Fig. 7) fall on the top of each other and reveal a St^-1 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 Re_p and the dust Stokes number St, but turbulence adds two additional dimensionless parameters specifying the ambient turbulent flow. They are the gas Reynolds number Re = v_LL/ν_mol and the turbulent intensity I = v_L/U_c 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 I and fix the planetesimal Reynolds number to Re_p = 400 (weakly turbulent wake) in order to limit computational costs. The gas Reynolds number Re is implicitly varied as it is coupled to I due to our specific forcing scheme that prescribes L 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 Figs. 2b 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 concentrations Shaw (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 E_I for various I. One observes that higher is the turbulent intensity I, 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 St ≈ 1 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 ΔE(St,I) = E_I/E₀ − 1 of the turbulent efficiency compared with the laminar one is shown in Fig. 9. It is represented as a function of St − St_c that is the distance from the critical Stokes number St_c as ΔE diverges when St → St_c. The large values of ΔE at small St decreases as a power law with an exponent close to −1. The dependence of ΔE on I is nearly quadratic. Indeed, the different curves almost fall on the top of each other when normalized by I^1.7 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 a and b. The dependence of these parameters on the turbulent intensity is shown in Fig. 10 together with simple fitting formulas. The parameter a is decreasing roughly exponentially as a function of I. b is close to zero for small intensities and strongly increasing starting from I ≈ 0.4.

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 I fall approximately on top of each other once shown as a function of St_L = t_sv_L/L. It is thus a turbulent time scale, namely the large-eddy turn-over time, that replaces the advection time d/U_c 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 I and St and allow for collisions of dust particles beyond the critical Stokes number St_c. In turbulent gas flows, the dust heads onto the planetesimal with a fluctuating velocity. The probability distribution of the headwind is given by the one-point 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 I. These headwind variations lead to a fluctuating planetesimal Reynolds number and to fluctuating dust Stokes numbers. Especially close to the laminar value St_c this results in higher collision probabilities and an absence of a critical Stokes number. Another consequence of a variable headwind are fluctuations of the stream-wise velocity gradient σ in the upstream boundary layer of the planetesimal. They modify the local Stokes number that can be defined by St = t_sσ. The probability for a collision of a dust particle with stopping time t_s is then P(St>St_c) = P(σ>St_c/t_s), 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 P(σ>St_c/t_s) ~ exp(−CSt_c/t_s) which is consistent with the formally observed exponential fall off of E_I(St) 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 large-scale turbulent fluctuations. Let us assume (in the reference frame of the sphere) a velocity of the form u = U_ce_x + I sin(2π/t_L,t) e_y. The volume swept by the sphere is evidently larger in the turbulent case (I> 0) than in the laminar (I = 0). 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 (I = 0) 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 I 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 et al. 2012).

In laminar disks the mean impact speed of small-size 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 time-scales that equal their response time. This coupling of inertial particles is known to create non-trivial phenomena such as the mentioned small-scale preferential concentrations that are the most effective for St_Kol = t_s/t_Kol ≈ 0.6 particles (Reade & Collins 2000; Bec et al. 2006).

We observe such a eddy-dust coupling also for the collision velocity v_c (the norm of the dust velocity vector at impact) of dust particles with a planetesimal (see Fig. 13). Asymptotically, small-St dust (small in the sense of particles with a small collision efficiency) still only mildly touches the surface, while large-St 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 (I = 1.18), 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 d (planetesimal diameter) t_d = t_L (d/L)^2/3 (t_d = 4.4,2.1,1.1 for I = 0.29,0.6,1.18) all maxima of the curves (located at St ≈ 3.2,1.6,0.8) in Fig. 13 align at St = t_s/t_d ≈ 0.6. 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 St ≈ 1 peak is even more important than for the mean impact speed. For a St = 0.4 dust particle in a I = 1.18 headwind, the velocity distribution is up to six times broader than in a non-turbulent disk.

The impact speed has, besides its average value and standard deviation, a probability distribution that varies with both I and St. It reveals that high-impact 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 I (see Fig. 15). The numerical data is there compared to the non-central chi-squared distribution that would be obtained if v_c were the norm of a three-dimensional 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 I at varying St. One observes deviations from the chi-squared 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 chi-squared 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 n. Small St dust preferentially touches the planetesimal with an mean impact angle close to 90°. High inertia dust heads straight onto the planetesimal and experiences an average impact angle of 45°. 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 Re_p = 50. Our simulations in the laminar case cover a range of Re_p from 100 to 1000. 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 Re_p (the critical Stokes number is a function of log (1 + Re_p/ 2) – 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 I, but a fixed planetesimal Reynolds number Re_p = 400. However, when turbulence becomes important, we expect the results to become very weakly dependent on Re_p. This is for two reasons: First as seen in Fig. 2, for values of I 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 d/ℓ_Kol ≪ 1, the situation becomes similar to the laminar case, but with a headwind that is increased by . We write (18)where f_hydro is the collision efficiency as defined by Guillot et al. (2014), .and E₀ and E_I the fitting formula (16) and (17). With this definition, 2R_pf_hydro 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 St^∗ 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 f_hydro varies as a function of particle size in an MMSN disk for a 1 km-radius planetesimal either at 7.1 au or at 1 au. The first case corresponds to a planetesimal Reynolds number Re_p = 400 equal to the one used in the hydrodynamical simulations with turbulence. The second case corresponds to a much higher Reynolds number Re_p = 5.4 × 10⁴ 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 cross-section approximately equal to the geometric one (i.e., f_hydro ≈ 1). In laminar disks or when turbulence is small, the cross section drops for small particles such that St< 1. 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 α = 10^-2.

The comparison of the laminar Re_p = 400 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 50. When turbulence is added, it is worth noticing that while the case with α = 10^-2 stands out and allows much smaller particles to collide, the cases with α = 10^-4 and 10^-6 are almost indistinguishable from the laminar case. This is a direct consequence from the fact that δ_Kol< 1 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 δ_Kol 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 f_hydro for a disk with low turbulence (e.g., α = 10^-6) 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 present-day 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 h_d and of the drift velocity of the dust particles. The latter depends on orbital distance r, orbital (Keplerian) frequency Ω and stopping time t_s. 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 f_hydro accounts for hydrodynamical effects discussed previously (the purely geometrical limit is recovered for f_hydro = 1).

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 three-body 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², but including f_hydro from Eq. (18). In doing so, we also adopt an important modification stemming from the work of Johansen et al. (2015) and Visser & Ormel (2016): 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 f_hydro 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 low-planetesimal 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 ~100μm and planetesimals between one and 100 km. The inclusion of turbulence effects alsoshifts 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 α = 10^-2 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 micron-sized 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 1% of the geometrical one (i.e., corresponding to f_hydro< 0.1 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 (α = 10^-2), but for α = 10^-4, these particles avoid collisions with planetesimals between about 0.3 and 30 km within a fraction of an au. One-micron particles are collected inefficiently inside a region extending from about 0.1 to 3 au depending on planetesimal size for α = 10^-2, but this region extends from 0.8 to 6 au for α = 10^-4. 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 small-size 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 α = 10^-2 and α = 10^-4 cases.

The change in behavior of the contour plots for α = 10^-2, dust sizes between 1 μm and 1 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 δ_Kol< 1, 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 identity 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 high-resolution hydrodynamical simulations of the flow around a spherical planetesimal of diameter d moving with a velocity U_c, 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 St< 1 (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 Re_p = 400, 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 α = 10^-2, the upper limit dust size in the hydrodynamical regime is decreased by a factor 100 and even sub-μm size particles collide efficiently with one-kilometer 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.

---

Acknowledgments

We thank Satoshi Okuzumi, Zoe Leinhart and Rico Visser for useful discussions. Most of the simulations were done using HPC resources from GENCI-IDRIS (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 ANR-12-BS09-011-04).

References

Appendix A: Definition of symbols

We summarize here the symbols used in this paper and their definitions.

All Tables

All Figures