Dust Coagulation Regulated by Turbulent Clustering in Protoplanetary Disks

The coagulation of dust particles is a key process in planetesimal formation. However, the radial drift and bouncing barriers are not completely resolved, especially for silicate dust. Since the collision velocities of dust particles are regulated by turbulence in a protoplanetary disk, the turbulent clustering should be properly treated. To that end, direct numerical simulations (DNSs) of the Navier Stokes equations are requisite. In a series of papers, Pan&Padoan used a DNS with the Reynolds number Re~1000. Here, we perform DNSs with up to Re=16100, which allow us to track the motion of particles with Stokes numbers of 0.01<~St<~0.2 in the inertial range. By the DNSs, we confirm that the rms relative velocity of particle pairs is smaller by more than a factor of two, compared to those by Ormel&Cuzzi (2007). The distributions of the radial relative velocities are highly non-Gaussian. The results are almost consistent with those by Pan&Padoan or Pan et al. at low-Re. Also, we find that the sticking rates for equal-sized particles are much higher than those for different-sized particles. Even in the strong-turbulence case with alpha-viscosity of 10^{-2}, the sticking rates are as high as>~50% and the bouncing probabilities are as low as ~10% for equal-sized particles of St<~0.01. Thus, the turbulent clustering plays a significant role for the growth of cm-sized compact aggregates (pebbles) and also enhances the solid abundance, which may lead to the streaming instability in a disk.


INTRODUCTION
Planetesimals are thought to be the precursors of both earth-like planets as well as the cores of gas giants and ice giants. The formation of planetesimals has been a longstanding and perplexing obstacle toward the full understanding of the origins of planets. Planetesimals are widely believed to form as a consequence of the hierarchical coagulation from submicronsize dust particles to kilometer-size bodies in protoplanetary disks (e.g. Lissauer 1993;Chiang & Youdin 2010;Johansen et al. 2014). The growth, however, faces several difficulties, such as bouncing, fragmentation, and radial drift barriers. As the gas disk is partially pressure-supported in the radial direction, the gas rotates with a sub-Keplerian velocity. The resultant friction by the gas forces dust particles to drift toward the central star. In particular, the centimeterto meter-size particles have large drift velocities and may rapidly fall into the central star on a timescale of a few 100 yr (Adachi et al. 1976;Weidenschilling 1977), which is often referred to as the meter-size problem. Therefore, the growth timescale is required to be shorter than the drift timescale for the successful formation of planetesimals. However, as the dust particles grow, they become less sticky, and the high-velocity collisions may lead to bouncing or fragmentation instead of coagulation (Blum & Wurm 2008;Brauer et al. 2008a,b;Birnstiel et al. 2009Birnstiel et al. , 2012Güttler et al. 2010;Zsom et al. 2010Zsom et al. , 2011Windmark et al. 2012).
The streaming instability is a potential mechanism to circumvent the radial drift barriers (Youdin & Goodman 2005).
In order for this instability to work successfully, the formation of cm-sized compact aggregates (pebbles) and the enhancement of solid abundance in a protoplanetary disc are crucial (Johansen et al. 2009;Bai & Stone 2010;Johansen et al. 2014;Carrera et al. 2015;Ida & Guillot 2016;Yang et al. 2017). On the other hand, in recent years, N -body molecular-dynamics simulations have shown that fluffy dust aggregates have the potential to overcome these bouncing or fragmentation barriers (Wada et al. 2009;Paszun & Dominik 2009;Wada et al. 2013;Meru et al. 2013;Seizinger & Kley 2013;Gunkelmann et al. 2016). Wada et al. (2013) derived the critical collision velocity u c , below which fluffy dust aggregates can coalesce: u c ≃ 60 − 80 m s −1 for ice dust and u c ≃ 6 − 8 m s −1 for silicate dust. This implies that icy aggregates overcome the fragmentation barriers more easily, compared to silicate aggregates.
Protoplanetary disks inevitably become turbulent due to the differential rotation, and therefore the collision velocities are regulated by the turbulent motion in a wide range of particle sizes, from 1 mm to 10 m (e.g., a review by Johansen et al. 2014). Völk et al. (1980) built up a framework set (the Völk-type model) based on a Langevin approach for the nonlinear response of dust particles to turbulent eddy motion, which was further developed by Markiewicz et al. (1991). Ormel & Cuzzi (2007)(OC07) provided closed-form analytical approximations to the Völk-type model. Okuzumi et al. (2012), based on an analytic formula by OC07, have simulated the growth of fluffy icy aggregates outside the snow line, taking into account the change in aggregate porosities, and found that the porosity evolution enables the icy aggregates to grow across the radial drift barrier. Kataoka et al. (2013) have explored the compression of fluffy aggregates and shown that the aggregates can evolve into compact icy planetesimals. The abovementioned works revealed that icy planetesimals can form in a wide range beyond the snowline in protoplanetary disks. However, the difficulties for silicate planetesimal formation are difficult to alleviate, because the critical collision velocity for silicate dust is smaller by an order of magnitude than that for icy dust.
The actual sticking rates of dust particles are strongly dependent on the probability distribution function (PDF) of the collision velocities. Even if the root mean square (rms) collision velocity exceeds the critical value, a subset of the dust particles can have collision velocities lower than the critical one and eventually grow, evading the fragmentation barrier. Employing the rms collision velocity derived by OC07, Windmark et al. (2012) and Garaud et al. (2013) have explored the effects of the PDF on the collisional dust growth barriers, under the assumption of a Gaussian (Maxwellian) distribution. However, the PDFs of turbulence-induced relative velocities are found to be highly non-Gaussian by numerical, experimental, and theoretical studies (Sundaram & Collins 1997;Gustavsson et al. 2008;Gustavsson & Mehlig 2011;Hubbard 2012). In addition, it is known that particles with small inertia preferentially concentrate in low-vorticity, high-strain regions during turbulence due to the centrifugal mechanism of the vorticity (Maxey 1987;Squires & Eaton 1991;Fessler et al. 1994). Furthermore, effects such as "caustics" (Wilkinson et al. 2006) and the "sling effect" (Falkovich & Pumir 2007) allow the particles with large inertia to become less coupled with the local fluid ve-locity field and assemble from different regions (e.g., see Bragg & Collins 2014). Such "turbulent clustering" effects are significant when considering the process of the collisional coagulation of dust particles in protoplanetary disks.
To treat the turbulent clustering properly, the direct numerical simulation (DNS) of the Navier-Stokes equations coupled with tracking dust particle motions are requisite. In the DNS, the smallest eddies in the turbulence are resolved without introducing numerical viscosity and turbulence models. Pan et al. (2011) handled the collision statistics with the turbulent clustering, using an Eulerian formulation instead of the Navier-Stokes equations. Then, in a series of the papers by Pan & Padoan (2013) Pan et al. (2014a,b); Pan & Padoan (2014), and Pan & Padoan (2015)(PP15), they used a DNS of the Navier-Stokes equations in the context of planetesimal formation. By analyzing the DNS data, they studied the statistics of colliding dust grains including the radial relative velocity, its probability distribution and the collision rate between dust grains. However, Reynolds number dependence of the results has not been investigated yet.
The motion of particles in turbulence is characterized by the Stokes number given by St = Ωτ p , where Ω is the Keplerian frequency at a radial distance and τ p is the stopping time by the gas friction. In addition, another Stokes number can be defined by the turnover timescale τ η of the smallest eddies of turbulence as St η = τ p /τ η . The St η is determined by the resolution of the simulations, and has a relation St η /St ∝ Re 1/2 , where Re is the Reynolds number. According to Kolmogorov theory, as the Reynolds number increases, the inertial range of turbulence that regulates the particle dynamics becomes wider. The scale ratio between the largest and the smallest eddies is known to increase in proportion to Re 3/4 . Considering the molecular viscosity of a protoplanetary disk, the Reynolds number is estimated to be as high as Re = 10 10 (α/10 −2 )(R/AU) −3/2 in the Minimum-Mass Solar Nebula (MMSN) Model (Hayashi 1981), where α is the turbulence parameter and R is the radial distance from the central star. However, even if the Reynolds number is smaller than this value, we can trace the particle behavior over a range of Stokes numbers according to the simulated inertial range of turbulence. If we focus on particle sizes from millimeters to meters, simulations of Re > O(10 4 ) are required. In the simulations by PP13, the resolution was Re ≃ 10 3 , which corresponds to St η /St = 23.5 and realizes the inertial range of turbulence over only one order of magnitude in linear dimensions. However, a recent development in supercomputers allows us to perform particle tracking simulations based on the DNS at Reynolds numbers as high as Re > O(10 4 )  In this paper, we perform high Reynolds number DNSs, where the number of grid points and the Reynolds number are up to N 3 = 2048 3 and Re = 16100, respectively, which corresponds to St η /St = 85 and can realize the inertial range over two order of magnitude in linear dimensions. By these DNSs, we obtain the rms relative velocity for particle pairs and the PDF of the collision velocities. We find that the rms relative velocity is smaller than that derived by the Völk-type model developed by OC07. In addition, we discuss the growth timescale of dust aggregates and the sticking rates in the context of overcoming the drift and fragmentation barriers.
In Section 2, we present the method of our particle tracking simulation, based on the DNS of forced incompressible homogeneous isotropic turbulence. The statistics of the motion of the particles obtained by the DNS are shown in Section 3. The statistics include the rms relative velocity, the collision kernel, and the PDF of the radial relative velocities. The results are compared to the Völk-type model by OC07. In Section 4, assuming the MMSN model, we assess the collision timescale for both compact and fluffy aggregates and the sticking rates of dust particles. Our conclusions are summarized in Section 5.

