Streaming Instability and Turbulence: Conditions for Planetesimal Formation

The streaming instability (SI) is a leading candidate for planetesimal formation, which can concentrate solids through two-way aerodynamic interactions with the gas. The resulting concentrations can become sufficiently dense to collapse under particle self-gravity, forming planetesimals. Previous studies have carried out large parameter surveys to establish the critical particle to gas surface density ratio ($Z$), above which SI-induced concentration triggers planetesimal formation. The threshold $Z$ depends on the dimensionless stopping time ($\tau_s$, a proxy for dust size). However, these studies neglected both particle self-gravity and external turbulence. Here, we perform 3D stratified shearing box simulations with both particle self-gravity and turbulent forcing, which we characterize via $\alpha_D$ that measures turbulent diffusion. We find that forced turbulence, at amplitudes plausibly present in some protoplanetary disks, can increase the threshold $Z$ by up to an order of magnitude. For example, for $\tau_s = 0.01$, planetesimal formation occurs when $Z \gtrsim 0.06$, $\gtrsim 0.1$, and $\gtrsim 0.2$ at $\alpha_D = 10^{-4}$, $10^{-3.5}$, and $10^{-3}$, respectively. We provide a single fit to the critical $Z$ as a function of $\alpha_D$ and $\tau_s$ required for the SI to work (though limited to the range $\tau_s = 0.01$--0.1). Our simulations also show that planetesimal formation requires a mid-plane particle-to-gas density ratio that exceeds unity, with the critical value being independent of $\alpha_D$. Finally, we provide the estimation of particle scale height that accounts for both particle feedback and external turbulence.


