Dust Coagulation Regulated by Turbulent Clustering in Protoplanetary Disks

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+({\boldsymbol{u}}\cdot {\rm{\nabla }}){\boldsymbol{u}}=-{\rm{\nabla }}p+\nu {{\rm{\nabla }}}^{2}{\boldsymbol{u}}+{\boldsymbol{f}}\end{eqnarray} \tag{ 1 }$$

and the continuity equation

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{u}}=0,\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{u}}$, p, $\nu$, and ${\boldsymbol{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), Kaneda et al. (2003) and Ishihara et al. (2016), which is briefly reviewed here for convenience. (The readers may refer to Yokokawa et al. 2002 and Ishihara et al. 2015 for details of the parallel computations.) In the DNS, the NS 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 }=\sqrt{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 ${\boldsymbol{f}}$ that generates turbulence is given by $\hat{{\boldsymbol{f}}}({\boldsymbol{k}})=-c\hat{{\boldsymbol{u}}}({\boldsymbol{k}})$ in the wave-vector space, where $\hat{{\boldsymbol{f}}}$ and $\hat{{\boldsymbol{u}}}$ are Fourier transforms of ${\boldsymbol{f}}$ and ${\boldsymbol{u}}$, respectively. The value of c is set to be 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\equiv 3u{{\prime} }^{2}/2={U}^{2}/2$, almost time independently (≈0.5). Here, $u^{\prime}$ is the rms value of the fluctuating velocity in one direction and $U={\langle {\boldsymbol{u}}\cdot {\boldsymbol{u}}\rangle }^{1/2}$ is the three-dimensional (3D) rms of the flow velocity ${\boldsymbol{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 as $L=\pi /(2{U}^{2}){\int }_{0}^{k}E(k)/k{dk}$ and T = L/U, respectively. T corresponds to Ω⁻¹ in 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 for 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 $\eta ={({\nu }^{3}/\langle \varepsilon \rangle )}^{1/4}$ with the mean energy dissipation rate $\langle \varepsilon \rangle$. The value of k_maxη represents a small-scale 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; Schumacher 2007; Watanabe & Gotoh 2007; Donzis et al. 2008; Yeung et al. 2015). It will be shown that the collision statistics are not that sensitive to the k_maxη values, provided that k_maxη ≳ 1.

Figure 1 shows the compensated energy spectra, ${k}^{5/3}E(k)/\langle \varepsilon {\rangle }^{2/3}$, of the generated turbulence for different values of Re and k_maxη. The spectrum of the velocity field used in PP13 and PP15 is 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 = 16,100, we realize L/η = 1.2 × 10³ and T/τ_η = 85, respectively, which are much larger than L/η = 1.4 × 10² and T/τ_η = 23.5 for Re = 1000.

Zoom In Zoom Out Reset image size

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 spectra of PP13 and PP15 near k = k_max comes from the difference in the numerical methods.

To see the effect of the wavenumber truncation on 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 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 using 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 small-scale resolution (see Figures 3 and 10).

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

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{X}}}{{dt}}={\boldsymbol{V}},\quad \displaystyle \frac{d{\boldsymbol{V}}}{{dt}}=\displaystyle \frac{1}{{\tau }_{{\rm{p}}}}({\boldsymbol{u}}-{\boldsymbol{V}}),\end{eqnarray} \tag{ 3 }$$

where ${\boldsymbol{X}}$, ${\boldsymbol{V}}$, and τ_p are the position, velocity, and stopping time of the particle, respectively, and ${\boldsymbol{u}}$ is the velocity of the fluid flow at ${\boldsymbol{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 ${\boldsymbol{u}}({\boldsymbol{x}},n{\rm{\Delta }}t/2)$ in this scheme). The velocity ${\boldsymbol{u}}$ at the particle position is evaluated by using an interpolation method. To obtain accurate statistics, we use 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 ${\tau }_{\eta }\equiv {(\nu /\langle \varepsilon \rangle )}^{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}}_{\eta }={\tau }_{{\rm{p}}}/{\tau }_{\eta }$. Note that ${{St}}_{\eta }/{St}=T/{\tau }_{\eta }\propto {{Re}}^{1/2}$. In the case of a typical protoplanetary disk of Re = O(10¹⁰), the ratio is O(10⁵). 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 = 16,100, the value of St_η/St is 85.