DNS of forced incompressible turbulence
In protoplanetary disks, the turbulence is known to be subsonic and thus essentially incompressible (Hayashi 1981). Therefore, in this paper, we consider threedimensional turbulence of an incompressible fluid of unit density that obeys the Navier-Stokes (NS) equations and the continuity equation where u, p, ν, and f are the velocity, pressure, kinematic viscosity, and external force, respectively. The numerical method used in the DNS is essentially identical to that used in Yokokawa et al. (2002) and Kaneda et al. (2003), which is briefly reviewed here for the convenience. (The readers may refer to Yokokawa et al. (2002) and Morishita et al. (2015) for the details of parallel computations.) In the DNS, the Navier-Stokes equations are solved by a fully alias-free Fourier spectral method, where aliasing errors are removed by the so-called phase-shift method. (Note that the Fourier spectral method gives the spatial derivatives to spectral accuracy and is not subject to the numerical viscosity. On the other hand, in finite difference schemes, the spa-tial derivatives are calculated to finite accuracy and the numerical viscosity is inevitably included.) The computation domain is assumed to be 2π-periodic in each Cartesian coordinate direction, so that both the minimum wavenumber and the wavenumber increment in the DNS are unity. The maximum wavenumber is given by where N is the number of grid points in each of the Cartesian coordinates in real space. Time integration is achieved using a fourth-order Runge-Kutta method with a constant time increment.
The forcing f that generates turbulence is given bŷ f (k) = −cû(k) in the wavevector space, wheref and u are Fourier transforms of f and u, respectively. The value of c is set as non-zero only in the wavenumber range of k < 2.5 and is adjusted at every time step so as to keep the total kinematic energy, E ≡ 3u ′2 /2 = U 2 /2, almost time-independent (≈ 0.5). Here, u ′ is the root mean square (rms) value of the fluctuating velocity in one direction and U = u · u 1/2 is the threedimensional (3D) rms of flow velocity u. A similar forcing was used in Kerr (1985), Vincent & Meneguzzi (1991), and Jimenez et al. (1993). The integral length scale L and the eddy turnover time T are defined by L = π/(2U 2 ) k 0 E(k)/kdk and T = L/U , respectively. T corresponds to Ω −1 in the actual dimensions. The variations of L and T in steady turbulence generated by the forcing are within 10% of the mean (Ishihara et al. 2007).
We use the velocity fields obtained by the DNS in Kaneda et al. (2003) and Ishihara et al. (2007) as the initial conditions for the present study and adopt the same values of the kinematic viscosity. Therefore, each velocity field at t = 0 is in a statistically steady state of turbulence and we do not have any initial transient period in the turbulence field used in the particle tracking simulation. The details of the turbulence characteristics of the DNS data are found in Ishihara et al. (2007).
In Table 1, we show a summary of the DNS parameters and turbulence characteristics. In the DNSs, the value of ν is chosen so that k max η = 1 or 2, where η is the Kolmogorov micro length scale given by η = (ν 3 / ε ) 1/4 with the mean energy dissipation rate ε . The value of k max η represents a smallscale resolution. Each "Run" in Table 1 is named by the combination of the values of N and k max η. It is known that the sensitivity/insensitivity to the values of k max η depends on the statistical quantity to be studied (Yamazaki et al. 2002;Watanabe & Gotoh 2007;Schumacher 2007;Yeung et al. 2015). It will be shown that the collision statistics are not so sensitive to the k max η values, provided that k max η 1. Figure 1 shows the compensated energy spectra,   k 5/3 E(k)/ ε 2/3 , of the generated turbulence for different values of Re and k max η. The spectrum of the velocity field used in PP13 and PP15is also shown for comparison. The horizontal range corresponds to the "inertial range", where E(k) ∝ k −5/3 . In Figure 1, we observe that the horizontal range becomes wider with increasing Re. Therefore, it is expected that the inertial range may be satisfactorily resolved by our DNS.
In the DNS of Re = 16100, we realize L/η = 1.2 × 10 3 and T /τ η = 85, respectively, which are much larger than L/η = 1.4 × 10 2 and T /τ η = 23.5 for Re = 1000. We recognize that the compensated energy spectra obtained by our DNSs have a pile-up near k = k max . It is known that such a pile-up is caused by the wavenumber truncation in the DNS based on a Fourier spectral method. It is also known that such a pile-up does not appear in the turbulence simulation using finite difference methods. Therefore, the difference between our spectra and the spectrum of PP13 and PP15 near k = k max comes from the difference in the numerical methods.
To see an effect of the wavenumber truncation in the Fourier spectral method, we compare the energy spectra between runs 512-1 and 1024-2 (and also between runs 256-1 and 512-2) in Figure 1. They are slightly different from each other in the low wavenum-ber range (k < 3) and in the high wavenumber range (kη > 0.7). The difference in the low wavenumber range is presumably caused by the difference in the energycontaining eddies at the forcing scales, while the difference at high wavenumbers is caused by the wavenumber truncation. Note that in contrast to the high and low wavenumber ranges, the difference between the spectra in the intermediate range for the two runs is very small. This comparison suggests that the collision statistics are not sensitive to the difference between k max η ∼ 1 and k max η ∼ 2, if they are insensitive to the detail of the fine-scale statistics in the energy dissipation range of the turbulence. In this study, we will mainly present the results obtained by the DNS data with k max η ∼ 1. However, we will also show the results of the DNS data with k max η ∼ 2 to confirm the insensitivity of the quantity to the small-scale resolution (see figures 3 and 10).

