Concentrating small particles in protoplanetary disks through the streaming instability
Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, 221 00 Lund, Sweden
email: ccyang@astro.lu.se
Received: 21 November 2016
Accepted: 13 June 2017
Laboratory experiments indicate that direct growth of silicate grains via mutual collisions can only produce particles up to roughly millimeters in size. On the other hand, recent simulations of the streaming instability have shown that mm/cmsized particles require an excessively high metallicity for dense filaments to emerge. Using a numerical algorithm for stiff mutual drag force, we perform simulations of small particles with significantly higher resolutions and longer simulation times than in previous investigations. We find that particles of dimensionless stopping time τ_{s} = 10^{2} and 10^{3} – representing cm and mmsized particles interior of the water ice line – concentrate themselves via the streaming instability at a solid abundance of a few percent. We thus revise a previously published critical solid abundance curve for the regime of τ_{s} ≪ 1. The solid density in the concentrated regions reaches values higher than the Roche density, indicating that direct collapse of particles down to mm sizes into planetesimals is possible. Our results hence bridge the gap in particle size between direct dust growth limited by bouncing and the streaming instability.
Key words: hydrodynamics / methods: numerical / minor planets, asteroids: general / planets and satellites: formation / protoplanetary disks / instabilities
© ESO, 2017
1. Introduction
In the coreaccretion scenario of planet formation, planetary cores are assembled beginning with interstellar μmsized dust grains (Safronov 1969). This process of growing planetary cores covers more than 30 orders of magnitude in mass and more than 13 orders of magnitude in size, required to be completed within the 1–5 Myr lifetime of their natal protoplanetary disks (e.g., Haisch et al. 2001; Mamajek 2009). The course of planet formation is usually divided into different stages according to the size of the solids or planetary objects involved, and each stage has its own major difficulties (e.g., Papaloizou & Terquem 2006; Baruteau et al. 2014; Testi et al. 2014). One of these stages yet to be understood is the formation of kilometerscale planetesimals (we refer to e.g., Johansen et al. 2014, and references therein).
A major obstacle to the formation of planetesimals is the “radialdrift barrier”. Solid particles marginally coupled to the gas via drag force drift radially inwards and are quickly removed from protoplanetary disks due to the gaseous head wind (Adachi et al. 1976; Weidenschilling 1977a). For example, the timescale for the radial drift of metersized boulders at ~1 au of the minimummass solar nebula (MMSN; Weidenschilling 1977b; Hayashi 1981) is ~100 yr, significantly shorter than the typical lifetime of protoplanetary disks. In fact, solid particles of a wide range of sizes at various locations in disks suffer from radial drift (e.g., Brauer et al. 2007; Youdin 2010). One possibility to circumvent this barrier is to increase the collisional crosssection of the particles by collecting porous icy dust aggregates (Kataoka et al. 2013). This process, however, may only operate outside the ice line and relies on the presence of large amounts of submicron monomers. Therefore, some mechanism(s) to efficiently concentrate solid materials into high density so that direct gravitational collapse can occur appears to be required to form planetesimals.
Two types of mechanism that concentrate solid materials exist; one passive, one active. The former includes longlived local pressure maxima (Whipple 1972; Johansen et al. 2009a; Bai & Stone 2014) or vortices (Barge & Sommeria 1995; Bracco et al. 1999; Klahr & Hubbard 2014; Lyra 2014), in which solids passively follow the underlying flow of gas by friction. By contrast, the streaming instability discovered by Youdin & Goodman (2005) is realized by the actionreaction pair of the drag force between solid particles and gas, with which the solids actively engage in the dynamics with the gas to spontaneously concentrate themselves into high densities (Johansen & Youdin 2007; Bai & Stone 2010a; Yang & Johansen 2014). As a result, not only does the dense filamentary structure of solids driven by the streaming instability have significantly reduced radial drift speeds and thus longer radial drift timescales, but also planetesimals can form by direct gravitational collapse within these dense filaments (Johansen et al. 2012; Johansen et al. 2015; Simon et al. 2016, 2017; Schäfer et al. 2017).
Nevertheless, it remains problematic for the streaming instability to drive planetesimal formation inside the ice line of protoplanetary disks (Dra¸żkowska & Dullemond 2014). Direct dust growth by coagulation of compact^{1} silicate grains is limited to mm sizes, due to bouncing and fragmentation at collisions (Zsom et al. 2010; Birnstiel et al. 2012). On the other hand, there exists a threshold in solid abundance only above which can the streaming instability drive solid concentration into high densities (Johansen et al. 2009b; Bai & Stone 2010a). Using a suite of twodimensional simulations, Carrera et al. (2015) found that this threshold abundance increases drastically with decreasing particle size. For mmsized particles inside the ice line of protoplanetary disks, they suggested that a solidtogas column density ratio of more than 10% is required, an exceedingly high value in usual disk conditions. Therefore, it seems that a significant gap between dust coagulation and the onset of planetesimal formation exists inside the ice line of young protoplanetary disks.
A few properties of the streaming instability are worth noticing, though, which may explain the exceedingly high critical solid abundance found for small particles. First, even though the growth rate of the linear streaming instability for small particles is relatively small compared to that for large particles when the local solidtogas density ratio is low, it significantly overtakes that for large particles when the density ratio is more than the order of unity (Youdin & Johansen 2007). Hence, the response of small particles to the streaming instability should be more prominent when the latter condition is reached. However, the wavelengths of the fastestgrowing modes are roughly proportional to the dimensionless stopping time τ_{s} ≡ Ω_{K}t_{s} and thus the size of the particles, where Ω_{K} is the local Keplerian frequency and t_{s} is the stopping time characterizing the mutual drag force between the gas and the particles (Youdin & Johansen 2007). These wavelengths are as small as ~10^{4}H for particles with τ_{s} ~ 10^{3}, where H is the local scale height of the gas (Bai & Stone 2010b). Resolving these faster growing modes of the streaming instability remains challenging in current numerical simulations. Although it remains unclear how these linear modes cooperate to determine the conditions for strong concentration of solid particles in the nonlinear saturation of the streaming instability (Youdin & Johansen 2007; Jacquet et al. 2011), systematic resolution studies considering higher resolutions seem warranted. Finally, in the nonlinear stage of the streaming instability, the timescales for sedimentation and radial drift of the particles should continue to be relevant for the system. These timescales are inversely proportional to τ_{s} when τ_{s} ≪ 1, and hence a proportionately longer time may be required for small particles to concentrate themselves than for their large counterparts.
In spite of that, due to stringent timestep constraint of the stiff mutual drag force, it is numerically challenging to simulate the particlegas system in question with a higher resolution and longer simulation time than was achieved by Carrera et al. (2015). To make such simulations feasible, we have devised a new algorithm in Yang & Johansen (2016) to relieve this timestep constraint, and have demonstrated that the algorithm manages small particles and/or strong local solid concentration with satisfactory numerical convergence and accuracy. Employing this technique in this work, we revisit Carrera et al. (2015) for the case of dimensionless stopping time τ_{s} ≪ 1 with significantly higher resolutions and longer simulation times.
We indeed find these small particles can spontaneously concentrate themselves into high density at a much reduced solid abundance: ~1–2% for particles of τ_{s} = 10^{2} while ~3–4% for particles of τ_{s} = 10^{3}. In the following section, we briefly describe our methodology. We then present the evolution of the system, focusing on the concentration of the solid particles, in both twodimensional (2D; Sect. 3) and threedimensional (3D; Sect. 4) models. The implications for the formation of planetesimals are discussed in Sect. 5, where we revise the critical solid abundance of Carrera et al. (2015) with the results of this work. We conclude with a short summary in Sect. 6.
2. Methodology
As in our previous publication (Yang & Johansen 2014), we use the localshearingbox approximation (Goldreich & LyndenBell 1965) to simulate a system of gas and solid particles. The computational domain is a small rectangular box at an arbitrary distance from the central star. Its origin rotates around the star at the local Keplerian angular speed Ω_{K} and its three axes constantly align in the radial, azimuthal, and vertical directions, respectively. We consider isothermal, nonmagnetized gas in a regular Eulerian grid with numerous Lagrangian solid particles in it. The gas interacts with each of the particles via their mutual drag force, which is characterized by the stopping time t_{s} (Whipple 1972; Weidenschilling 1977a) or its dimensionless counterpart τ_{s} ≡ Ω_{K}t_{s} (Youdin & Goodman 2005). We include the linearized vertical gravity due to the central star on the particles but ignore it for the gas, since the computational domain considered in this work is so small compared to the vertical scale height of the gas that there is no appreciable vertical density stratification in the gas over the domain. For simplicity, we assume that all particles have the same stopping time (we refer to the discussion in Sect. 5). We also ignore collisions between and selfgravity of the particles in order to isolate the effects driven by the streaming instability. The standard sheared periodic boundary conditions are imposed (Brandenburg et al. 1995; Hawley et al. 1995), and we assume the vertical dimension is also periodic.
To evolve this system of gas and particles, we use the Pencil Code^{2}, a highorder finitedifference simulation code for astrophysical fluids and particles (Brandenburg & Dobler 2002). It employs sixthorder centered differences for all spatial derivatives and a thirdorder RungeKutta method to integrate the system of equations. To maintain numerical stability, hyperdiffusion operators are needed in each dynamical field variable for the gas, and we fix the mesh Reynolds number for these operators so that noise close to the Nyquist frequency is properly damped while the power over large dynamical range is preserved (Yang & Krumholz 2012). Furthermore, we use a sixthorder Bspline interpolation to integrate the shear advection terms (cf., Johansen et al. 2009a). As mentioned above, the solid particles are modeled as Lagrangian (super)particles, which have individual positions and velocities that are integrated in unison with the Runge–Kutta steps. Their interaction with the Eulerian gas is achieved by the standard particlemesh method (e.g., Hockney & Eastwood 1988). To obtain high spatial accuracy, we choose the TriangularShapedCloud scheme for the particlemesh interpolation and assignment (Youdin & Johansen 2007).
Since we focus our attention in this work on small particles with potentially strong local solid concentrations, the mutual drag force between the gas and the particles is stiff. In this regard, we employ the numerical algorithm developed in Yang & Johansen (2016) to relieve the stringent timestep constraint due to this stiffness of the drag force. This algorithm intricately decomposes the globally coupled system of equations and makes it possible to integrate the equations on a cellbycell basis. It uses the analytical solutions in each cell to achieve numerical stability with an arbitrary time step. The momentum feedback from the particles to the gas cells is also an essential ingredient to expedite numerical convergence.
For comparison purposes, we investigate both 2D (radialvertical) and 3D models. The computational domain has a size of either 0.2H or 0.4H in each direction, where H is the vertical scale height of the gas. We consider solid particles of τ_{s} = 10^{2} and τ_{s} = 10^{3}, which are approximately cm and mm in size (when their porosity is low) at ~2.5 au in the primordial MMSN, respectively (e.g., Johansen et al. 2014). At later stages of disk evolution, these stopping times represent particles approximately ten times smaller (Bitsch et al. 2015). An external radial pressure gradient to the gas is imposed so that the azimuthal velocity of the gas is reduced by Δv ≡ Πc_{s} = 0.05c_{s}, a typical value in the inner region of the MMSN, where c_{s} is the isothermal speed of sound (Bai & Stone 2010a; Bitsch et al. 2015). This external radial pressure gradient remains constant throughout the entire duration of each simulation.
The initial conditions are set as follows. Given that it experiences no vertical gravity, the gas has a uniform density of ρ_{0}. On the other hand, we allocate, on average, one particle per gas cell but randomly position the particles, such that the particle density distribution is vertically Gaussian with a constant scale height of 0.02H. The initial scale height of particles is arbitrarily chosen in order to reduce the computing time for the initial sedimentation process, especially of small particles. All the (super)particles have the same mass, which is determined by the solidtogas mass ratio Z ≡ Σ_{p,0}/ Σ_{g,0}, where Σ_{p,0} and are the initial uniform column densities of the particles and the gas, respectively (Yang & Johansen 2014). With the density fields of the gas and the particles fixed, we apply the Nakagawa et al. (1986) equilibrium solutions to the horizontal velocities of the gas and the particles. Their initial vertical velocities are set zero.
Model specifications.
The atomic nature of the particles with a finite, constant mass limits the regions these particles can probe. The mass of each particle is given by (1)where Δx and Δy are the sizes of each cell in x and ydirections, respectively, is the average number of particles per cell, and N_{z} is the number of cells in zdirection. This implies that the mass density each particle contributes to a cell is on the order of(2)where Δz and L_{z} are the sizes of each cell and the computational domain in zdirection, respectively. Since each cell is sampled discretely by the particles, ρ_{p,1} is also the minimum density of particles a simulation model can represent. For L_{z} = 0.2H and Z = 0.02, . At equilibrium state, the distribution of particles resembles a Gaussian function, and the maximum altitude a simulation model can reach is approximately (3)where H_{p} is the scale height of the particles. For , L_{z} = 0.2H, and H_{p} = 0.015H, z_{max} ≃ 1.8H_{p}, at which the density of particles is ~19% of that in the midplane.
Even though the layer of particles our simulation models can sample is restricted to  z  ≲ z_{max}, the dynamics of the particlegas system near the midplane remains representative. The particles in the highaltitude regions ( z  ≳ z_{max}) only constitutes ~7% of the total mass of the particle layer in the above estimate. The scale height of the particles is hence well approximated by the standard deviation of the vertical particle position z_{p}: (4)which we use to identify H_{p} throughout this work. We also note that near the midplane, the number of particles per cell is on the order of (5)For , L_{z} = 0.2H, and H_{p} = 0.015H, n_{p}(0) ≃ 5. In Appendix A, we vary the number of particles in our setup for comparison purposes.
Fig. 1 Evolution of the particle density distribution for 2D models with particles of dimensionless stopping time τ_{s} = 10^{2}. The top, middle, and bottom panels show the models with an increasing solid abundance of Z = 0.01, 0.02, and 0.04, respectively, and the time t in terms of the local orbital period P increases from left to right. The particle density ρ_{p} is measured with respect to the initial gas density in the midplane ρ_{0}, and the radial and vertical positions are expressed in terms of the vertical scale height of the gas H. The models have a computational domain of 0.2H × 0.2H and a resolution of 2560H^{1}. 

