Effects of the Compressibility of Turbulence on the Dust Coagulation Process in Protoplanetary Disks

Planetesimals are believed to be formed by the coagulation of dust grains in the protoplanetary disk turbulence. However, the bouncing and fragmentation barriers have not been completely solved, particularly for silicate dust. To circumvent these barriers, the turbulent clustering of dust particles must be properly treated. According to the minimum-mass solar nebula (MMSN) model, the Mach number of the turbulence ranges from M rms ≃ 0.01–0.32, and thus the turbulence is often regarded as essentially incompressible. However, it has not been quantitatively investigated whether the incompressible limit is adequate for protoplanetary disk simulations. We therefore compare in this study the motions of inertial particles in direct numerical simulations (DNSs) of the Navier–Stokes equation between weakly compressible turbulence and incompressible turbulence. In the DNSs of compressible turbulence, we use an external force to set the total dissipation and the dilatational-to-solenoidal dissipation ratio. The DNSs reveal that despite the small Mach number M rms( ≲ 0.3), the compressible turbulence field notably differs from the incompressible field in terms of the density fluctuations, pressure fluctuations, and shocklet generation, depending on the ratio of the dilatational forcing. However, we quantitatively confirmed that these effects on the particle collision statistics are weak and that the motion of inertial particles in weakly compressible turbulence is dominated by the solenoidal velocity components. Therefore we can conclude that the incompressible assumption is appropriate for an investigation of the dust coagulation process in protoplanetary disk turbulence, as assumed in the MMSN model.