Particle tracking simulation
We consider the motion of small solid particles with density ρ s in a gas flow with density ρ g . The ratio β = ρ s /ρ g is assumed to be much larger than unity. Then, the equation of motion of each inertial particle is given by where X, V , and τ p are the position, velocity, and stopping time of the particle, respectively, and u is the velocity of the fluid flow at X (see, e.g., Davila & Hunt 2001). Equation (3) is solved with the fourth-order Runge-Kutta method, where the time step (∆t) is set to be twice as large as that used to solve Equation (1) (because we need u(x, n∆t/2) in this scheme). The velocity u at the particle position is evaluated by using an interpolation method. To obtain accurate statistics, we use the cubic spline interpolation (see Yeung & Pope 1988). The interpolation scheme is implemented by solving tridiagonal matrix problems in parallel with the method developed by Mattor et al. (1995). See Ishihara et al. (2015) for the actual parallel implementation.
In this paper, we normalize Equation (3) in two ways. One is by using L and T and the other is by using η and τ η ≡ (ν/ ε ) 1/2 . In the former, the motion of the particles is characterized by the Stokes number given by St = τ p /T , and in the latter, the motion is characterized by another Stokes number, given by St η = τ p /τ η . Note that St η /St = T /τ η ∝ Re 1/2 . In the case of the typical protoplanetary disk of Re = O(10 10 ), the ratio is O(10 5 ). However, even if the Reynolds number is smaller than this value, we can trace the particle behavior over a range of Stokes numbers according to the simulated inertial range of the turbulence. In our largest DNS of Re = 16100, the value of St η /St is 85.
In our particle tracking simulation, we set seven different values of St as follows.
The corresponding values of St η = St(T /τ η ) are 0.1, 0.2, 0.5, 1, 2, 5, and 10 for N = 512, 1024 and 0.1, 0.5, 1, 2, 5, 10, and 20 for N = 2048. For each value of St, we maximally track 512 3 particles. The particles are distributed randomly in the whole computational domain at t = 0 with initial velocity V 0 = 0. The statistics related to the particle motions are taken at t = 3T in the following analyses. By performing several preliminary runs, we have confirmed that the particle statistics become almost timeindependent after t = 3T .
In this paper, we focus on the particles with the Stokes number less than 0.3. The collision statistics of particles with St = O(0.1) are expected to be mainly affected by the eddies in the inertial range which are able to simulate properly by high-resolution DNSs of homogeneous isotropic turbulence using appropriate forcing methods. Therefore, the aim of this paper is to provide reliable results of collision statistics of the particles in the inertial range by performing a series of the high-resolution DNSs. Since the relative velocity attains its maximum at St ∼ 1, such collisions may be important to discuss the fragmentation barrier. However, it is also true that such collisions occur via the processes of collision and sticking of the dusts with smaller values of St which are dominated by eddies in the inertial range. Therefore, we took a bottom-up approach.

Relative velocity
The relative velocity between a pair of particles with a separation r and its rms value are calculated at locations X 1 and X 2 as w = V 1 − V 2 and respectively, where V 1 and V 2 are the velocities of the particles and N p is the number of pairs. Figure 2 shows the r-dependence of the rms relative velocity between equal-sized particles. Actually, collision velocities should be measured at particle sizes, which are in general much smaller than the Kolmogorov length scale in turbulence in protoplanetary disks. To estimate the rms relative velocity at a separation of less than η/4 in the DNSs, one may need a much larger number of particles than that used in our DNSs. However, Figure  2 demonstrates that the relative velocity for the particles of St > 0.05 is not significantly affected by eddies smaller than η(∼ 10 −3 L) and therefore almost constant at a separation r 10 −3 L. Figure 3 presents the St-dependence of the rms relative velocity at a fixed small separation of r = 10 −3 L for different Reynolds numbers. Here, r/L = 10 −3 approximately corresponds to r/η = 1/8, 1/4, 1/2, and 1 for runs at Re = 936, 2100 (and 2310), 6700, and 16100, respectively. (The other comparison for different Re values using a fixed separation of r/η would be also possible. However, here we are interested in St(= τ p /T ) dependence of the rms relative velocity. So, we measure the relative velocity using a fixed small separation r normalized by L.) The data from PP13 (Re = 1000, r = η/4) agree well with our data (Re = 936, r = η/8) in the range of large St ( 0.1). They deviate from our results in the range St 10 −2 . This deviation is presumably due to the difference in the value of r/η, and comes from the range of St η < 1. However, in this paper, we focus on the collision statistics in the range St η > 1. In addition, the data for Re = 2100(k max η = 1) agree well with the data for Re = 2300(k max η = 2). This St=0.117 St=0.235 Figure 2. rms relative velocity as a function of the separation (r) between a pair of particles. The data are plotted for each St number for 512 3 particles at t = 3T in the DNS of Run2048-1 (Re = 16100) and normalized by u 2 1/2 ≡ u · u 1/2 or uη.  agreement indicates that the relative velocities at small separations are not so sensitive to the difference in the value of k max η.
In Figure 3, we perceive that the relative velocity at the small separation of r/L = 10 −3 is an increasing function of Re for a fixed value of St. Also, we observe that the curves tend to approach a line with slope 1/2 at the portion of larger values of St( 0.1) for Re > 10 4 . This dependence ( w 2 1/2 ∝ St 1/2 ) is consistent with the inertial range scaling of the relative velocity in Völk-type models; see PP15.