Open with DEXTER 
All the models investigated here are listed in Table 1. For each set of model parameters, we evolve the system for at least 1000P, where P = 2π/ Ω_{K} is the local orbital period. This is a significantly long simulation time compared to previous works in the literature, which is necessary to capture the timescale for the streaming instability to operate on and concentrate sedimented small particles, as discussed in detail in the following sections.
3. Twodimensional models
In this section, we discuss the evolution of the density distribution of the particle layer found in our 2D models. The system is axisymmetric in this case, with only radial and vertical variations. We divide the discussion of the models by the dimensionless stopping time τ_{s} of the particles.
3.1. Particles of τ_{s} = 10^{2}
Figure 1 shows the evolution of the particle layer from three of our 2D models with particles of τ_{s} = 10^{2}. They have the same computational domain of 0.2H × 0.2H and resolution of 2560H^{1}, but differ in their solid abundance Z: Z = 0.01, 0.02, and 0.04.
Fig. 2 Scale height of the particle layer H_{p} (top panels) and maximum of the local particle density ρ_{p} (bottom panels) as a function of time for all our 2D models that trigger strong concentration of solids at high resolutions. Each column represents one set of dimensionless stopping time τ_{s} and solid abundance Z. The scale height H_{p}, particle density ρ_{p}, and time t are normalized by the vertical scale height of the gas H, the initial midplane gas density ρ_{0}, and the local orbital period P, respectively. The solid and dotted lines are for a computational domain of 0.2H × 0.2H and 0.4H × 0.4H, respectively, and the resolutions are differentiated by different colors. 

Open with DEXTER 
We first compare the cases of Z = 0.01 and Z = 0.02. For t ≲ 100P, both models follow similar evolution. The particles continue their sedimentation toward the midplane during the first ~20P, from the initial scale height of 0.02H down to ~0.013H and ~0.012H for Z = 0.01 and 0.02, respectively (Fig. 2). We note that the (efolding) timescale for sedimentation is approximately (2πτ_{s})^{1}P, which is ~16P for particles of τ_{s} = 10^{2}. In the meantime, random voids of various sizes driven by the streaming instability are developed and they move around throughout the particle layer, which is the same phenomenon found in the unstratified simulations of saturated streaming turbulence of tightlycoupled particles (Johansen & Youdin 2007; Yang & Johansen 2016). After t ~ 20P, the particles do not sediment anymore due to their random motion. Furthermore, vertical oscillation of the particle layer starts to develop and the layer becomes increasingly more corrugated (Fig. 1). In this manner, a small but appreciable fraction of the particles can be slung away from the midplane and move close to the vertical boundary. As a result, the scale height of the particle layer gradually increases with time and reaches a level of ~0.020H and ~0.016H for Z = 0.01 and 0.02, respectively, when t ~ 100P (Fig. 2). A similar corrugation mode was also reported by LorénAguilar & Bate (2015), who suggested that it is driven by largescale toroidal vortices. However, we have not found evidence for such largescale vortices in the particle layer in our models, due to complicated smallscale particle flow in the streaming turbulence. Moreover, it remains unclear how boundary conditions affect this structure, which needs to be further investigated (R. Li et al., in prep.; see also the discussion on the effect of other aspects of the model below).
Fig. 3 Evolution of the radial concentration of the particle layer for our 2D models that trigger strong concentration of solids. Three cases of different dimensionless stopping time τ_{s}, solid abundance Z, and model resolution are presented. The colors show the column density of the particles Σ_{p} as a function of radial position x and time t, where Σ_{p}, x, and t are normalized by the initial gas column density Σ_{g,0}, the vertical scale height of the gas H, and the local orbital period P, respectively. In each case, the left and right panels are from models with a computational domain of 0.2H× 0.2H and 0.4H× 0.4H, respectively. The timespan for case (c) is 2.5 times that for cases (a) and (b). 