INTRODUCTION
Planet formation involves a variety of physical processes ranging from collisions of tiny dust grains to migration of massive planets within a circumstellar disk.It must be efficient and rapid; it requires that growth from micron-sized dust into fully-fledged planets of up to 10 5 km in scale be completed within several Myr be-fore gas dissipates, possibly even earlier as suggested by some observations (e.g., Manara et al. 2018, Segura-Cox et al. 2020).The intermediate stage of the process, where millimeter-to centimeter-sized pebbles coalesce into kilometer-scale planetesimals, entails one of the most challenging questions to answer: how do the planetesimals form out of their much smaller pebble constituents (see Simon et al. 2022 for a recent review)?
The challenge in answering this question lies with the growth barriers that must be overcome in order to form planetesimals. First, while collisional coagulation of dust grains is effective at producing grains of size ∼mm-cm, growth beyond this scale is stalled by the fragmentation and/or bouncing, the latter of which may even prevent growth even to the fragmentation limit (Zsom et al. 2010, see also Dominik & Dullemond 2023 for more recent work). of these larger grains in the inner regions (≲ 10 au) of protoplanetry disks (PPDs; e.g., Blum & Wurm 2008;Güttler et al. 2010;Zsom et al. 2010;Birnstiel et al. 2012).Second, as particles experience orbit at the Keplerian speed, they feel a headwind caused by the more slowly rotating, sub-Keplerian gas; this headwind removes angular momentum from the particles, causing them to drift radially inward.Since the drift timescale is short compared to the disk lifetime, (e.g., ∼ 300 orbital periods for cm-sized particles at 50 AU; Adachi et al. 1976), the rapid radial drift imposes limits on the growth of particles in the outer regions (Birnstiel et al. 2012).
These growth barriers can be bypassed by the streaming instability (SI, Youdin & Goodman 2005), which arises from angular momentum exchange between gas and solid particles via aerodynamic coupling.Linear studies of the SI (Youdin & Goodman 2005;Youdin & Johansen 2007) have demonstrated that the coupled gas-solid system in disks is unstable, with the exponential growth rate depending mainly on the dimensionless stopping time of particles (τ s ) and particle-to-gas density ratio.Beyond the linear regime, numerical simulations have shown that nonlinear evolution of the SI leads to particles concentrating in narrow filaments, and that under some circumstances, these filaments can reach sufficiently high densities that gravitational collapse ensues, forming planetesimals (Johansen et al. 2007(Johansen et al. , 2009;;Johansen et al. 2015;Simon et al. 2016;Schäfer et al. 2017;Abod et al. 2019).
Several numerical studies have delved into the nature of the SI clumping by examining the critical ratio of pebble surface density to gas surface density (Z), beyond which the SI triggers strong clumping that facilitates planetesimal formation.These studies have explored large regions in (τ s , Z) parameter space (Carrera et al. 2015;Yang et al. 2017;Li & Youdin 2021, hereafter LY21), quantifying a boundary (the "clumping boundary") that determines which combinations of these two parameters could result in strong clumping.They have identified that the lowest critical Z values occur around τ s ∼ 0.1 -0.3, with the exact critical Z values depending on the specific numerical setup (e.g., vertical boundary conditions and domain size; see Li & Youdin 2021, Section 4.1.4).While thresholds for the SI clumping have been established for much of parameter space, there is still work to be done to understand this boundary un-der the presence of more realistic physics, e.g., including particle self-gravity and external turbulence.
First, the previous SI simulations did not take (driven) external turbulence 1 into account.Observations show that, although PPDs are weakly turbulent compared to predictions for fully-ionized accretion disks, there is evidence for a non-zero level of turbulence that varies substantially between disks.Observations of turbulent line broadening derive an upper limit of turbulent velocity to sound speed ratio δv/c s ∼ 0.01 for HD 163296 (Flaherty et al. 2017; see also Flaherty et al. 2018 andFlaherty et al. 2020 for other disks), whereas other disks (DM Tau - Flaherty et al. 2020 andIM Lup -Paneque-Carreño et al. 2023) show δv/c s ∼ 0.3 above 1-2H (H being vertical scale height of gas) away from the midplane.Observations sensitive to the vertical settling of dust grains suggest even weaker turbulence close to mid-plane in Class II disks (i.e., δv/c s ≲ 0.003 in the outer regions of Oph163131 disk, assuming a turnover timescale of turbulent eddies similar to the local orbital timescale; Villenave et al. 2022), or more modest turbulence in Class I disks (e.g., δv/c s ≳ 0.01 in the outer regions of IRAS04302 disk under the same assumption for the turnover timescale and assuming τ s = 0.01; Villenave et al. 2023).Clearly, the strength of turbulence (at least as inferred from these observations) varies between disks.It is therefore crucial that we understand the effect of turbulence on the SI clumping criterion and planetesimal formation, across a range of turbulence levels from the minimum set by the SI itself up to at least the levels inferred in DM Tau and IM Lup.
Second, in previous studies of the clumping boundary, particle self-gravity has been neglected, with the threshold for planetesimal formation being instead defined by whether or not the maximum particle density would exceed the Hill density (a condition often referred to as "strong clumping"; LY21).This choice is largely justified since in the presence of very weak turbulence (such as that produced by the SI itself), the Hill criterion is likely a sufficient condition for collapse; a particle cloud (i.e., local particle overdensity) within a narrow filament needs to be sufficiently dense to be stable against tidal shear.However, if turbulence becomes important, gravity has to overpower not only tidal shear but also the dif-fusion to collapse particle clouds under their own weight.Thus, as both turbulent diffusion and tidal shear can counteract self-gravity, a Toomre-like criterion for gravitational instability becomes a more appropriate criterion for collapse, with diffusion playing a similar role to pressure (Gerbig et al. 2020;Klahr & Schreiber 2020;Gerbig & Li 2023).
These considerations strongly motivate an inclusion of both external turbulence and particle self-gravity in determining the clumping boundary.Thus, in this paper, we incorporate both the self-gravity of particles and forced turbulence into 3D gas-plus-particle numerical simulations.We carry out 41 simulations in order to quantify the clumping boundary, here defined as the boundary above which gravitationally bound clumps form in our simulations.
We note that a number of previous works have addressed particle concentration and planetesimal formation under the influence of turbulence.For example, Chen & Lin (2020) and Umurhan et al. (2020) analytically show that even moderate level of turbulence (i.e., α SS ≳ 10 −4 ; where α SS is standard Shakura-Sunyaev turbulent viscosity parameter; Shakura & Sunyaev 1973) can suppress linear growth of the SI.The notion that turbulence may hinder the SI was numerically confirmed by nonlinear SI simulations with externally driven hydrodynamic turbulence (Gole et al. 2020).Moreover, there have been number of studies that investigate the effect of turbulence driven by (magneto-)hydrodynamical instabilities present in PPDs, such as the MRI (Johansen et al. 2007;Yang et al. 2018;Xu & Bai 2022a) or the vertical shear instability (VSI; Nelson et al. 2013; see Schäfer & Johansen 2022 for the VSI and SI study), on the particle concentration.The main takeaway from these works is particle concentration can occur in large-scale features (such as localized concentrations of gas pressure) that emerge naturally from the (magneto-)hydrodynamic proccessat at work.We will return to a discussion of these results later in this manuscript, but for now it is worth mentioning that these numerical studies explored a relatively narrow region in parameter space in terms of τ s and Z.Thus, the influence of turbulence has yet to be unveiled in a broader parameter space.Exploring such a parameter space with turbulence (and particle self-gravity) is the goal of this paper.
This paper is organized as follows.We describe our numerical approach in Section 2. In Section 3.1, we examine the effect of turbulence on particle concentration.Section 3.2 is dedicated to the simulations with the particle self-gravity included.Thus, it is in this section where we present a new clumping boundary.We demonstrate the effect of particle feedback in Section 3.3 and its influence on vertical profiles of particle density in Section 3.4.Finally, we discuss and summarize our results in Sections 4 and 5, respectively.

Numerical Method
We use the Athena hydrodynamics+particle code (Stone et al. 2008;Stone & Gardiner 2010;Bai & Stone 2010b) to perform 3D simulations of the SI.We artificially force turbulence (see Section 2.3 for details) in a vertically stratified shearing box, which is a co-rotating patch of a disk sufficiently small such that the domain can be treated in Cartesian coordinates without curvature effects (see e.g., Hawley et al. 1995).More specifically, we use a local reference frame at a fiducial radius R 0 that rotates at angular frequency Ω.The equations of gas and particles are written in Cartesian coordinates (x, ŷ, ẑ), where x, ŷ, and ẑ denote unit vectors pointing to radial, azimuthal, and vertical directions, respectively, and the local Cartesian frame is related to the cylindrical coordinate of the disk (R, ϕ, z ′ ) by Solid particles are included and treated as Lagrangian super-particles, each of which is a statistical representation of a much larger number of particles with the same physical properties.The particles are coupled to the gas via aerodynamic forces.We adopt periodic boundary conditions in the azimuthal and vertical directions, while shearing periodic boundary conditions (Hawley et al. 1995) are used in the radial direction.
To solve the hydrodynamic fluid equations, we use the unsplit Corner Transport Upwind (CTU) integrator (Colella 1990), third-order spatial reconstruction (Colella & Woodward 1984), and an HLLC Riemann solver (Toro 2006).The hydrodynamic equations in the shearing box approximation are where ρ g is the gas mass density, u is the gas velocity, and c s is the sound speed.In the momentum equation (Equation 2), I is the identity matrix, and P is the gas pressure defined in Equation (3).On the right Lim et al.
hand side of Equation (2), the first three terms are radial tidal forces (gravity and centrifugal), vertical gravity, and Coriolis force, respectively.The penultimate term denotes the back-reaction from the particles to the gas where ρ p , v v v, and t stop are the particle mass density, particle velocity, and stopping time of particles, respectively.The last term f turb is the forcing term for turbulence (see Section 2.3).Note that we use an orbital advection scheme (Masset 2000), in which the Keplerian shear u K ≡ −(3/2)Ωxŷ is subtracted and we only evolve the deviation u ≡ u tot − u K where u tot is the total gas velocity.The last of the above equations is our isothermal equation of state, which we assume here for simplicity; thus, c s is constant in space and time.
The equation of motion for particle i (out of N par total particles) is written as and solved with a semi-implicit integrator (Bai & Stone 2010b).The orbital advection scheme, which is applied to the gas dynamics, is implemented for the particles as well.In Equation ( 4), the first, the second, and the third terms on the right hand side are radial tidal forces (again, gravity and centrifugal), vertical gravity, and the Coriolis force, respectively.The fourth term is the drag force on individual particles.The second-tolast term represents a constant inward acceleration of particles that allow them to drift radially inward in the presence of the radial gas pressure gradient (see Bai & Stone 2010b).In our treatment (and again following Bai & Stone 2010b), the parameter η is the fraction of the Keplerian velocity by which the orbital velocity of particles is effectively increased (i.e., made super-Keplerian), while the gas stays at Keplerian velocities; the frame is effectively boosted by an amount ηu K , which won't change the essential physics.The last term is the particle self-gravity, which is calculated by solving the following Poisson equation: where G is the gravitational constant.We use fast Fourier transforms (FFT) to solve the Poisson equation for the gravitational potential (see Simon et al. 2016 for more details on the implementation of the Poisson solver) from which we then calculate the selfgravitational acceleration (a g = −∇Φ).Gas self-gravity is ignored in our simulations.We use the triangular-shaped cloud (TSC) method outlined in Bai & Stone (2010b) to interpolate the gas velocity at the grid cell centers to the particle locations for the calculation of the gas drag term in Equation ( 4).The same interpolation scheme is used to map momenta from the particle locations to the grid cell centers to calculate the back-reaction term on the gas in Equation (2).The TSC method is also used to calculate a g ; the mass density of particles is mapped to grid cell centers to solve Equation ( 5) and calculate a g , which is then interpolated back to the location of particles.

Initial Conditions and Parameters
Due to the typical length scales of the SI being much smaller than the vertical gas scale height (H = c s /Ω), our domain size is necessarily smaller than shearing box simulations of other processes, such as magnetically driven turbulence (see, e.g., Hawley et al. 1995).More specifically, our domain size is (L x , L y , L z ) = (0.4,0.2, 0.8)H, where L x , L y , and L z are radial (length), azimuthal (width), and vertical (height) sizes of the boxes, respectively.For the runs with smaller τ s and/or stronger forced turbulence, we increase L z to reduce the number of particles crossing the vertical boundaries (see Table 1).The number of grid cells in each direction is (N x , N y , N z ) = (256,128,512), equating to 640 grid cells per H in each direction, respectively.We adjust the vertical resolution for the taller boxes so as to maintain an equivalent resolution to 640 cells per H.
The dynamics of gas and particles in our simulations are determined by five dimensionless parameters.We focus on the four parameters relevant to the particles here and describe the other one characterizing the forced turbulence in the next subsection.
First, the aerodynamic coupling of the particles to the gas is controlled by the dimensionless stopping time which also represents particle size (t stop is proportional to the grain size; see, e.g., equation 1.48a in Youdin & Kenyon 2013 for the formula of t stop used in this work.).
In each simulation, all particles are the same size, but we vary τ s in different simulations, exploring values of τ s ranging from 0.01 to 0.1.For reference, these τ s values correspond to particle sizes of millimeters to centimeters at 50 AU in a protoplanetary disk with reasonable choices for disk mass, disk size, etc. (e.g., Carrera et al. 2021), though there is some variation in these numbers depending on the disk model employed.
Second, the abundance of the particles relative to the gas is characterized by the ratio of initial surface density of particles (Σ p0 ) to that of gas (Σ g0 ≡ √ 2πρ g0 H): where ρ g0 is the initial gas density.The surface density ratio sets the strength of particle feedback on the gas as well (Equation 2).We consider multiple Z values spanning from 0.015 to 0.4.Fourth, we control the strength of the particle selfgravity by using the dimensionless parameter Here, Q is the Toomre parameter (Toomre 1964) for the gas disk.We fix G = 0.05 (Q ≃ 32), which allows us to compare our results to previous numerical studies of the SI.The assumed Q value gives the Hill density2 We set the number of particles as follows: where N cell is the total number of cells.The prefactor on the right hand side results from setting n p = 1 as the number of particles per cell.Bai & Stone (2010b) show that setting n p to 1 is necessary to accurately capture density distributions of particles in unstratified SI simulations.However, as LY21 pointed out, the effective particle resolution will be > 1 for vertically stratified simulations since particle settling leads to a particle layer thickness H p < L z .Since our simulations include forced turbulence, the effective particle resolution is not known in advance.Nonetheless, we expect that the effective resolution is higher than 1 for all of our simulations since H p < L z is always satisfied in our simulations.Furthermore, we test how the effective particle resolution impacts our results in Section 4.3.2.
The gas is initialized as a Gaussian density profile with vertical thickness H in hydrostatic balance.The particles are initially distributed with a Gaussian profile about the mid-plane in the vertical direction; the scale height of this Gaussian is 0.02H.In the horizontal directions, initial positions of particles are randomly chosen from uniform distribution.The initial velocities of the gas are zero (with Keplerian shear subtracted), whereas the particles have initial azimuthal velocity of ηu K = 0.05c s .Importantly, the particles are not initially in an equilibrium state.This is because particles have an initial azimuthal velocity causing them to drift radially inward (and we do not assume the Nakagawa-Sekiya-Hayashi equilibrium condition; Nakagawa et al. 1986).Moreover, forced turbulence already exists from the initialization of particles (next section), and the turbulent diffusion will not immediately counterbalance the vertical gravity from a host star.All of our simulations are shown (with associated relevant parameters) in Table 1.

Turbulence Forcing
As mentioned earlier, we inject turbulence into our simulation domain.Unlike Gole et al. (2020) in which turbulence is driven in Fourier space and added to real space via an inverse Fourier transform, we drive turbulence in real space with a cadence of t drive = 0.001Ω −1 .We compared one of our runs to that in Gole et al. (2020), both of which have the same τ s , Z and comparable kinetic energy of gas (Run T30Z2A3.5), and found no significant differences.We use a vector potential method to force turbulence, which guarantees that the velocity perturbations introduced into the domain are incompressible (i.e., divergence-free): where A is the vector potential, and Λ is a forcing amplitude that determines the velocity magnitude of the forced turbulence.The amplitude does not change with time for a constant energy input at every t drive .The velocity perturbations are obtained by taking a curl of A numerically in a way that guarantees that f turb is a cellcentered quantity.A and thus f turb are sinusoidal with a phase that varies randomly every time the forcing is done.We elaborate on the equations for A and the way we handle the shearing-periodic boundary conditions in (1): run name (the numbers after T and Z are in units of one hundredth, while those after A are the absolute values of power indices; (2): dimensionless stopping time of particles (see Equation 6); (3): initial surface density ratio of particle to gas (see Equation 7); (4): dimensionless diffusion parameter of turbulence (see Equation 13); (5); dimensions of the simulation domain in unit of gas scale height; (6): the number of grid cells in each direction; (7): number of particles; (8): when particle self-gravity is switched on; (9): whether or not gravitational collapse of particles occurs (Y for Yes and N for No) (10); time-averaged standard deviation of particles' vertical positions (see Equation 20); (11): time-averaged particle density at the disk mid-plane; (12): time-averaged maximum particle density; (13): time interval over which quantities in column ( 10)-( 12) are averaged.We do not report the quantities for runs that have relatively short pre-clumping phases.All runs have the global radial pressure gradient of Π = 0.05 and the self-gravity parameter of G = 0.05.Time evolution, vertical density profile, and final 2D snapshots of particle density for each run listed in the table are available at https://github.com/simon-research-group/JeonghoonPublic.* We decided not to turn on particle self-gravity in these runs.Nevertheless, we conclude that they are not capable of forming planetesimals even if the gravity is turned on because a run with one step higher Z value but identical τs and αD (e.g., T1Z10A3 vs T1Z20A3) does not lead to planetesimal formation with the self-gravity on.** Since this run is to study the effect of the number of particles on the particle concentration (see Section 4.3), we do not include particle self-gravity here.
Table 2. Summary of our forcing parameter and relevant quantities from τs = 0.1, Z = 10 −5 simulations Note-We report αD, αD,z, and δu ′ values that result from turbulence forced with the interpolated Λ (Figure 1).For the SI simulations listed in Table 1, we use these Λ values as the initial condition for the forcing amplitude for corresponding αD.
Appendix A. We stress that while our perturbations are initially incompressible, there is no guarantee that the injected turbulence maintains a divergenceless (i.e., incompressible) velocity field.However, we have verified that in the absence of particles, the divergenceless component of the velocity field accounts for ∼ 99% of the total velocity field (see Appendix A for details).
We use a parameter α D to quantify the level of the forced turbulence throughout this paper.The α D is by definition a spatial diffusion of turbulent gas in a dimensionless form as follows: where D g ≃ (δu ′ ) 2 t eddy is a dimensional version of α D representing gaseous diffusion due to velocity fluctuations (Fromang & Papaloizou 2006;Youdin & Lithwick 2007): , and τ eddy ≡ t eddy Ω is the dimensionless eddy turnover time of turbulence.We adjust Λ to obtain our target values of α D , which are 10 −4 , 10 −3.5 , and 10 −3 .In the following, we explain how we obtain the forcing amplitude (Λ) for the three α D values.
First of all, we assume that in the limit ρ p ≪ ρ g , the particle scale height H p can be determined by the balance between the particle settling and vertical diffusion by turbulence (Youdin & Lithwick 2007): where α D,z corresponds to the α D for vertical diffusion only.Second, to obtain H p /H and the resulting α D,z , we perform four simulations, each of which has a different Λ value, the dimension of (L x , L y , L z ) = (0.4,0.2, 0.8)H, and the resolution of (N x , N y , N z ) = (256,128,512).We set τ s = 0.1 and Z = 10 −5 in these simulations; the very  12), while the latter is the parameter (Equation 13) we use to denote the level of turbulence in our SI simulations.The Λ values shown in this plot are obtained by a linear interpolation (see the text for details).Each red horizontal line indicates each αD we choose: 10 −4 , 10 −3.5 , and 10 −3 .The error bar shows ±1 standard deviation.We confirm that the interpolated Λ values result in the desired αD values with high accuracy.
small Z guarantees ρ p ≪ ρ g and that the Equation ( 14) holds true (self-gravity of particles is deactivated).Once each simulation reaches a saturated state, we average H p /H over time for 250Ω −1 and obtain α D,z .Third, since α D,z accounts for the diffusion in the vertical direction only, we relate α D,z to α D by where l ≡ δu ′ x /δu ′ z and m ≡ δu ′ y /δu ′ z .In other words, l 2 +m 2 +1 = 3 means isotropic turbulence 3 .We timeaverage l and m within the same time interval that we use for time-averaging H p /H.The eddy turnover time (τ eddy ) is assumed to be isotropic so that τ eddy = τ eddy,i ; i = x, y, z.In this manner, we can calculate α D in each of the four simulations without needing to know 3 We note that the forced turbulence is not perfectly isotropic.The anisotropy is mainly caused by δu ′ y which is systematically lower than the other two components.We believe this is because eddies with large τ eddy are more easily destroyed by orbital shear which acts along y direction.If this is true, δu ′ y can be lower than the other two components since larger eddies contribute to the kinetic energy of turbulence more according to Kolmogorov turbulence theory.It is also possible that since Ly < Lx, the narrower side limits sizes of eddies along the y direction and in turn lowers δu ′ y .
Lim et al.
τ eddy so that we establish the relation between Λ and α D .Using the relation, we perform linear interpolation between the data points (i.e., α D at each Λ) to find the appropriate Λ that is expected to produce the desired α D values.We perform three additional simulations, each of which has Λ obtained from the interpolation.Then, we measure α D via the same procedure we describe above to confirm whether the interpolated Λ results in the desired values of α D . Figure 1 shows the resulting α D for each interpolated Λ value, which is summarized in Table 2.The error bar shows ±1 standard deviation due to the temporal fluctuation of H p and gas velocities.Clearly, we achieve our desired α D values, denoted by red vertical lines, with very high accuracy.These interpolated Λ values are adopted as the initial condition for the forcing amplitude and kept as constant in our SI simulations.However, instead of Λ, we use α D to denote the level of the forced turbulence to make contact with standard disk dynamics notation.Hence, each of our SI simulations has an unique combination of (τ s , Z, α D ) as listed in Table 1; "No forcing" run refers to a simulation where f turb = 0.
We clarify that the α D (α D,z ) measures bulk (vertical) diffusion in the gas set by (vertical) velocity fluctuations and the eddy turnover time of the forced turbulence.As noted by previous studies (e.g., Youdin & Lithwick 2007;Yang et al. 2018), one should not interpret the parameter as α SS , which is responsible for angular momentum transport due to turbulent shear stress in a disk.
Before initializing particles in the SI simulations, we force turbulence with the interpolated Λ values up to t par = 300Ω −1 in order to let turbulence fully develop without being affected by particles and also to reach a statistically steady state.At t = t par , we initialize particles as described in Section 2.2.Before the initialization of particles, we verified that the Fourier spectrum of (δu ′ ) 2 (see Appendix A) is reasonable with power across a range of scales, and with the most power at the driving scale (∼ 0.1H).In the no-forcing run, particles are initialized at the beginning.We use the notation of t−t par (in units of Ω −1 ) for our temporal dimension; particle initialization is done at t−t par = 0.

RESULTS
We present a summary of statistics of our simulations in Table 1.We report temporal averages of particle scale height, maximum particle density, and mid-plane density ratio of the particle and the gas.We adopt a measurement strategy similar to that of LY21 to facilitate comparison with their results: quantities are averaged after the particles have settled from their ini-tial positions and are in a statistical equilibrium (vertically) against turbulent diffusion.The time average is done either before the maximum density first exceeds (2/3)ρ H (termed as pre-clumping phase, LY21), or before the particle self-gravity is turned on (hereafter, t sg ), whichever occurs first.Table 1 presents t sg of simulations at which the self-gravity is implemented.We opt not to include the self-gravity in three of our runs as they seem incapable of forming planetesimals.This choice is made when the run with one step higher Z value, but identical τ s and α D , does not result in planetesimal formation.Since the maximum density of these three runs never reaches (2/3)ρ H , the time-averaging is done all the way to the end of the simulations.We exclude from our reported statistics any simulation where strong particle concentrations happens too rapidly.More specifically, we exclude simulations where a quasi-steady state only lasts for a few tens of Ω −1 or is never achieved before the maximum density reaches the threshold.

Effect of Turbulence on the Particle Concentration
We show the effect of the turbulence on particle dynamics without particle self-gravity before presenting more detailed analysis of our results.Figure 2 shows the results of simulations with τ s = 0.1, Z = 0.02 but for different α D values.From top to bottom, gas rms velocity (black) and maximum density of particles (orange) as a function of time, azimuthally averaged, and vertically integrated particles densities (Σ p = ρ p dz) are shown.The horizontal lines in the top panels denote Hill density, ρ H ≃ 180ρ g0 .
As can be seen from the top panel, gas rms velocity for α D = 10 −4 , 10 −3.5 , and 10 −3 levels off at ∼ 0.01c s , ∼ 0.03c s , and ∼ 0.05c s , respectively, while that of the no-forcing run stays only at ∼ 0.006c s .These values were calculated via time averaging between t s and t e , where this time-averaging interval is determined based on the criterion outlined above (also see Table 1).This directly affects the evolution of maximum density of particles; the peak density is very close to or exceeds the Hill density in the no-forcing and α D = 10 −4 runs, while it saturates at only ∼ 10ρ g0 in the other two runs due to the stronger turbulence.Interestingly, weak turbulence (i.e., α D = 10 −4 ) seems to enhance the particle concentration more than the no-forcing run does.However, given the stochastic nature of clumping, it is hard to draw any firm conclusions about whether α D = 10 −4 actually produces stronger clumping than the same system without forced turbulence.
The middle and the bottom panels in Figure 2 presents snapshots for each run at t−t par = 350Ω −1 that reveal the degrees of particle concentration.First, the  10).
Particle self-gravity is disabled in the runs presented here.As turbulence becomes stronger from left to right, the particle layer become thicker, and the filaments become more diffused and less well-defined.
side view of the particle layer (i.e., middle panels) clearly shows that the turbulence vertically stirs particles and thus thickens the particle layer.The particle scale heights at t−t par = 350Ω −1 are ∼ 0.008H, ∼ 0.01H, ∼ 0.02H, and ∼ 0.05H for the no-forcing, α D = 10 −4 , 10 −3.5 , and 10 −3 runs, respectively.Second, as can be seen from the bottom panels, the noforcing and α D = 10 −4 runs have azimuthally elongated particle filaments, which implies that the SI is active.
Even though the filaments form in both runs, the latter run has fewer filaments than the former run.The α D = 10 −3.5 run shows marginal filament formation, and the α D = 10 −3 run shows no evidence for filament formation, instead showing a very diffused particle medium.
The two-dimensional distributions of particle density in Figure 2 suggest that turbulence can weaken or even completely suppress the SI.When turbulence is weak or moderate (i.e., α D = 10 −4 and 10 −3.5 ), the SI forms elongated filaments.However, the resulting filaments are fewer and less dense than the no-forcing case.Conversely, the α D = 10 −3 case shows no filamentary structures at all (see the bottom right panel).This could be attributed to very strong vertical stirring that prevents particle settling.It is also possible that the SI is active but that filament formation is overpowered by destructive diffusion.The latter scenario is evident in Runs T1Z10A3.5 and T2Z4A4.Despite these runs showing ρ p ≳ ρ g at the mid-plane, as indicated in Table 1 and Figure 4, the presence of the turbulent diffusion prevents the SI from forming dense particle filaments, consequently inhibiting planetesimal formation.