Bidisperse case
In the above, we considered the monodisperse case of the identical particles. Here, we present the results for bidisperse cases, i.e., the relative velocities of particles with different St numbers (St 1 and St 2 ). Figure  4 shows the St(= St 1 )-dependence of the rms relative velocity between the particles with different ratios of f ≡ St 2 /St 1 . To compare the relative velocities at a scale of smallest eddies, we show the results at r = η/4. We find that the rms relative velocity is higher for the different-sized particles than for the equal-sized particles. This trend is prominent for St ≈ 10 −2 . The result can be understood by considering that the equal-sized particles with small separations move in the same way in a turbulent flow, and therefore the relative velocities remain low.
It is known that collisions between particles at high speeds may lead to bouncing or fragmentation. Therefore, the sticking rate of collision particles depends on the statistics of relative velocities. The present results suggest that the collisions between the identical particles are more appropriate for sticking than the collisions between particles of different sizes. We will discuss the sticking rate quantitatively in Section 4.

Comparison with the Völk-type model
As for the relative velocities of particles induced by turbulence, a closed-form expression developed by OC07 has been widely adopted (e.g., a review by Johansen et al. 2014). The original formalism (the Völk-type model) was developed by Völk et al. (1980) and Markiewicz et al. (1991). Cuzzi & Hogan (2003) obtained closed-form expressions for the Völk-type model. OC07 generalized the approach and results of Cuzzi & Hogan (2003) to obtain closed-form expressions for relative velocities between particles of arbitrary, and unequal, sizes.
For equal-size particles with a stopping time of τ p , the closed-form expression (OC07) is given by where , and the critical wavenumber k * for the particle with stopping time τ p is determined by where the relative velocity v rel between a particle and an eddy is given by and k L is the smallest wavenumber. OC07 assumed the energy spectrum given by for which the total energy is V 2 g /2. As shown in PP15, in turbulent flows with a wide inertial range, the above model predicts a τ 1/2 p scaling for w 2 1/2 , if τ p is in the inertial range.
In Figure 5, we compare the DNS results of the relative velocity with the Völk-type model given by Equation (8) in which Equation (11) is used for E(k) by setting k η /k L = Re 4/3 and Re = 10 8 . For comparison, the Völk-type models which tentatively employ the DNS data for E(k) in Equation (8) are also plotted. In the range St = O(0.1), the latter models agree well with the Völk-type model that assumes the model spectrum (11).
We can recognize that the DNS results for Re = 16100 has a slope 1/2 for 0.05 St 0.2. The slope is consistent with the Völk-type model for high-Re turbulence. However, the DNS values are smaller than the Völktype model by a factor of two. PP15 did not show the clear slope of the relative velocity, but their result suggested that the Völk-type models typically overestimate the rms of the particle relative velocity. Our results suggest that in the inertial range of high-Re turbulence the rms of the particle relative velocity obeys the scaling law, but the values are smaller than the Völk-type model by a factor of two. PP15 proposed an improved closed expression of the relative velocity (their Eq. 25) and an improved formula of Völk-type models (at their Fig.10). We plot these expressions in Figure 5 for comparison. The difference between the former and the Völk-type model by OC07 is small. On the other hand, the latter formula seem to work well at St 0.1.
In Figure 6, we compare the DNS results of the relative velocities to those of the closed-form expression derived by OC07. The DNS results for Re = 16100 are shown by colored squares with rms values as functions of St 1 = τ 1 /t L and St 2 = τ 2 /t L . For unequal-sized particles with stopping times of τ 1 and τ 2 , the closed-form expression (OC07) is given by where t L = (V L k L ) −1 , V 2 L = (2/3)V 2 g , and t η = Re −1/2 t L . τ * 1 and τ * 2 are the solutions of Equation (9) for τ p = τ 1 and τ p = τ 2 , respectively. The relative velocity variance w 2 1/2 by this formula is shown by contours in Figure 6 for the case Re = 10 8 . The comparison with the relative velocities obtained by the DNS clarifies that the DNS results are smaller than half of the values of the closed-form expression given by OC07, regardless of the St ratios. In particular, for equal-sized particles with St 10 −2 , the DNS results are smaller by an order of magnitude. These reduced relative velocities have a significant impact on avoiding the fragmentation barrier, which is discussed in details in Section 4.