Open with DEXTER 
For t> 100P, the model with solid abundance Z = 0.01 maintains the same state up to t = 4000P, without any sign of further concentration of solid particles. (We only show the state up to t = 1000P in Fig. 1.) By drastic contrast, however, the model with Z = 0.02 shows appreciable radial concentration of solids in the midplane over time, as shown in Fig. 1, the left panels of Fig. 2, and Fig. 3a. A first dominant concentration begins at t ~ 100P and reaches its peak at t ~ 230P. A second even stronger concentration closely upstream appears soon afterward while the first is dispersed. This second concentration reaches its saturation state at t ~ 400P and then maintains its dominance and equilibrium up to the end of the simulation (t = 1000P). At this state, the solid concentration attains a peak local density of approximately 30 times the background gas density (Fig. 2).
The left panels of Fig. 2 show the evolution of the scale height and maximum local density of the particle layer for all of our 2D models with dimensionless stopping time τ_{s} = 10^{2} and solid abundance Z = 0.02. For the first ~50P, during which the streaming turbulence is established, it is seen that the higher the resolution, the earlier the establishment and the larger the local maximum density of solids result. This is consistent with the fact that more and more faster growing modes of the linear streaming instability as well as smallerscale structures in the particle layer are resolved (Youdin & Goodman 2005; Youdin & Johansen 2007). We note that even at our highest resolution, 5120H^{1}, the fastest growing mode of the linear streaming instability remains well underresolved, the wavelength of which is on the order of 10^{4}H (cf., Bai & Stone 2010b; Yang & Johansen 2016). It is also seen that the higher the resolution, the higher the scale height of the particle layer in this initial stage. On the other hand, the timescale to establish the corrugation of the particle layer does not seem to depend on resolution, but sensitively on computational domain. As shown in Fig. 2, these timescales for our 0.2H × 0.2H and 0.4H × 0.4H models are ~100P and ~400P, respectively. Nevertheless, the degree of the vertical excitation of the particle layer is rather similar between the two computational domains and between different resolutions.
Saturation state of the 2D models.
Fig. 4 Evolution of the particle density distribution for 2D models with particles of dimensionless stopping time τ_{s} = 10^{3}. The top and bottom panels show the models with a solid abundance of Z = 0.03 and Z = 0.04, respectively. Both of which have a computational domain of 0.2H × 0.2H and a resolution of 1280 H^{1}. The notations and the layout are otherwise the same as those in Fig. 1. 

Open with DEXTER 
The corresponding time averages of the properties shown in Fig. 2 at the saturation state are summarized in Table 2. There exists some critical model resolution above which significantly stronger solid concentration occurs; this critical resolution is around 512 grid points per dimension (~2560H^{1} and ~1280H^{1} for the 0.2H × 0.2H and 0.4H × 0.4H models, respectively). Convergence in the maximum local solid density at the saturated state for each computational domain is seen by comparing the models with the highest two resolutions. The axisymmetric filaments of solids in the lowresolution models are less efficient in collecting particles upstream and easier to become dispersed, and new filaments can sporadically form in between the existing ones, competing for solid materials. On the other hand, those in the highresolution models are well concentrated and separated, leaving little material in between for any more filaments to form, as exemplified by Fig. 3a. This dichotomy may be related with how well the scale height of the particle layer is resolved, which determines how many unstable modes of the linear streaming instability can operate in the simulation models.
Also listed in Table 2 are the estimated time to reach saturation and the final number of dominant axisymmetric solid filaments for each computational domain and resolution. There exists no obvious dependence of the saturation time on either dimensions or resolution, indicating the nonlinear and stochastic nature of the system. It depends on the contentious interaction between the first generation of filaments before a stable configuration is secured (Fig. 3). Nevertheless, we establish that the timescale for particles of τ_{s} = 10^{2} to spontaneously concentrate themselves is on the order of 400–950P. Moreover, one and three major solid filaments tend to form in our 0.2H × 0.2H and 0.4H × 0.4H models, respectively. This is consistent with the approximate 3:2 ratio in the maximum local solid density at the saturation state and indicates that the mass budget in the solid filaments in a 0.2H × 0.2H model can be overestimated (Schäfer et al. 2017). It has been shown that for τ_{s} ≃ 0.3 particles, a model with at least 0.4H on each side is needed to generate consistent density distribution of solids in the saturated state of the streaming instability, where multiple solid filaments can form (Yang & Johansen 2014; R. Li et al., in prep.). Therefore, our measurement for the 0.4H × 0.4H models is likely to be more accurate.
Our models demonstrate that the critical solid abundance Z_{c} above which spontaneous strong concentration of particles of τ_{s} = 10^{2} is triggered lies in 0.01 < Z_{c} < 0.02; this is a somewhat smaller value than the range 0.02 < Z_{c} < 0.03 reported in Carrera et al. (2015). The reason for this discrepancy is that their models did not have enough resolution and simulation time. Carrera et al. (2015) used a 128 × 128 grid for a computational domain of 0.2H × 0.2H. With an initial solid abundance of Z = 0.005, they evolved the models for ~50P then artificially reduced the background gas density exponentially with time so that Z = 0.08 was attained at t ~ 200P, when the simulations ended. However, as discussed above, we find that a grid of at least 512 × 512 with a simulation time of at least 400P is necessary to properly establish the saturation state of the solid concentration for particles of τ_{s} = 10^{2}, a condition which was not satisfied by the models of Carrera et al. (2015). We show in the following subsection that this difference in measured values of critical solid abundance is significantly more drastic for particles of τ_{s} = 10^{3}, which further indicates the importance of simulating small particles with high resolutions and long simulation times.
For comparison purposes, we also investigate the case of a solid abundance of Z = 0.04 for particles of τ_{s} = 10^{2}. The evolution of the particle scale height and the maximum local solid density for various resolutions and computational domains is shown in the middle panels of Fig. 2, and the properties at the saturation stage for these models are listed in Table 2. Similar to the case of Z = 0.02, the streaming turbulence is established for the first ~50P, and then the particle layer is slightly stirred up by the turbulence up to t ~ 100–200P for the 0.2H × 0.2H models and t ~ 400–600P for the 0.4H × 0.4H models. However, the equilibrium scale height of the particle layer is appreciably less, at a level of ~0.012H, and the corrugation mode is much less pronounced than those observed in the models for Z = 0.02, as shown in Fig. 1. In general, higher solid abundance gives a more sedimented layer of particles.
On the other hand, conspicuous radial concentrations of solids occur as early as t ~ 50P, as shown in Fig. 1, the middle panels of Fig. 2, and Fig. 3b. At this point, multiple concentrated filaments of solids with a density of a few tens of initial midplane gas density ρ_{0} are already present, in contrast to the models for Z = 0.02. In the following few hundreds of orbital periods, these filaments undergo a few merging events, forming denser filaments. At the saturation stage when t ≳ 500–900P, the 0.2H × 0.2H models obtain a peak solid density of ~100–200ρ_{0} with two dense filaments, while the 0.4H × 0.4H models reach ~200–300ρ_{0} with four dense filaments (Table 2). Therefore, more solid filaments form with higher solid abundance Z, and the concentration of solids in these filaments appears to be a super linear function of Z.
3.2. Particles of τ_{s} = 10^{3}
The evolution of the particle layer in our 2D models for particles of τ_{s} = 10^{3} is qualitatively similar to that for particles of τ_{s} = 10^{2} described in Sect. 3.1, but there exist essential quantitative differences. Figure 4 shows such an evolution for two models with solid abundance Z = 0.03 and Z = 0.04, respectively, both of which have the same computational domain of 0.2H × 0.2H and resolution of 1280H^{1}. Similar to particles of τ_{s} = 10^{2}, particles of τ_{s} = 10^{3} initially continue their sedimentation process until the streaming turbulence is well saturated and supporting the particle layer. Given that the sedimentation timescale for particles of τ_{s} = 10^{3} is ~160P, being one order of magnitude longer than that for particles of τ_{s} = 10^{2}, the balance between sedimentation and turbulent excitation is not established until t ~ 100–200P, as shown in the right panels of Fig. 2. At this point, the scale height of the particle layer is ~0.016H, somewhat higher than that for particles of τ_{s} = 10^{2}. We note that the largest voids of particles in the streaming turbulence are significantly smaller than those found in the models of particles of τ_{s} = 10^{2}, which is consistent with the fact that the critical wavelength of the linear streaming instability decreases with decreasing τ_{s} (Youdin & Goodman 2005; Youdin & Johansen 2007). The development of the corrugation of the particle layer can also be seen afterward. However, the magnitude of this effect is not as strong as that for particles of τ_{s} = 10^{2}. The scale height of the particles slowly increases with time and reaches its peak at t ~ 450P, but the resulting excitation of the particle layer is barely noticeable (Fig. 4).
We find that particles of τ_{s} = 10^{3} with a solid abundance of either Z = 0.02 or Z = 0.03 remain in this statistically equilibrium state without strong concentration up to t = 1000P and t = 5000P, respectively, in any of the model specifications listed in Table 1. On the other hand, particles of τ_{s} = 10^{3} with Z = 0.04 can strongly concentrate themselves in the midplane later on (Fig. 4). As shown in Fig. 3c, one axisymmetric solid filament begins to develop at t ~ 700P and increases its concentration over time. The concentration reaches its peak at t ~ 1200P with a solid density of ~160ρ_{0}, where ρ_{0} is the initial gas density in the midplane (Fig. 2). The filament is saturated and maintains its equilibrium state up to the end of the simulation (t = 2500P).
The right panels of Fig. 2 shows the evolution of the scale height and maximum local density of the particle layer for all of our 2D models with dimensionless stopping time τ_{s} = 10^{3} and solid abundance Z = 0.04. By comparing with the other models in Fig. 2, it is evident that although particles of either size follow a similar pattern of evolution, all the timescales involved with particles of τ_{s} = 10^{3} are significantly longer than those with particles of τ_{s} = 10^{2}. As is discussed above, these include the timescales for the initial balance between sedimentation and streaming turbulence to establish, for the corrugation mode to develop, and for the radial concentration of the particles in the midplane to reach saturation, the latter requiring more than 1000P. We note also that in contrast to particles of τ_{s} = 10^{2}, the initial scale height of the particles of τ_{s} = 10^{3} slightly decreases with increasing resolution.
Table 2 lists the timeaveraged properties of the particle layer at the final saturation stage for all of our 2D models. Similar to particles of τ_{s} = 10^{2}, there exists a critical resolution above which strong concentration of particles of τ_{s} = 10^{3} occurs; this resolution is ~1280H^{1} for both 0.2H × 0.2H and 0.4H × 0.4H models. We observe that the timescale for particles of τ_{s} = 10^{3} to spontaneously concentrate themselves into dense axisymmetric filaments is on the order of 1000–2000P. One and two major solid filaments tend to form in our 0.2H × 0.2H and 0.4H × 0.4H models, respectively. The maximum local solid density in these filaments can be as high as ~300ρ_{0}. This is one order of magnitude higher than what is achieved by particles of τ_{s} = 10^{2} with an abundance of Z = 0.02. It becomes comparable when either type of particles has Z = 0.04, but the former have relatively stronger concentration due to the smaller number of filaments formed. Finally, the scale height of the particles in these models is also noticeably reduced when strong concentration occurs, which is due to further sedimentation of particles in the dense filaments.
Our models indicate that the critical solid abundance Z_{c} needed to trigger strong concentration of particles of τ_{s} = 10^{3} lies in the range 0.03 <Z_{c}< 0.04. By contrast, Carrera et al. (2015) did not find any sign of significant concentration for solid abundances up to Z ~ 0.2. It should be clear by now that the reason for this discrepancy is due to the long timescale (≳1000P) and high model resolution (≳1280H^{1}) required to realize the process. The 2D models in Carrera et al. (2015) had a resolution of 640H^{1} and a simulation time of ~200P. As shown in Fig. 2, only an equilibrated streaming turbulence can be observed during this period.
Fig. 5 Side and top views of the particle layer for a 3D model with particles of τ_{s} = 10^{2} and a solid abundance of Z = 0.02. The time t in terms of the local orbital period P increases from left to right. The top panels show the azimuthal average of the particle density ρ_{p} with respect to the initial gas density in the midplane ρ_{0}, while the bottom panels show the column density of solids Σ_{p} with respect to the initial column density of gas Σ_{0}. The coordinates are normalized by the vertical scale height of the gas H. The model has a computational domain of 0.2H on each side and a resolution of 640H^{1}. 