In our particle tracking simulation, we set seven different values for St as follows:

$$\begin{eqnarray}{St} & = & 0.00158,0.00632,0.0158,0.0316,\\ & & \ 0.0632,\ 0.158,0.316,\ \mathrm{for}\ N=512,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}{St} & = & 0.00188,0.0038,\ 0.0094,0.0188,\\ & & \ 0.038,0.094,0.188,\ \mathrm{for}\ N=1024,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}{St} & = & 0.00117,0.00587,\ 0.0117,0.0234,\\ & & \ 0.0587,0.117,0.234,\ \mathrm{for}\ N=2048.\end{eqnarray} \tag{ 6 }$$

The corresponding values of ${{St}}_{\eta }={St}(T/{\tau }_{\eta })$ 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 track a maximum of 512³ particles. The particles are distributed randomly in the whole computational domain at t = 0 with initial velocity ${{\boldsymbol{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 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 that are able to be simulated 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 the collision statistics of the particles in the inertial range by performing a series of 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.

The relative velocity between a pair of particles with a separation r and its rms value are calculated at locations ${{\boldsymbol{X}}}_{1}$ and ${{\boldsymbol{X}}}_{2}$ as ${\boldsymbol{w}}={{\boldsymbol{V}}}_{1}-{{\boldsymbol{V}}}_{2}$ and

$$\begin{eqnarray}&&\langle {w}^{2}{\rangle }^{1/2}={\left\{\displaystyle \frac{1}{{N}_{{\rm{p}}}}\displaystyle \sum _{| {{\boldsymbol{X}}}_{1}-{{\boldsymbol{X}}}_{2}| =r}{{\boldsymbol{w}}}^{2}\right\}}^{1/2},\end{eqnarray} \tag{ 7 }$$

respectively, where ${{\boldsymbol{V}}}_{1}$ and ${{\boldsymbol{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 $\eta (\sim {10}^{-3}L)$ and therefore almost constant at a separation r ≲ 10⁻³L.

Zoom In Zoom Out Reset image size

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⁻³ approximately corresponds to r/η = 1/8, 1/4, 1/2, and 1 for runs at Re = 936, 2100 (and 2310), 6700, and 16,100, respectively. (The other comparison for different Re values using a fixed separation of r/η would also be possible. However, here we are interested in the 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⁻². 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η.

Zoom In Zoom Out Reset image size

In Figure 3, we perceive that the relative velocity at the small separation of r/L = 10⁻³ 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 with larger values of St (≳0.1) for Re > 10⁴. This dependence ($\langle {w}^{2}{\rangle }^{1/2}\propto {{St}}^{1/2}$) is consistent with the inertial range scaling of the relative velocity in Völk-type models; see PP15.

In the above, we considered the monodisperse case of identical particles. Here, we present the results for bidisperse cases, i.e., the relative velocities of particles with different St numbers (St₁ and St₂). Figure 4 shows the St( = St₁) dependence of the rms relative velocity between particles with different ratios of f ≡ St₂/St₁. To compare the relative velocities at the scale of the smallest eddies, we show the results at $r=\eta /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⁻². 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.

Zoom In Zoom Out Reset image size

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 collisions between identical particles are more appropriate for sticking than collisions between particles of different sizes. We will discuss the sticking rate quantitatively in Section 4.

As for the relative velocities of particles induced by turbulence, a closed-form expression developed by OC07 has been widely adopted (e.g., see the 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

$$\begin{eqnarray}&&\langle {w}^{2}\rangle ={\int }_{{k}^{* }}^{\infty }4E(k)(1-F{(k)}^{2}){dk},\end{eqnarray} \tag{ 8 }$$

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

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\tau }_{{\rm{p}}}}=\displaystyle \frac{1}{{\tau }_{k}({k}^{* })}+{k}^{* }{v}_{\mathrm{rel}}({k}^{* }),\end{eqnarray} \tag{ 9 }$$

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

$$\begin{eqnarray}&&{v}_{\mathrm{rel}}^{2}({k}^{* })={\int }_{{k}_{L}}^{{k}^{* }}2E(k^{\prime} ){[{\tau }_{{\rm{p}}}/({\tau }_{{\rm{p}}}+{\tau }_{k}(k^{\prime} ))]}^{2}{dk}^{\prime} \end{eqnarray} \tag{ 10 }$$

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

$$\begin{eqnarray}E(k)=\left\{\begin{array}{ll}{V}_{{\rm{g}}}^{2}/(3{k}_{L}){(k/{k}_{L})}^{-5/3} & \ \mathrm{for}\ {k}_{L}\lt k\lt {k}_{\eta },\\ 0\ & \ \mathrm{for}\ k\lt {k}_{L},k\gt {k}_{\eta },\end{array}\right.\end{eqnarray} \tag{ 11 }$$

for which the total energy is ${V}_{{\rm{g}}}^{2}/2$. As shown in PP15, in turbulent flows with a wide inertial range, the above model predicts a ${\tau }_{{\rm{p}}}^{1/2}$ scaling for $\langle {w}^{2}{\rangle }^{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}_{\eta }/{k}_{L}={{Re}}^{4/3}$ and Re = 10⁸. For comparison, the Völk-type models that 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 recognize that the DNS results for Re = 16,100 has a slope of 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 those of the Völk-type model by a factor of two. PP15 did not show a 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 those of the Völk-type model by a factor of two. PP15 proposed an improved closed expression of the relative velocity (their Equation (25)) and an improved formula of Völk-type models (at their Figure 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 seems to work well at St ≳ 0.1.

Zoom In Zoom Out Reset image size

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 = 16,100 are shown by colored squares with rms values as functions of ${{St}}_{1}={\tau }_{1}/{t}_{L}$ and ${{St}}_{2}={\tau }_{2}/{t}_{L}$. For unequal-sized particles with stopping times of τ₁ and τ₂, the closed-form expression (OC07) is given by

$$\begin{eqnarray}\langle {w}^{2}\rangle & = & \displaystyle \frac{{V}_{{\rm{g}}}^{2}}{{t}_{L}}\left({\left[{\tau }_{k}+\displaystyle \frac{{\tau }_{1}^{2}}{{\tau }_{1}+{\tau }_{k}}\right]}_{{\tau }_{k}={t}_{\eta }}^{{\tau }_{k}=\max ({\tau }_{1}^{* },{\tau }_{2}^{* })}\right.\\ & & +\displaystyle \frac{{\tau }_{2}-{\tau }_{1}}{{\tau }_{1}+{\tau }_{2}}{\left[\displaystyle \frac{{\tau }_{1}^{2}}{{\tau }_{1}+{\tau }_{k}}\right]}_{{\tau }_{k}=\max ({\tau }_{1}^{* },{\tau }_{2}^{* })}^{{\tau }_{k}={t}_{L}}\\ & & +(1\leftrightarrow 2)),\end{eqnarray} \tag{ 12 }$$

where ${t}_{L}={({V}_{L}{k}_{L})}^{-1}$, ${V}_{L}^{2}=(2/3){V}_{{\rm{g}}}^{2}$, and ${t}_{\eta }={{Re}}^{-1/2}{t}_{L}$. ${\tau }_{1}^{* }$ and ${\tau }_{2}^{* }$ are the solutions of Equation (9) for τ_p = τ₁ and τ_p = τ₂, respectively. The relative velocity variance $\langle {w}^{2}{\rangle }^{1/2}$ of this formula is shown by contours in Figure 6 for the case Re = 10⁸. 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⁻², 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 detail in Section 4.

Zoom In Zoom Out Reset image size

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₁ and a₂ can be expressed as

$$\begin{eqnarray}&&{n}_{1}({a}_{1}){n}_{2}({a}_{2}){\rm{\Gamma }}({a}_{1},{a}_{2}),\end{eqnarray} \tag{ 13 }$$

where ${n}_{1}({a}_{1})$ and n₂(a₂) are the average number densities and Γ(a₁, a₂) is the collision kernel. The kernel formula hitherto used in dust coagulation models is ${{\rm{\Gamma }}}^{\mathrm{com}}=\pi {d}^{2}\langle {w}^{2}{\rangle }^{1/2}$ with d = a₁ + a₂, where the effect of turbulent clustering is not taken into account and the rms relative velocity, $\langle {w}^{2}{\rangle }^{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 ${\rm{\Gamma }}=2\pi {d}^{2}\langle | {w}_{{\rm{r}}}| \rangle$ with the radial relative velocity, ${w}_{{\rm{r}}}={\boldsymbol{w}}\cdot ({{\boldsymbol{X}}}_{2}-{{\boldsymbol{X}}}_{1})/| {{\boldsymbol{X}}}_{2}-{{\boldsymbol{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),

$$\begin{eqnarray}&&{{\rm{\Gamma }}}^{\mathrm{sph}}=2\pi {d}^{2}\langle | {w}_{{\rm{r}}}| \rangle g(d),\end{eqnarray} \tag{ 14 }$$

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 $\xi (r)$ by g(r) = 1 + ξ(r), and evaluated by

$$\begin{eqnarray}&&g(r)=\displaystyle \frac{N(r)}{4\pi {r}^{2}\bar{n}{dr}},\end{eqnarray} \tag{ 15 }$$

where N(r) is the average number of particles in a spherical shell of volume 4πr²dr at a distance r from a reference particle and $\bar{n}$ is the average particle number density. For a uniform distribution of particles, we have g(r) = 1, i.e., ξ(r) = 0. For a non-uniform distribution of 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³) to study the effect of turbulent clustering on the collision kernel.

Here, we use our DNS data with Reynolds numbers up to Re = 16,100 to obtain $\langle | {w}_{{\rm{r}}}| \rangle$ and g(r), and to elucidate the Re dependence of the collision kernel. Figure 7 shows the normalized radial relative velocity, $2\langle | {w}_{{\rm{r}}}| \rangle /\langle {w}^{2}{\rangle }^{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 $\langle {w}^{2}{\rangle }^{1/2}$ for the collision kernel leads to the overestimation. Additionally, since the values are a decreasing or increasing function of St, $\langle | {w}_{{\rm{r}}}| \rangle$ does not obey the power law (∝ St^1/2) like $\langle {w}^{2}{\rangle }^{1/2}$ does (Figure 3). Therefore, it is concluded that there are quantitative and qualitative discrepancies from $\langle {w}^{2}{\rangle }^{1/2}$ for evaluating the collision kernel Γ. The low value of $\langle | {w}_{{\rm{r}}}| \rangle$ 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\langle | {w}_{{\rm{r}}}| \rangle /\langle {w}^{2}{\rangle }^{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.

Zoom In Zoom Out Reset image size

Figure 8 shows the St dependence of g(r) in the DNS of Run 2048-1 (Re = 16,100). The behaviors of g(r) are quantitatively consistent with the recent DNS results by Ireland et al. (2016). In Figure 8, we see that g(r) is approximately constant at $r/L\lesssim {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).

Zoom In Zoom Out Reset image size

Figure 9 shows the St dependence of g(r) at r/L = 10⁻³ and r/L = 10⁻² for three different values of Re. Here, r/L = 10⁻³ (10⁻²) corresponds to r/η = 1/4 (5/2), 1/2 (5), and 1 (10) for Re = 2100, 6700, and 16,100, 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 also observe 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 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⁴).

Zoom In Zoom Out Reset image size

Ireland et al. (2016) showed that, for fixed values of St_η ≳ 1 and r/η, g(r) is an increasing function of Re. Their Re dependence does not contradict the results in Figure 9. In fact, our DNS data also give almost the same Re dependence as that shown in Figure 21 of Ireland et al. (2016). However, it should be noted that the collision statistics for fixed values of St_η and r/η are completely different from those for fixed values of St and r/L.

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⁻² for Re = 16,100 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⁴. From this finding, it is expected that the collision kernel at St ≈ 0.02–0.2 for a much higher Reynolds number (${Re}\,\gg \,{10}^{4}$) is similar to that obtained by our DNS of Re = O(10⁴).

Zoom In Zoom Out Reset image size

As for bidisperse cases, the St(=St₁) 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 identical particles than for different-sized particles. The comparison with the data from PP15 shows that the DNS results for St ≳ 0.1 at Re = 16,100 are not far from those at Re = 1000.

Zoom In Zoom Out Reset image size

In the collision kernel, we have adopted the average value of the radial relative velocities, $\langle | {w}_{{\rm{r}}}| \rangle$. However, to assess the fraction of particles that 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\equiv {{St}}_{2}/{{St}}_{1}$.

Based on the DNS data of Re = 16,100, 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 as

$$\begin{eqnarray}&&V\equiv \langle {(x-\langle x\rangle )}^{2}\rangle \end{eqnarray} \tag{ 16 }$$

and

$$\begin{eqnarray}&&K\equiv \langle {(x-\langle x\rangle )}^{4}\rangle /\langle {(x-\langle x\rangle )}^{2}{\rangle }^{2},\end{eqnarray} \tag{ 17 }$$

respectively, and listed in Table 2. It is shown that when ${St}\geqslant 1.17\times {10}^{-2}({{St}}_{\eta }\geqslant 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 as

$$\begin{eqnarray}&&V\approx {St}^{1.3}\,{\rm{a}}{\rm{n}}{\rm{d}}\,K\approx 3.0+0.79\times {St}^{-1.3}\end{eqnarray} \tag{ 18 }$$

in the range 5.87 × 10⁻² ≤ St ≤ 2.35 × 10⁻¹.

Zoom In Zoom Out Reset image size

Table 2.  Variance and Kurtosis of the Normalized Relative Velocity (w_r/U) Measured at t = 3T for Equal-sized Particles with a Separation r = η/4 in the DNS of Re = 16,100

${St}({{St}}_{\eta })$	0.0117(1.0)	0.0235(2.0)	0.0587(5.0)	0.117(10.0)	0.235(20.0)
Variance (V)	8.04E–4	2.58E–3	2.21E–2	6.15E–2	1.39E–1
Kurtosis (K)	283	85.4	34.1	13.8	8.33
μ	0.273	0.352	0.453	0.632	0.813
β	1.72E–5	3.17E–4	4.75E–3	3.32E–2	1.01E–1

Note. The parameters μ and β in Equation (19), determined using K_se = K and V_se = V, are also listed.

Download table as:  ASCIITypeset image

The PDFs for different-sized particles, ${P}^{(\mathrm{diff})}$, are shown for $f\equiv {{St}}_{2}/{{St}}_{1}=1$, 1/2, and 1/4, fixing St₁ = 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). A similar trend has already been observed in low-Re simulations (see, e.g., Figures 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 different-sized particles ($f\ne 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.

Zoom In Zoom Out Reset image size

Table 3.  Same as Table 2 but for Pairs of Different-sized Particles

f	1/2	1/4	1/2
${{St}}_{1}({{St}}_{\eta 1})$	0.117(10.0)	0.235(20.0)	0.235(20.0)
${{St}}_{2}({{St}}_{\eta 2})$	0.0587(5.0)	0.0587(5.0)	0.117(10.0)
Variance (V)	8.91E–2	1.74E–1	1.50E–1
Kurtosis (K)	9.40	5.13	7.03
μ	0.760	1.12	0.899
β	0.069	0.203	0.130

Note. The value of f denotes the ratio ${{St}}_{2}/{{St}}_{1}$.

Download table as:  ASCIITypeset image

Following Sundaram & Collins (1997), Wang et al. (2000), PP13, and Pan et al. (2014b), we attempt to fit the PDF of x = w_r/U with a stretched exponential function given by

$$\begin{eqnarray}&&{P}_{\mathrm{se}}(x)=\displaystyle \frac{\mu }{2\beta {\rm{\Gamma }}(1/\mu )}\exp \left[-{\left(\displaystyle \frac{| x| }{\beta }\right)}^{\mu }\right],\end{eqnarray} \tag{ 19 }$$

where ${P}_{\mathrm{se}}(x)$ is normalized to satisfy ${\int }_{-\infty }^{\infty }{P}_{\mathrm{se}}(x){dx}=1$ as in PP13 and Pan et al. (2014b). For this function, the variance and kurtosis are analytically given by

$$\begin{eqnarray}&&{V}_{\mathrm{se}}={\beta }^{2}{\rm{\Gamma }}(3/\mu )/{\rm{\Gamma }}(1/\mu )\end{eqnarray} \tag{ 20 }$$

and

$$\begin{eqnarray}&&{K}_{\mathrm{se}}=\displaystyle \frac{{\rm{\Gamma }}(1/\mu ){\rm{\Gamma }}(5/\mu )}{{\rm{\Gamma }}{(3/\mu )}^{2}},\end{eqnarray} \tag{ 21 }$$

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

$$\begin{eqnarray}&&{K}_{\mathrm{se}}=K,\quad \beta ={\left[\displaystyle \frac{{\rm{\Gamma }}(1/\mu )V}{{\rm{\Gamma }}(3/\mu )}\right]}^{1/2},\end{eqnarray} \tag{ 22 }$$

where V and K are given by the approximated formulas (Equation (18)). The values of μ and β determined by the DNS for the case Re = 16,100 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 the non-Gaussian behavior of the flow velocity affects 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)\equiv {\int }_{-x}^{x}{P}^{(\mathrm{eq})}(x){dx}$.

Zoom In Zoom Out Reset image size

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 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⁻⁴ in panels (a) and (c), and for α = 10⁻² in panels (b) and (d). The upper panels are the results for compact aggregates, and the lower ones 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.

Zoom In Zoom Out Reset image size

We estimate the "sticking rates," which are defined to be the probabilities of colliding pairs sticking 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 $| {\boldsymbol{w}}|$, and the collision-rate-weighted distribution of the 3D amplitude is obtained from the unweighted PDF, $P(| {\boldsymbol{w}}| )$, simply as

$$\begin{eqnarray}&&{P}_{{\rm{cyl}}}(| {\boldsymbol{w}}| )=| {\boldsymbol{w}}| P(| {\boldsymbol{w}}| )/\langle | {\boldsymbol{w}}| \rangle .\end{eqnarray} \tag{ 23 }$$

As shown in Pan & Padoan (2014), the spherical kernel formulation of the weighted PDF is also possible but is more complicated. Therefore, we use Equation (23) for estimating the sticking rates, which are given by

$$\begin{eqnarray}&&{P}_{\mathrm{stick}}^{(\mathrm{eq},\ \mathrm{diff})}={\int }_{0}^{{u}_{{\rm{c}}}^{(\mathrm{eq},\mathrm{diff})}}{P}_{{\rm{cyl}}}^{(\mathrm{eq},\mathrm{diff})}(x){dx}\end{eqnarray} \tag{ 24 }$$

for equal-sized and different-sized collisions, respectively, where $x=| {\boldsymbol{w}}|$ and u_c is the critical velocity above which the collisions lead to 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}_{{\rm{c}}}=1.1{(a/0.76\mu {\rm{m}})}^{-5/6}\,{\rm{m}}$ s⁻¹, where a is the size of monomers in aggregates. It agrees well with 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

$$\begin{eqnarray}&&{u}_{{\rm{c}}}^{(\mathrm{eq})}=60{(a/\mu {\rm{m}})}^{-5/6}\,{\rm{m}}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{u}_{{\rm{c}}}^{(\mathrm{diff})}=80{(a/\mu {\rm{m}})}^{-5/6}\,{\rm{m}}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 26 }$$

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

$$\begin{eqnarray}&&{u}_{{\rm{c}}}^{(\mathrm{eq})}=6{(a/\mu {\rm{m}})}^{-5/6}\,{\rm{m}}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{u}_{{\rm{c}}}^{(\mathrm{diff})}=8{(a/\mu {\rm{m}})}^{-5/6}\,{\rm{m}}\ {{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 28 }$$

Recently, Gunkelmann et al. (2016) studied the porosity dependence of the fragmentation of fluffy aggregates composed of silicate grains with radius of 0.76 μm, 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}_{{\rm{c}}}=1\,{\rm{m}}$ s⁻¹ for compact aggregates and Equations (27) and (28) for fluffy aggregates.

We evaluate the sticking rates of ${P}_{\mathrm{stick}}^{(\mathrm{eq})}$ and ${P}_{\mathrm{stick}}^{(\mathrm{diff})}$ assuming R = 1 au for the cases with α = 10⁻⁴ and 10⁻², based on the PDFs obtained by the DNS with Re = 16,100. 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.

Zoom In Zoom Out Reset image size

Panels (a) and (b) show, respectively, the DNS results for α = 10⁻⁴ and 10⁻², assuming ${u}_{{\rm{c}}}=1\,{\rm{m}}$ s⁻¹. Quite interestingly, the sticking rates are not strongly dependent on α if St ≲ 0.01 and remain at a high level of ≳50%, although the theoretical prediction from a Gaussian distribution declines steeply for α = 10⁻². Besides, the decrease of the rates at St ≳ 0.01 is much more gradual compared to the theoretical prediction, especially in the case of α = 10⁻². 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\ne 1$). In the case of α = 10⁻², the difference is more than an order of magnitude. This comes from the fact that the variance of the relative velocities is smaller for equal-sized 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 $\alpha ={10}^{-4}$ and 10⁻², 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 that the sticking rates are slightly higher for different-sized particles ($f\ne 1$). As shown in Figure 13, the PDFs have small differences for different f values. However, ${u}_{{\rm{c}}}^{(\mathrm{diff})}$ is larger than ${u}_{{\rm{c}}}^{(\mathrm{eq})}$, and therefore ${P}_{\mathrm{stick}}^{(\mathrm{diff})}$ is a bit higher than ${P}_{\mathrm{stick}}^{(\mathrm{eq})}$. But, in the case of α = 10⁻², the rates for f = 1 are higher than those for $f\ne 1$ and decline more gradually compared to the theoretical prediction at St ≳ 0.01.

Using the PDFs, we can evaluate the bouncing/fragmentation fraction of colliding particles, if we specify the critical collision velocity ${u}_{{\rm{c}}}$ as in the previous subsection. The bouncing/fragmentation probabilities are estimated as

$$\begin{eqnarray}&&{p}_{{\rm{B}}/{\rm{F}}}^{(\mathrm{eq},\mathrm{diff})}={\int }_{-\infty }^{-{u}_{{\rm{c}}}^{(\mathrm{eq},\mathrm{diff})}/U}{P}^{(\mathrm{eq},\mathrm{diff})}(x){dx}\end{eqnarray} \tag{ 29 }$$

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

We evaluate the B/F (bouncing/fragmentation) probabilities of ${p}_{{\rm{B}}/{\rm{F}}}^{(\mathrm{eq})}$ and ${p}_{{\rm{B}}/{\rm{F}}}^{(\mathrm{diff})}$ assuming R = 1 au for the cases with α = 10⁻⁴ and 10⁻², based on the PDFs obtained by the DNS with Re = 16,100. 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\ne 1$). For St ≲ 0.01, the B/F probabilities for f = 1 are lower than 10% even in the case of α = 10⁻². As discussed above, the rms relative velocity in the DNSs is smaller than that in the closed-form expression, and also, the average 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⁻². Hence, equal-sized particles preferentially survive even in a highly turbulent flow.

Zoom In Zoom Out Reset image size

As shown above, the fluffy silicate aggregates may coalesce with a high probability in a weakly turbulent disk (α = 10⁻⁴). However, the coagulation of compact aggregates is much less effective for St ≳ 0.1 in a highly turbulent disk (α = 10⁻²), and hence there is still the obstacle of the radial drift barrier. For the growth from centimeter-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 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 a spontaneous strong concentration of solids occurs, where Z is the solid-to-gas surface mass density ratio, with Z ∼ 0.01 being the solar abundance. Using 2D simulations, Carrera et al. (2015) found that the critical abundance is supersolar and increases drastically with decreasing particle size for St < 0.1. Very recently, using 2D and 3D high-resolution simulations, Yang et al. (2017) 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⁻², and 0.03 < Z_c < 0.04 for particles of St = 10⁻³. Although these are slightly supersolar, some mechanisms are still required to enhance the solid abundance to allow the formation of planetesimals via streaming instability.

In our simulations, as shown in Figure 8, the solid abundance is enhanced by turbulent clustering, which is dependent on St. Since the clustering is dependent on scales, the density contrast of dust particles, ${\rho }_{d}/{\bar{\rho }}_{d}$, depends on a coarse-grain scale. In Table 5, 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 ${\rm{\Delta }}{L}^{3}$. As seen in Table 5, the enhancement of solid abundance becomes larger as St increases and L decreases. If α = 10⁻² 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 critical abundance shown in Figure 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⁻². Furthermore, Figure 17 shows that the B/F 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 timescale of dust particles, which is given by ${t}_{\mathrm{coll}}\equiv 1/[{n}_{{\rm{d}}}g(r)\langle | {w}_{{\rm{r}}}| \rangle \sigma ]$, where ${n}_{{\rm{d}}}$ is the number density of the dust particle of size a, g(r) is the RDF of the particle, $\langle | {w}_{{\rm{r}}}| \rangle$ is the averaged radial relative velocity, and $\sigma (=\pi {(2a)}^{2})$ is the cross-section area. Note that the collision timescale is independent of the value of α, since ${n}_{{\rm{d}}}\propto {\alpha }^{-1/2}$ and $\langle | {w}_{{\rm{r}}}| \rangle \propto {\alpha }^{1/2}$. As for g(r) and $\langle | {w}_{{\rm{r}}}| \rangle$, we use the values at r/η = 1/4 for Re = 16,100. 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 6, the evaluated collision times are listed for compact aggregates (${\rho }_{{\rm{s}}}=1{\rm{g}}$ cm⁻³) and fluffy aggregates (${\rho }_{{\rm{s}}}={10}^{-5}\,{\rm{g}}$ cm⁻³). Since g(r) is a decreasing function of St and $\langle | {w}_{{\rm{r}}}| \rangle$ is an increasing function, the collision timescale is insensitive to St and much shorter than the drift timescale even for compact aggregates.

Table 6.  Collision Timescale of Dust Particles Defined as ${t}_{\mathrm{coll}}\equiv 1/[{n}_{{\rm{d}}}g(r)\langle | {w}_{{\rm{r}}}| \rangle \sigma ]$, where ${n}_{{\rm{d}}}(={\rho }_{{\rm{d}}}/{m}_{{\rm{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, $\langle | {w}_{{\rm{r}}}| \rangle$ is the Average of the Radial Relative Velocities, and σ(=π (2a)²) is the Cross-section Area

St	Compact (${\rho }_{{\rm{s}}}=1\,{\rm{g}}$ cm−3)					Fluffy (${\rho }_{{\rm{s}}}={10}^{-5}\,{\rm{g}}$ cm−3)
	0.0117	0.0235	0.0587	0.117	0.235	0.0117	0.0235	0.0587	0.117	0.235
a [cm]	4.36	6.17	9.75	13.8	19.5	1379	1950	3083	4361	6167
md [kg]	0.347	0.982	3.88	11.0	31.1	110	311	1228	3473	9823
${t}_{\mathrm{coll}}$[yr]	1.07	1.01	1.07	1.12	1.32	0.00338	0.00319	0.00337	0.0035	0.00418

Note. The timescale is shown for compact aggregates (ρ_s = 1 g cm⁻³) and fluffy aggregates (${\rho }_{{\rm{s}}}={10}^{-5}\,{\rm{g}}$ cm⁻³).

Download table as:  ASCIITypeset image

Considering the above assessments, we can expect that mostly equal-sized particles of 0.01 ≲ St ≲ 0.1 grow in a timescale of ∼Ω⁻¹. Simultaneously, 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 on a timescale of ≳100 Ω⁻¹.