Collision kernel
So far, we have focused on the variance of relative velocities, which gives a measure of the collision rate. Here, we consider the collision kernel to incorporate the turbulent clustering effect properly. The collision rate per unit volume between two particles with radii a 1 and a 2 can be expressed as . Turbulence-induced rms relative velocities w 2 1/2 normalized to u 2 1/2 between two particles characterized by St1 and St2. Colored squares (with rms values) represent the DNS results for a particle separation of r = η/4 at t = 3T (Re = 16100). Gray contours (with attached rms values) denote the prediction of the Völk-type model (12), which is equivalent to Figure 4C of OC07. The diagonal line (St1 = St2) corresponds to Figure 5. Note that the DNS results for St1 (or St2) < 0.05 are more or less affected by viscosity.
where n 1 (a 1 ) and n 2 (a 2 ) are the average number densities and Γ(a 1 , a 2 ) is the collision kernel. The kernel formula hitherto used in dust coagulation models is Γ com = πd 2 w 2 1/2 with d = a 1 + a 2 , where the effect of turbulent clustering is not taken into account and the rms relative velocity, w 2 1/2 , is usually taken from the model of Völk et al. (1980). As a statistical mechanical description of Γ for zero-inertia particles, Saffman & Turner (1995) proposed a spherical formulation, in which the kernel is given by Γ = 2πd 2 |w r | with the radial rela- tive velocity, w r = w · (X 2 − X 1 )/|X 2 − X 1 |. In addition, Sundaram & Collins (1997) considered the turbulent clustering effect on Γ and derived an expression for finite-inertia particles. Wang et al. (2000) suggested that a formula (based on the spherical formulation), is more accurate than the formula considered by Sundaram & Collins (1997), where g(d) is the radial distribution function (RDF) at r = d. The RDF is related to the two-point correlation function ξ(r) by g(r) = 1 + ξ(r), and evaluated by where N (r) is the average number of particles in a spherical shell of volume 4πr 2 dr at a distance r from a reference particle andn is the average particle number density. For a uniform distribution of the particles, we have g(r) = 1, i.e., ξ(r) = 0. For a non-uniform distribution of the particles due to turbulent clustering, we have g(r) > 1. Recently, Pan & Padoan (2014) and PP15 evaluated Γ sph by performing the DNS with Re = O(10 3 ) to study the effect of turbulent clustering on the collision kernel.
Here, we use our DNS data with the Reynolds numbers up to Re = 16100 to obtain |w r | and g(r), and elucidate the Re-dependence of the collision kernel.  Figure 7 shows the normalized radial relative velocity, 2 |w r | / w 2 1/2 , as a function of St for different Re numbers. As seen in this figure, the values are always less than unity, and therefore the use of w 2 1/2 for the collision kernel leads to the overestimation. Additionally, since the values are a decreasing or increasing function of St, |w r | does not obey the power law (∝ St 1/2 ) like w 2 1/2 does (Figure 3). Therefore, it is concluded that there are quantitative and qualitative discrepancies from w 2 1/2 for evaluating the collision kernel Γ. The low value of |w r | leads to the enhancement of the particle sticking rates as discussed in Section 4. Pan & Padoan (2014) studied the St-dependence of the ratio 2 |w r | / w 2 1/2 and showed that the values of the ratio approach its Gaussian value of (8/3π) 1/2 ≃ 0.92 for St > 1. They argued that the non-Gaussianity peaks at St η ≃ 1, resulting in a dip in the ratio at St η ≃ 1. Figure  7 shows that the bottom position of the dip is around St η ≃ 1 − 2, and therefore the dip becomes shallower accordingly as Re increases. Figure 8 shows the St-dependence of g(r) in the DNS of Run 2048-1 (Re = 16100). The behaviors of g(r) are quantitatively consistent with the recent DNS results by Ireland et al. (2016a). In Figure 8, we see that g(r) is approximately constant at r/L 10 −2 for St 0.1. This fact allows us to use the constant value when considering the collision statistics (as discussed in Section 4). Figure 9 shows the St-dependence of g(r) at r/L = 10 −3 and r/L = 10 −2 for three different values of Re. Here, r/L = 10 −3 (10 −2 ) corresponds to r/η = 1/4 (5/2), 1/2 (5), and 1 (10) for Re = 2100, 6700, and 16100, respectively. Figure 9 demonstrates that g(r) has a peak for a fixed value of r. Also, Figure 8 suggests that the value of St at the peak depends on r, and g(r)  has a peak near St η = 1 for small values of r = O(η). In Figure 9, we notice that, for fixed values of St = O(0.1) and r/L, g(r) is a decreasing function of Re provided that St η ≥ 2.0; this Re-dependence is weaker for the larger values of St. We observe also that, for fixed values of St = O(0.1) and Re, g(r) becomes larger as r/L decreases; this r/L-dependence is weaker for the larger values of St (as shown in Figure 8). The Re-dependence and r/L-dependence of g(r) suggest that an asymptotic behavior of g(r) (in the range of St = O(0.1)) for much smaller r and for much higher Re can be surmised from the DNS data of Re = O(10 4 ). Ireland et al. (2016a) showed that, for fixed values of St η 1 and r/η, g(r) is an increasing function of Re. Their Re-dependence does not contradict with the results in Figure 9. In fact, our DNS data also gives almost the same Re-dependence as that shown in Fig St f=1 f=1/2 f=1/5 f=0 f=1(PP15) f=1/2(PP15) f=1/4(PP15) f=0(PP15) Figure 11. Collision kernel per unit cross section in the spherical formulation at the distance of η/4 for particle pairs with fixed Stokes ratios, f ≡ St2/St1 = 1, 1/2, 1/5, and 0. The data are obtained at t = 3T in the DNS of Re = 16100. The data from PP15 (Re = 1000, r = η/4 ∼ 2 × 10 −3 L) are also plotted for comparison.
noted that their collision statistics are completely different from each other. Figure 10 presents the Re-dependence of the resulting collision kernel obtained by our DNS. In Figure 10, we see a noticeable increase in the collision kernel in the small St range with increasing Re. For example, the value of the collision kernel at St = 2 × 10 −2 for Re = 16100 is six times larger than that for Re = 1000. However, the Reynolds number dependence of the collision kernel at St ≈ 0.02 − 0.2 seems to converge provided that Re 10 4 . From this finding, it is expected that the collision kernel at St ≈ 0.02 − 0.2 for a much higher Reynolds number (Re ≫ 10 4 ) is similar to those obtained by our DNS of Re = O(10 4 ).
As for bidisperse cases, the St(= St 1 )-dependence of the collision kernel (Equation (14)) is presented for different f = St 2 /St 1 in Figure 11. To compare the collision kernel at the scale of the smallest eddies, we show the results at the distance of η/4. As seen in this figure, the values of Γ for different-sized particles are smaller than those for the identical particles. This result is caused by the fact that the concentration of particles due to turbulent clustering occurs more effectively for the identical particles than for the different-sized particles. The comparison with the data from PP15 shows that the DNS results for St 0.1 at Re = 16100 are not far from those at Re = 1000.