Open with DEXTER 
As a final remark, Table 2 also lists the timeaveraged radial and vertical velocity dispersions of the solid particles in the saturated stage of our 2D models with significant concentration. They are in general about 0.2–0.3% of the speed of sound, irrespective of the stopping time and solid abundance. However, we note that being an ensemble average of all the particles, the measured velocity dispersions do not reveal their spatial variation, and tend to be biased towards less dense regions. Moreover, as presented in Appendix B, dynamically significant perturbation in the gas density may occur when the dimensionless stopping time τ_{s} ≪ 1 and the local solidtogas density ratio ϵ ≫ 1. Detailed analysis of the density and velocity structure in these models and the corresponding particlegas dynamics will be the target of a future investigation.
4. Threedimensional models
In this section, we use 3D models to study the effect of the azimuthal dimension on the spontaneous concentration of small solid particles. We focus on a computational domain of 0.2 times gas scale height H on each side and a solid abundance that can trigger significant solid concentration in the 2D models reported in Sect. 3. We note that Johansen & Youdin (2007) found that the morphological features of the streaming turbulence for tightly coupled particles in unstratified models are not axisymmetric, which can lead to stronger density fluctuations. By comparing each 3D model with its 2D counterpart, we find how the azimuthal dimension can change the properties of a strictly axisymmetric distribution of sedimented solids.
Fig. 6 Scale height of the particle layer H_{p} (top panels) and maximum of the local particle density ρ_{p} (bottom panels) as a function of time for all our 3D models. Each column represents one set of dimensionless stopping time τ_{s} and solid abundance Z. All the models have a computational domain of 0.2H on each side. Different line colors represent different resolutions. The dotted lines are obtained by first averaging over the azimuthal direction before taking the maximum of ρ_{p}. 

Open with DEXTER 
4.1. Particles of τ_{s} = 10^{2}
Figure 5 shows the evolution of the particle layer for particles of dimensionless stopping time τ_{s} = 10^{2} and a solid abundance of Z = 0.02, at a resolution of 640H^{1}. Similar to its 2D counterpart, the particles obtain their balance between sedimentation and streaming turbulence within ~20P, where P is the local orbital period, as also shown in the left panels of Fig. 6. The resulting scale height of particles is ~0.010H, slightly lower than that in the 2D model. From the top view of the particle disk, it can be seen that the streaming turbulence is indeed nonaxisymmetric, with smallscale filamentary structures driven by background Keplerian shear. This is consistent with what was found in unstratified models for tightly coupled particles (Johansen & Youdin 2007). From t ~ 20P to ~200P, the particles are gradually excited in the vertical direction, and the scale height slightly increases to the level of ~0.012H. In contrast to the corrugation mode found in 2D models, however, the excitation of the particles in the 3D models is not regular in the radial direction, and the level of the excitation is appreciably less.
Fig. 7 Evolution of the radial concentration of the particle layer for two 3D models. The colors show the azimuthal average of the column density of the particles Σ_{p} as a function of radial position x and time t, where Σ_{p}, x, and t are normalized by the initial gas column density Σ_{g,0}, the vertical scale height of the gas H, and the local orbital period P, respectively. The left panel has particles of dimensionless stopping time τ_{s} = 10^{2} and a solid abundance of Z = 0.02, while the right panel has τ_{s} = 10^{3} and Z = 0.04. Both models have the same computational domain of 0.2H on each side and resolution of 640H^{1}. 