Planetesimal Formation
In this section, we present results from the simulations that incorporate both forced turbulence and particle self-gravity.However, before discussing the results, we describe how we differentiate between "Collapse" and "No collapse" runs.The left panel of Figure 3 (which corresponds to the run with τ s = 0.1, α D = 10 −3.5 ) shows the maximum particle density as a function of time for Z = 0.02 (green) and Z = 0.032 (blue) runs, the former being "No collapse" and the latter being "Collapse".The horizontal line corresponds to the Hill density, and the vertical line indicates t sg in both simulations.The "Collapse" run shows the maximum density sharply increasing by more than a factor of 10 right after t sg , which is evidence for the gravitational collapse of particles.On the other hand, the "No collapse" run shows the steady evolution of the maximum density even though the density slightly increases upon turning on the self-gravity.Furthermore, we investigate the spatial distribution of the particle density as well.The middle and the right panels of the figure show the final snapshots of the surface density of the "No collapse" and "Collapse" runs, respectively.The right panel reveals gravitation-  ally bound objects, which we call planetesimals4 , while the middle panel has no such objects but shows weak filaments have formed.In summary, we consider both the temporal evolution and the spatial distribution of particle density to categorize runs into "Collapse" and "No collapse".

Density at Midplane
In unstratified SI, the mid-plane density ratio of particle and gas, ϵ ≡ ρ p /ρ g (z = 0), is a crucial parameter.More specifically, when τ s ≪ 1, the linear growth rate of the SI increases with ϵ, with a sharp increase as ϵ approaches and surpasses unity (Youdin & Goodman 2005).Therefore, ϵ ≳ 1 is often assumed as a condition for the SI to produce dense clumps of particles.While this could be true in the linear regime of unstratified SI, this is not necessarily true in the fully nonlinear stratified SI, as pointed out by Yang et al. (2018).Numerical simulations have shown that critical ϵ values can deviate from unity depending on τ s and other factors.More specifically, LY21 reported a critical ϵ (for strong clumping to occur) from ≈ 0.3 to ≈ 3 depending on the value of τ s .This is consistent with previous work by Gole et al. (2020) in which the critical ϵ (for planetesimal formation to occur) is ≈ 0.5 in the presence of external turbulence.Here, we follow up on this work and fur-ther examine the critical ϵ for planetesimal formation to occur (since we include particle self-gravity) but for a wider range of parameters than Gole et al. (2020).
Figure 4 shows temporally and horizontally averaged values of ϵ from the simulations for which we could take time sufficiently long time averages during the preclumping phase (i.e., not every simulation in Table 1 is shown).We calculate ϵ by taking particle and gas densities at ±1 one grid cell above and below the midplane."Collapse" and "No collapse" runs are denoted by circles and triangles, respectively.Red, green, and blue colors denote α D = 10 −4 , 10 −3.5 , and 10 −3 , respectively; we did not color-code Z values.The black curve is the best fit by least squares (see below) to the data marking the approximate location of the critical ϵ value, above which collapse occurs.In choosing this fit, we assume a quadratic function in log-log space as in LY21 but with different coefficients.More precisely, we find a critical ϵ of: where A = 0.41, B = 0.71, C = 0.37 for all α D values.
We also include a fit to the critical ϵ in LY21 shown as the skyblue curve (ϵ crit,LY21 ).We calculate this least squares fit as follows.We assume that at a given τ s and α D , a critical ϵ falls in the middle between the adjacent no-collapse data point (triangle) and the collapse data point (circle).Moreover, we take the range between the adjacent no-collapse and collapse points as the 95% confidence interval for the location of ϵ crit .Therefore, we have a data point (i.e., a critical ϵ) and its associated error at each τ s and α D .These are then input into a weighted least squares algorithm that accounts for the varying error sizes (het-  Temporally and horizontally averaged ratio of particle density to gas density at the mid-plane as a function of τs in our simulations.We do not include simulations where the pre-clumping phase is too short to take time-averages.Runs where the collapse occurs are shown as circles, whereas those where the collapse does not occur are shown as triangles.Red, green, and blue represent αD = 10 −4 , 10 −3.5 , and 10 −3 , respectively, while we did not denote Z values.The black curve denotes the least squares best fit to the critical value of ϵ as described in Equation ( 16).The thin grey curves are random fits drawn from multivariate normal distribution of A, B, C (Equation 16) that show uncertainties in the best fit curves.We also include the critical curve from LY21, depicted as a sky blue curve (ϵcrit,LY21).The figure shows that ϵcrit does not vary with αD (i.e., even changing αD, every no-collapse case is below every collapse case), which results in a single critical curve that fits all ϵ values regardless of αD.
eroscedastic errors).The thin, grey curves in Figure 4 are ten random sample fits, each of whose coefficients are drawn from a multivariate normal distribution of (A, B, C).The width of the distribution in each "dimension" (i.e., variable) is taken from the covariance matrix from the fit.Those random samples thus provide a sense of the uncertainty in the best fit curve.Figure 4 has several implications.First, as α D changes, every collapse run remains above every nocollapse run.Thus, the ϵ crit curve precisely cuts between the no-collapse and collapse points regardless of α D values.This implies that ϵ crit varies with τ s but not significantly with α D .This could potentially be at-tributed to the use of higher Z values for larger α D at a given τ s .For instance, at τ s = 0.1, the Z values of the runs just above the curve are 0.02, 0.024, and 0.045 for α D = 10 −4 , 10 −3.5 , and 10 −3 , respectively.In other words, although larger α D makes a particle layer thicker, using higher Z adds more mass of particles within the layer, resulting in similar particle densities around the mid-plane (i.e., similar ϵ) in the three α D cases.However, we again emphasize that the accuracy of our fit is limited due to the sparsity of data points across the τ s range considered in this study.Consequently, the fit could potentially exhibit variation with α D if additional  19shows the formula), whereas the skyblue one uses the ϵcrit from Equation ( 16) to calculate the critical Z assuming that the vertical profile for particle density is Gaussian.
numerical simulations were to be conducted to further populate parameter space.Second, based on Equation ( 16), the values of ϵ crit approximate to 3.98, 1.75, and 1.18 for τ s = 0.01, 0.03, and 0.1, respectively.These values are several times larger than those from ϵ crit,LY21 5 .Moreover, since we have only one simulation at τ s = 0.3, ϵ crit at this τ s value remains unknown in our work.Moreover, the majority of the data points, regardless of whether or not a corresponding run shows planetesimal formation, are well above ϵ crit,LY21 .A potential explanation for this is that even when particles are able to settle to the mid-plane and form a layer to have ϵ ≳ 1, further concentration may be required for them to withstand the turbulence that disperses particles (through aerodynamic coupling with the gas) in all directions.
Lastly, since ϵ crit we report in this paper is the critical value for gravitational collapse of particle clumps 5 LY21 found similar ϵ crit with Gole et al. (2020) at τs = 0.3, whereas our results are inconsistent with LY21.While this may imply that the two different forcing methods produce inconsistent results, we note that Gole et al. ( 2020) used ϵ crit ≈ Z/(Hp/H) ≈ Z(τs/α crit ) 1/2 ∼ 0.5 with α crit = 10 −3.25 .This calculation for ϵ crit assumes ρp ≪ ρg, which is not always the case in our simulations, especially close to the mid-plane where ϵ ≳ 1; see Figure 4 rather than just for the SI-induced concentration, the condition that ϵ simply be greater than unity does not necessarily guarantee gravitational collapse (see Gerbig et al. 2020;Klahr & Schreiber 2020, 2021;Gerbig & Li 2023 for details of the collapse criterion).This is because conditions for gravitational collapse of a particle cloud should be dependent on internal (to the pebble cloud) turbulent diffusion of particles within the cloud as well as its density.

