Dust Coagulation Regulated by Turbulent Clustering in Protoplanetary Disks
Takashi Ishihara Naoki Kobayashi Kei Enohata Masayuki Umemura Kenji Shiraishi
astro-ph.EP







Dust coagulation regulated by turbulent clustering
























1Graduate School of Environmental and Life Science, Okayama University, Okayama 700-8530, Japan; ishihara@ems.okayama-u.ac.jp
2Department of Computational Science and Engineering, Graduate School of Engineering, Nagoya University, Nagoya 464-8603, Japan
3Center for Computational Sciences, University of Tsukuba, Tsukuba 305-8577, Japan; umemura@ccs.tsukuba.ac.jp
4Institute of Materials and Systems for Sustainability, Nagoya University, Nagoya 464-8601, Japan




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 α-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 
submicron-size dust particles to kilometer-size bodies in protoplanetary disks
<cit.>.
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 centimeter- to meter-size particles have large drift velocities and may rapidly fall into 
the central star on a timescale of a few 100 yr 
<cit.>,
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 
<cit.>.

The streaming instability is a potential mechanism to circumvent the radial drift barriers
<cit.>. 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
<cit.>.
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
<cit.>.
<cit.> 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 
<cit.>.
<cit.> 
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 
<cit.>.
<cit.>(OC07) 
provided closed-form analytical approximations 
to the Völk-type model. 
<cit.>,
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. 
<cit.> 
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,
<cit.> and <cit.> 
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 
<cit.>.
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 
<cit.>.
Furthermore, effects such as “caustics” <cit.> and the “sling effect” <cit.>
allow the particles with large inertia to become less coupled with the local fluid velocity field and 
assemble from different regions <cit.>. 
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. 
<cit.> 
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 
<cit.>(PP13), <cit.>, and <cit.>(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 <cit.>, 
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) <cit.>

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.




§ PARTICLE TRACKING USING DNS OF TURBULENCE
 



 §.§ DNS of forced incompressible turbulence




In protoplanetary disks, the turbulence is known to be subsonic and 
thus essentially incompressible <cit.>. 
Therefore, in this paper, we consider three-dimensional turbulence of an incompressible fluid of unit density
that obeys the Navier-Stokes (NS) equations

    ∂ u/∂ t + ( u ·∇)  u 
     = -∇ p + ν∇ ^2  u +  f