Open with DEXTER 
In the meantime, multiple radial concentrations of solids near the midplane begin to appear, as shown in Figs. 5 and 7. As the solids continue to accumulate, the filamentary structures become more and more aligned in the azimuthal direction. Interestingly, though, these filaments migrate radially outwards, contrary to what is seen in the 2D models. Similar behavior was also seen in the 3D models conducted by Carrera et al. (2015), who used an explicit method to integrate the mutual drag force. We present a more detailed analysis of this behavior in Appendix C. In spite of this, the maximum local density of solids in these filaments reaches approximately ten times the initial gas density in the midplane ρ_{0} at t ~ 600P (Fig. 6), a similar level and timescale observed in the 2D model of the same resolution and dimensions (Fig. 2). At the end of the simulation when t = 1000P, the system has not yet obtained a statistically steady state. We also note that the resolution we have achieved here remains lower than the critical resolution required to trigger significant concentration in the corresponding 2D models (Sect. 3.1).
Also shown in the left panels of Fig. 6 are the evolution of the scale height of the particle layer and the maximum local solid density for otherwise the same models but with two lower resolutions. Similar to the trend found in the 2D models for particles of τ_{s} = 10^{2} (Sect. 3.1), the equilibrium scale height increases slightly with increasing resolution at t ~ 20P, and all the models undergo moderate vertical excitation of particles during t ~ 20P to ~150P. The two lowerresolution models, however, only drive one axisymmetric filamentary concentration of particles during t ~ 200P to 1000P. Because of this, the model with a resolution of 320H^{1} reaches a higher local solid density of ~20ρ_{0} in the filament than the model with a resolution of 640H^{1}. It remains to be investigated whether or not a convergent saturation state can be achieved with even higherresolution models.
4.2. Particles of τ_{s} = 10^{3}
Figure 8 shows the evolution of the particle layer for particles of dimensionless stopping time τ_{s} = 10^{3} and a solid abundance of Z = 0.04, at a resolution of 640H^{1}. Similar to its 2D counterpart at the same resolution, the balance between sedimentation and streaming turbulence is attained at t ~ 50P, as also shown in the right panels of Fig. 6. The scale height of the particle layer at this point, though, is ~0.011H, which is appreciably less than ~0.017H maintained by the 2D model (Fig. 2). Nevertheless, the particles are gradually excited to a level of scale height ~0.013H over the period of 50P ≲ t ≲ 200P.
As shown in the right panel of Fig. 7 as well as Fig. 8, three filaments of concentrated solids begin to develop at t ~ 100P. Similar to the 3D model for particles of τ_{s} = 10^{2} at the same resolution of 640H^{1} (Sect. 4.1), two of the three filaments migrate outwards. They eventually merge with the third at t ~ 250P and t ~ 550P, respectively. As a result, the merged filament of solids becomes so dense that it virtually stops any further radial drift, and maintains its state up to the end of the simulation at t = 1000P. The maximum local solid density reaches the order of 10^{2}ρ_{0} and the average scale height of the particles slightly decreases due to further sedimentation in the dense filament, as shown in the right panels of Fig. 6. This concentration of solids is at the same level of what is obtained by the corresponding 2D models, albeit at a lower resolution than the critical one required by the latter (Sect. 3.2).
The right panels of Fig. 6 show the scale height of the particle layer and the maximum local solid density as a function of time for our 3D models with particles of τ_{s} = 10^{3} and Z = 0.04 at three different resolutions. The model with the resolution of 160H^{1} does not show any sign of significant concentration of solids over the course of the simulation. At the resolution of 320H^{1}, accumulation of solids into one broad filament appears, and the density of the filament gradually and steadily increases with time, reaching ~30ρ_{0} at t = 1000P. As discussed above, the model with the highest resolution of 640H^{1} shows strong concentration of solids and the system reaches its final saturated state at t ~ 650P, which is considerably less than the timescale of ~1000P required by its 2D counterpart (Fig. 2 and Table 2).
5. Implications for planetesimal formation
The 2D and 3D models we present in Sects. 3 and 4 indicate that significant spontaneous concentration of solids in protoplanetary disks can occur at a considerably less solid abundance than reported in Carrera et al. (2015), especially for particles as small as τ_{s} = 10^{3}. It is a matter of timescale and resolution for the process to operate. For particles of τ_{s} = 10^{2} and those of τ_{s} = 10^{3}, we identify that timescales of 400–1000P and 600–2000P are required, respectively. We note that even though a considerably longer timescale is necessary for small particles to concentrate themselves, their radial drift timescale is also proportionately longer due to their small stopping times. The radial drift timescale is given by (Adachi et al. 1976): for τ_{s} ≪ 1, where Π ≡ Δv/c_{s} and Δv is the reduction in gas velocity due to radial pressure gradient. The timescales for τ_{s} = 10^{2} and τ_{s} = 10^{3} are thus one order of magnitude longer than those required for driving radial concentration of solids by the streaming turbulence. Therefore, these small particles do not suffer from the radialdrift barrier any worse than large particles. Moreover, the radialdrift barrier is further alleviated once significant concentration of solids is established due to the back reaction of the drag force from the solids to the gas (Nakagawa et al. 1986; Johansen & Youdin 2007).
Fig. 8 Side and top views of the particle layer for a 3D model with particles of τ_{s} = 10^{3} and a solid abundance of Z = 0.04. The notations, layout, and model specification are otherwise the same as those in Fig. 5. 

Open with DEXTER 
Given our findings of a much reduced critical solid abundance Z_{c} for which particles of τ_{s} ≪ 1 can spontaneously concentrate themselves, we hereby revise the Z_{c} curve as a function of τ_{s} reported in Carrera et al. (2015). The critical curve for τ_{s} ≳ 0.1 is the same as in Carrera et al. (2015) and is given by (8)where log is the logarithm with base 10. On the other hand, we arbitrarily choose a quadratic function of log τ_{s} for log Z_{c} such that it passes Z_{c} = 0.035 at τ_{s} = 10^{3} and smoothly joins the other curve at τ_{s} = 0.1 with a slope of zero. The resulting function is given by (9)The updated critical curve is shown in Fig. 9.
Moreover, we note that even though particles of τ_{s} = 10^{2} can concentrate at a lower solid abundance of Z = 0.02, their concentration is one order of magnitude less dense than particles of τ_{s} = 10^{3} at Z = 0.04. The maximum local solid density reached by particles of τ_{s} = 10^{3} is about 300 times the background gas density in the midplane ρ_{0}, while that by particles of τ_{s} = 10^{2} is only about 20–30ρ_{0}. By comparing Z = 0.02 and Z = 0.04 for particles of τ_{s} = 10^{2} (Sect. 3.1), we note also that with only a difference of twice the solid abundance, the response of the saturated state for particles of τ_{s} = 10^{2} can be more than one order of magnitude stronger. At Z = 0.04, they reach a peak density of ~100–300ρ_{0}, on the same order of that by particles of τ_{s} = 10^{3}. However, due to the formation of relatively more dense filaments of solids at smaller separations for particles of τ_{s} = 10^{2}, the saturated peak density remains comparatively lower than that of particles of τ_{s} = 10^{3}.
This observation leads to an interesting scenario for planetesimal formation. Solid particles of smaller sizes might be more predisposed towards gravitational collapse inside selfinduced dense filaments of solids than those of larger sizes. To see this, the Roche density ρ_{R} is given by where M_{⋆} and a are the mass of and the distance to the central star, respectively, and ρ_{0} is the gas density in the midplane. The normalization is applied at the location of the Asteroid Belt in the MMSN. We note that the observed gas densities in the terrestrial region of a young protoplanetary disk can be one or two magnitudes higher (we refer to e.g., Dutrey et al. 2014, and references therein), and a gas density of 10^{8}–10^{10} g cm^{3} for a< 5 au can also be seen in numerical simulations of starirradiated accretion disk models at their earliest evolution stages (Bitsch et al. 2015; Baillié et al. 2016). In any case, the solid concentration achieved in our models with particles of τ_{s} = 10^{3} is well above the Roche density in the terrestrial region of a wide range of disks when the solid abundance Z ~ 0.04, while it may require particularly dense disks or further outward locations for particles of τ_{s} = 10^{2} to drive planetesimal formation unless Z ≳ 0.04 also. Most importantly, we have demonstrated that inside the ice line, strong concentration of mmsized solid particles, whose growth is limited by bouncing and fragmentation (Zsom et al. 2010; Birnstiel et al. 2012), and hence planetesimal formation therein is possible, bridging the problematic gap between dust coagulation and planetesimal formation (Dra¸żkowska & Dullemond 2014).
Despite that, for a disk with an average solid abundance below the critical one, some mechanisms are still required to enhance the abundance in a local but somewhat larger scale for planetesimals to be able to form via the streaming instability. Some interesting possibilities include photoevaporation of the gaseous disk (Carrera et al. 2017), radial pileup of solids (Dra¸żkowska et al. 2016; Gonzalez et al. 2017), and ice condensation and sublimation (Ros & Johansen 2013; Ida & Guillot 2016; Schoonenberg & Ormel 2017).
Fig. 9 Revised critical solid abundance as a function of the dimensionless stopping time τ_{s} from Carrera et al. (2015). The filled and open circles show our models with and without significant radial concentration of solids, respectively. The solid black line for τ_{s} ≳ 10^{1} and the dashed blue line for τ_{s} ≲ 10^{1} are the unmodified critical curve from Carrera et al. (2015). The black solid line for τ_{s} ≲ 10^{1} shows the revised part of the critical curve, where the black dashed line is the extrapolation of it. For solid abundances Z above the black line (green region), spontaneous concentration of solid particles by the streaming instability can occur. 