Threshold for Planetesimal Formation: Critical Z
Our main results are shown in Figure 5.It demonstrates for which parameters (τ s , Z, α D ) planetesimals form.In the figure, filled and open circles correspond to "Collapse" and "No collapse" runs, respectively.All simulations in the figure maintain a constant Π and G, both of which are 0.05.We do not show Run T10Z2 in the figure, which includes self-gravity but not forced turbulence (see Table 1 for the details of this simulation).
The figure demonstrates that planetesimal formation via the SI may be very difficult in the presence of external turbulence.Taking the τ s = 0.01 cases for example, the critical Z values are ≳ 0.06, ≳ 0.1, and ≳ 0.2 for α D = 10 −4 , 10 −3.5 , and 10 −3 , respectively.On the other hand, Z ∼ 0.02 is enough for the SI to produce dense clumps in the absence of external turbulence for this value of τ s (Carrera et al. 2015;Yang et al. 2017, LY21), which would lead to planetesimal formation if self-gravity of particles was activated (none of these studies included particle self-gravity).
We also plot the critical Z curve assuming a Gaussian scale height for the particles (Z crit,Gaussian Hp , LY21) and using ϵ crit as fit by Equation ( 16) in sky-blue.The curve is given by where where H p,SI = h η ηr, and H p,α /H = α D,z /(α D,z +τ s ).
The former represents the contribution of the SI-driven turbulence to the particle scale height, whereas the latter represents that of externally driven turbulence.We adopt the same value for h η as in LY21, which is ≃ 0.2.This value is an approximation for the particle scale height within the range of τ s (and in units of ηr) as measured from their SI simulations.As a result, H p,SI /H ≃ Π/5.To obtain H p,α for each α D , we use α D,z in Table 2.
As can be seen from Figure 5, the critical curve (skyblue, Equation 17) and the data are not consistent at all; we find systematically lower critical Z values than those calculated from Equation ( 17), implying that Equation (18) does not accurately predict the actual scale height of particles in simulations.The explanation for the inconsistency is given in Sections 3.3 and 3.4, where we focus on the indirect impact of particle feedback on the particle scale height and the vertical profiles of particle density.
Since Equation ( 17) does not match our simulation results, we attempt to provide a new fit to critical Z values.Assuming log Z crit is a function of both log τ s and log α D : with conditions of α D,min = 10 −3 τ s and 0.01 ≤ τ s ≤ 0.1, where A ′ = 0.15, B ′ = −0.24,C ′ = −1.48,and D ′ = 1.18.To find these coefficients, we performed a multivariate least squares fit, assuming (as we did for the ϵ crit ) that at a given τ s and α D , the critical Z lies in the middle between adjacent empty (no collapse) and filled (collapse) circles.The resulting fits are shown as red curves in each panel of Figure 5.The grey curves in each panel are random sample fits whose coefficients are randomly drawn from a multivariate normal distribution of (A ′ , B ′ , C ′ , D ′ ).As we did for ϵ crit , we used a covariance matrix of (A ′ , B ′ , C ′ , D ′ ) to produce the normal distribution.We emphasize that Equation ( 19) is valid only in the range of τ s and α D provided above and should not be extrapolated beyond the range.This is because τ s significantly affects particle dynamics; the simple form of Z crit would not be able to encompass a wider range of τ s than that we explore.In addition, we found that when α D ≲ α D,min at a given τ s , Equation ( 19) has a turning point and ends up showing Z crit increasing with decreasing α D , which is unlikely.
It is also worth pointing out that we do not include Π in Equation ( 19).This is because while Sekiya & Onishi ( 2018) demonstrate that Z/Π is the fundamental parameter combination (instead of Π and Z separately) for stratified SI, the role of Z/Π in SI clumping when external turbulence is included has not yet been explored.
One of the interesting findings of LY21 is a sharp transition in the critical Z value around τ s ∼ 0.015.Although our parameter space data are too sparsely populated to examine this feature, the α D = 10 −4 result implies that the sharp transition may disappear.As seen from the left panel of Figure 5, Z crit lies between 0.04 and 0.065 for τ s = 0.01 and between 0.04 and 0.05 for τ s = 0.02.This is a much shallower transition than was seen in LY21 in which the critical Z is ∼ 0.016 and ∼ 0.007 for τ s = 0.01 and 0.02, respectively.To investigate this more, we plot particle density that is integrated in z and averaged in y versus x and time (i.e., space-time plots) in Figure 6.The three panels correspond to different τ s values, which are 0.01, 0.02, and 0.03, for Z = 0.04 and α D = 10 −4 .We did not turn on particle self-gravity during the time span considered in the figure.
Figure 6 reveals a stochastic evolution of the filaments; weak filaments are disrupted in all three cases, and only a few strong filaments survive in τ s = 0.02 and 0.03 cases.This may be due to the external turbulence, which can contribute to unevenly distributed filaments by providing additional diffusion.On the contrary, LY21 found that for τ s = 0.01 and Z = 0.0133, filaments are so evenly spaced that they do not interact with each other, while those for τ s = 0.02 and Z = 0.01 are less uniform, with the filaments merging with each other to form a few dense ones.Since we have only looked at the the jump in Z crit for α D = 10 −4 , and even here, our data points around τ s = 0.01 are very sparse, further studies are needed to delve into this issue more.However, our results do imply that turbulence may disrupt the evenly spaced filaments that were seen in LY21, ultimately leading to a smoother transition in Z crit be- tween τ s = 0.01 and τ s = 0.02.It is also possible that the three-dimensional nature of the problem allows for this behavior, as the corresponding simulations in LY21 were all 2D.