PDF of radial relative velocity
In the collision kernel, we have adopted the averaged value of the radial relative velocities, |w r | . However, to assess the fraction of particles which have velocities lower than the critical collision velocity, we should derive the PDF of w r . Since the critical collision velocity may be different between equal-sized collisions and different- sized ones, we need to obtain the PDF depending on the Stokes number ratio f ≡ St 2 /St 1 . Based on the DNS data of Re = 16100, we acquire the PDF for equal-sized particles, P (eq) . Figure 12 shows the resultant PDF of the normalized radial relative velocity w r /U at a separation of r = η/4 for each St number, where U is the rms value of the fluctuating velocity in one direction. The negative and positive values of w r /U represent the approaching and receding pairs, respectively. The variance and kurtosis of x ≡ w r /U are defined by and respectively, and listed in Table 2. It is shown that when St ≥ 1.17 × 10 −2 (St η ≥ 1.0), the variance is an increasing function of St and the kurtosis is a decreasing function of St. The data in Table 2 indicate that the values of V and K are respectively approximated by in the range 5.87 × 10 −2 ≤ St ≤ 2.35 × 10 −1 . The PDFs for different-sized particles, P (diff) , are shown for f ≡ St 2 /St 1 = 1, 1/2, and 1/4, fixing St 1 = 0.235, in Figure 13. As expected, Figure 13 suggests that the non-Gaussianity of the PDFs is weakened as f decreases (as the size difference becomes large). The similar trend has already been observed in low-Re simulations (see, e.g., Figs. 2, 4, and 8 in Pan et al. (2014b)). Our DNS results confirm that the trend is also true for higher Re turbulence. The variance and the kurtosis of the PDFs for the different-sized particles (f = 1) are listed in Table 3. The data in Tables 2 and 3 imply that, as f decreases, the value of K monotonically decreases while the variance monotonically increases. Table 2. Variance (V ) and kurtosis (K) of the normalized relative velocity (wr/U ) measured at t = 3T for the equal-sized particles with a separation r = η/4 in the DNS of Re = 16100. Parameters, µ and β, in Equation (19)   Following Sundaram & Collins (1997); Wang et al. (2000); PP13; Pan et al. (2014b), we attempt to fit the PDF of x = w r /U with a stretched exponential function given by where P se (x) is normalized to satisfy ∞ −∞ P se (x)dx = 1 as in PP13 and Pan et al. (2014b). For this function, the variance and kurtosis are analytically given by and respectively. Since the value of K se does not depend on the value of β, we can first determine the value of µ with K se = K and then the value of β with V se = V , where V and K are the DNS values. Specifically, the formula to determine the values of µ and β is given by where V and K are given by the approximated formulas (Equation (18)). The values of µ and β determined by the DNS for the case Re = 16100 are shown in Tables 2  and 3. Theoretically, a stretched exponential PDF with µ = 4/3 was predicted for inertial-range particles under the assumption of exactly Gaussian flow velocity and Kolmogorov scaling (Gustavsson et al. 2008;PP13). The values of µ for the inertial-range particles of St = O(0.1) in our DNS are less than the theoretical value and suggest that non-Gaussian behavior of the flow velocity affects on the PDF of the radial relative velocity of the inertial range particles. Hence, the fitting function (Equation (19)) and the values of µ and β in Tables 2 and 3 are useful for modeling the PDF of the radial relative velocities for the inertial range particles. Note that the values of µ for the particles of 0.1 < St < 0.3 in our DNS are not far from those for the corresponding St values reported for low-Re simulations in Pan et al. (2014b). This fact suggests that not only the variance but also the PDF of the normalized relative velocity for St 0.1 is not so sensitive to the values of the Reynolds number.
An example of the stretched exponential fit to the PDF is plotted in Figure 14(a) for the case of identical particles. We confirm that the stretched exponential function qualitatively approximates well the PDFs of the relative velocities. In Figure 14(b), we show the accuracy of the fitting quantitatively by evaluating the integral Q(x) ≡ x −x P (eq) (x)dx.  Table 2 are used.