Open with DEXTER 
Some discussion on our simplifying assumption of particles of the same size is in order here. It has been suggested by Bai & Stone (2010a) that the particlegas dynamics for a population of particles of various sizes is predominantly driven by the largest ones. They found a critical solid abundance of Z_{c} ~ 0.05 for particles of seven species ranged from τ_{s} = 10^{4} to 10^{1}, and conjectured that it is the total abundance of the particles with τ_{s} ≳ 10^{2} that forms the criterion for triggering strong concentration of particles; that is, Z(τ_{s} ≳ 10^{2}) ≳ 0.02. However, if the dust coagulation is limited by the bouncing barrier, the largest particles in the distribution may only have τ_{s} ~ 10^{3} (Zsom et al. 2010; Dra¸żkowska & Dullemond 2014). In this case, the conjecture by Bai & Stone (2010a) breaks down and it is not clear which size range of the particles would dictate the particlegas dynamics in the system. On the other hand, when the dust coagulation is limited by the bouncing barrier, the distribution of the particle size tends to concentrate within about one order of magnitude near the barrier (Zsom et al. 2010). Therefore, our assumption of particles of the same size might still assimilate this scenario.
Finally, several open questions remain in concentrating small particles by the streaming instability. We have observed further sedimentation in dense filaments in our models, indicating that the velocity dispersion of the particles decreases in the dense regions. It remains to be seen if further dust growth could proceed inside such an environment, given the reduced relative velocities, and if the dust growth in turn could drive more spontaneous concentration of solid materials via the streaming instability. Moreover, the presence of the magnetic fields could render rich structure in the gaseous protoplanetary disks, and the disk midplane can be turbulent with various strengths at different locations (Turner et al. 2014), which significantly affects the ability of solid materials to sediment onto the midplane (e.g., Okuzumi & Hirose 2011; Zhu et al. 2015). This has nontrivial effects on the dust coagulation and the radial concentration of solids for planetesimal formation to occur (Johansen et al. 2014; Testi et al. 2014).
6. Summary
In this work, we revisit the condition for the spontaneous concentration of solid particles driven by the streaming instability in the regime of dimensionless stopping time τ_{s} ≪ 1, as originally investigated by Carrera et al. (2015). Using the localshearingbox approximation, we simulate a sedimented layer of solid particles of the same size in a laminar gaseous environment, developing streaming turbulence and possibly the ensuing radial concentration of particles. With the assistance of the numerical algorithm for the stiff mutual drag force by Yang & Johansen (2016), we are able to conduct the same simulation models of Carrera et al. (2015) with significantly higher resolutions and longer simulation times. We focus on two stopping times τ_{s} = 10^{2} and 10^{3}, and systematically vary the solidtogas column density ratio Z and the resolution, with two computational domains of 0.2H × 0.2H and 0.4H × 0.4H in the radialvertical plane, where H is the local scale height of the gas. We also conduct a few 3D counterparts of the same models and find similar results.
We find that small solid particles can indeed spontaneously concentrate themselves into high density via the streaming instability, given enough time and resolution. For particles of τ_{s} = 10^{2} and 10^{3}, the critical solid abundance Z_{c} above which strong concentration of solids occurs is in the range of 0.01 <Z_{c}< 0.02 and 0.03 <Z_{c}< 0.04, respectively. The timescales required for this process to reach saturation are approximately 400–1000P and 600–2000P, respectively, where P is the local orbital period. For τ_{s} = 10^{2} at a solid abundance of Z = 0.02, a maximum local solid density of 20–30ρ_{0} is observed, where ρ_{0} is the initial midplane density of the gas, while τ_{s} = 10^{3} at a solid abundance of Z = 0.04, results in a density of ~300ρ_{0}. We also find a superlinear dependence of the solid concentration on the solid abundance. With these new measurements, we hereby revise the critical solid abundance curve of Carrera et al. (2015) in the regime of τ_{s} ≪ 1 in Fig. 9. For a solid density of ≳10, direct gravitational collapse of mmsized particles can occur and form planetesimals in the terrestrial region of typical protoplanetary disks, as long as the solid abundance reaches the level of a few percent. This may help alleviate the problematic size gap between the largest particles from dust coagulation and the smallest particles that can spontaneously concentrate themselves via the streaming instability, especially inside the ice line of a protoplanetary disk.
The behavior of bouncing might only occur for dust aggregates of volume filling factor greater than ~0.3 (Wada et al. 2011) or ~0.5 (Seizinger & Kley 2013), as determined by numerical simulations.
The Pencil Code is publicly available at http://pencilcode.nordita.org/.
Acknowledgments
We thank the anonymous referee, Tristan Guillot, and Joanna Dra¸żkowska for their constructive comments, which helped to significantly improve the manuscript. All the simulation models presented in this work were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC Center for High Performance Computing in KTH Royal Institute of Technology and at Lunarc in Lund University, both of which are located in Sweden. This research was supported by the European Research Council under ERC Starting Grant agreement 278675PEBBLE2PLANET. A.J. is grateful for financial support from the Knut and Alice Wallenberg Foundation and from the Swedish Research Council (grant 20103710).
References
 Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Prog. Th. Phys., 56, 1756 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010a, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010b, ApJS, 190, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2014, ApJ, 796, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Baillié, K., Charnoz, S., & Pantin, E. 2016, A&A, 590, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] [Google Scholar]
 Baruteau, C., Crida, A., Paardekooper, S.J., et al. 2014, in Protostars and Planets VI, eds. T. Montmerle, D. Ehrenreich, & A.M. Lagrange (Tucson: Univ. Arizona Press), 667 [Google Scholar]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bracco, A., Chavanis, P. H., Provenzale, A., & Spiegel, E. A. 1999, Phys. Fluids, 11, 2280 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., & Dobler, W. 2002, Comput. Phys. Commun., 147, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741 [NASA ADS] [CrossRef] [Google Scholar]
 Brauer, F., Dullemond, C. P., Johansen, A., et al. 2007, A&A, 469, 1169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Drążkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Drążkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dutrey, A., Semenov, D., Chapillon, E., et al. 2014, in Protostars and Planets VI, eds. T. Montmerle, D. Ehrenreich, & A.M. Lagrange (Tucson: Univ. Arizona Press), 317 [Google Scholar]
 Goldreich, P., & LyndenBell, D. 1965, MNRAS, 130, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, J.F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984 [NASA ADS] [Google Scholar]
 Haisch, Jr., K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Hayashi, C. 1981, Prog. Th. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1988, Computer Simulation Using Particles (New York: CRC Press) [Google Scholar]
 Ida, S., & Guillot, T. 2016, A&A, 596, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A., & Klahr, H. 2009a, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A., & Mac Low, M.M. 2009b, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, A&A, 537, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, eds. T. Montmerle, D. Ehrenreich, & A.M. Lagrange (Tucson: Univ. Arizona Press), 547 [Google Scholar]
 Johansen, A., Mac Low, M.M., Lacerda, P., & Bizzarro, M. 2015, Sci. Adv., 1, 1500109 [NASA ADS] [CrossRef] [Google Scholar]
 Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 557, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21 [NASA ADS] [CrossRef] [Google Scholar]
 LorénAguilar, P., & Bate, M. R. 2015, MNRAS, 453, L78 [NASA ADS] [CrossRef] [Google Scholar]
 Lyra, W. 2014, ApJ, 789, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Mamajek, E. E. 2009, in AIP Conf. Ser. 1158, eds. T. Usuda, M. Tamura, & M. Ishii, 3 [Google Scholar]
 Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Terquem, C. 2006, Rep. Prog. Phys., 69, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Ros, K., & Johansen, A. 2013, A&A, 552, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Safronov, V. S. 1969, Evoliutsiia doplanetnogo oblaka (Moscow: Nauka) [Google Scholar]
 Schäfer, U., Yang, C.C., & Johansen, A. 2017, A&A, 597, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schoonenberg, D., & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seizinger, A., & Kley, W. 2013, A&A, 551, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, ApJ, 847, L12 [NASA ADS] [CrossRef] [Google Scholar]
 Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, eds. T. Montmerle, D. Ehrenreich, & A.M. Lagrange (Tucson: Univ. Arizona Press), 339 [Google Scholar]
 Turner, N. J., Fromang, S., Gammie, C., et al. 2014, in Protostars and Planets VI, eds. T. Montmerle, D. Ehrenreich, & A.M. Lagrange (Tucson: Univ. Arizona Press), 411 [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1977a, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1977b, Ap&SS, 51, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Whipple, F. L. 1972, in TwentyFirst Nobel Symposium, From Plasma to Planet, ed. A. Elvius (New York: Wiley Interscience Division), 211 [Google Scholar]
 Yang, C.C., & Johansen, A. 2014, ApJ, 792, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, C.C., & Johansen, A. 2016, ApJS, 224, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, C.C., & Krumholz, M. 2012, ApJ, 758, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N. 2010, in Physics and Astrophysics of Planetary Systems, eds. T. Montmerle, D. Ehrenreich, & A.M. Lagrange (Les Ulis: EDP), EAS Pub. Ser., 41, 187 [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Stone, J. M., & Bai, X.N. 2015, ApJ, 801, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Particle resolution
The solid particles in a local patch of the protoplanetary disk are too numerous to be evolved individually in current numerical simulations. Therefore, the approach of the socalled superparticles or representative particles is often adopted. Each superparticle models an arbitrarily high number of identical physical particles and assumes the collective dynamical properties of them. The more superparticles are used, the more the system of physical particles is resolved.
On the other hand, under the particlemesh construct, the numerical accuracy of treating any interaction between the gas and the particles is limited by the resolution of the fixed grid (e.g., Youdin & Johansen 2007), and hence advantage may not necessarily be gained by a great number of particles within one grid cell. In the context of the nonlinear evolution of the streaming instability without vertical gravity, it has been shown that the results are insensitive to the number of superparticles used (Johansen & Youdin 2007; Bai & Stone 2010b). In particular, Bai & Stone (2010b) showed that the resolution of the grid is appreciably more important in obtaining a convergent density distribution of particles, and one particle per cell on average is sufficient in reproducing the distribution.
To investigate the potential effects of the particle resolution on our results, we select one of our 2D models presented in Sect. 3.1 and simulate the same model with ten times more superparticles, that is, with an average number of particles per cell of (see Sect. 2). The model has particles of dimensionless stopping time τ_{s} = 10^{2}, a solid abundance of Z = 0.02, a computational domain of 0.2H × 0.2H, and a resolution of 2560H^{1}. Figure A.1 compares the scale height of the particle layer and the maximum local particle density as a function of time for this model with and 10 (cf., Fig. 2). The simulation with follows similar evolution as the one with : initial sedimentation, saturation of the streaming turbulence, excitation of the particle layer, and finally radial concentration of the solid particles. The former reaches the saturated state somewhat earlier than the latter (t ~ 200P versus t ~ 400P) and has the same scale height of (0.0140 ± 0.003)H. On the other hand, the simulation with obtains a saturated maximum local solid density of (12 ± 1)ρ_{0}, roughly a factor of two smaller than that of the simulation with . The radial and vertical velocity dispersions are also slightly lower at Δv_{p,x} = (0.0027 ± 0.0001)c_{s} and Δv_{p,z} = (0.0028 ± 0.0001)c_{s}, respectively (cf., Table 2).
Fig. A.1 Scale height of the particle layer H_{p} (top panel) and maximum of the local particle density ρ_{p} (bottom panel) as a function of time for two 2D models with different numbers of superparticles. Both models have particles of dimensionless stopping time τ_{s} = 10^{2}, a solid abundance of Z = 0.02, a computational domain of 0.2H × 0.2H, and a resolution of 2560H^{1}. The blue and green lines represent an average number of particles per cell of and 10, respectively. 

Open with DEXTER 
Fig. A.2 Sideview of the particle density at time t = 300P from the 2D model of Fig. A.1 with an average number of particles per cell of . 

Open with DEXTER 
The reason for the smaller maximum solid density in the simulation with is that two wellseparated axisymmetric filaments of solids are formed in the computational domain, as shown in Fig. A.2 (cf., Fig. 1). We note that the simulation with also exhibits a transient stage of two coexisting filaments around t ~ 240P before one final dominant filament emerges (Fig. 3a). As demonstrated in Figs. 3, nonlinear interaction between adjacent filaments may or may not occur in different systems, and the final number of dominant filaments may be uncertain by one or two. Therefore, the maximum local density of particles can be uncertain within a factor of a few in the saturated state of radial concentration of solid materials.
Appendix B: Radial pressure gradient in the midplane
As far as the back reaction from the solid particles to the gas via the mutual drag force is concerned, interesting dynamics may occur when the local solidtogas density ratio ϵ ≡ ρ_{p}/ρ_{g} ≳ 1. In this case, while solids tend to radially drift towards the star, they have significant collective momentum to drive the gas to move radially outwards by the conservation of angular momentum. It has been suggested that a local pressure bump in the gas may form via this mechanism when ϵ ~ 1 and τ_{s} ~ 1 (Gonzalez et al. 2017).
Fig. B.1 Final density profiles of gas (blue, left axis) and particles (green, right axis) in the midplane for various 2D models with the dimensions of 0.2H × 0.2H. The background density gradient is added to the gas profile for easy recognition of any pressure bump. 

Open with DEXTER 
In Figs. B.1, we plot the final density profiles of the gas and the particles in the midplane for several of our 2D models, in which strong local concentration of solids occurs. The background density gradient −(2ΠΩ_{K}/c_{s})ρ_{0}x (see Sect. 2) is added to the gas profile for easy recognition of any pressure bump in our systems. In all our models, the solidtogas density ratio in the midplane reaches ϵ ~ 1 simply due to the process of sedimentation. For particles of τ_{s} = 10^{2}, we find that the perturbation in the gas density is in general small with respect to the background density gradient unless ϵ ≳ 10. Even for ϵ ~ 10–100, a pressure bump barely forms near a local density peak of solids (Figs. B.1a and b). On the other hand, it appears that a local pressure bump that coincides with the local density peak of solids does readily occur for particles of τ_{s} = 10^{3} when ϵ ≫ 1 (Fig. B.1c). Therefore, it seems that in our systems, τ_{s} ≪ 1 and ϵ ≫ 1 is required to generate dynamically significant perturbation in the gas via the back reaction of the drag force.
Appendix C: Nonaxisymmetric particlegas dynamics in three dimensions
The 3D models presented in Sect. 4 exhibit an interesting phenomenon that does not occur in the 2D models described in Sect. 3. The dense filaments of solid materials appear to migrate radially outwards, as shown in Fig. 7. In this appendix, we discuss the particlegas dynamics associated with this phenomenon in more detail.
Fig. C.1 Azimuthaveraged radial profiles of the particle/gas properties in the midplane at t = 600P for the 3D model with particles of τ_{s} = 10^{2} and a resolution of 640H^{1}. The top, middle, and bottom panels show densities, radial velocities, and azimuthal velocities with respect to the background Keplerian shear, respectively. The blue (left axis) and the green (right axis) lines represent the gas and the solid particles, respectively. The background density gradient is added to the gas profile for easy recognition of any pressure bump (see Appendix B). The gas and the particle densities are normalized by the initial gas density in the midplane ρ_{0}, while the velocities are weighted by the densities and normalized by the velocity reduction of the gas Δv = Πc_{s} due to external radial pressure gradient (Sect. 2). We note that a zero velocity is equivalent to moving at the Keplerian velocity, irrespective of the location. 

Open with DEXTER 
Figure C.1 shows the azimuthaveraged radial profiles of the particle/gas properties in the midplane at t = 600P for the 3D model with particles of dimensionless stopping time τ_{s} = 10^{2} and a resolution of 640H^{1} (cf., Figs. 5 and 7). The dense filament of solids at this moment is located at x ~ +0.05H. Due to the tight coupling between the gas and the particles, the relative velocity between them are small with respect to the speed of sound c_{s}. The azimuthal velocities of the gas and the particles near the filament remain subKeplerian. More importantly, the radial velocities of the particles in the filament are predominantly inward, on the order of −0.02Πc_{s}. For comparison, the outward radial velocity of the filament shown in the left panel of Fig. 7 is about + 4 × 10^{4}Πc_{s}. On the other hand, the gas density profile near the filament indicates a weak density bump. However, the magnitude of the radial pressure gradient this bump generates is appreciably less than the external one. In other words, this density bump is not a dynamically significant one. Therefore, we conclude that the apparent outward radial migration of the dense filament of solids is a pattern movement, instead of a material one.
Fig. C.2 Streamlines in the midplane at t = 600P for the 3D model with particles of τ_{s} = 10^{2} and a resolution of 640H^{1} for a) the particles and b) the gas. The background shows the densities of the particles and the gas, respectively. 