Particle Feedback and the Particle Scale Height
Massless particles settle toward the mid-plane while competing with turbulent stirring.This competition establishes a Gaussian distribution (Dubrulle et al. 1995) of the particle density with scale height H p that is related to the strength of turbulent stirring (i.e., α D ) and τ s as in Equation ( 14).However, for particles with mass, their mass (i.e, Z) can affect the vertical stirring as well by imposing mass loading on the gas, reducing H p ; the magnitude of this effect naturally depends on Z (Yang et al. 2017, 2018, LY21, Xu & Bai 2022a).Moreover, the particle feedback can alter the vertical profile of the particle density in other ways.For example, Xu & Bai (2022a) carried out MHD simulations of the MRI in the low ionization limit and found that particle feedback enhances vertical settling by reducing the eddy correlation time and not by changing the vertical velocity of the gas.Furthermore, they found a non-Gaussian vertical particle density profile with a cusp around the mid-plane.
In order to examine the effect of the particle feedback on the particle scale height (H p ), we calculate the standard deviation of the vertical particle position (which  20) in units of the gas scale height.The dashed black line in each panel denotes the prediction for the particle scale height described in Equation ( 18).As in Figure 4, we only present the simulations where we can do the time-average within a sufficiently long pre-clumping phase.The standard deviation measured in the simulations is always lower than the prediction.21).We do not report results from simulations with a pre-clumping phase that is too short (as in Fig. 7).The figure shows that except for a few runs at τs = 0.01, this ratio is in excellent agreement with the dashed line (Equation 18) for all τs at a given αD.
for a Gaussian profile would equal H p ): where z i is the vertical position of the ith particle and ⟨z i ⟩ is the mean vertical position.Figure 7 shows timeaveraged σ p,z values from all runs in Figure 4 as a function of τ s and at each α D .The color scale denotes Z values.The circle and triangle markers show whether or not a run results in planetesimal formation via gravitational collapse (see Section 3.2 for details).The dashed line in each panel denotes the prediction for the Gaussian particle scale height when external turbulence is considered (H p,SI Turb ; Equation 18).The time-averaging is done during the pre-clumping phase to prevent the scale height from being skewed to the high density regions.
It is evident from Figure 7 that the measured σ p,z is always smaller than H p,SI Turb .Furthermore, at a given τ s and α D , larger Z corresponds to smaller σ p,z , indicating the particle feedback effect on the particle layer thickness.
In order to understand this result in greater detail, we consider the effective scale height of a dust-gas mixture (Yang & Zhu 2020) Lim et al.
is the effective sound speed of the mixture (Shi & Chiang 2013;Laibe & Price 2014;Lin & Youdin 2017;Chen & Lin 2018); although the cited papers used ρ p /ρ g to characterize particle density in Equation ( 22), we decide to use ϵ, which is the density ratio at the midplane, because our simulations are vertically stratified.
As evident from the two equations above, increasing ϵ decreases the effective sound speed ( c s ) and the effective scale height ( H).This can be interpreted as the effect of the mass loading of particles on the gas, which increases the mixture's inertia but does not contribute to thermal pressure of the gas.With this in mind, we write the particle scale height in the presence of both external turbulence and the particle feedback as H p / H = H p,SI Turb /H instead of H p = H p,SI Turb .This results in: In other words, the particle scale height is reduced by a factor of √ 1+ϵ from H p,SI Turb in the presence of particle feedback.To compare this expected particle scale height with our results, we compute time-averages of σ p,z / H = σ p,z √ 1+ϵ/H in the runs shown in Figure 7, which should be close or equal to H p,SI Turb /H (dashed lines in Figure 7) if Equation ( 23) accurately predicts σ p,z .Figure 8 shows the comparison between σ p,z / H and H p,SI Turb /H.The figure clearly shows that most of the data points are on the dashed line (i.e., H p,SI Turb /H), which demonstrates that Equation ( 23) predicts the scale height of particles from our simulations very accurately.This implies that the effect of the mass loading is likely the reason why σ p,z decreases with increasing Z at a given τ s and α D as found in Figure 7.However, runs with τ s = 0.01 and high Z values are still above the dashed lines, deviating from the prediction.We delve into these discrepancies in the next section.

Vertical Distribution of Particle Density
As we just showed, there are a few outliers at τ s = 0.01 that deviate from Equation ( 23), namely Runs T10Z6.5A4,T1Z12.5A3.5, and T1Z25A3 all with Z ≳ Z crit .In what follows, we examine vertical profiles of the particle density in these simulations to explain the discrepancy with the prediction.The profiles are calculated with particle self-gravity turned off.
In Figure 9, we show the space-time plots of the horizontally averaged particle density vs. z and time (left) and the time-averaged vertical profiles (right) for the three simulations.In all panels, we zoom in to z = ±0.4Hfrom the mid-plane.In the right panels, we present the Gaussian (dashed) and Voigt (dotted) fits to the simulation data (solid).
First, it is evident from the left panels that the particles build up a very thin layer close to the mid-plane, with an additional extended distribution vertically away from this layer.The vertical extent of both the thin layer and the extended distribution increases with increasing α D .This is not surprising as we know that stronger vertical diffusion leads to a thicker layer.
Second, the plots on the right clearly reveal non-Gaussian particle density profiles.The profiles have a cusp near the mid-plane and extended wings on the outskirts of the cusp (see also Xu & Bai 2022a as they see a similar shape to the density distribution).The Gaussian fit (dashed) in each panel matches the data (solid) only very close to the mid-plane, whereas the Voigt fit (dotted) very approximately traces the data up to larger height (i.e., ∼ ±0.1H to ∼ ±0.2H).The cusp develops due to the particle-loading on the gas that increases the inertia of the particle-gas fluid; this leads to the reduction of the vertical velocities of the gas at the mid-plane, with particles being more easily settled to form thin layers.Indeed, we found that ⟨u ′ z 2 ⟩/c s at the mid-plane is only ∼ 30% of the same quantity averaged over the entire volume in the three simulations.On the other hand, some particles are still diffused away from the mid-plane and produce the extended region described above.We find that a Voigt profile, which has a Gaussian shape near the mid-plane and Lorentzian wings above and below the mid-plane, fits the simulation data much better than a Gaussian.Nonetheless, the data do deviate from the Voigt profile at sufficiently large |z|.
The discrepancy between σ p,z and Equation ( 23) in the simulations with very large Z is simply the result of the density profile deviating significantly from a Gaussian.Since Equations ( 18) and ( 23) are based on a Gaussian profile, neither equation is appropriate for the scale height of particles that cause such strong feedback onto the gas.
We caution that a Voigt profile may not be an actual solution for the particle density distribution, and there is no physical motivation behind the fitting.Furthermore, Lyra & Kuchner (2013) analytically derived a Gaussian particle density in the presence of particle feedback, which seems to contradict the non-Gaussian profiles in Figure 9.This disagreement may stem from the fact that they assumed a constant diffusion coefficient.This is likely not true in our simulations since we find that vertical velocity of gas is reduced around the mid-plane (Xu & Bai 2022b also found that the gas velocity is reduced by particle mass-loading.).Future analytical and numerical studies are needed to examine the vertical profile of particle density in more depth.23).Particle self-gravity is off during the time span considered here.From top to bottom, αD = 10 −4 , 10 −3.5 , and 10 −3 .In the left panels, the horizontal dashed lines show the FWHM of the Voigt profiles shown in the right panels.The solid, dashed, and dotted lines in the right panels denote simulation data, Gaussian fit, and Voigt fit to the data, respectively.The title of each panel on the right shows the time interval over which each profile is averaged.Note that we only show the distributions between z = ±0.4H,while the vertical extents of the actual computational boxes are beyond this region.Particles form a thin, dense layer around the mid-plane and have a vertical density profile that deviates significantly from Gaussian.
In an attempt to quantify the characteristic width of the thin layer, we calculate the Full Width at Half Maximum (FWHM) of the Voigt profile for each run presented in Figure 9.The FWHM (f V ) is given by where f L and f G are FWHMs of Lorentzian and Gaussian profiles, respectively.The values of the f V for the three runs from top to bottom are ∼ 0.019H, ∼ 0.032H, and ∼ 0.097H, respectively.We mark the FWHMs in the left panels as two horizontal dashed lines for each run; the FWHMs constrain the thickness of the thin particle layers very well.
To demonstrate how the profiles change with larger τ s values, we present plots similar to those in Figure 9, but for τ s = 0.1, Z = 0.02 in Figure 10.Here, we zoom in  further and present the distributions within z = ±0.2H,unlike than in Figure 9.
These runs show negligible discrepancy between σ p,z and Equation ( 23) in Figure 8. Furthermore, as expected, the profiles are significantly changed compared to the τ s = 0.01 profiles.First, the vertical extent of the distributions is much narrower as a result of stronger settling of the τ s = 0.1 particles compared with the τ s = 0.01 particles.Second, thin particle layers that are separated from a more extended particle region do not form, and the time-averaged vertical profiles (right panels) do not show cusps at the mid-plane.This results from the Z values being too small for particles to add considerable mass-loading on the gas.Third, as a result of this lack of the significant mass-loading, the density profiles are well approximated by a Gaussian but with reduced width according to Equation ( 23).Overall, because of the absence of a cusp and Gaussian-like density profiles, σ p,z and Equation ( 23) well matches each other than the simulations shown in Figure 9.
In summary, the results in this subsection indicate two related considerations.First, due to mass-loading on the gas, the vertical particle density profile may still be Gaussian, but with a reduced width compared to the massless particle case.In this case, the criterion (Equation 17) overestimates critical Z values.However, for even higher Z values, the mass loading becomes so significant that the vertical particle density profile significantly deviates from Gaussian.This deviation is particularly noticeable in simulations with τ s = 0.01, where Z reaches such high levels that particles form a thin, dense layer at the mid-plane.Since an analytical expression for such a profile does not yet exist, we argue that caution is warranted when characterizing the width of particle layers under the influence of turbulence at large Z.However, for more moderate Z values Equation ( 23) can be used to estimate particle scale height in the presence of both external turbulence and particle feedback.