and the continuity equation

    ∇· u = 0,

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 <cit.>
and <cit.>, which is briefly reviewed here for the convenience.
(The readers may refer to <cit.> and <cit.> 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 spatial 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 k_ max=√(2)N/3, 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 by f̂( k)=-cû( k) in the wavevector space, 
where f̂ and û 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 ≡ 3 u'^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 three-dimensional (3D) rms  of flow velocity u. 
A similar forcing was used in <cit.>, <cit.>, and <cit.>.
The integral length scale L and the eddy turnover time T are defined by 
L=π/(2U^2)∫_0^k E(k)/k dk 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 <cit.>. 




We use the velocity fields obtained by the DNS in <cit.> and <cit.> 
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 <cit.>.

In Table <ref>, 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 small-scale resolution.
Each “Run” in Table <ref> 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 <cit.>.
It will be shown that the collision statistics are not so sensitive to 
the k_ maxη values, provided that k_ maxη≳ 1.

Figure <ref> shows the compensated energy spectra, 
k^5/3E(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 <ref>, 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 <ref>. 
They are slightly different from each other 
in the low wavenumber 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 energy-containing 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 <ref> and <ref>).




 §.§ 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

    d  X/dt =  V,   d  V/dt = 1/τ_ p( u -  V),

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 <cit.>.
Equation (<ref>) 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 (<ref>) (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 <cit.>.
The interpolation scheme is implemented by solving tridiagonal matrix problems in parallel with the method developed by <cit.>. See <cit.> for the actual parallel implementation.

In this paper, we normalize Equation (<ref>) 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.

    St   =    0.00158, 0.00632, 0.0158, 0.0316, 
           0.0632,   0.158,  0.316,  for N=512, 
    
    St   =    0.00188, 0.0038,  0.0094, 0.0188, 
           0.038, 0.094, 0.188,   for N=1024, 
    
    St   =    0.00117, 0.00587,   0.0117, 0.0234, 
           0.0587,  0.117,  0.234,  for N=2048.

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 time-independent 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.



§ NUMERICAL RESULTS




 §.§ 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

    ⟨ w^2 ⟩^1/2={1/N_ p∑_| X_1- X_2|=r w^2}^1/2,

respectively, where V_1 and V_2 are the velocities of the particles 
and N_ p is the number of pairs.
Figure <ref> 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 <ref> demonstrates that 
the relative velocity for the particles of St>0.05 is not significantly affected 
by eddies smaller than η(∼ 10^-3L)
and therefore almost constant at a separation r ≲ 10^-3L.






Figure <ref> presents the St-dependence of the rms relative velocity 
at a fixed small separation of r=10^-3L 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 agreement indicates that the relative velocities at small separations are not so sensitive to the difference in the value of k_ maxη.

In Figure <ref>, 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 <ref> 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 <ref>.





 §.§ 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 <cit.>.
The original formalism (the Völk-type model) was developed by <cit.> and <cit.>.
<cit.> obtained closed-form expressions for the Völk-type model.
OC07 generalized the approach and results of <cit.> 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

    ⟨ w^2 ⟩=∫^∞_k^*4E(k)(1-F(k)^2)dk,

where F(k)=τ_ p/(τ_ p+τ_k(k)), τ_k(k) = [2E(k)k^3]^-1/2, and 
the critical wavenumber k^* for the particle with stopping time τ_ p is determined by

    1τ_ p=1τ_k(k^*)+k^*v_ rel(k^*),

where the relative velocity v_ rel between a particle and an eddy is given by

    v_ rel^2(k^*)=∫_k_L^k^*2E(k')[τ_ p/(τ_ p+τ_k(k'))]^2dk'

and k_L is the smallest wavenumber.
OC07 assumed the energy spectrum given by 

    E(k) = {[ V_ g^2/(3k_L)(k/k_L)^-5/3            for k_L<k<k_η,;                         0        for k<k_L,  k>k_η, ].

for which the total energy is V_ g^2/2.
As shown in PP15, in turbulent flows with a wide inertial range, the above model predicts a τ_ p^1/2 scaling for ⟨ w^2⟩^1/2, if τ_ p is in the inertial range.


In Figure <ref>, we compare the DNS results of the relative velocity with the Völk-type model given by Equation (<ref>) in which Equation (<ref>) 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 (<ref>) 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 (<ref>).
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ölk-type 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  <ref>  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 <ref>, 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

    ⟨ w^2 ⟩=   V_ g^2 t_L([τ_k+τ_1^2τ_1+τ_k]_τ_k=t_η^τ_k= max(τ_1^*,τ_2^*)
        +τ_2-τ_1τ_1+τ_2[τ_1^2τ_1+τ_k]_τ_k= max(τ_1^*,τ_2^*)^τ_k=t_L
        +(1↔ 2)),

where t_L=(V_Lk_L)^-1, V_L^2=(2/3)V_ g^2, and t_η=Re^-1/2t_L. τ_1^* and τ_2^* are the solutions of Equation (<ref>) for τ_ p=τ_1 and τ_ p=τ_2, respectively. 
The relative velocity variance ⟨ w^2 ⟩ ^1/2 by this formula is shown by contours 
in Figure <ref> 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 <ref>. 





 §.§ 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 

    n_1(a_1)n_2(a_2)Γ(a_1, a_2),

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 <cit.>.
As a statistical mechanical description of Γ for zero-inertia particles, <cit.> proposed 
a spherical formulation, in which
the kernel is given by Γ=2π d^2<|w_ r|> with 
the radial relative velocity, w_ r= w·( X_2- X_1)/| X_2- X_1|.
In addition, <cit.> 
considered the turbulent clustering effect on Γ and derived an expression for finite-inertia particles. 
<cit.> suggested that a formula (based on the spherical formulation),

    Γ^ sph=2π d^2⟨|w_ r|⟩ g(d),

is more accurate than the formula considered by <cit.>, 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

    g(r)=N(r) 4π r^2 n̅dr,

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 and n̅ 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, <cit.> 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 <ref> 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 <ref>).
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 <ref>.
<cit.> 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 <ref> shows that 
the bottom position of the dip is around St_η≃1-2, and therefore
the dip becomes shallower accordingly as Re increases.