Open with DEXTER 
The outward pattern movement of the filament of solids seems to result from the largescale nonaxisymmetric structure of the particlegas system, which is absent in our 2D models. As shown in Fig. 5, the solid particles concentrate not only radially, but also azimuthally. Figure C.2 shows a slice through the midplane at the same time from the same model as in Fig. C.1. Also shown are the streamlines of the gas and the particles. We note that the streamlines are slightly bent radially outwards near the concentration of the particles as well as the gas. This indicates that the particles azimuthally above the concentration tend to move outwards while those below tend to move inwards. The less dense region at the opposite azimuthal phase display an opposite trend, that is, the streamlines are slightly bent radially inwards. We also note that the radially depleted region between the filaments of solids signals smallscale vortex motions. The mechanism that leads to this type of nonaxisymmetric structure and movement remains unclear and requires future dedicated investigation.
Similar pattern of the particlegas dynamics described above is also manifested in the 3D models with particles of τ_{s} = 10^{3}. However, as shown in Fig. 8 and the right panel of Fig. 7, once the apparent radially outward movement of the two weak filaments joins the third stronger one, an axisymmetric dense filament of solids is produced. The filament appears to be almost stationary afterwards, which is consistent with the particlegas drag equilibrium solution of Nakagawa et al. (1986). This further supports the idea that the apparent radially outward movement of filaments is exclusively a nonaxisymmetric phenomenon.
All Tables
All Figures
Fig. 1 Evolution of the particle density distribution for 2D models with particles of dimensionless stopping time τ_{s} = 10^{2}. The top, middle, and bottom panels show the models with an increasing solid abundance of Z = 0.01, 0.02, and 0.04, respectively, and the time t in terms of the local orbital period P increases from left to right. The particle density ρ_{p} is measured with respect to the initial gas density in the midplane ρ_{0}, and the radial and vertical positions are expressed in terms of the vertical scale height of the gas H. The models have a computational domain of 0.2H × 0.2H and a resolution of 2560H^{1}. 