DISCUSSION
This section is dedicated to providing a better understanding of the robustness and impact of our results.In Section 4.1, we consider the influence of turbulence on planetesimal formation, drawing comparisons with previous studies.We then discuss the potential implications of our findings on observations of protoplanetary disks in Section 4.2.Finally, we highlight potential limitations and uncertainties in our work in Section 4.3.

Does Turbulence Hinder or Help Planetesimal
Formation?
While our results suggest that turbulence acts as a hindrance to planetesimal formation, there are a number of other possible routes to forming planetesimals in turbulent disks.In this subsection, we discuss these other routes and the connection with this work.

Self-consistently driven turbulence
The results we present here stand in contrast to other work in which turbulence is included self-consistently (as opposed to forcing isotropic, incompressible turbulence as we do here) and gives rise to structures and behaviors that can concentrate particles.For example, the MRI is known to produce localized reductions in the radial pressure gradient (generally referred to as "zonal flows"; see e.g., Simon & Armitage 2014).As radial drift slows in these regions (and can become trapped if the local pressure profile has a maximum), these zonal flows serve as natural sites for particle concentration and thus planetesimal formation as shown in, e.g., Johansen et al. (2007) and Xu & Bai (2022a).Similarly, Schäfer & Johansen (2022) demonstrate that the VSI can also produce localized changes to the pressure gradient that then triggers the SI.6 On the contrary, there is no evidence of significant pressure variations to trap particles in our simulations, meaning that particle concentration is driven by the SI in this work.
At face value, these works suggest that turbulence can act to help planetesimal formation.More concretely, in the work by Xu & Bai (2022a), the maximum density of particles surpassed ρ H for τ s = 0.1, Z = 0.02, even with α D ≳ 10 −3 .For comparison, our Run T10Z2A3 (τ s = 0.1, Z = 0.02, α D = 10 −3 ) shows no indication of significant particle clumping at all (see the rightmost column of Figure 2).
However, the potential discrepancy with our work can be elucidated by a more in-depth examination of Z.In Xu & Bai (2022a), the quoted Z = 0.02 value is a "global" value (i.e., it was the average over the entire domain), and toward the pressure maximum (i.e., where Π = 0), Z is enhanced.While we cannot quantify the precise value of Z at this location without more data, these considerations are in qualitative agreement with our results: the local value of Z at the pressure maxima is (very possibly significantly) higher than the background value of Z = 0.02.
Such a direct comparison with Xu & Bai (2022a) should be treated with caution, however, as the particle concentrations in their work occur at the pressure maxima (i.e., where Π = 0); the SI does not operate in such regions.In fact, that they did not see particle concentration in regions outside of the pressure bump, where Z is locally smaller and where Π ̸ = 0 aligns with our findings that planetesimal formation is hindered when τ s = 0.1, α D = 10 −3 , and Z ≲ 0.04.
The key point here is that it is not the turbulence itself that is aiding in planetesimal formation (though, see arguments about the Turbulent Concentration mechanism below) but rather localized increases in Z due to enhancements in the gas pressure.While this argument may appear to be one of semantics, it is important to distinguish between long-lived coherent structures, such as zonal flows, and random turbulent fluctuations, such as eddies, that are very short-lived by comparison.Thus, while these pressure bumps are a side-effect of the turbulence, they are arguably distinct from the turbulence itself.Our simulations provide a way to interpret these simulations by means of using Z as a control parameter rather than one that changes with location based on the dynamics at work.That all being said, Yang et al. (2018) demonstrate that particle concentration in a (at least somewhat) turbulent disk environment is possible even in the absence of such pressure bumps.In particular, they find modest to strong particle concentration in an Ohmic dead zone (see e.g., Gammie 1996 for a description of the dead zone model) that results from an anisotropy in turbulent diffusion.That is, there was no pressure bump to trigger the SI but rather anisotropic turbulence.A direct comparison with these results is difficult since our turbulence is forced isotropically, and thus we leave a study of the effect of anisotropic turbulence to future work.

Lim et al.
Overall, our simulations should be treated as more controlled experiments; i.e., Z is an input parameter and not something that arises from the simulation as a result of particle concentration.Furthermore, with the exception of Schäfer & Johansen (2022), all of the works involving turbulence driven by (magneto-)hydrodynamical instabilities in PPDs use τ s = 0.1, whereas our higher resolution simulations (i.e., more grid zones per H) allow us to resolve the SI for τ s = 0.01.7 Thus, our work serves as an important counterpart to the numerous papers that include turbulence driven self-consistently and that (in most cases) are limited to larger τ s by resolution requirements.

Turbulent Concentration
Beyond the processes we just described, planetesimal formation may be aided by turbulence through the turbulent concentration mechanism (hereafter, TC; Cuzzi et al. 2001Cuzzi et al. , 2008;;Hartlep et al. 2017;Hartlep & Cuzzi 2020).Specifically, particles with certain τ s are preferentially concentrated by eddies whose eddy turnover time is comparable to τ s .In other words, there is an optimal τ s value that makes the turbulent Stokes number at scale ℓ St ℓ ≡ τ s /τ eddy,ℓ ∼ 1, where τ eddy,ℓ , is the dimensionless eddy turnover time at scale ℓ (τ eddy,ℓ ≡ t eddy,ℓ Ω, and t eddy,ℓ is the dimensional turnover time).While we did not address the TC in this study, we now investigate whether TC might be present in our simulations by doing a simple scaling calculation.
Assuming the forced turbulence follows the Kolmogorov relations, u ℓ ∝ ℓ 1/3 and ℓ ∼ u ℓ t eddy,ℓ ∝ t 3/2 eddy,ℓ , where u ℓ and t eddy,ℓ are a characteristic velocity and a turnover time of an eddy at scale ℓ, respectively.Next, the forcing is done between ∼ 0.1H and ∼ 0.2H (see Appendix A for details), and we set 0.2H as the outer scale (L) of the forced turbulence.The eddy turnover time for the outer scale is obtained from τ eddy,L ∼ α D /(δu/c s ) 2 (Equation 13), which is ∼ 0.51, ∼ 0.44, and ∼ 0.37 for α D = 10 −4 , 10 −3.5 , and 10 −3 , respectively; (δu/c s ) 2 is computed without particles (see Table 2).For the purposes of this analysis, we set τ eddy,L = 0.44.We thus obtain ℓ at which St ℓ ∼ 1.0 (i.e., τ s = τ eddy,ℓ ) where TC becomes efficient: Particles with τ s = 0.01 would be concentrated by eddies at scale ℓ ∼ 7×10 −4 H, which is well below the width of a grid cell ∆ = H/640 ∼ 0.0016H.The τ s = 0.03 and τ s = 0.1 particles on the other hand would be concentrated at ℓ ∼ 0.004H (equating to ℓ/∆ ∼ 2) and ℓ ∼ 0.02H (ℓ/∆ ∼ 14), respectively.These scales are above the grid scale (though in the case of τ s = 0.03, this is only marginally true).However, a very approximate estimate for the dissipation scale (inferred from kinetic energy power spectra, see Appendix A) in our simulations is 0.01H.The τ s = 0.03 particles would be concentrated on a scale less than this, whereas the τ s = 0.1 particles would be concentrated on a scale only twice that of the dissipation scale.
It is worth noting that while the most recent results suggest that St ℓ ≈ 0.3 is the most optimal value (Hartlep et al. 2017;Hartlep & Cuzzi 2020), making it easier to resolve TC, there remains enough uncertainty in both this value and our approximate analysis that the absence of the TC should not be weighed too heavily.A much deeper dive into the TC as a possible mechanism for planetesimal formation is certainly required, but is beyond the scope of this paper.
In any case, our work clearly demonstrates that over a large parameter space of τ s , Z, and α D values, turbulence can significantly weaken the SI and prevent planetesimals from forming.This happens through turbulent diffusion counteracting vertical settling and/or the radial concentration of filaments within the disk plane.