Figure <ref> 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 <cit.>.
In Figure <ref>, 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 <ref>). 

Figure <ref> 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 <ref> demonstrates that g(r) has a peak for a fixed value of r. 
Also, Figure <ref> 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 <ref>, 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 <ref>).
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).

<cit.> 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 <ref>. 
In fact, our DNS data also gives almost the same Re-dependence as that shown in Figure 21 of <cit.>. 
However, it should be noted that their collision statistics are completely different from each other.

Figure <ref> presents the Re-dependence of the resulting collision kernel obtained by our DNS.
In Figure <ref>, 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 (<ref>)) is presented
for different f= St_2/St_1 in Figure <ref>. 
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 <ref> 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 

    V≡⟨ (x-⟨ x⟩)^2 ⟩

and

    K≡⟨ (x-⟨ x⟩)^4 ⟩/⟨ (x-⟨ x⟩)^2⟩^2,

respectively, and listed in Table <ref>.
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 <ref> indicate that 
the values of V and K are respectively approximated by 

    V≈ St^1.3 and K≈ 3.0+0.79× St^-1.3,

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 <ref>. 
As expected, Figure <ref> 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 <cit.>).
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 <ref>.
The data in Tables <ref> and <ref> imply that, as f decreases, the value of K monotonically decreases 
while the variance monotonically increases. 








Following <cit.>; PP13; <cit.>,
we attempt to fit the PDF of x=w_r/U with a stretched exponential function given by

    P_ se(x)=μ 2βΓ(1/μ)exp[-(|x|β)^μ],

where P_ se(x) is normalized to satisfy ∫_-∞^∞ P_ se(x) dx =1 as in PP13 and <cit.>.
For this function, the variance and kurtosis are analytically given by 

    V_ se=β^2Γ(3/μ)/Γ(1/μ)

and

    K_ se=Γ(1/μ)Γ(5/μ)Γ(3/μ)^2,

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

    K_ se=K,   β=[Γ(1/μ)VΓ(3/μ)]^1/2,