IMPLICATIONS FOR PLANETESIMAL FORMATION
The density of the compact dust aggregates is virtually equal to that of the monomer (ρ s ≈ 1g cm −3 ). The formation of cm-sized (St = 0.01 − 0.1) compact aggregates (pebbles) is a key for the growth to planetesimals via the streaming instability (e.g. Johansen et al. 2014). On the other hand, dust aggregates can have fluffy structures with much lower bulk densities if the compression by impacts is sufficiently weak (Wada et al. 2009;Okuzumi et al. 2009Okuzumi et al. , 2012Zsom et al. 2010Zsom et al. , 2011. In recent years, the collisions of fluffy dust aggregates have been explored by performing N -body molecular-dynamics simulations (Wada et al. 2009;Paszun & Dominik 2009;Wada et al. 2013;Meru et al. 2013;Seizinger & Kley 2013;Gunkelmann et al. 2016). These studies have revealed that fluffy aggregates, depending on the breaking energy, may resolve the difficulties due to the bouncing and fragmentation barriers. So far, fluffy aggregates are thought to likely form as a result of the coalescence of icy grains. However, quite recently, Arakawa & Nakamoto (2016) have shown that fluffy aggregates can result from the collisions of nanometer-sized silicate grains. Here, according to Johansen et al. (2014), we consider two cases of compact aggregates with ρ s = 1g cm −3 and extremely fluffy aggregates with ρ s = 10 −5 g cm −3 .

Collision velocity in the MMSN model
Based on the MMSN model, we evaluate the rms relative velocities induced by turbulence. In Figure 6, we showed the rms relative velocities as a function of St number pairs. Translating the Stokes numbers into the particle sizes according to Johansen et al. (2014), we can obtain the rms collision velocities for compact and fluffy dust aggregates. In Figure 15, we show the resultant rms collision velocities at 1 AU from the central star, for α = 10 −4 in panels (a, c) and for α = 10 −2 in panels   Figure 15. rms relative velocities (in m s −1 ) between two particles of different sizes for compact aggregates (ρs = 1g cm −3 ) (top panels) and fluffy aggregates (ρs = 10 −5 g cm −3 ) (bottom panels). Here, the MMSN model is employed assuming 1 AU from the central star for (a, c) α = 10 −4 and (b, d) α = 10 −2 . Filled squares (with rms values) represent the DNS results, while the contours (with rms values) are the theoretical estimates obtained from OC07, using the particle sizes corresponding to the Stokes numbers given by Johansen et al. (2014).
(b, d). The upper panels are the results for compact aggregates, and the lower are those for fluffy aggregates. For comparison, the estimates derived with the closedform analytic formula by OC07 are also depicted. We find that the rms relative velocities by the DNSs are smaller by more than a factor of two, compared to the prediction of the closed-form expression. In particular, for particles of equal size, the DNS results are much lower. These results have a considerable impact on the sticking rates and the collision timescale as discussed below.

Sticking rates of colliding pairs
We estimate the "sticking rates", which is defined to be the probabilities of the sticking of colliding pairs per unit time. Hubbard (2012) pointed out that to estimate the sticking rates of colliding pairs the weighting factor proportional to the collision rate should be taken into account. In the cylindrical kernel formulation, the collision rate is proportional to |w|, and the collision-rate weighted distribution of the 3D amplitude is obtained from the unweighted PDF, P (|w|), simply by = 8m s −1 for different-sized collisions, which are given for silicate dust by Wada et al. (2013). Colored solid curves represent the DNS results calculated from the PDFs, using pairs of particles at a separation r = η/4 in the DNS of Re = 16100. Colored dotted curves denote the theoretical prediction assuming a Gaussian (Maxwell) distribution whose variance is given by OC07.
As shown in Pan & Padoan (2014), the spherical kernel formulation of the weighted PDF is also possible but is more complicated. Therefore, we use (23) for estimating the sticking rates, which are given by for equal-sized and different-sized collisions, respectively, where x = |w| and u c is the critical velocity above which the collisions lead to the bouncing or fragmentation (see e.g. Wada et al. (2013)). Wada et al. (2009), using N -body molecular dynamics simulations, have obtained the critical velocity for silicate aggregates to be u c = 1.1(a/0.76µm) −5/6 m s −1 , where a is the size of monomers in aggregates. It agrees well with the laboratory experiments (Blum & Wurm 2000). Additionally, Wada et al. (2013) have derived a scaling relation of the critical collision velocity as a function of the breaking energy. They have considered fluffy aggregates composed of ballistic particle-cluster aggregation clusters, which are fairly compact (fractal dimension D ∼ 3). Even in such compact aggregates, all surface interactions between monomers in contact in the aggregates determine the critical collision velocity, which depends on the size of the monomers. For icy aggregates, they have derived the critical collision velocity as for collisions between equal-sized and different-sized particles, respectively, and for silicate aggregates u (eq) c = 6(a/µm) −5/6 m s −1 , Recently, Gunkelmann et al. (2016) studied the porosity-dependence of the fragmentation of fluffy aggregates composed of silicate grains of 0.76 µm radius, and found that the critical velocity for agglomerate fragmentation decreases with the porosity of the aggregates. Although the critical collision velocities are still under debate, those for silicate aggregates are smaller by an order of magnitude than those for icy aggregates. Therefore, the sticking rates are expected to be much lower in silicate aggregates. Here, we consider two cases of u c for silicate aggregates; u c = 1m s −1 for compact aggregates and equations (27) and (28) for fluffy aggregates.
We evaluate the sticking rates of P (eq) stick and P (diff) stick assuming R = 1 AU for the cases with α = 10 −4 and 10 −2 , based on the PDFs obtained by the DNS with Re = 16100. In Figure 16, we plot the results as a function of St. For comparison, we also plot the theoretical prediction based on the collision-rate weighted distribution of the 3D amplitude assuming a Gaussian (Maxwell) distribution with the variance given by OC07.
Panels (a) and (b) show respectively the DNS results for α = 10 −4 and 10 −2 , assuming u c = 1m s −1 . Quite interestingly, the sticking rates are not strongly dependent on α if St 0.01 and keep a high level of 50%, although the theoretical prediction from a Gaussian distribution declines steeply for α = 10 −2 . Besides, the declination of the rates at St 0.01 is much more gradual compared to the theoretical prediction, especially in the case of α = 10 −2 . It shows that the non-Gaussianity of the velocity distribution function due to turbulent clustering makes the sticking rates for St 0.01 remarkably higher than those theoretically expected. Also, it is worth noting that the rates for equal-sized particles (f = 1) are higher than those for different-sized particles (f = 1). In the case of α = 10 −2 , the difference is more than an order. This comes from the fact that the variance of the relative velocities is smaller for equalsized particles (f = 1), as shown in Figure 15. These results imply that equal-sized particles grow much faster than different-sized particles, and therefore the fraction of equal-sized particles increases with time.
Panels (c) and (d) show respectively the DNS results for α = 10 −4 and 10 −2 , assuming the critical collision velocity for fluffy aggregates given by Wada et al. (2013). In panel (c), the DNS results show that the sticking rates are approximately unity at St 0.03, and also the sticking rates are slightly higher for different-sized particles (f = 1). As shown in Figure 13, the PDFs have small differences for different f values. However, u (diff) c is larger than u (eq) c , and therefore P (diff) stick is a bit higher than P (eq) stick . But, in the case of α = 10 −2 , the rates for f = 1 are higher than those for f = 1 and decline more gradually compared to the theoretical prediction at St 0.01.

Bouncing/fragmentation probabilities
Using the PDFs, we can evaluate the bouncing/fragmentation fraction of colliding particles, if we specify the critical collision velocity u c as in the previous subsection. The bouncing/fragmentation probabili-ties are estimated as for equal-sized and different-sized collisions, respectively, where x = w r /U . We evaluate the B/F probabilities of p (eq) B/F and p (diff) B/F assuming R = 1 AU for the cases with α = 10 −4 and 10 −2 , based on the PDFs obtained by the DNS with Re = 16100. In Figure 17, we plot the results as a function of St. In any case, the obtained B/F probabilities are significantly lower than the theoretical prediction from a Gaussian distribution. Also, importantly, the B/F probabilities for equal-sized particles (f = 1) are much lower than those for different-sized particles (f = 1). For St 0.01, the B/F probabilities for f = 1 are lower than 10% even in the case of α = 10 −2 . As discussed in the above, the rms relative velocity in the DNSs is smaller than that in the closed-form expression, and also the averaged radial relative velocity is even lower than the rms relative velocity, as shown in Section 3.4. Owing to such low values of relative velocities and the non-Gaussianity of the PDF in the DNSs, the B/F probabilities are not dramatically enhanced even for α = 10 −2 . Hence, equal-sized particles preferentially survive even in a highly turbulent flow.

Streaming instability
As shown above, the fluffy silicate aggregates may coalesce with a high probability in a weakly turbulent disk (α = 10 −4 ). However, the coagulation of compact aggregates are much less effective for St 0.1 in a highly turbulent disk (α = 10 −2 ) , and hence there is still an obstacle of the radial drift barrier. For the growth from cm-sized (St = 0.01 − 0.1) compact aggregates (pebbles) to planetesimals, the streaming instability may be a possible route to circumvent this obstacle. Youdin & Goodman (2005) discovered the streaming instability promoted by the action-reaction pair of the drag force between solid particles and gas, and the growth of pebbles via the streaming instability has been extensively explored (Johansen et al. 2014, and references therein). Importantly, Johansen et al. (2009) and Bai & Stone (2010) have pointed out that there is a critical solid abundance Z c , above which spontaneous strong concentration of solids occurs, where Z is the solidto-gas surface mass density ratio, Z ∼ 0.01 being the solar abundance. Using 2D simulations, Carrera et al. (2015) found that the critical abundance is super-solar and increases drastically with decreasing particle size for St < 0.1. Very recently, using 2D and 3D highresolution simulations, Yang et al. (2017) Table 5. Collision timescale of dust particles defined by t coll ≡ 1/[n d g(r) |wr| σ], where n d (= ρ d /m d ) is the number density of the dust particles with the radius a (ρ d is the dust density in the disk and m d the mass of the dust aggregate), g(r) is the RDF of the particle, |wr| is the average of the radial relative velocities, and σ(= π(2a) 2 ) is the cross-section area. The timescale is shown for compact aggregates (ρs = 1g cm −3 ) and fluffy aggregates (ρs = 10 −5 g cm −3 ). 0.03 < Z c < 0.04 for particles of St = 10 −3 . Although these are slightly super-solar, some mechanisms are still required to enhance the solid abundance to allow the formation of planetesimals via the streaming instability. In our simulations, as shown in Figure 8, the solid abundance is enhanced by the turbulent clustering, dependent on St. Since the clustering is dependent on scales, the density contrast of dust particles, ρ d /ρ d , depends on a coarse-grain scale. In Table 6, we show the maximal density contrast of dust particles as a function of St, where ∆L is a coarse-grain scale in units of α 1/2 H and the density is calculated in the volume of ∆L 3 . As seen in Table 6, the enhancement of solid abundance becomes larger according as St increases and ∆L lessens. If α = 10 −2 and ∆L ≤ 0.03H, the enhancement is a factor of ∼ 1.5 for St = 0.01 and ∼ 3.5 for St = 0.1. This satisfies the condition for the critical abundance shown in Fig. 9 of Youdin & Goodman (2005).
As shown in Figure 16, the sticking of equal-sized particles is faster than that of different-sized ones. This supports the assumption of particles of the same size employed by Yang et al. (2017). The sticking rate of equal-sized particles is as high as 50% at St 0.01, and reduces to 10% at St 0.1, in the case of α = 10 −2 . Furthermore, Figure 17 shows that the bouncing/fragmentation probabilities are as low as 10% at St 0.01, and increase steeply toward St 0.1. Therefore, equal-sized particles of St ∼ 0.01 are expected to selectively grow.
Also, we evaluate the collision time-scale of dust particles, which is given by t coll ≡ 1/[n d g(r) |w r | σ], where n d is the number density of the dust particle of size a, g(r) is the RDF of the particle, |w r | is the averaged radial relative velocity, and σ(= π(2a) 2 ) is the cross-section area. Note that the collision timescale is independent of the value of α, since n d ∝ α −1/2 and |w r | ∝ α 1/2 . As for g(r) and |w r | , we use the values at r/η = 1/4 for Re = 16100. If we specify the values of St and ρ s , we have n d and a in the Stokes regime, using the vertical scale-height for dust particles (Okuzumi et al. 2012). In Table 5, the evaluated collision times are listed for compact aggregates (ρ s = 1g cm −3 ) and fluffy aggregates (ρ s = 10 −5 g cm −3 ). Since g(r) is a decreasing function of St and |w r | is an increasing function, the collision timescale is insensitive to St and much shorter than the drift timescale even for compact aggregates.
Considering the above assessments, we can expect that mostly equal-sized particles of 0.01 St 0.1 grow in a timescale of ∼ Ω −1 . Simultaneously, the turbulent clustering enhances the solid abundance. According to Yang et al. (2017), the critical solid abundance tends to be minimal toward St ∼ 0.1. Hence, the streaming instability may be triggered at 0.01 < St < 0.1, and the strong concentration of solids proceeds in a timescale of 100Ω −1 .

CONCLUSIONS
In order to investigate the dynamical statistics of dust particles through turbulent clustering in a protoplanetary disk, we have performed the high-resolution DNSs of the Navier-Stokes equations. The number of grid points and the Reynolds number are up to 2048 3 and Re = 16100, respectively, which are of the highest resolution ever in astrophysical DNSs. These large-scale DNSs have allowed us to track the motion of dust particles with Stokes numbers of 0.01 St 0.2 in the inertial range for the first time. As results of these simulations, we have found the following: • As the Reynolds number of the turbulence in-creases (or the inertial range extends), the rms relative velocity, w 2 1/2 , of particle pairs at a fixed small separation (normalized by the integral length scale) is augmented for small St number particles, and is asymptotically proportional to St 1/2 in the inertial range.
• The rms relative velocities by the DNSs are smaller by more than a factor of two, compared to those from the closed-form expression derived by Ormel & Cuzzi (2007), irrespective of the St number ratios of the particle pairs. Also, the averaged radial relative velocity is even lower than the rms relative velocity. Hence, the findings by Pan & Padoan (2013 have been confirmed by high-Re DNSs. • The PDFs of the radial relative velocities are highly non-Gaussian, and are well fitted by a stretched exponential function like Equation (19). The PDF of the normalized relative velocity for St 0.1 is not so sensitive to the values of the Reynolds number. Hence, the results are consistent with those at low-Re by Pan et al. (2014b).
• Almost independently of α, the sticking rates of colliding particles are as high as 50% for particles of St 0.01 and declines gradually at St 0.01, although the theoretical prediction from a Gaussian distribution declines steeply for α = 10 −2 . This comes from the non-Gaussianity of the radial relative velocities and the smaller variance of the relative velocities as a result of the turbulent clustering.
• Since the variance of the relative velocities for equal-sized particles (f = 1) is smaller than that for different-sized particles (f = 1), the sticking rates for f = 1 are higher than those for f = 1. The difference is larger than an order of magnitude in the case of α = 10 −2 . It implies that equal-sized particles grow much faster than different-sized particles, and therefore the fraction of equal-sized particles increases with time.
• The bouncing/fragmentation probabilities are significantly lower than the theoretical prediction from a Gaussian distribution. The probabilities for equal-sized particles (f = 1) are much lower than those for different-sized particles (f = 1), in the case of α = 10 −2 . Hence, equal-sized particles preferentially survive even in a highly turbulent flow.
• The turbulent clustering enhances the solid abundance. The enhancement in a scale of 0.03 scaleheight of the disk increases from a factor of ∼ 1.5 at St = 0.01 to ∼ 3.5 at St = 0.1, in the case of α = 10 −2 . Therefore, the streaming instability may be triggered at 0.01 < St < 0.1.
In the present DNSs, we have assessed sticking rates, but the actual coagulation of dust particles has not been incorporated. Hence, we cannot predict the mass function of dust aggregates as a result of the hierarchical coagulation. This is a significant issue for planetesimal formation. In DNSs in the near future, we will explore the growth of dust aggregates in turbulence.
We are grateful to K. Furuya, E. Kokubo, S. Michikoshi, T. Nakamoto, S. Okuzumi, and K. Yoshida for their valuable discussions. The computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Projects ID: hp150174, ID: hp160102, and ID: hp170087) was partially used in this study. It also used the FX100 system at the Information Technology Center, Nagoya University. This research was supported in part by the Interdisciplinary Computational Science Program of the Center for Computational Sciences, University of Tsukuba, Grant-in-Aid for Scientific Research (B) by JSPS (15H03603,15H03638), and MEXT as "Exploratory Challenge on Post-K computer" (Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental Variations of Planets in the Solar System).