Implication for Observations
Several recent observations have quantified the scale height of millimeter-emitting dust in Class II disks, such as those surrounding HL Tau (Pinte et al. 2016) or Oph 163131 (Villenave et al. 2022).The findings from these studies generally suggest that the dust particles are well settled in the outer regions of the disks, with their scale heights being less than 1 au at a radial distance of 100 au.This behavior could imply the presence of very weak turbulence in the outer regions as indicated by the models used in the cited observations.However, those models neglected particle feedback.Thus, it is valuable to examine the potential implications for these observaiotns of particle feedback and the damping of turbulence.Villenave et al. (2022) use a radiative transfer model to estimate the scale height of dust that radiates 1.3 mm continuum emission in the disk around Oph 163131.The resulting scale height is ∼ 0.5 au at 100 au, while the scale height of the gas is estimated to be 9.7 ± 3.5 au from scattered-light data (Wolff et al. 2021).If we use the mean value for the gas scale height, this equates to H p /H ∼ 0.05.Assuming that the millimeter-emitting dust in their observation has τ s = 0.01 at 100 au and ρ p ≪ ρ g , α D,z ∼ 2.5×10 −5 (Equation 14).This α D,z value is roughly consistent with the upper limit of ∼ 10 −5 calculated in Villenave et al. (2022) by adopting the dust settling model of Fromang & Nelson (2009).
However, as we discussed in Sections 3.3 and 3.4, the mass of the particle can increase the inertia of the dustgas mixture, resulting in a very thin particle layer even in the presence of stronger turbulence.For example, the dust thickness in our simulations, σ p,z /H, spans from ∼ 0.04 to ∼ 0.1 for τ s = 0.01 depending on the values of α D,z and Z (see Figure 7).The observed dust scale height in Oph 163131 equates to H p /H ranging from ∼ 0.04 to ∼ 0.08 if we account for the uncertainty in the gas scale estimation (i.e., H = 9.7±3.5 au), which is consistent with our results but only for α D,z ≫ 10 −5 .
To summarize, well-settled dust disks (i.e., small H p ) do not necessarily have very weak turbulence, especially if the dust-to-gas ratio is ≫ 1.Thus, the level of turbulence inferred by ignoring feedback should be regarded as a lower limit.

Strength of Particle Self-Gravity
Throughout this paper, we fix G (Equation 9) to 0.05 for all simulations we perform in order to compare our results to previous studies that assume (LY21) or use the same value (Gole et al. 2020).Moreover, our choice of the G value should be viewed as a conservative way to establish the collapse threshold since planetesimal formation could be triggered in even weaker concentrations when a higher G (or equivalently, a lower Q) is used (see Equation 10or Gerbig et al. 2020).In other words, critical Z values will be lower in young massive disks or the outer regions of a Class II disk (e.g., G ∼ 0.2 at r = 45 au based on the disk model in Carrera et al. (2021)).Therefore, our Z crit in Section 3.2 is subject to change when different values of G are used.

Number of Particles
Every simulation considered in this work has the same average number of particles per grid cell (i.e., n p = 1).However, the effective particle resolution, which is the number of particles per 2σ p,z (i.e., n p,eff = n p [L z /2σ p,z ]), varies with α D and the other parameters.For example, Run T10Z2A4 has σ p,z ∼ 0.013 and n p,eff ∼ 30, whereas Run T10Z2A3 has σ p,z ∼ 0.05 and n p,eff ∼ 8.In order to test how changing the effective particle resolution changes our results, we run Run T10Z2A4 again but with the same effective particle resolution of Run T10Z2A3, named T10Z2A4-n p .The Run T10Z2A4-n p has N par ∼ 4.61×10 6 , which is approximately a factor of 4 smaller than that of Run T10Z2A4.
We compare the two models in Figure 11, which presents the maximum density of particles as a function of time in the left panel and the radial concentration over time in the middle and the right panels.In the left panel, Runs T10Z2A4 and T10Z2A4-n p are denoted as black and grey, respectively, and the Hill density at G = 0.05 is shown as the orange horizontal line.The particle self-gravity is turned off in the data shown here.From the maximum density plots, the two runs exhibit almost identical evolution until t−t par ∼ 200Ω −1 after which they diverge.The time-averaged maximum densities from t−t par = 200Ω −1 to 500Ω −1 are ∼ 138ρ g0 and ∼ 160ρ g0 for Runs T10Z2A4 and T10Z2A4-n p , respectively.Given that the maximum density is inherently stochastic, the discrepancy does not seem to be significant.The radial concentration of particles of the two runs (middle and right panels) looks very similar as well.Although Run T10Z2A4-n p has one more filament at the end of the simulation, we do not believe this is a significant difference because the interaction between filaments is highly nonlinear, with the final number of dominant filaments and the maximum density values be-  26) as a function of δ for τs = 0.01 (blue), 0.03 (green), and 0.1 (red).The yaxes on the left and on the right show rcrit/∆ and rcrit/H, respectively, where ∆ is the width of a grid cell.We assume that δ ≈ αD,z, which ranges from 10 −6 to ∼ 4×10 −4 .The larger number in this range corresponds to the value for αD = 10 −3 , while the smaller number is an assumed value based on the fact that particles become less diffusive around dense clumps.The black horizontal line denotes ∆.Except for the case where δ is very low, the critical radius is larger than the cell size.
ing uncertain to some degree.Overall, we conclude that the effective particle resolutions we employ is not likely to significantly affect our results.

Numerical Resolution
Other than G and n p , we also fix the grid resolution to 640/H in all simulations.However, it is possible that increasing this resolution would affect the ϵ crit and Z crit curves.First considering 2D simulations, Yang et al. (2017) found that at τ s = 0.01, the critical Z above which the SI produces filaments lies between 0.02 and 0.04 at 640/H resolution, whereas it lies below 0.02 at ≳ 1280/H resolution; these results suggest that increased resolution lowers the critical value.On the other hand, LY21 found that the critical Z value above which strong clumping occurs (again, their definition of strong clumping requires the maximum particle density to exceed the Hill density) was 0.0133 (at τ s = 0.01) for all resolutions up to 5120/H (though smaller boxes were used for higher resolution simulations).In their 3D simulations, Yang et al. (2017) found that for τ s = 10 −3 and Z = 0.04 at 160/H resolution, there was no sign of significant particle concentration or filament formation and ρ p,max < 10ρ g0 (i.e., a critical Z ≳ 0.04).However, at 320/H and 640/H resolutions, they observed a significant increase in concentration (ρ p,max ∼ 30ρ g0 and ∼ 200ρ g0 in the lower and the higher resolutions, respectively) and the formation of dense, persistent fil-aments (i.e., a critical Z ≲ 0.04).The differences in whether a resolution dependence for the critical Z was observed may be a result of the different codes used (including other factors, such as boundary conditions).However, taken together, these results suggest that while resolution may affect the critical Z values, these values will likely only be changed by a factor of order unity.Therefore, while more investigation is certainly needed to demonstrate if increasing grid resolution indeed affects ϵ crit and Z crit , we expect that it would produce at most a modest change in these critical values.
In addition to SI concentration, whether gravitational collapse of particles occurs depends on the numerical resolution.More specifically, numerical simulations should resolve the critical length r crit derived in Klahr & Schreiber (2020) at which gravitational contraction balances internal (to a collapsing cloud of pebbles) turbulent diffusion to accurately capture the collapse and planetesimal formation.
Here, δ is a dimensionless parameter for internal diffusion within a particle cloud.Particle clouds whose sizes are greater than r crit are too massive to be held up by diffusion, and thus, they gravitationally collapse.Here, we examine whether or not the scale defined by r crit is resolved in our simulations.Previous works measure the radial diffusion of particles, assuming they undergo a random walk in the radial direction, to quantify δ (see, e.g., Baehr et al. 2022;Gerbig & Li 2023).Instead of directly measuring the radial diffusion in our simulations, we let δ be a free parameter ranging from 10 −6 to 4×10 −4 to cover δ values that likely result from the turbulence in our simulations.The lower limit is to account for the fact that denser regions are less diffusive (Gerbig & Li 2023), and the upper limit corresponds to α D,z for α D = 10 −3 .
Figure 12 shows r crit as a function of δ for three selected τ s values, which are 0.01 (blue), 0.03 (green), and 0.1 (red).The first (left) and the second (right) vertical axes show r crit in the unit of the size of a grid cell (∆) and in the unit of H, respectively.We denote ∆ as the black horizontal line.As can be seen, r crit is generally larger than the grid cell size, meaning that our simulations should largely resolve gravitational collapse if a local particle clump is gravitationally unstable.For τ s = 0.1 (red) and very small δ, the critical length becomes smaller than the cell size.Thus, it is possible that Runs T10Z1.5A4and/or T10Z2A3.5 would produce collapsed regions if the numerical resolution was higher.However, it is less likely that an increased reso-lution would change the results of the α D = 10 −3 runs because the radial diffusion in these runs is much larger than 10 −6 .
Overall, the ϵ crit and Z crit we report should be viewed as upper limits since both SI concentration and gravitational collapse typically benefit from higher resolutions.Although our findings may be adjusted based on the caveats highlighted in this subsection, it is crucial to note that our simulations are 3D and incorporate essential physics, such as self-gravity of particles and external turbulence.Thus, our simulations represent a significant advancement in our quantification of the critical planetesimal formation curves.