Introduction
Planetesimals are believed to form through the hierarchical coagulation of dust grains, from submicron-sized dust particles to kilometer-sized bodies during the turbulence occurring in protoplanetary disks (e.g., Lissauer 1993;Chiang & Youdin 2010;Johansen et al. 2014). However, in this scenario, several barriers such as radial drift (Adachi et al. 1976;Weidenschilling 1977), bouncing, and fragmentation barriers (Blum & Wurm 2008;Brauer et al. 2008aBrauer et al. , 2008bBirnstiel et al. 2009Birnstiel et al. , 2012Güttler et al. 2010;Zsom et al. 2010Zsom et al. , 2011Windmark et al. 2012) remain incompletely resolved. The radial drift barrier is a problem because meter-sized particles have high radial drift velocities owing to gas friction force and may rapidly fall into the central star on a timescale of 10 2 or 10 3 yr (Adachi et al. 1976). To overcome this barrier, dusts should be effectively coagulated in turbulent gas flows. A bouncing or fragmentation barrier is a problem that occurs when the collision energies of particles with a relatively high velocity in turbulence cause particle bouncing or destruction, and the growth of dust aggregates is halted at the centimeter to meter scale (Windmark et al. 2012). To overcome this barrier, dust should collide simultaneously in turbulent gas flows at speeds slower than certain critical collision velocities. Because the critical collision velocity for silicate dust is significantly lower than that of ice dust, a bouncing or fragmentation barrier is more severe for silicate dust (Wada et al. 2013). To consider a possible scenario for dust particles to grow and overcome these barriers, it is important to quantitatively understand the role of turbulence in the coagulation processes.
According to the minimum-mass solar nebula (MMSN) model (Hayashi 1981), the Mach number M rms , which is the ratio of the root mean square (rms) of the velocity fluctuations to the speed of sound, of turbulence in protoplanetary disks ranges from M rms = 0.01 to M rms = 0.32 for α = 10 −4 to α = 10 −1 , where α is the turbulence parameter . In protoplanetary disks, the turbulence may therefore be regarded as subsonic and thus essentially incompressible. To study the effect of turbulence on the coagulation process, Pan et al. (2011) used an approximation method to solve the Euler equations to simulate compressible turbulence with M rms ≈ 1. They assumed that the clustering statistics in compressible turbulence would be close to that in incompressible turbulence, and demonstrated in their simulation that turbulent clustering may considerably enhance the particle collision rate. They emphasized that a definite result on the Reynolds number dependence of the clustering statistics requires simulations with a higher numerical resolution. In a series of papers (Pan & Padoan 2013Pan et al. 2014aPan et al. , 2014b, they investigated the effects of turbulent clustering using a direct numerical simulation (DNS) of compressible Navier-Stokes equations. Using the DNS, they obtained reliable results for turbulence and showed, for example, that the rms relative velocity of particle pairs is twice that of the conventional theoretical estimate based on the Völk-type model (Völk et al. 1980) developed by Ormel & Cuzzi (2007). However, the Reynolds number attained by their DNS was low. Moreover, the Reynolds number dependence of DNS results had not yet been studied in the context of planetesimal formation in protoplanetary disks. Recently, Ishihara et al. (2018) conducted DNSs of turbulence at high Reynolds numbers of up to Re = 16,100. The DNS results in Ishihara et al. (2018) are consistent with those at a low Re by Pan & Padoan (2013 and Pan et al. (2014aPan et al. ( , 2014b and show that the rms relative velocity of the inertial particles obeys a power law in the inertial range of turbulence, as suggested by Ormel & Cuzzi (2007). It was demonstrated that the clustering and collision statistics of particles with high inertia are not sensitive to the increase in the Reynolds number. However, the effects of the compressibility of turbulence in protoplanetary disks on the particle collision statistics have not been quantitatively investigated.
For the inertial particles in incompressible turbulent flows, there have been many studies from various perspectives in the past few decades, particularly in the fields of fluid mechanics and atmospheric science (e.g., Maxey 1987;Wang & Maxey 1993;Eaton & Fessler 1994;Balkovsky et al. 2001;Shaw 2003;Falkovich & Pumir 2004). In general, particles in flow fields have inertia depending on their size and density. Studies based on numerical simulations of turbulence show that particles with low inertia are ejected from strong vortices and are concentrated preferentially in regions of low vorticity and high strain rate (Maxey 1987;Squires & Eaton 1991;Eaton & Fessler 1994). In addition, particles with high inertia are known to be insensitive to the motions of small-scale eddies and are concentrated in other regions (Bec et al. 2007;Yoshimoto & Goto 2007). The latter studies based on the DNSs of incompressible turbulence suggest that the motions of the particles are determined by the eddies that have a timescale comparable to the characteristic timescale of the particles. Therefore, to study the motions of inertial particles in the inertial range of turbulence, a large-scale DNS should be conducted to attain a high Reynolds number.
By contrast, for the inertial particles in compressible turbulence, studies have been conducted specifically in the field of astrophysics. In addition to the studies on particle clustering in protoplanetary disks (e.g., Johansen et al. 2004Johansen et al. , 2006Pan & Padoan 2013Pan et al. 2014aPan et al. , 2014b, the clustering of particles in molecular clouds has also been studied (e.g., Hopkins & Lee 2016;Lee et al. 2017;Tricco et al. 2017;Mattsson et al. 2019aMattsson et al. , 2019bSeligman et al. 2019;Li & Mattsson 2021). In the former, turbulence is subsonic, whereas in the latter, turbulence is highly compressible and is dominated by shock waves. According to studies on inertial particles in high Mach number turbulence (e.g., Hopkins & Lee 2016;Mattsson et al. 2019a), owing to inertia, large grains (with a typical radius larger than 1.0 μm) decouple from the gas flow, whereas small grains (smaller than 0.1 μm) will tend to trace the motions of the gas better. However, the Reynolds number in these studies is much lower than that in reality; therefore, further studies are expected.
In the field of computational fluid dynamics, some studies have been conducted on particles in compressible turbulence (e.g., Mashayek & Jaberi 1999;Yang et al. 2014;Zhang et al. 2016;Dai et al. 2017;Zhang & Xiao 2018). Yang et al. (2014) investigated the dynamics of inertial particles in transonic homogeneous turbulence and found that heavy particles develop large clouds with high number densities downstream of the shocklets. Zhang et al. (2016) found that particles with intermediate and high inertia are collected not only in lowvorticity regions, but also in high-vorticity regions behind the randomly formed shocklets. These results are consistent with those in the astrophysical literature. By contrast, Mashayek & Jaberi (1999) studied the behavior of solid particles in forced isotropic low Mach number turbulent flows (M rms < 0.2) and showed that if the flow is subsonic and completely free of shocks, no significant M rms dependence can be observed in the particle behavior. This suggests that the motion of dust particles in protoplanetary disks may be regarded as the same as in incompressible turbulence. However, according to previous studies (Samtaney et al. 2001;Benzi et al. 2008) based on the DNS of compressible turbulence, a locally supersonic flow occurs at moderate Mach numbers, such as 0.3 and 0.4. Therefore, once shocklets are generated even in low M rms turbulence in protoplanetary disks, the dust particle statistics in compressible turbulence are expected to differ from those in incompressible turbulence. In addition, according to a recent study by Donzis & John (2020), the statistical state of compressible turbulence cannot be described by the Reynolds and Mach numbers alone, and a characteristic dilatational velocity needs to be incorporated into the governing parameters instead.
This study aims to quantitatively investigate the effect of the compressibility on the statistics related to the clustering and collisions of dust particles at scales in the universal equilibrium range of turbulence. Considering the conditions in protoplanetary disks, we compare the particle statistics in the DNSs of compressible turbulence (M rms  0.3) with those in the DNSs of incompressible turbulence. In Section 2 we describe the methods of our DNSs of compressible and incompressible turbulence and the methods of particle tracking. To maintain the turbulence at a statistically steady state, we adopted appropriate forcing methods in the DNSs. In the case of the DNS of compressible turbulence, both the solenoidal and dilatational components of velocity are forced by specifying the values of its ratio. Therefore, in the series of our DNSs, the results (statistics of turbulence and those of particles) depend on the Reynolds number, Mach number, and the ratio of the dilatational forcing to the solenoidal forcing. In Section 3 we first describe the characteristic features of the compressible turbulence field obtained by the DNSs. We then compare the particle statistics in compressible turbulence with those in incompressible turbulence and study how the results of particle statistics in incompressible turbulence ) are modified when the turbulence is compressible. In Section 4 we summarize the results.

DNS of Forced Turbulence
In this study, we consider the three-dimensional homogeneous isotropic turbulence of a compressible fluid that obeys the continuity equation, Navier-Stokes equations, and total energy equation. These equations can be written as follows: where ρ, u i , p, and E are the density, velocity, pressure, and total energy, respectively. Furthermore, f i is the forcing term on the momentum equations (see below). The total energy E, viscous stress tensor τ ij , and heat flux q i are given by where ε and T are the internal energy and temperature, respectively. The dynamic viscosity and thermal conductivity are assumed to be given bym m = T 0 0.76 andk k = T 0 0.76 , respectively, where μ 0 and κ 0 are the reference (dimensional) constants,˜= T T T 0 is the nondimensional temperature, and T 0 is the reference temperature. The equation of state is given by where R is the gas constant. The sound velocity is given by g = c RT , where γ = C p /C v is the ratio of the specific heat at constant pressure C p to that at constant volume C v . Equations (1)-(7) can be nondimensionalized using the characteristic velocity ( = á + + ñ U u u u  (7) with that in incompressible turbulence, we nondimensionalize p as Ma 2 from Equation (7), wherer r r = 0 . In our study, we consider two cases of the reference Mach number, Ma = 0.1 and Ma = 0.3, referring to the MMSN model (Hayashi 1981). In addition, we set γ = 1.4 and a Prandtl number Pr ≡ μ 0 C p /κ 0 = 0.7 considering hydrogen gas of approximately 280 K, which is the temperature at 1 au from the central star.
In our DNS, the turbulence field is assumed to be periodic in each direction of the Cartesian coordinates with a fundamental periodic box of size 2π, and thus it can be expressed as a Fourier series with both the minimum wavenumber and the wavenumber increment being one. To maintain a statistically quasi-stationary state of turbulence, we need to force the field. For this purpose, we force the field in a deterministic manner, as in Petersen & Livescu (2010). That is, the Fourier coefficient of f i in Equation (2) is not zero only at low wavenumbers |k| < 3 and is given by¯¯( where u i s and u i d are the solenoidal and dilatational components of u i , respectively, andū i is the spectrally filtered modified velocity obtained by setting the Fourier coefficient of u i for |k| 3 to zero. Because u i s and u i d satisfy = u div 0 s and rot u d = 0, the first and second terms on the right-hand side of Equation (9) drive the rotational and compressive motions of the turbulent field at large scales. The coefficients c s and c d are determined to keep the mean rate of energy dissipation per unit mass ò = 〈τ ij (∂u i /∂x j )〉/ρ 0 approximately constant and to set the ratio ò r = ò d /ò s as the target values (ò r target = 0, 0.01, and 0.1). Here, the solenoidal and compressible dissipations per unit mass are defined as (Sarkar et al. 1991), where ω = rot u is the vorticity. The relation ò = ò s + ò d is exact for constant viscosity and homogeneous turbulence. In our DNSs, we set ò r target = 0 or 0.01 for Ma = 0.1, and ò r target = 0.01 and 0.1 for Ma = 0.3. Note that the higher the values of Ma and ò r target , the stronger the effect of the compressibility.
In the simulations of turbulence aimed at describing the conditions in a planetary disk, shearing box type boundary conditions are often used (e.g., see Johansen et al. 2009). In such simulations, the large-scale shear drives the turbulence at scales in a range of universal equilibrium (from inertial to dissipation). The purpose of our study is to understand the effect of compressibility on particle motion in the inertial range of turbulence in a planetary disk. Therefore, we use the largescale forcing f i in the Equation (9) to generate the universal equilibrium range of turbulence in a statistically steady state. The term f i in Equation (2) generates the term f j u j in the righthand side of Equation (3). However, we have introduced a function F e = − f j u j as in Petersen & Livescu (2010) to ensure that the total energy is conserved and the temperature can reach a statistically stationary state. Other choices of F e , for example, F e = − 〈f j u j 〉, may also be possible. However, it is expected that the properties of particle clustering are insensitive to the details of the cooling function when the temperature field is statistically stationary (e.g., see Wang et al. 2010). We confirm this expectation in the appendix by comparing the particle statistics obtained using F e = − f j u j with those obtained in the DNSs of "isothermal" compressible turbulence, which are often used in astrophysics simulations, including Pan & Padoan (2013).
To numerically solve the Navier-Stokes equations of a compressible fluid, we use an eighth-order compact difference (CD) scheme (Lele 1992) for the spatial derivatives in the advection terms and use the eighth-order central finite difference (FD) scheme for the second derivatives in the viscous terms. The reason for the latter is to save computational costs because the viscous terms are typically smaller than the advection terms by three orders of magnitude (e.g., see Wang et al. 2010). The CD scheme is implemented using the method developed by Mattor et al. (1995) for solving tridiagonal matrix problems in parallel. Time marching is conducted using the third-order TVD Runge-Kutta method (Gottlieb & Shu 1998) with a constant time step (Δt). We used an eighth-order lowpass filtering scheme (Gaitonde & Visbal 1998) to remove unphysical numerical oscillations at high wavenumbers. To minimize the damping in low-wavenumber modes, the free parameter α f (|α f | < 0.5) in the filtering scheme (Gaitonde & Visbal 1998) was set to α f = 0.49.
For comparison, we consider the three-dimensional turbulence of an incompressible fluid of unit density (ρ = ρ 0 = 1) that obeys the continuity equation and the Navier-Stokes equations, i , and ν = ν 0 ( = μ 0 /ρ 0 ) is a constant kinematic viscosity and f i I is a solenoidal forcing term that maintains the total energy constant (see Ishihara et al. 2018 for details). In the DNS of incompressible turbulence, Equations (10) and (11) are solved through a completely alias-free Fourier spectral method, where aliasing errors are removed by the phase-shift method. The DNS method is the same as that used in Ishihara et al. (2018).
We used the incompressible velocity fields obtained by the DNS in Kaneda et al. (2003) and Ishihara et al. (2007) as the initial conditions for the DNSs of incompressible turbulence. For compressible turbulence, we first performed preliminary DNSs of forced compressible turbulence to obtain statistically steady states. We used the statistically steady states of compressible turbulence thus obtained as the initial conditions for the simulations of particle tracking in compressible turbulence.
The DNS parameters and the turbulence characteristics of the simulated turbulent flows (both compressible and incompressible cases) are listed in Table 1. The instantaneous kinetic energy (K ), mean energy dissipation rate (〈ò〉), and integral length scale (L) are computed in terms of the three-dimensional energy spectrum , respectively. Note that the values of K for the runs of incompressible turbulence are almost constant (0.5), while those for the runs of compressible turbulence fluctuate and are roughly constant (K ; 0.52 ± 0.05). The values of Re, L, 〈ν〉, and 〈ò〉 in Table 1 were obtained through time averaging. The run names consist of the initial character ("I" or "C") representing an incompressible or compressible flow and the number of the grid points (N) in each of the Cartesian coordinates in real space. Compressible flows are followed by the reference Mach number (Ma) and the character representing the types of forcing used in our DNS of compressible turbulence i.e., "n" stands for ò r target = 0, "w" indicates ò r target = 0.01, and "s" is for ò r target = 0.1. The time step was set to Δt = 1.0 × 10 −3 for N = 256 and Δt = 0.5 × 10 −3 for N = 512 and 1024. The reference Mach number Ma, target ratio ò r target , and value of ν 0 that determines the kinematic viscosity˜ñ n m r = 0 are the main parameters in our DNS. In the case of incompressible turbulence, we set the kinematic viscosity as a constant and used the same value of ν as that used in Ishihara et al. (2018) for each value of N.
The actual values of M rms ,˜ñ n m r á ñ = á ñ 0 , and ò r at a statistically steady state are close to the values of Ma, ν 0 , and ò r target , respectively (see Table 1). In our DNSs, we confirmed that regardless of the values of Ma and ò r target , the relation approximately holds for the nonconstant viscosity (˜˜) n n m r = 0 . Figure 1 shows that the compensated energy spectra approximately overlap within the range of 0.02 < kη < 0.7, which suggests that small-scale fluid motions within this range agree well with each other. We recognize that the compensated energy spectra obtained by the DNSs of incompressible turbulence have a pile-up near = k k max , where = k N 2 3 max is the maximum wavenumber in the DNSs. It is known that such a pile-up is caused by the wavenumber truncation in the DNS based on a Fourier spectral method. However, no such pile-up appears in the DNS of compressible turbulence based on FD methods. Ishihara et al. (2018) compared the collision statistics of dust particles obtained through the DNSs of different small-scale resolutions ) and confirmed that the collision statistics are insensitive to the difference in small-scale resolutions. This fact suggests that the existence of a pile-up near = k k max may have no crucial effect on the collision statistics.

Particle Tracking
We consider the motion of small solid particles with density ρ s in a gas flow field of fluid with density ρ. The ratio ρ s /ρ is assumed to be significantly higher than unity. Therefore the Table 1 Simulation (DNS) Parameters and Turbulence Characteristics equation of motion of each inertial particle is given by where X, V, and τ p are the position, velocity, and stopping time of the particle, respectively, and u is the velocity of a fluid flow at X. The stopping time is given by τ p = (ρ s /ρ)(a p /c) in the Epstein regime (when the particle size a p is smaller than the mean free path of the gas particles λ g ), whereas it is given by in the Stokes regime (when a p is larger than λ g ; e.g., see Pan & Padoan 2013;Johansen et al. 2014). For a typical gas density in protoplanetary disks, ρ ; 10 −9 g cm −3 at 1 au, λ p is approximately 1 cm, as discussed in Pan & Padoan (2013). In the case of a p  1 cm, as considered in Ishihara et al. (2018), the Stokes regime may therefore be more appropriate. Assuming that the relative velocity |u − V| is much lower than c, we ignore the modification of the Epstein drag, as discussed in Hopkins & Lee (2016) for the case of particles in molecular clouds. In addition, for simplicity, we ignore particle-particle interactions and assume that the turbulence is not modified by the motion of the particles. Equation (12) is solved using the fourthorder Runge-Kutta method, where the time step is set to twice as large as that used to solve the fluid equations. The velocity u at the particle position is evaluated using cubic spline interpolation (see Yeung & Pope 1988). The interpolation (a tridiagonal matrix problem) is calculated in parallel in the same way as that used for the CD scheme.
Equation (12) can be normalized in two ways, i.e., using L and T eddy , and using η and τ η . In the former, the motion of the particles is characterized by the Stokes number given by St = τ p /T eddy , and in the latter, it is characterized by the Stokes number given by St η = τ p /τ η . In this study, we mainly use St.
The corresponding values of St η are obtained from Table 1 for the values of T eddy and τ η ). In the DNS, we tracked eight sets of particles with = St 0, 0.003, 0.006, 0.015, 0.03, 0.06, 0.15, and 0.3, where the particles with St = 0 are fluid particles that obey dX/dt = u instead of Equation (12). For each value of St, 256 3 particles were tracked for the DNS with N = 256, 512, and 512 3 particles for the DNS with N = 1024. The particles are distributed randomly within the whole computational domain at t = 0. The initial velocity of the inertial particles (St ≠ 0) is set to V 0 = 0. The particles are tracked until t = 6T eddy , which is sufficient for the particle statistics to reach their statistically steady states.

Statistics of Compressible Turbulence
To investigate the effect of compressibility on the turbulent field, we focus on density (ρ, Figure 2), velocity divergence (q º u div , Figure 3), and the local Mach number (M loc , Figure 4), which are basic physical quantities in compressible turbulence. We examined their probability density functions (PDFs) and studied their typical spatial fluctuations using visualization. In the visualization of each physical quantity, we compare the two cases of C512-0.1w (Ma = 0.1, ò r target = 0.01; weakly compressible) and C512-0.3 s (Ma = 0.3, ò r target = 0.1; weakly but relatively strongly compressible). In Figures 2 and 3 the values of the density and velocity divergence are normalized based on their average and rms values, respectively.
Figure 2(a) shows that the density fluctuation increases with an increase in M rms and ò r , and that the PDFs of density become positively skewed as M rms and ò r increase. According to the previous studies, inertial particles tend to accumulate in highdensity regions of the background flow in both one-way (Zhang et al. 2016) and two-way (Dai et al. 2017) interactions (where a one-way interaction is a case in which the turbulent flow field is not modified by the presence of the particles, and a two-way interaction refers to a case in which the flow is affected by the particles). According to Zhang et al. (2016), there are two reasons why particles in compressible turbulence gather in high-density regions. As the first reason, the particles tend to accumulate in low-vorticity regions, which correspond to high-density regions. The accumulation of the inertial particles in low-vorticity regions is also observed in an incompressible turbulent flow. As the second reason, the particles concentrate in high-density regions formed downstream of the shocklets (see also Yang et al. 2014). Note that the latter is purely due to the compressible fluid motion. In Figures 2(b) and (c) we observe discontinuous regions where the density ρ rapidly changes. Such regions, which are remarkable, particularly in C512-0.3 s, correspond to the shocklet regions (see Figure 3). It was also confirmed in our weakly compressible (Ma = 0.3) case that the shocklet regions form high-density regions downstream, as observed in Yang et al. (2014). To evaluate the relationships between the particle positions and density, we investigate the density at the position of the inertial particles (see Figure 10). Figure 3(a) indicates that the PDFs of velocity divergence are negatively skewed. This is consistent with the previous results (e.g., Lee et al. 1991;Pirozzoli & Grasso 2004;Wang et al. 2012). In Figure 3(a) we observe that the flatness of the PDF increases as M rms , ò r , and Re increase. The M rms and Re dependence of the PDF is consistent with the results in Jagannathan & Donzis (2016). Our PDF with Ma = 0.1 and ò r target = 0 follows approximately a Gaussian distribution: However, the PDF obtained under the same conditions of Ma and ò r target in Jagannathan & Donzis (2016) is more intermittent. The discrepancy in the results is considered to be , as a function of kη for different runs. The value in parentheses in the legend represents that of ò r target . The arrow is in the direction of increasing Re. due to differences in the external forces. The external force used in Jagannathan & Donzis (2016) is random; thus, it significantly changes the large-scale turbulence structure. Conversely, our external force is deterministic and retains the turbulence structure for our PDF, which agrees well with the result of decaying compressible turbulence with Ma = 0.1, shown in Figure 7 of Pirozzoli & Grasso (2004). Our results indicate that the PDF of the velocity divergence similarly depends on ò r . The visualizations in Figures 2 and 3 demonstrate that the sheet-like narrow (blue) regions (called shocklet regions), where the values of the velocity divergence are high and negative, correspond to the discontinuous regions of the density. Jagannathan & Donzis (2016) demonstrated that strong compressions for which θ < − 3θ rms and strong expansions for which θ > 3θ rms occur equally at a low M rms ( < 0.3), where θ rms is the rms value of θ. Following Jagannathan & Donzis (2016), we examined the percentage volume in which the normalized velocity divergence θ/θ rms is within the indicated range (see Table 2). Our results reveal that the volume of strong compressions is larger than that of the strong expansions, even   at lower M rms ( ≈ 0.1) cases (except for C256-0.1n). Note that only solenoidal forcing was used in the DNSs in Jagannathan & Donzis (2016), whereas compressible mode was forced in the cases of ò r target ≠ 0 in our DNS. Therefore the results in Table 2 indicate that the percentage volume of θ/θ rms depends on both the Mach number and the ratio of dilatational to solenoidal forcing. As shown in Table 2, the percentage volume of the compression region for which θ < − 3θ rms increases as the values of M rms and ò r increase; however, it is still smaller than 2%. Samtaney et al. (2001) defined shocklet regions as the regions where θ < − 3θ rms and Δρ = 0. Therefore θ < − 3θ rms is a necessary condition for the region to be a shocklet. However, Figures 2 and 3 demonstrate that such compression regions coincide well with regions in which the density changes sharply. Figure 4(a) shows the PDFs of the local Mach number defined in Equation (8). According to the PDFs, the local Mach number (M loc ) takes values within a certain range that depends on the value of Ma (whereas its ò r dependence is weak). The range is wider for Ma = 0.3 than for Ma = 0.1. The same is also the case for higher Mach numbers (e.g., see Wang et al. 2017). When the Mach number is almost equal to 0.3, the turbulent flow field has regions in which the local Mach number exceeds one, and we can observe that shocklet regions in which θ < − 3θ rms and the fluid density sharply change are spontaneously formed within the turbulent flow field (see Figures 2 and 3).
To understand the difference between incompressible turbulence and compressible turbulence better, we compare their pressure fields.  Table 1 show that the enstrophy decreases with the Mach number and ò r . It is believed that hierarchical multiscale vortex structures occur in turbulence. However, the relationship between large-scale vortex motion and compressibility has not been well studied. To understand this relationship, we study the energy dissipation spectra defined by D(k) = k 2 E s (k) for compressible turbulence, where E s (k) is the energy spectra of the solenoidal component of the velocity, and E s (k) = E(k) in incompressible turbulence. The energy dissipation spectra D(k) satisfy the relationship as á . These relationships reveal that the energy dissipation spectrum D(k) can be regarded as the energy spectrum of the swirling motion in turbulence. In Figure 6 we plot the values of D(k), which are normalized by large-scale characteristics (á ñ u L 2 ). The figure shows that the values of D(k) at small wavenumbers (kL < 50) obtained from the compressible turbulence of Ma = 0.3 are slightly smaller than those obtained from the incompressible turbulence of approximately the same total energy. This result indicates a tendency under the same total energy in which the larger the compressible component of the energy, the weaker the incompressible component of the energy (the swirling motions of turbulence).
This section is summarized as follows. (i) Depending on the values of M rms and ò r , there are noticeable differences between compressible and incompressible turbulent fields, particularly in the density ρ, velocity divergence θ, pressure fluctuation ¢ p , and the dilatational component p d ; however, all of these differences come from the dilatational component of the velocity. (ii) The differences in solenoidal motions, that is, those in p s and D(k), are small (see Figures 5(b) and 6; note that when the total kinematic energies are almost the same, the enstrophy as well as the energy dissipation spectrum at small wavenumbers including the inertial range is slightly smaller in compressible turbulence than in incompressible turbulence.) (iii) Depending on the values of ò r , shocklets are observed in Table 2 Percentage of Normalized Velocity Divergence (θ/θ rms ) in Different Ranges

Particle Distribution
Next, we study the effect of compressibility of turbulence on the statistics of the particle distribution. As the first step, we calculate the radial distribution function (RDF) to study the degree of agglomeration of particles in turbulence. The RDF is defined by where m(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, whereas for a nonuniform distribution of the particles owing to turbulent clustering, we have g(r) > 1. Figure 7 shows the St dependence of the RDF obtained from I1024 and C1024-0.3 s. The result of I1024 is qualitatively consistent with the result shown in Ishihara et al. (2018). The curves g(r) obtained from compressible turbulence are nearly identical to those from incompressible turbulence. However, a close inspection shows that the curves of C1024-0.3 s agree well with those of I1024 for St  0.006 and St  0.06; however, they are slightly smaller than the incompressible case for St = 0.015 and 0.03. The values of St = 0.015 and 0.03 correspond to St η ≈ 1. Thus, these inertial particles are mostly affected by eddies at the Kolmogorov scale. According to Miura (2004) and Jagannathan & Donzis (2016), strong eddies at the Kolmogorov scale become less prominent as the effect of the compressibility increases. This effect is consistent with our DNS result in which the enstrophy, Ω, in C1024-0.3 s is smaller than that in I1024 (see Table 1). Owing to this effect, it is therefore regarded that the aggregation rate of the inertial particles of St η ≈ 1 (St = 0.015 and 0.03) in compressible turbulence becomes slightly lower than that in incompressible turbulence. Note that because µ h St St Re 1 2 , the St values that correspond to St η ≈ 1 decrease as Re increases.
The Stokes number dependence of g(r) at r = 10 −3 L and r = 10 −2 L for different Re values is plotted in Figure 8. For fixed values of St = O(0.1) and r/L, g(r) is a decreasing function of Re (this result is consistent with the result of Ishihara et al. 2018). For a fixed value of r/L, the difference between the values of g(r) for inertial particles in compressible turbulence and those in incompressible turbulence is largest at St η ≈ 1. That is, owing to the decrease in the enstrophy, the values of g(r) for particles around St η = 1 in compressible  turbulence are effectively lower than those in incompressible turbulence. Furthermore, Figure 8 shows that the difference is a decreasing function of Re.
As illustrated in Figure 8, St where the RDF is peaking shifts depending on the value of r/L for which the RDF is evaluated. To avoid this dependence, we also measure in Figure 9 the correlation dimension D 2 (see Martínez & Saar 2002;Chun et al. 2005), which is a fractal dimension uniquely determined by the particle distribution. Here, D 2 is related to the threedimensional RDF as ( ) for r/L = 1, and D 2 is asymptotic to 3 as St → 0 and St → ∞ . Figure 9 demonstrates that St, where D 2 reaches the minimum value, is given by that corresponding to St η ≈ 0.5 − 1, and the minimum value of D 2 obtained for incompressible turbulence is approximately 2.3. These results are in agreement with previous studies Bec et al. 2007;van Aartrijk & Clercx 2008). In the case of compressible turbulence, D 2 around St η = 1 is slightly higher than that of incompressible turbulence. This tendency is consistent with the results shown in Figure 8.
To study the relationships between the compressible turbulence field and spatial positions of the inertial particles, we investigate the ensemble-averaged fluid vorticity w á ñ p and density r á ñ p at the positions of the inertial particles. Figure 10 shows the values of w w á ñ á ñ p and r r á ñ á ñ p with respect to St η , where w á ñ and r á ñ are the spatial averages of ω and ρ, respectively. For comparison, the data by Zhang et al. (2016), who studied the inertial particle statistics in a transonic turbulent flow (M rms = 1.01), are also plotted. Figure 10(a) shows that the values of w w á ñ á ñ p take the minimum value at St η ≈ 0.5 in incompressible cases and have a tendency to increase with the values of M rms and ò r . The former result in which w w á ñ á ñ p takes the minimum value at around St η = 0.5 implies that the inertial particles with St η ≈ 0.5 are most effectively ejected from the center of strong eddies at the Kolmogorov length scales. This is consistent with the results of Sakurai & Ishihara (2018), who studied the average number density of inertial particles in high-vorticity regions. As mentioned previously, high-vorticity regions become less prominent with increasing levels of compressibility of the turbulent field. In addition, in a low Mach number flow, particles generally tend to concentrate in regions of low vorticity owing to the centrifuge effect of the vortices, and the particle concentration decreases monotonically with increasing magnitude of the vorticity (Zhang et al. 2016). Therefore the latter result in which w w á ñ á ñ p increases with the values of M rms and ò r can be explained as a consequence of the reduction of the enstrophy in compressible turbulence because the inertial particles can stay at relatively high-vorticity regions owing to the reduction of the centrifuge effect in compressible turbulence. From this perspective, Figure 10(a) further suggests that increasing the level of compressibility reduces not only the centrifugal forces acting on the dust particles with St η ≈ 0.5, but also those acting on the particles with St η  1. This is consistent with the slight reduction of the swirling motions in compressible turbulence (see Figure 6). Figure 10(b) shows that the values of r á ñ p are increasing functions of M rms and ò r , which is consistent with the results in Yang et al. (2014) and Zhang et al. (2016). In other words, the inertial particles tend to concentrate in high-density regions in compressible turbulence. However, compared to the cases of M rms ≈ 1.0 (Yang et al. 2014;Zhang et al. 2016), for M rms  0.3, the tendency is weak and thus r á ñ p is only slightly larger than r á ñ at a rate lower than 1.02. As for the averaged fluid vorticity at the positions of the inertial particles, we observed systematic differences between weakly compressible turbulence (M rms  0.3) and incompressible   turbulence. These differences can be attributed to the small reduction in enstrophy in compressible turbulence. However, differences in the RDF as well as in the correlation dimension D 2 are quantitatively small and observed only for particles with St η = O(1). As for the averaged fluid density at the positions of the inertial particles, we observed that the differences between weakly compressible turbulence (M rms  0.3) and incompressible turbulence are negligibly small. Therefore it can be concluded that when M rms  0.3, as in the case of the MMSN model, the distribution of particles with  St 1 Re 1 2 is not overly sensitive to the effect of compressibility.

Relative Velocity of Colliding Particles
The relative velocity between a pair of particles at locations X 1 and X 2 with distance r and its rms value is calculated as respectively, where N p is the number of particle pairs at the specified distance r. In Figure 11 we plot the rms relative velocity between the particles at small distances r = 10 −3 L as functions of St. The rms relative velocity (at a small distance specified in terms of r/L) is an increasing function of Re, which is consistent with the result obtained by Ishihara et al. (2018).
At the same value of Re, the rms relative velocity of the particles in compressible turbulence agrees well with that in incompressible turbulence and is lower by more than a factor of two than the conventional theoretical estimate derived by Ormel & Cuzzi (2007). However, a close inspection reveals that the rms relative velocity of particles with high inertia slightly decreases with the increase in the values of M rms and ò r , whereas the rms relative velocity of particles with St η ≈ 1 increases. The former appears to be due to the reduction of large-scale swirling motions in compressible turbulence, as shown in Figure 6. As described in Yang et al. (2014), the velocity of particles of St η ≈ 1 encountering shocklets at t s from upstream (the lower pressure side) slowly increases when t < t s and starts to decrease at t s . Therefore the latter may be explained by the motion of the particles encountering shocklets one after the other because the relative velocity between such particle pairs may increase.

Collision Kernel
Next, we consider a collision kernel that properly incorporates the turbulence clustering effects. The Reynolds number dependence of the collision kernel was studied in Ishihara et al. (2018). Here, we consider the effect of compressibility on the collision kernel. The collision rate per unit volume between two particles with radii a 1 and a 2 is 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 commonly used in dust coagulation models is Γ com = πd 2 〈w 2 〉 1/2 , where d = a 1 + a 2 and 〈w 2 〉 1/2 are the rms relative velocities, which is taken from the model of Völk et al. (1980). According to Saffman & Turner (1956), a cylindrical formulation of the kernel is given by Γ cyl = πd 2 〈|w|〉, while a spherical formulation is given by Γ sph = 2πd 2 〈|w r |〉 with the radial velocity · ( ) | | = -w X X X X w r 1 2 1 2 . In these formulas, the effect of turbulence clustering is not considered. The turbulence clustering effects were considered by Sundaram & Collins (1997) and Wang et al. (2000), and the resulting formulas are 2 . For a low fixed value of d/L, the Reynolds number dependence of the ratio | | ( ) á ñ á ñ =G G w w 2 r 2 1 2 sph com as a function of St was studied in Ishihara et al. (2018). The authors showed that the values of the ratio are constantly lower than unity; therefore the use of 〈w 2 〉 1/2 for the collision kernel leads to an overestimation. Figure 12 shows that the values of the ratio are lower than ( ) p » 8 3 0.92 1 2 , which is obtained if w is Gaussian (e.g., Pan & Padoan 2014). Similarly, the values of the ratio for compressible turbulence are slightly lower than the corresponding values for incompressible turbulence. The latter result implies that the degree of non-Gaussianity of w is slightly stronger in compressible turbulence than in incompressible turbulence. Figure 13 represents the collision kernel G =  is a noticeable increase in the collision kernel within the small St range with an increasing Re, although the Reynolds number dependence of the collision kernel at St  0.1 is weak. Figure 13 shows that the collision kernel is not overly sensitive to the effect of the compressibility. However, a detailed view of Figure 13 reveals that within the large St range (St ≈ 0.1) the collision kernel in compressible turbulence is slightly smaller than that in incompressible turbulence. As observed in Figure 8, the effect of compressibility on g(r) is small within the large St range. However, within the large St range, 〈|w r |〉 in compressible turbulence is slightly smaller than in incompressible turbulence (as shown in Figure 11 in the case of 〈w 2 〉).
Owing to the effect of compressibility that appears in 〈|w r |〉, the collision kernel in compressible turbulence within the large St range is therefore slightly smaller than that in incompressible turbulence.

PDF of the Relative Velocities
To investigate the relative velocities between particles with high inertia in detail, we calculate the PDF of the relative velocity considering a weighting factor that is proportional to the collision rate (see Hubbard 2012 for the weighting factor). In the cylindrical kernel formulation, the collision-rateweighted distribution of the three-dimensional amplitude of the relative velocity is given by u 2 1 2 and ( ) P x is the unweighted PDF. In Figure 14 we plot the collision-rate-weighted PDFs for particles with St = 0.06, 0.15, and 0.3 obtained from I1024 and C1024-0.3 s. The curve of the collision-rate-weighted distribution for particles with St = 0.06 obtained in compressible turbulence agrees well with the corresponding curve obtained in compressible turbulence. However, there is a difference between the curves for the particles with high inertia (St  0.15) in both compressible and incompressible turbulence. The collision-rate-weighted distributions obtained from compressible turbulence are slightly larger than those from incompressible turbulence within the range of | | á ñ  w u 0.5 2 1 2 . This result indicates that the probability of low-speed collisions occurring between particles with high inertia in compressible turbulence is slightly higher than that in incompressible turbulence.
In addition, we studied the differences between compressible and incompressible turbulence in the PDF of the radial relative velocity of the particle pairs. Figure 15 reveals that the PDF obtained in compressible turbulence agrees well with that obtained in incompressible turbulence. This result indicates that the distribution of the radial relative velocity of particle pairs in compressible turbulence can be well fitted by a stretched exponential function (Sundaram & Collins 1997;Wang et al. 2000;Pan & Padoan 2013;Pan et al. 2014b), as confirmed in incompressible turbulence by Ishihara et al. (2018). Figure 15 shows that the PDF of St  0.06 (St  0.03) obtained in compressible turbulence is narrower (wider) than that in incompressible turbulence. This causes a difference between the bouncing and fragmentation probabilities of colliding particles with high inertia (see Figure 17). Figure 11. rms relative velocity between particles at small fixed distance r = 10 −3 L. The values are normalized by á ñ u 2 1 2 for each run. The arrow is in the direction of increasing Re. The black line represents St 1/2 . The gray line shows the Völk-type model by Ormel & Cuzzi (2007). The lower and upper insets are enlarged views at St 0.1 for the runs of N 3 = 256 3 and 512 3 , respectively.

Sticking Rates and Bouncing and Fragmentation Probabilities of Colliding Pairs in the MMSN Model
Next, we consider the effect of the compressibility on the sticking rates and bouncing and fragmentation probabilities of colliding dust particles in protoplanetary disks. For this purpose, we estimate them on the basis of the DNS of compressible turbulence and compare the results with those of incompressible turbulence. According to the MMSN model, the characteristic velocity in a protoplanetary disk is given by U = α 1/2 c s , where α is the turbulence parameter, 1 is the speed of sound, and R is the distance from the central star. The typical value of α is 10 −4 -10 −1 ). When we assume α = 10 −2 and R = 1 au, we then have U = 110 m s −1 . According to Wada et al. (2013), the critical velocity for collisions between equal-sized silicate aggregates is u c = 6 m s −1 . Note that collisions above the critical velocity u c lead to bouncing or fragmentation. Taking U = 110 m s −1 and u c = 6 m s −1 into account, we set the critical velocity of the colliding particles in the DNS as u c /U = 6/110 ≈ 0.0545. By using the collision-rateweighted distribution of the three-dimensional amplitude of the relative velocity given by Equation (16), the sticking rates P stick can be defined as follows: We evaluate the sticking rate P stick by assuming u c = (6/110)U in Equation (17) and by using the DNS data of P cyl . Figure 16 shows the results with different Reynolds numbers as a function of St. This reveals that the sticking rate is a decreasing function of Re and that the decreasing rate is low within the large St range (St  0.1), suggesting that the rate asymptotically approaches a certain value within the large St range. The latter observation is consistent with the observation in Ishihara et al. (2018), that is, the average of the relative velocity as a function of St tends to approach a power law proportional to St 1/2 within the large St range with an increasing Re. We observed that the values of the sticking rate within the large St range in compressible turbulence with Ma = 0.3 are slightly higher than those in incompressible turbulence. This result can be explained by our observation of the collision-rate-weighted distribution of the three-dimensional amplitude of the relative velocity shown in Figure 14.
In Ishihara et al. (2018), the bouncing and fragmentation probabilities of colliding particles were estimated using the PDFs of the normalized radial component of the relative velocity as Pan et al. (2014a), and Ishihara et al. (2018) that the shape of ( ) P w r is approximately symmetric and highly non-Gaussian. The DNSs of incompressible turbulence by Ishihara et al. (2018) illustrated that the values of p B/F are significantly lower than the theoretical prediction based on a Gaussian distribution. Moreover, the values of p B/F for equal-sized particles are much lower than those for different-sized particles. Here, we compare p B/F for equal-sized particles in compressible turbulence with those in incompressible turbulence. As shown in Figure 17, the values of p B/F for the particles with St = 0.15-0.3 in compressible turbulence are slightly lower than those in incompressible turbulence, by a factor of 0.04-0.07. This result can be explained by the difference in the PDF of the radial relative velocity obtained between compressible and incompressible turbulence (see Figure 15).

Conclusions and Discussion
According to the minimum-mass solar nebula (MMSN) model, gas turbulence in a protoplanetary disk is weakly compressible (M rms ; 0.01 − 0.32) and is often regarded as essentially incompressible. However, this treatment has not been justified quantitatively so far. To quantitatively elucidate whether the incompressible limit is adequate for protoplanetary disk simulations, we have compared the motions of inertial particles in DNSs between compressible turbulence and incompressible turbulence. In the DNSs of compressible turbulence, we varied the Reynolds number ( ) = -Re 1000 6000 , the values of the Mach number M rms ( = 0.1 − 0.3), and the ratio ò r ( = 0 − 0.1) of dilatational to solenoidal dissipation.
In the comparison between compressible and incompressible turbulence, we found notable differences in some physical   On the other hand, the following items (i) to (vii) denote small differences between compressible and incompressible turbulence.
(i) The differences in the turbulent fields owing to the solenoidal velocity component, e.g., those in the solenoidal component of the pressure and in the energy dissipation spectrum, are small. (ii) The PDFs of the solenoidal pressure in compressible turbulence agree well with those of the pressure fluctuation in incompressible turbulence. (iii) When the total kinematic energies are almost the same, the enstrophy as well as the energy dissipation spectrum at small wavenumbers including the inertial range is slightly smaller in compressible turbulence than in incompressible turbulence. (iv) For the averaged fluid vorticity at the positions of the inertial particles, the value for weakly compressible turbulence (M rms  0.3) is slightly higher than that for incompressible turbulence. The differences can be regarded as being caused by (iii). As for the averaged fluid density at the positions of the inertial particles, the differences between weakly compressible turbulence (M rms  0.3) and incompressible turbulence are negligibly small. (v) As for the statistics of particle distribution, i.e., the RDF and the correlation dimension, the differences between weakly compressible turbulence and incompressible turbulence are small (despite the difference in the averaged fluid vorticity at the particle positions and the remarkable differences in turbulent fields, which come from dilatational component of the velocity). The differences, which are largest for the particles with St η = O(1), are small (less than 10% at their largest) and become much smaller for the particles with  St 1 Re 1 2 . Note that the RDF is slightly smaller and the correlation dimension is slightly higher in compressible turbulence than in incompressible turbulence. This may be explained as a consequence of (iii), that is, owing to the small reduction in enstrophy in compressible turbulence. (vi) As for the relative mean velocity and the collision kernel, the differences between weakly compressible turbulence and incompressible turbulence are small (less than 7% at their largest). The rms relative velocity of particles at fixed small distance (r = 10 −3 L) in compressible turbulence is an increasing function of Re, which is consistent with the incompressible results by Ishihara et al. (2018). The result indicates that the rms relative velocities from the compressible DNSs are lower by more than a factor of two than those from the closed-form expression derived by Ormel & Cuzzi (2007), regardless of the St number of the particle pairs. The rms relative velocity of particles with high inertia slightly decreases with increasing values of M rms and ò r , whereas the rms relative velocity of particles with St η ≈ 1 increases. The former may be due to (iii). In addition, the latter may be explained by the motions of the particles of St η ≈ 1, which move against Figure 17. Same as Figure 16, but the plot is for bouncing and fragmentation probabilities of colliding particles. shocklets from upstream (the lower pressure side) in compressible turbulence (see the explanation in Figure 11). The collision kernel within the large St range (St 0.1) in compressible turbulence is slightly smaller than that in incompressible turbulence. This may be due to (iii) because (iii) causes the reduction of radial relative velocity of particles with high inertia. (vii) The PDFs of the relative velocities of colliding particles obtained in compressible turbulence agree well with those obtained in incompressible turbulence, and the differences between weakly compressible turbulence and incompressible turbulence in the sticking rates and bouncing and fragmentation probabilities of colliding pairs are small. In the case of Re ≈ 6000, the difference is smaller than 14% for the sticking rates and smaller than 7% for the bouncing and fragmentation probabilities.
Note that the sticking rates of colliding particles are a decreasing function of Re and asymptotically approach a certain value within the large St range. The sticking rates of particles with St = 0.15-0.3 in compressible turbulence are slightly higher than those in incompressible turbulence. Moreover, the bouncing and fragmentation probabilities of colliding particles with St = 0.15-0.3 in compressible turbulence are lower than those in incompressible turbulence by approximately 0.04-0.07. These results suggest that the compressibility of the turbulence has a small positive effect for the colliding particles to preferentially survive without bouncing or fragmentation.
Our results demonstrate that within the range of the Mach number assumed in the MMSN model, compressible turbulence fields have properties that differ from those of incompressible turbulent fields, mainly owing to the existence of the dilatational component of the turbulence. For the behavioral properties of the dust particles, the difference between compressible and incompressible turbulence was found to be small, i.e., (iv) to (vii). This can be explained as follows: (1) When M rms < 0.3, the properties of the solenoidal component of the velocity in compressible turbulence agree well with those of the incompressible turbulence.
(2) When M rms < 0.3, the clustering and collision statistics of the dust particles are mainly determined based on the solenoidal component of the turbulent velocity fluctuation.
Therefore we can conclude that the incompressible assumption is appropriate for the investigation of dust coagulation process in the protoplanetary disk turbulence as assumed in the MMSN model and an appropriate simplifying assumption in analytical theories. Moreover, we observed small but systematic differences in the clustering and collision statistics of dust particles in compressible turbulence from those in incompressible turbulence. The differences owing to the compressibility include the decrease in the degree of agglomeration of particles around St η = 1 as well as the increase and decrease of the rms relative velocity for the particles with low and high inertia, respectively. Within the Mach number range considered in this study, the effects of the compressibility were small. However, it is expected that as the compressibility increases, these effects become large and may dominate the dynamical properties of colliding particles in general turbulent flows.
Our DNSs showed that shocklets can occur in protoplanetary disk turbulence. The relationships between shock waves (strong shocklets) and chondrule formation in protoplanetary disks have been discussed for several years (see, e.g., Hood & Horanyi 1991, 1993Iida et al. 2001;Arakawa & Nakamoto 2019). Because chondrules are spherical, it is believed that once the solid rock components have been melted, they are then solidified again. The melting point of the component is approximately 2000 K. Therefore it is necessary to heat it to a temperature of higher than 2000 K. Furthermore, because volatile components are contained in chondrules, it is assumed that they occur suddenly within a short time, i.e., approximately within minutes. Iida et al. (2001) considered the influential shock-wave heating model as a candidate for the chondrule formation process. That is, chondrules are thought to be rapidly heated by shock waves. Their results show that the optimal shock velocities for chondrule formation are mostly similar to the Keplerian velocities within the asteroidal and Jupiter orbital regions of the MMSN model (Hayashi 1981). The mechanism of shock-wave generation is believed to be due to the eccentric planetesimals excited by Jovian resonances (Weidenschilling et al. 1998) or to the gravitational instability of the protoplanetary disk (Boss & Durisen 2005;Boley & Durisen 2008). Our DNS results suggest that regions with stronger shocks could appear for a realistic Reynolds number. However, the volume fraction of such regions is expected to be significantly too small to account for the observed chondrules.
incompressible DNS, but agrees well with that of compressible DNS. This indicates that the particle statistics analyzed in this study do not significantly depend on the thermodynamics process, including the cooling process. Note that if the Mach number is even higher, e.g., in an interstellar medium, the statistics in nonisothermal turbulence can be different from those in isothermal turbulence (Nolan et al. 2015).