where V and K are given by the approximated formulas (Equation (<ref>)).
The values of μ and β determined by the DNS for the case Re=16100 
are shown in Tables <ref> and <ref>.



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 (<ref>)) and the values of μ and β in Tables <ref> and <ref> 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 <cit.>. 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 <ref>(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 <ref>(b), we show the accuracy of the fitting quantitatively by evaluating the integral Q(x)≡∫_-x^x P^( eq)(x) dx.










§ 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 <cit.>.
On the other hand, dust aggregates can have fluffy structures with much lower bulk densities 
if the compression by impacts is sufficiently weak 
<cit.>.
In recent years, the collisions of fluffy dust aggregates have been explored 
by performing N-body molecular-dynamics simulations 
<cit.>.
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, <cit.> have shown that 
fluffy aggregates can result from the collisions of nanometer-sized silicate grains. 
Here, according to <cit.>, we consider two cases of
compact aggregates with ρ_ s=1g cm^-3 and extremely fluffy aggregates with ρ_ s=10^-5g cm^-3.


As for a protoplanetary disk, we assume an α-model, in which the turbulent viscosity is given as ν_ t=α c_ sH, 
c_ s is the sound speed, H is the vertical scale height, and the typical value of α is between ∼ 10^-4 to 10^-2. 
In addition, we employ the Minimum-Mass Solar Nebula (MMSN) model (Hayashi 1981), 
which provides the gas temperature
T_ g = 280(R/AU)^ - 1/2K, 
the gas density 
ρ_ g=1.2 × 10^-9 (R/ AU)^-11/4g cm^-3,
the sound speed
c_ s=1.1× 10^3(R/ AU)^-1/4 m  s^-1, 
the surface mass density of gas
Σ_ g = 1.7 × 10^3 (R/AU)^-3/2g cm^-2,
and the vertical scale height of the gas
H=4.7× 10^-2 (R/ AU)^5/4 AU,
where R is the distance from the central star.
Assuming a cross section of 2.5× 10^-15 cm^2 for hydrogen molecules, we have the kinematic viscosity given by ν=6.0× 10^4(R/ AU)^5/2 cm^2  s^-1.
In the α-model, the characteristic velocity and the characteristic length scale in the turbulence are given by U=α^1/2c_ s and L=α^1/2H, respectively, i.e., 
U=1.1×10^2(α/10^-2)^1/2(R/ AU)^-1/4  m  s^-1,
L=5.5×10^8(α/10^-2)^1/2(R/ AU)^5/4  m.
Thus, the Reynolds number Re is given by
Re=UL/ν∼10^10(α/10^-2)(R/ AU)^-3/2
and the Kolmogorov length and time scales are estimated using η=Re^-3/4L and τ_η=Re^-1/2(L/U).
Turbulence characteristics for α=10^-2 and 10^-4 are listed in Table <ref>.












 §.§ Collision velocity in the MMSN model

Based on the MMSN model, we evaluate the rms relative velocities induced by turbulence. 
In Figure <ref>, we showed the rms relative velocities as a function of St number pairs. 
Translating the Stokes numbers into the particle sizes according to <cit.>,
we can obtain the rms collision velocities for compact and fluffy dust aggregates.
In Figure <ref>, 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 (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 closed-form 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.






0.5cm


 §.§ 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

    P_cyl(| w|)=| w|P(| w|)/<| w|>.

As shown in <cit.>,  the spherical kernel formulation of the weighted PDF is also possible but is more complicated. 
Therefore, we use (<ref>) for estimating the sticking rates, which are given by

    P_ stick^( eq, diff)= ∫_0^u_ c^( eq,diff) P^( eq,diff)_cyl(x) dx

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. <cit.>).

<cit.>, 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/6m s^-1, where a is the size of monomers in aggregates.
It agrees well with the laboratory experiments <cit.>.
Additionally, <cit.> 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

    u_ c^( eq)   =   60 (a/μ m)^-5/6 m s^-1, 
    
    u_ c^( diff)   =   80 (a/μ m)^-5/6 m s^-1,

for collisions between equal-sized and different-sized particles, respectively, 
and for silicate aggregates

    u_ c^( eq)   =   6 (a/μ m)^-5/6 m s^-1, 
     
    u_ c^( diff)   =   8 (a/μ m)^-5/6 m s^-1.

Recently, <cit.> 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 (<ref>) and (<ref>)
for fluffy aggregates.

We evaluate the sticking rates of P_ stick^( eq) and P_ stick^( diff) 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 <ref>, 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 equal-sized particles (f=1), as shown in Figure <ref>. 
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 <cit.>. 
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 <ref>, the PDFs have small differences for different f values. 
However, u_ c^( diff) is larger than u_ c^( eq), and therefore 
P_ stick^( diff) is a bit higher than P_ stick^( eq). 
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 probabilities are estimated as 

    p_ B/F^( eq,diff)=∫_-∞^-u_ c^( eq,diff)/U P^( eq,diff)(x) dx

for equal-sized and different-sized collisions, respectively, where x=w_r/U.

We evaluate the B/F probabilities of p_ B/F^( eq) and p_ B/F^( diff) 
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 <ref>, 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 <ref>.
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.

<cit.> 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 <cit.>.
Importantly, <cit.> and <cit.>  have pointed out that
there is a critical solid abundance Z_c, above which spontaneous strong concentration of solids occurs,
where Z is the solid-to-gas surface mass density ratio, Z∼ 0.01 being the solar abundance. 
Using 2D simulations, <cit.> 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 high-resolution simulations, <cit.> have shown that
the critical abundance is not a steep function of St, and 
0.01 <Z_c < 0.02 for particles of St= 10^-2 and 
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 <ref>, 
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 <ref>, 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/2H 
and the density is calculated in the volume of Δ L^3.
As seen in Table <ref>, 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 <cit.>.

As shown in Figure <ref>,
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 <cit.>.
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 <ref> 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 <cit.>. 
In Table <ref>, the evaluated collision times are listed 
for compact aggregates (ρ_ s=1g cm^-3) and fluffy aggregates (ρ_ s=10^-5g 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 <cit.>, 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 increases (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 <cit.>, 
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
<cit.> 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 (<ref>).
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 <cit.>.



  * 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 scale-height 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). 




[Adachi et al.(1976)]Adachi1976 Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756 
[Arakawa & Nakamoto(2016)]Arakawa2016 Arakawa, S., & Nakamoto, T. 2016, , 832, L19 
[Bai & Stone(2010)]BS2010 Bai, X.-N., & Stone, J. M. 2010, , 722, 1437 
[Birnstiel et al.(2009)]Birnstiel2009 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, , 503, L5 
[Birnstiel et al.(2012)]Birnstiel2012 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, , 539, A148 
[Blum & Wurm(2000)]Blum2000 Blum, J., & Wurm, G. 2000, , 143, 138
[Blum & Wurm(2008)]Blum2008 Blum, J., & Wurm, G. 2008, , 46, 21 
[Bragg & Collins(2014)]Bragg2014 Bragg, A. D., & Collins, L. R. 2014, New Journal of Physics, 16, 055014
[Brauer et al.(2008a)]Brauer2008a Brauer, F., Henning, T., & Dullemond, C. P. 2008, , 487, L1 
[Brauer et al.(2008b)]Brauer2008b Brauer, F., Dullemond, C. P., & Henning, T. 2008, , 480, 859 
[Carrera et al.(2015)]Carrera2015 Carrera, D., Johansen, A., & Davies, M. B. 2015, , 579, A43 
[Chiang & Youdin(2010)]Chiang2010 Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493 
[Cuzzi & Hogan(2003)]Cuzzi2003 Cuzzi, J. N., & Hogan, R. C. 2003, , 164, 127 
[Davila & Hunt(2001)]Davila2001 Davila, J., & Hunt, J. 2001, J. Fluid Mech., 440, 117
[Donzis et al.(2008)]Donzis2008 Donzis, D. A., Yeung, P. K., & Sreenivasan, K. R. 2008, Physics of Fluids, 20, 045108 
[Donzis et al.(2008)Donzis, Yeung, & Sreenivasan]Donzis_etal2008 Donzis, D. A., Yeung, P. K., & Sreenivasan, K. R. 2008, Phys. Fluids, 20, 045108
[Falkovich & Pumir(2007)]Falkovich2007 Falkovich, G., & Pumir, A. 2007, Journal of Atmospheric Sciences, 64, 4497 
[Fessler et al.(1994)]Fessler1994 Fessler, J. R., Kulick, J. D., & Eaton, J. K. 1994, Physics of Fluids, 6, 3742 
[Garaud et al.(2013)]Garaud2013 Garaud, P., Meru, F., Galvagni, M., & Olczak, C. 2013, , 764, 146 
[Gunkelmann et al.(2016)]Gunkelmann2016 Gunkelmann, N., Ringl, C., & Urbassek, H. M. 2016, , 589, A30 
[Gustavsson & Mehlig(2011)]Gustavsson2011 Gustavsson, K., & Mehlig, B. 2011, , 84, 045304 
[Gustavsson et al.(2008)]Gustavsson2008 Gustavsson, K., Mehlig, B., Wilkinson, M., & Uski, V. 2008, Physical Review Letters, 101, 174503 
[Güttler et al.(2010)]Gttler2010 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, , 513, A56 
[Hayashi(1981)]Hayashi1981 Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35 
[Hubbard(2012)]Hubbard2012 Hubbard, A. 2012, , 426, 784 
[Ida & Guillot(2016)]Ida2016 Ida, S., & Guillot, T. 2016, , 596, L3 
[Ireland et al.(2016a)]Ireland2016a Ireland, P. J., Bragg, A. D., & Collins, L. R. 2016, Journal of Fluid Mechanics, 796, 617 
[Ireland et al.(2016b)]Ireland2016b Ireland, P. J., Bragg, A. D., & Collins, L. R. 2016, Journal of Fluid Mechanics, 796, 659 
[Ishihara et al.(2007)]Ishihara2007 Ishihara, T., Kaneda, Y., Yokokawa, M., Itakura, K., & Uno, A. 2007, Journal of Fluid Mechanics,592, 335 
[Ishihara et al.(2009)]Ishihara2009 Ishihara, T., Gotoh, T., & Kaneda, Y. 2009, Annual Review of Fluid Mechanics, 41, 165 
[Ishihara et al.(2015)Ishihara, Enohata, Morishita, Yokokawa, &Ishii]Ishihara2015 Ishihara, T., Enohata, K., Morishita, K., Yokokawa, M., & Ishii, K. 2015, Parallel Computing Technologies, Vol. 9251 (Springer International Publishing), 522–527
[Ishihara et al.(2016)]Ishihara2016 Ishihara, T., Morishita, K., Yokokawa, M., Uno, A., & Kaneda, Y. 2016, Physical Review Fluids, 1, 082403 [Jimenez et al.(1993)Jimenez, Wray, Saffman, & Rogallo]Jimenez1993 Jimenez, J., Wray, A., Saffman, P., & Rogallo, R. 1993, J. Fluid Mech., 255
[Johansen et al.(2009)]Johansen2009 Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, , 704, L75 
[Johansen et al.(2014)]Johansen2014 Johansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI, 547 
[Kaneda et al.(2003)]Kaneda2003 Kaneda, Y., Ishihara, T., Yokokawa, M., Itakura, K., & Uno, A. 2003, Physics of Fluids, 15, L21 
[Kataoka et al.(2013)]Kataoka2013 Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, , 557, L4 
[Kerr(1985)]Kerr1985 Kerr, R. M. 1985, J. Fluid Mech., 153, 31
[Lissauer(1993)]Lissauer1993 Lissauer, J. J. 1993, , 31, 129 
[Markiewicz et al.(1991)]Markiewicz1991 Markiewicz, W. J., Mizuno, H., & Voelk, H. J. 1991, , 242, 286 
[Mattor et al.(1995)Mattor, Williams, & Hewett]Mattor1995Mattor, N., Williams, T. J., & Hewett, D. 1995, Parallel Computing, 21, 1769
[Maxey(1987)]Maxey1987 Maxey, M. R. 1987, Physics of Fluids, 30, 1915 
[Meru et al.(2013)]Meru2013 Meru, F., Geretshauser, R. J., Schäfer, C., Speith, R., & Kley, W. 2013, , 435, 2371 
[Morishita et al.(2015)Morishita, Yokokawa, Uno, Ishihara, & Kaneda]Morishita2015 Morishita, K., Yokokawa, M., Uno, A., Ishihara, T., & Kaneda, Y. 2015, in Parallel CFD 2015
[Okuzumi et al.(2009)]Okuzumi2009 Okuzumi, S., Tanaka, H., & Sakagami, M.-a. 2009, , 707, 1247 
[Okuzumi et al.(2012)]Okuzumi2012 Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, , 752, 106 
[Ormel & Cuzzi(2007)]Ormel2007 Ormel, C. W., & Cuzzi, J. N. 2007, , 466, 413 (OC07)
[Pan & Padoan(2013)]Pan2013 Pan, L., & Padoan, P. 2013, , 776, 12 (PP13)
[Pan & Padoan(2014)]Pan2014_IV Pan, L., & Padoan, P. 2014, , 797, 101 
[Pan & Padoan(2015)]Pan2015 Pan, L., & Padoan, P. 2015, , 812, 10 (PP15)
[Pan et al.(2014a)]Pan2014_II Pan, L., Padoan, P., & Scalo, J. 2014a, , 791, 48 
[Pan et al.(2014b)]Pan2014_III Pan, L., Padoan, P., & Scalo, J. 2014b, , 792, 69
[Pan et al.(2011)]Pan2011 Pan, L., Padoan, P., Scalo, J., Kritsuk, A. G., & Norman, M. L. 2011, , 740, 6 
[Paszun & Dominik(2009)]Paszun2009 Paszun, D., & Dominik, C. 2009, , 507, 1023 
[Saffman & Turner(1995)]Saffman1995 Saffman, P. G., & Turner, J. S. 1995, J. Fluid Mech., 1, 16
[Schumacher(2007)]Schumacher2007 Schumacher, J. 2007, Europhysics Lett., 80, 54001
[Seizinger & Kley(2013)]Seizinger2013 Seizinger, A., & Kley, W. 2013, , 551, A65 
[Squires & Eaton(1991)]Squires1991 Squires, K. D., & Eaton, J. K. 1991, Physics of Fluids A, 3, 1169 
[Sundaram & Collins(1997)]Sundaram1997 Sundaram, S., & Collins, L. R. 1997, Journal of Fluid Mechanics, 335, 75 
[Völk et al.(1980)Volk, Jone, Morfill, & Roeser]Volk1980 Volk, H. J., Jone, F. C., Morfill, G. E., & Roeser, S. 1980, Astron. Astrophys., 85, 316
[Vincent & Meneguzzi(1991)]Vincent1991 Vincent, A., & Meneguzzi, M. 1991, J. Fluid Mech., 225, 1 
[Voelk et al.(1980)]Voelk1980 Voelk, H. J., Jones, F. C., Morfill, G. E., & Roeser, S. 1980, , 85, 316 
[Wada et al.(2009)]Wada2009 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, , 702, 1490
[Wada et al.(2013)]Wada2013 Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, , 559, A62 
[Wang et al.(2000)Wang, Wexler, & Zhou]Wang2000 Wang, L.-P., Wexler, A. S. Q., & Zhou, Y. 2000, J. Fluid Mech., 415
[Watanabe & Gotoh(2007)]Watanabe2007 Watanabe, T., & Gotoh, T. 2007, J. Fluid Mech., 590, 117
[Weidenschilling(1977)]Weidenschilling1977 Weidenschilling, S. J. 1977, , 180, 57 
[Wilkinson et al.(2006)]Wilkinson2006 Wilkinson, M., Mehlig, B., & Bezuglyy, V. 2006, Physical Review Letters, 97, 048501 
[Windmark et al.(2012)]Windmark2012 Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012, , 544, L16 
[Yamazaki et al.(2002)]Yamazaki2002 Yamazaki, Y., Ishihara, T., & Kaneda, Y. 2002, Journal of the Physical Society of Japan, 71, 777
[Yang et al.(2017)]Yang2017 Yang, C.-C., Johansen, A., & Carrera, D. 2017, , 606, A80 
[Yeung & Pope(1988)]Yeung1988 Yeung, P. K., & Pope, S. B. 1988, Journal of Computational Physics, 79, 373
[Yeung et al.(2015)Yeung, Zhai, & Sreenivasan]Yeung2015 Yeung, P. K., Zhai, X. M., & Sreenivasan, K. R. 2015, Proc. Natl. Acad. Sci. USA, 112, 12633
[Yokokawa et al.(2002)Yokokawa, Itakura, Uno, Ishihara, & Kaneda]Yokokawa2002Yokokawa, M., Itakura, K., Uno, A., Ishihara, T., & Kaneda, Y. 2002, in Proceedings of the IEEE/ACM SC2002 Conference, 50
[Youdin & Goodman(2005)]YG05 Youdin, A. N., & Goodman, J. 2005, , 620, 459 
[Zsom et al.(2010)]Zsom2010 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, , 513, A57 
[Zsom et al.(2011)]Zsom2011 Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, , 534, A73 