SUMMARY
In this paper we have presented results from stratified shearing box simulations in which gas and particles are aerodynamically coupled to each other.In order to study the effect of turbulence on the streaming instability (SI) and the formation of planetesimals, we include both the self-gravity of particles and externally driven incompressible turbulence.Our simulations explore a relatively broad range of parameter space, namely different dimensionless stopping times τ s , particle-to-gas surface density ratios Z, and forcing amplitudes α D .We summarize our main results as follows: 1. Incompressible turbulence can impede SI-driven concentration of particles via turbulent diffusion in two possible ways.First, this diffusion can prevent particles from settling, thereby preventing the mid-plane dust-to-gas density ratio ϵ from exceeding the critical value for filament formation.Second, even if the particles do settle and form a layer around mid-plane, the formation of filaments can be counterbalanced by turbulent diffusion acting in the plane of the disk (in addition to vertically).
2. The critical ϵ, at or above which planetesimal formation occurs, is ϵ ≳ 1.This is a factor of a few larger than the corresponding values in the absence of externally driven turbulence, but is still of order unity.
3. To balance the stronger diffusion associated with larger α D , more total mass in the particles (i.e., through the Z parameter) is needed.As such, the critical Z values (Z crit ) at or above which planetesimal formation occurs are much higher than those obtained without external turbulence.Quantitatively, when τ s = 0.01, Z crit ∼ 0.06 and ∼ 0.2 for α D = 10 −4 and 10 −3 , respectively, whereas Z ∼ 0.02 is sufficient for planetesimal formation in the absence of turbulence (e.g., LY21).4. Due to particle feedback, the characteristic particle height in our simulations (σ p,z ) is always lower than the particle scale height with negligible feedback.This behavior is the direct result of enhanced particle mass loading on the gas.
5. As a result of the strong influence of particle feedback on the dust scale height, observational measurements of the turbulent velocity should be regarded as a lower limit.It is possible to have stronger turbulence and a small dust scale height if the dust-to-gas ratio is sufficiently large.
6.For sufficiently large Z, the vertical particle density profiles can be significantly modified from a Gaussian.For τ s = 0.01 and Z ≥ Z crit , our simulations exhibit a cusp near the mid-plane resulting in a thin, dense layer of particles, and extended wings outside the layer, resembling a Voigt profile out to |z| ≲ 0.2H.
In closing, while there remain a number of uncertainties to be addressed in future work, our results demonstrate the crucial role that gas turbulence plays in limiting where in the disk and under what conditions planetesimals can form.).The power spectra have a peak at kH/(2π) ∼ 10 followed by a power-law toward larger wavenumbers, which may be an indication of energy cascade.
Finally, at even larger k, the spectra become steeper due to numerical dissipation.in our simulations.The color-coding denotes α D values.Additionally, the black dashed line shows a -5/3 power-law relation, which is a characteristic of Kolmogorov turbulence.The spectra are truncated at kH/(2π) = 320, which corresponds to the Nyquist scale (2∆).We summarize key observations from the figure as follows: 1. Simulations with larger α D consistently show larger power across all considered scales.
2. The power spectra exhibit common features of turbulence in numerical simulations.That is to say, they peak at kH/(2π) ∼ 10, which is a scale of ∼ 0.1H.At k larger than where the peak occurs, there is an indication of a (small) inertial range where a cascade of energy is likely taking place.At even larger wavenumbers (or, smaller scales), the power-law is broken, and the slope becomes steeper due to numerical dissipation being dominant over inertia of the turbulence.
3. Finally, the shape of the power spectra's (small) inertial range is in rough agreement with a -5/3 power law, which suggests a Kolmogorov type behavior.However, we emphasize that our simulations, which are both rotating and stratified, will not strictly follow Kolmogorov phenomenology, and thus we cannot firmly conclude that the turbulence in our simulations is behaving as in the classical Kolmogorov picture.

Figure 1 .
Figure1.Relationship between Λ and αD in simulations with τs = 0.1, Z = 10 −5 .The former is an initial condition for a forcing amplitude (Equation12), while the latter is the parameter (Equation13) we use to denote the level of turbulence in our SI simulations.The Λ values shown in this plot are obtained by a linear interpolation (see the text for details).Each red horizontal line indicates each αD we choose: 10 −4 , 10 −3.5 , and 10 −3 .The error bar shows ±1 standard deviation.We confirm that the interpolated Λ values result in the desired αD values with high accuracy.

Figure 2 .
Figure 2. Top: Time evolution of gas rms velocity (black) and maximum density of particles (orange) in simulations with τs = 0.1, Z = 0.02.Middle: Azimuthally averaged particle density (side view) zoomed in to the region z/H ∈ [−0.1, 0.1].Bottom: Vertically integrated particle mass density (i.e., surface density of particles; Σp; top view).The snapshots are taken at t−tpar = 350Ω −1 .Note that the middle and the bottom panels have different colorbar scales.From left to right, no-forcing, αD = 10 −4 , 10 −3.5 , 10 −3 runs are shown.The horizontal line in each panel on the top marks Hill density (see Equation10).Particle self-gravity is disabled in the runs presented here.As turbulence becomes stronger from left to right, the particle layer become thicker, and the filaments become more diffused and less well-defined.

Figure 3 .
Figure 3. Proof of concept for categorizing runs as "Collapse" or "No Collapse".Left: Time evolution of maximum particle density for T10Z2A35 (green) and T10Z3.2A35(blue) runs.The horizontal and vertical lines indicate Hill density and tsg, respectively.Middle: Final snapshot of Σp divided by its spatial average for T10Z2A3.5 run.Right: Same as the middle panel but for T10Z3.2A3.5 run.We define a run as "Collapse" if the maximum density drastically increases as shown by the blue curve in the left panel and if the 2D snapshot shows bound objects (i.e., strong overdensities on very small scales as shown in the right panel).

Figure 4 .
Figure4.Temporally and horizontally averaged ratio of particle density to gas density at the mid-plane as a function of τs in our simulations.We do not include simulations where the pre-clumping phase is too short to take time-averages.Runs where the collapse occurs are shown as circles, whereas those where the collapse does not occur are shown as triangles.Red, green, and blue represent αD = 10 −4 , 10 −3.5 , and 10 −3 , respectively, while we did not denote Z values.The black curve denotes the least squares best fit to the critical value of ϵ as described in Equation (16).The thin grey curves are random fits drawn from multivariate normal distribution of A, B, C (Equation16) that show uncertainties in the best fit curves.We also include the critical curve from LY21, depicted as a sky blue curve (ϵcrit,LY21).The figure shows that ϵcrit does not vary with αD (i.e., even changing αD, every no-collapse case is below every collapse case), which results in a single critical curve that fits all ϵ values regardless of αD.

Figure 6 .
Figure6.Space-time plots of the particle density that is integrated over z (Σp) and averaged in y (⟨Σp⟩y).We plot the quantity vs. x and time.From top to bottom, τs = 0.01, 0.02, and 0.03, respectively.The runs have the same Z and αD, which are 0.04 and 10 −4 , respectively.Particle self-gravity is turned off during the time span considered here.The filaments are frequently disrupted unless they are sufficiently dense, resulting in stochastic evolution for all three τs values.

Figure 7 .
Figure 7. Similar to Figure 4 but for the time-averaged standard deviation of particles' vertical positions (Equation20) in units of the gas scale height.The dashed black line in each panel denotes the prediction for the particle scale height described in Equation (18).As in Figure4, we only present the simulations where we can do the time-average within a sufficiently long pre-clumping phase.The standard deviation measured in the simulations is always lower than the prediction.

Figure 8 .
Figure 8. Similar to Figure 7 but for the time-averaged ratio of σp,z to H (Equation21).We do not report results from simulations with a pre-clumping phase that is too short (as in Fig.7).The figure shows that except for a few runs at τs = 0.01, this ratio is in excellent agreement with the dashed line (Equation18) for all τs at a given αD.

Figure 9 .
Figure9.The space-time plots of particle density averaged in x and y (⟨ρp⟩xy) vs. z and time (left) and time-averaged vertical profiles of the particle density (right) for three runs at τs = 0.01 that show σp,z deviating from Equation (23).Particle self-gravity is off during the time span considered here.From top to bottom, αD = 10 −4 , 10 −3.5 , and 10 −3 .In the left panels, the horizontal dashed lines show the FWHM of the Voigt profiles shown in the right panels.The solid, dashed, and dotted lines in the right panels denote simulation data, Gaussian fit, and Voigt fit to the data, respectively.The title of each panel on the right shows the time interval over which each profile is averaged.Note that we only show the distributions between z = ±0.4H,while the vertical extents of the actual computational boxes are beyond this region.Particles form a thin, dense layer around the mid-plane and have a vertical density profile that deviates significantly from Gaussian.

Figure 10 .
Figure 10.Similar to Figure 9 but for τs = 0.1, Z = 0.02.Note that we show the distributions between z = ±0.2Hand do not show Voigt fits in the right panels.. Particle self-gravity is off during the time span presented here.Unlike τs = 0.01 cases, a thin, dense layer doesn't form, and the vertical profiles are closer to Gaussian.

Figure 11 .
Figure11.Comparison of two simulations with all parameters the same apart from the effective particle resolution (n p,eff ).The left panel shows the time evolution of the maximum particle density, where the black and the grey curves are for the fiducial run (Run T10Z2A4, n p,eff ∼ 8) and the run with the effective resolution of Run T10Z2A3 (Run T10Z2A4-np, n p,eff ∼ 8), respectively.The orange horizontal line is the Hill density (ρH ).The middle and the right panels present the space-time plots of particle density averaged in y and integrated over z vs. x and time for the two runs.Particle self-gravity is disabled here.Despite slight differences in the maximum density and the radial concentration, the two runs with different n p,eff are very similar.

Figure 12 .
Figure12.The critical radius of a particle cloud above which collapse occurs (rcrit; Equation26) as a function of δ for τs = 0.01 (blue), 0.03 (green), and 0.1 (red).The yaxes on the left and on the right show rcrit/∆ and rcrit/H, respectively, where ∆ is the width of a grid cell.We assume that δ ≈ αD,z, which ranges from 10 −6 to ∼ 4×10 −4 .The larger number in this range corresponds to the value for αD = 10 −3 , while the smaller number is an assumed value based on the fact that particles become less diffusive around dense clumps.The black horizontal line denotes ∆.Except for the case where δ is very low, the critical radius is larger than the cell size.

Figure A. 1 .
Figure A.1.Time-averaged power spectra of turbulent velocity squared before initializing particles in simulations with (Lx, Ly, Lz) = (0.4,0.2, 0.8)H.Three αD values are denoted by different colors as shown in the legend.The black dashed line represents Kolmogorov dependence (∝ k −5/3).The power spectra have a peak at kH/(2π) ∼ 10 followed by a power-law toward larger wavenumbers, which may be an indication of energy cascade.

Table 1 .
List of simulations and time-averaged quantities Figure5.Overview of gravitational collapse in the SI simulations listed in Table1except Run T10Z2 in which no forcing is applied.Runs where the collapse occurs are shown as filled circles, whereas those where the collapse does not occur are shown as open circles.From left to right, αD = 10 −4 , 10 −3.5 , and 10 −3 , respectively.In each panel, we plot two different Zcrit, skyblue being Equation (17) and red being Equation (19).Both of the curves are fits to the critical Z values but by using different approaches; the red curve is the best fit by multivariate least squares (Equation