Open with DEXTER  
In the text 
Fig. 2 Scale height of the particle layer H_{p} (top panels) and maximum of the local particle density ρ_{p} (bottom panels) as a function of time for all our 2D models that trigger strong concentration of solids at high resolutions. Each column represents one set of dimensionless stopping time τ_{s} and solid abundance Z. The scale height H_{p}, particle density ρ_{p}, and time t are normalized by the vertical scale height of the gas H, the initial midplane gas density ρ_{0}, and the local orbital period P, respectively. The solid and dotted lines are for a computational domain of 0.2H × 0.2H and 0.4H × 0.4H, respectively, and the resolutions are differentiated by different colors. 

Open with DEXTER  
In the text 
Fig. 3 Evolution of the radial concentration of the particle layer for our 2D models that trigger strong concentration of solids. Three cases of different dimensionless stopping time τ_{s}, solid abundance Z, and model resolution are presented. The colors show the column density of the particles Σ_{p} as a function of radial position x and time t, where Σ_{p}, x, and t are normalized by the initial gas column density Σ_{g,0}, the vertical scale height of the gas H, and the local orbital period P, respectively. In each case, the left and right panels are from models with a computational domain of 0.2H× 0.2H and 0.4H× 0.4H, respectively. The timespan for case (c) is 2.5 times that for cases (a) and (b). 

Open with DEXTER  
In the text 
Fig. 4 Evolution of the particle density distribution for 2D models with particles of dimensionless stopping time τ_{s} = 10^{3}. The top and bottom panels show the models with a solid abundance of Z = 0.03 and Z = 0.04, respectively. Both of which have a computational domain of 0.2H × 0.2H and a resolution of 1280 H^{1}. The notations and the layout are otherwise the same as those in Fig. 1. 

Open with DEXTER  
In the text 
Fig. 5 Side and top views of the particle layer for a 3D model with particles of τ_{s} = 10^{2} and a solid abundance of Z = 0.02. The time t in terms of the local orbital period P increases from left to right. The top panels show the azimuthal average of the particle density ρ_{p} with respect to the initial gas density in the midplane ρ_{0}, while the bottom panels show the column density of solids Σ_{p} with respect to the initial column density of gas Σ_{0}. The coordinates are normalized by the vertical scale height of the gas H. The model has a computational domain of 0.2H on each side and a resolution of 640H^{1}. 

Open with DEXTER  
In the text 
Fig. 6 Scale height of the particle layer H_{p} (top panels) and maximum of the local particle density ρ_{p} (bottom panels) as a function of time for all our 3D models. Each column represents one set of dimensionless stopping time τ_{s} and solid abundance Z. All the models have a computational domain of 0.2H on each side. Different line colors represent different resolutions. The dotted lines are obtained by first averaging over the azimuthal direction before taking the maximum of ρ_{p}. 

Open with DEXTER  
In the text 
Fig. 7 Evolution of the radial concentration of the particle layer for two 3D models. The colors show the azimuthal average of the column density of the particles Σ_{p} as a function of radial position x and time t, where Σ_{p}, x, and t are normalized by the initial gas column density Σ_{g,0}, the vertical scale height of the gas H, and the local orbital period P, respectively. The left panel has particles of dimensionless stopping time τ_{s} = 10^{2} and a solid abundance of Z = 0.02, while the right panel has τ_{s} = 10^{3} and Z = 0.04. Both models have the same computational domain of 0.2H on each side and resolution of 640H^{1}. 

Open with DEXTER  
In the text 
Fig. 8 Side and top views of the particle layer for a 3D model with particles of τ_{s} = 10^{3} and a solid abundance of Z = 0.04. The notations, layout, and model specification are otherwise the same as those in Fig. 5. 

Open with DEXTER  
In the text 
Fig. 9 Revised critical solid abundance as a function of the dimensionless stopping time τ_{s} from Carrera et al. (2015). The filled and open circles show our models with and without significant radial concentration of solids, respectively. The solid black line for τ_{s} ≳ 10^{1} and the dashed blue line for τ_{s} ≲ 10^{1} are the unmodified critical curve from Carrera et al. (2015). The black solid line for τ_{s} ≲ 10^{1} shows the revised part of the critical curve, where the black dashed line is the extrapolation of it. For solid abundances Z above the black line (green region), spontaneous concentration of solid particles by the streaming instability can occur. 

Open with DEXTER  
In the text 
Fig. A.1 Scale height of the particle layer H_{p} (top panel) and maximum of the local particle density ρ_{p} (bottom panel) as a function of time for two 2D models with different numbers of superparticles. Both models have particles of dimensionless stopping time τ_{s} = 10^{2}, a solid abundance of Z = 0.02, a computational domain of 0.2H × 0.2H, and a resolution of 2560H^{1}. The blue and green lines represent an average number of particles per cell of and 10, respectively. 

Open with DEXTER  
In the text 
Fig. A.2 Sideview of the particle density at time t = 300P from the 2D model of Fig. A.1 with an average number of particles per cell of . 

Open with DEXTER  
In the text 
Fig. B.1 Final density profiles of gas (blue, left axis) and particles (green, right axis) in the midplane for various 2D models with the dimensions of 0.2H × 0.2H. The background density gradient is added to the gas profile for easy recognition of any pressure bump. 

Open with DEXTER  
In the text 
Fig. C.1 Azimuthaveraged radial profiles of the particle/gas properties in the midplane at t = 600P for the 3D model with particles of τ_{s} = 10^{2} and a resolution of 640H^{1}. The top, middle, and bottom panels show densities, radial velocities, and azimuthal velocities with respect to the background Keplerian shear, respectively. The blue (left axis) and the green (right axis) lines represent the gas and the solid particles, respectively. The background density gradient is added to the gas profile for easy recognition of any pressure bump (see Appendix B). The gas and the particle densities are normalized by the initial gas density in the midplane ρ_{0}, while the velocities are weighted by the densities and normalized by the velocity reduction of the gas Δv = Πc_{s} due to external radial pressure gradient (Sect. 2). We note that a zero velocity is equivalent to moving at the Keplerian velocity, irrespective of the location. 

Open with DEXTER  
In the text 
Fig. C.2 Streamlines in the midplane at t = 600P for the 3D model with particles of τ_{s} = 10^{2} and a resolution of 640H^{1} for a) the particles and b) the gas. The background shows the densities of the particles and the gas, respectively. 

Open with DEXTER  
In the text 