Concentrating small particles in protoplanetary disksthrough the streaming instability

Concentrating small particles in protoplanetary disks
through the streaming instability

C.-C. Yang (楊朝欽) Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-221 00 Lund, Sweden
11email: ccyang@astro.lu.se
   A. Johansen Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-221 00 Lund, Sweden
11email: ccyang@astro.lu.se
   D. Carrera Present address: Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, U.S.A.Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-221 00 Lund, Sweden
11email: ccyang@astro.lu.se
Received 21 November 2016 / Accepted 14 June 2017
Key Words.:
hydrodynamics – instabilities – methods: numerical – minor planets, asteroids: general – planets and satellites: formation – protoplanetary disks
{CJK}

UTF8bkai

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/cm-sized 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 and – representing mm- and cm-sized 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 . 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.

1 Introduction

In the core-accretion scenario of planet formation, planetary cores are assembled beginning with interstellar m-sized 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 kilometer-scale planetesimals (we refer to e.g., Johansen et al. 2014, and references therein).

A major obstacle to the formation of planetesimals is the “radial-drift 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 meter-sized boulders at 1 au of the minimum-mass 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 cross-section 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 sub-micron 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 long-lived 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 action-reaction 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 time-scales, but also planetesimals can form by direct gravitational collapse within these dense filaments (Johansen et al. 2012, 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 compact111The 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. 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 two-dimensional simulations, Carrera et al. (2015) found that this threshold abundance increases drastically with decreasing particle size. For mm-sized particles inside the ice line of protoplanetary disks, they suggested that a solid-to-gas 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 solid-to-gas 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 fastest-growing modes are roughly proportional to the dimensionless stopping time and thus the size of the particles, where is the local Keplerian frequency and 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 for particles with , where 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 when , 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 time-step constraint of the stiff mutual drag force, it is numerically challenging to simulate the particle-gas 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 time-step 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 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 while 3–4% for particles of . 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 two-dimensional (2D; Sect. 3) and three-dimensional (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 local-shearing-box approximation (Goldreich & Lynden-Bell 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 and its three axes constantly align in the radial, azimuthal, and vertical directions, respectively. We consider isothermal, non-magnetized 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 (Whipple 1972; Weidenschilling 1977a) or its dimensionless counterpart (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 self-gravity 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 Code222The Pencil Code is publicly available at http://pencil-code.nordita.org/, a high-order finite-difference simulation code for astrophysical fluids and particles (Brandenburg & Dobler 2002). It employs sixth-order centered differences for all spatial derivatives and a third-order Runge-Kutta method to integrate the system of equations. To maintain numerical stability, hyper-diffusion 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 sixth-order B-spline interpolation to integrate the shear advection terms (c.f., 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 particle-mesh method (e.g., Hockney & Eastwood 1988). To obtain high spatial accuracy, we choose the Triangular-Shaped-Cloud scheme for the particle-mesh 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 time-step 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 cell-by-cell 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.

Dimensionless Solid Computational Resolutions Simulation Time
Stopping Time Abundance Domain () ()
2D (radial-vertical) Models
0.01 640, 1280, 2560 4000
320, 640, 1280 4000
0.02 640, 1280, 2560, 5120 1000
320, 640, 1280, 2560 1000
0.04 640, 1280, 2560 1000
320, 640, 1280 1000
0.02 640, 1280, 2560 1000
320, 640, 1280 1000
0.03 640, 1280, 2560 5000
320, 640, 1280 5000
0.04 640, 1280, 2560 2500
320, 640, 1280 2500
3D Models
0.02 160, 320, 640 1000
0.04 160, 320, 640 1000
Table 1: Model specifications

For comparison purposes, we investigate both 2D (radial-vertical) and 3D models. The computational domain has a size of either 0.2 or 0.4 in each direction, where is the vertical scale height of the gas. We consider solid particles of and , 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 , a typical value in the inner region of the MMSN, where 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 . 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.02. 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 solid-to-gas mass ratio , where 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–Sekiya–Hayashi (1986) equilibrium solutions to the horizontal velocities of the gas and the particles. Their initial vertical velocities are set zero.

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 and are the sizes of each cell in - and -directions, respectively, is the average number of particles per cell, and is the number of cells in -direction. This implies that the mass density each particle contributes to a cell is on the order of

(2)

where and are the sizes of each cell and the computational domain in -direction, respectively. Since each cell is sampled discretely by the particles, is also the minimum density of particles a simulation model can represent. For and , . At equilibrium state, the distribution of particles resembles a Gaussian function, and the maximum altitude a simulation model can reach is approximately

(3)

where is the scale height of the particles. For , , and , , at which the density of particles is 19% of that in the mid-plane.

Even though the layer of particles our simulation models can sample is restricted to , the dynamics of the particle-gas system near the mid-plane remains representative. The particles in the high-altitude regions () 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 :

(4)

which we use to identify throughout this work. We also note that near the mid-plane, the number of particles per cell is on the order of

(5)

For , , and , . In Appendix A, we vary the number of particles in our setup for comparison purposes.

Figure 1: Evolution of the particle density distribution for 2D models with particles of dimensionless stopping time . The top, middle, and bottom panels show the models with an increasing solid abundance of , 0.02, and 0.04, respectively, and the time in terms of the local orbital period increases from left to right. The particle density is measured with respect to the initial gas density in the mid-plane , and the radial and vertical positions are expressed in terms of the vertical scale height of the gas . The models have a computational domain of 0.2 0.2 and a resolution of 2560.

All the models investigated here are listed in Table 1. For each set of model parameters, we evolve the system for at least 1000, where 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 Two-dimensional 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 of the particles.

3.1 Particles of

Figure 1 shows the evolution of the particle layer from three of our 2D models with particles of . They have the same computational domain of and resolution of 2560, but differ in their solid abundance : , 0.02, and 0.04.

Figure 2: Scale height of the particle layer (top panels) and maximum of the local particle density (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 and solid abundance . The scale height , particle density , and time are normalized by the vertical scale height of the gas , the initial mid-plane gas density , and the local orbital period , respectively. The solid and dotted lines are for a computational domain of 0.2 0.2 and 0.4 0.4, respectively, and the resolutions are differentiated by different colors.

We first compare the cases of and . For , both models follow similar evolution. The particles continue their sedimentation toward the mid-plane during the first 20, from the initial scale height of 0.02 down to 0.013 and 0.012 for and 0.02, respectively (Fig. 2). We note that the (-folding) timescale for sedimentation is approximately , which is 16 for particles of . 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 tightly-coupled particles (Johansen & Youdin 2007; Yang & Johansen 2016). After , 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 mid-plane 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.020 and 0.016 for and 0.02, respectively, when (Fig. 2). A similar corrugation mode was also reported by Lorén-Aguilar & Bate (2015), who suggested that it is driven by large-scale toroidal vortices. However, we have not found evidence for such large-scale vortices in the particle layer in our models, due to complicated small-scale 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. 2016, in preparation; see also the discussion on the effect of other aspects of the model below).

(a)
(b)
(c)
Figure 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 , solid abundance , and model resolution are presented. The colors show the column density of the particles as a function of radial position and time , where , , and are normalized by the initial gas column density , the vertical scale height of the gas , and the local orbital period , respectively. In each case, the left and right panels are from models with a computational domain of 0.2 0.2 and 0.4 0.4, respectively. The timespan for case (c) is 2.5 times that for cases (a) and (b).

For , the model with solid abundance maintains the same state up to , without any sign of further concentration of solid particles. (We only show the state up to in Fig. 1.) By drastic contrast, however, the model with shows appreciable radial concentration of solids in the mid-plane over time, as shown in Fig. 1, the left panels of Fig. 2, and Fig. 2(a). A first dominant concentration begins at and reaches its peak at . A second even stronger concentration closely upstream appears soon afterward while the first is dispersed. This second concentration reaches its saturation state at and then maintains its dominance and equilibrium up to the end of the simulation (). 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 and solid abundance . For the first 50, 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 smaller-scale structures in the particle layer are resolved (Youdin & Goodman 2005; Youdin & Johansen 2007). We note that even at our highest resolution, 5120, the fastest growing mode of the linear streaming instability remains well under-resolved, the wavelength of which is on the order of (c.f., 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 and models are 100 and 400, respectively. Nevertheless, the degree of the vertical excitation of the particle layer is rather similar between the two computational domains and between different resolutions.

Dimensions Resolution
() () () () () ()
(1) (2) (3) (4) (5) (6) (7) (8)
,
0.2 0.2 0640 550 1–2 0.0155(4) 09(2) 0.0028(1) 0.0028(1)
1280 800 1–2 0.0150(4) 07(1) 0.0027(1) 0.0026(1)
2560 400 1 0.0141(3) 28(4) 0.0032(1) 0.0033(2)
5120 450 1 0.0140(2) 26(3) 0.0038(2) 0.0040(2)
0.4 0.4 0320 500 1–5 0.0172(5) 05(1) 0.0024(1) 0.0027(2)
0640 950 2–3 0.0156(2) 07(1) 0.0022(1) 0.0022(1)
1280 800 3 0.0144(3) 18(4) 0.0023(1) 0.0023(1)
2560 500 3 0.0153(4) 19(2) 0.0029(1) 0.0030(1)
,
0.2 0.2 0640 350 2 0.0120(2) 079(15) 0.0029(1) 0.0027(1)
1280 500 2 0.0119(2) 176(73) 0.0027(1) 0.0027(1)
2560 950 2 0.0110(1) 114(9) 0.0027(1) 0.0027(1)
0.4 0.4 0320 600 5–6 0.0130(2) 021(3) 0.0026(1) 0.0025(1)
0640 750 4 0.0113(2) 225(45) 0.0023(1) 0.0021(1)
1280 650 4 0.0110(4) 339(70) 0.0021(1) 0.0020(1)
,
0.2 0.2 0640 1200 1 0.0159(2) 024(6) 0.0034(1) 0.0033(1)
1280 1200 1 0.0144(1) 156(18) 0.0028(1) 0.0027(1)
2560 1500 1 0.0141(1) 297(35) 0.0023(0) 0.0023(0)
0.4 0.4 0320 0.0172(1) 0.0032(1) 0.0029(1)
0640 2000 3–4 0.0159(0) 007(1) 0.0026(1) 0.0024(1)
1280 2200 2 0.0137(0) 266(23) 0.0021(0) 0.0020(0)
333Columns: (1) Dimensions of the computational domain in terms of the gas scale height ; (2) Resolution in number of cells per ; (3) Estimated time to reach saturation in terms of the local orbital period ; (4) Number of major axisymmetric solid filaments; (5) Time-averaged scale height of the particles in ; (6) Time-averaged maximum local particle density in terms of the initial mid-plane gas density ; (7) Time-averaged radial velocity dispersion in terms of the isothermal speed of sound ; (8) Time-averaged vertical velocity dispersion in . The time averages are taken from to the end of the simulation with their standard deviation shown in parentheses.
Table 2: Saturation state of the 2D models

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 (2560 and 1280 for the and 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 low-resolution 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 high-resolution models are well concentrated and separated, leaving little material in between for any more filaments to form, as exemplified by Fig. 2(a). 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 to spontaneously concentrate themselves is on the order of 400–950. Moreover, one and three major solid filaments tend to form in our and 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 model can be overestimated (Schäfer et al. 2017). It has been shown that for particles, a model with at least 0.4 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. 2016, in preparation). Therefore, our measurement for the models is likely to be more accurate.

Figure 4: Evolution of the particle density distribution for 2D models with particles of dimensionless stopping time . The top and bottom panels show the models with a solid abundance of and , respectively. Both of which have a computational domain of and a resolution of 1280 . The notations and the layout are otherwise the same as those in Fig. 1.

Our models demonstrate that the critical solid abundance above which spontaneous strong concentration of particles of is triggered lies in ; this is a somewhat smaller value than the range 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 128128 grid for a computational domain of . With an initial solid abundance of , they evolved the models for 50 then artificially reduced the background gas density exponentially with time so that was attained at , when the simulations ended. However, as discussed above, we find that a grid of at least 512512 with a simulation time of at least 400 is necessary to properly establish the saturation state of the solid concentration for particles of , 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 , 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 for particles of . 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 , the streaming turbulence is established for the first 50, and then the particle layer is slightly stirred up by the turbulence up to –200 for the models and –600 for the models. However, the equilibrium scale height of the particle layer is appreciably less, at a level of 0.012, and the corrugation mode is much less pronounced than those observed in the models for , 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 , as shown in Fig. 1, the middle panels of Fig. 2, and Fig. 2(b). At this point, multiple concentrated filaments of solids with a density of a few tens of initial mid-plane gas density are already present, in contrast to the models for . In the following few hundreds of orbital periods, these filaments undergo a few merging events, forming denser filaments. At the saturation stage when –900, the models obtain a peak solid density of 100–200 with two dense filaments, while the models reach 200–300 with four dense filaments (Table 2). Therefore, more solid filaments form with higher solid abundance , and the concentration of solids in these filaments appears to be a super linear function of .

3.2 Particles of

The evolution of the particle layer in our 2D models for particles of is qualitatively similar to that for particles of described in Sect. 3.1, but there exist essential quantitative differences. Figure 4 shows such an evolution for two models with solid abundance and , respectively, both of which have the same computational domain of and resolution of 1280. Similar to particles of , particles of 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 is 160, being one order of magnitude longer than that for particles of , the balance between sedimentation and turbulent excitation is not established until , as shown in the right panels of Fig. 2. At this point, the scale height of the particle layer is 0.016, somewhat higher than that for particles of . We note that the largest voids of particles in the streaming turbulence are significantly smaller than those found in the models of particles of , which is consistent with the fact that the critical wavelength of the linear streaming instability decreases with decreasing (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 . The scale height of the particles slowly increases with time and reaches its peak at , but the resulting excitation of the particle layer is barely noticeable (Fig. 4).

We find that particles of with a solid abundance of either or remain in this statistically equilibrium state without strong concentration up to and , respectively, in any of the model specifications listed in Table 1. On the other hand, particles of with can strongly concentrate themselves in the mid-plane later on (Fig. 4). As shown in Fig. 2(c), one axisymmetric solid filament begins to develop at and increases its concentration over time. The concentration reaches its peak at with a solid density of 160, where is the initial gas density in the mid-plane (Fig. 2). The filament is saturated and maintains its equilibrium state up to the end of the simulation ().

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 and solid abundance . 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 are significantly longer than those with particles of . 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 mid-plane to reach saturation, the latter requiring more than 1000. We note also that in contrast to particles of , the initial scale height of the particles of slightly decreases with increasing resolution.

Table 2 lists the time-averaged properties of the particle layer at the final saturation stage for all of our 2D models. Similar to particles of , there exists a critical resolution above which strong concentration of particles of occurs; this resolution is 1280 for both and models. We observe that the timescale for particles of to spontaneously concentrate themselves into dense axisymmetric filaments is on the order of 1000–2000. One and two major solid filaments tend to form in our and models, respectively. The maximum local solid density in these filaments can be as high as 300. This is one order of magnitude higher than what is achieved by particles of with an abundance of . It becomes comparable when either type of particles has , 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 needed to trigger strong concentration of particles of lies in the range . By contrast, Carrera et al. (2015) did not find any sign of significant concentration for solid abundances up to . It should be clear by now that the reason for this discrepancy is due to the long timescale (1000) and high model resolution (1280) required to realize the process. The 2D models in Carrera et al. (2015) had a resolution of 640 and a simulation time of 200. As shown in Fig. 2, only an equilibrated streaming turbulence can be observed during this period.

Figure 5: Side and top views of the particle layer for a 3D model with particles of and a solid abundance of . The time in terms of the local orbital period increases from left to right. The top panels show the azimuthal average of the particle density with respect to the initial gas density in the mid-plane , while the bottom panels show the column density of solids with respect to the initial column density of gas . The coordinates are normalized by the vertical scale height of the gas . The model has a computational domain of 0.2 on each side and a resolution of 640.
Figure 6: Scale height of the particle layer (top panels) and maximum of the local particle density (bottom panels) as a function of time for all our 3D models. Each column represents one set of dimensionless stopping time and solid abundance . All the models have a computational domain of 0.2 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 .

As a final remark, Table 2 also lists the time-averaged 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 and the local solid-to-gas density ratio . Detailed analysis of the density and velocity structure in these models and the corresponding particle-gas dynamics will be the target of a future investigation.

4 Three-dimensional 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 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.

4.1 Particles of

Figure 5 shows the evolution of the particle layer for particles of dimensionless stopping time and a solid abundance of , at a resolution of 640. Similar to its 2D counterpart, the particles obtain their balance between sedimentation and streaming turbulence within 20, where is the local orbital period, as also shown in the left panels of Fig. 6. The resulting scale height of particles is 0.010, 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 non-axisymmetric, with small-scale 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 to 200, the particles are gradually excited in the vertical direction, and the scale height slightly increases to the level of 0.012. 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.

Figure 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 as a function of radial position and time , where , , and are normalized by the initial gas column density , the vertical scale height of the gas , and the local orbital period , respectively. The left panel has particles of dimensionless stopping time and a solid abundance of , while the right panel has and . Both models have the same computational domain of 0.2 on each side and resolution of 640.

In the meantime, multiple radial concentrations of solids near the mid-plane 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 mid-plane at (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 , 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 (Sect. 3.1), the equilibrium scale height increases slightly with increasing resolution at , and all the models undergo moderate vertical excitation of particles during to 150. The two lower-resolution models, however, only drive one axisymmetric filamentary concentration of particles during to 1000. Because of this, the model with a resolution of 320 reaches a higher local solid density of 20 in the filament than the model with a resolution of 640. It remains to be investigated whether or not a convergent saturation state can be achieved with even higher-resolution models.

4.2 Particles of

Figure 8: Side and top views of the particle layer for a 3D model with particles of and a solid abundance of . The notations, layout, and model specification are otherwise the same as those in Fig. 5.

Figure 8 shows the evolution of the particle layer for particles of dimensionless stopping time and a solid abundance of , at a resolution of 640. Similar to its 2D counterpart at the same resolution, the balance between sedimentation and streaming turbulence is attained at , as also shown in the right panels of Fig. 6. The scale height of the particle layer at this point, though, is 0.011, which is appreciably less than 0.017 maintained by the 2D model (Fig. 2). Nevertheless, the particles are gradually excited to a level of scale height 0.013 over the period of .

As shown in the right panel of Fig. 7 as well as Fig. 8, three filaments of concentrated solids begin to develop at . Similar to the 3D model for particles of at the same resolution of 640 (Sect. 4.1), two of the three filaments migrate outwards. They eventually merge with the third at and , 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 . The maximum local solid density reaches the order of 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 and at three different resolutions. The model with the resolution of 160 does not show any sign of significant concentration of solids over the course of the simulation. At the resolution of 320, accumulation of solids into one broad filament appears, and the density of the filament gradually and steadily increases with time, reaching 30 at . As discussed above, the model with the highest resolution of 640 shows strong concentration of solids and the system reaches its final saturated state at , which is considerably less than the timescale of 1000 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 . It is a matter of timescale and resolution for the process to operate. For particles of and those of , we identify that timescales of 400–1000 and 600–2000 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): {align} t_drift &= (1 + τs24πΠτs) (HR)^-1P
&≃3×10^4 (τs10-3)^-1 (Π0.05)^-1 (H / R0.05)^-1 P, for , where and is the reduction in gas velocity due to radial pressure gradient. The timescales for and 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 radial-drift barrier any worse than large particles. Moreover, the radial-drift 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).

Figure 9: Revised critical solid abundance as a function of the dimensionless stopping time 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 and the dashed blue line for are the unmodified critical curve from Carrera et al. (2015). The black solid line for shows the revised part of the critical curve, where the black dashed line is the extrapolation of it. For solid abundances above the black line (green region), spontaneous concentration of solid particles by the streaming instability can occur.

Given our findings of a much reduced critical solid abundance for which particles of can spontaneously concentrate themselves, we hereby revise the curve as a function of reported in Carrera et al. (2015). The critical curve for is the same as in Carrera et al. (2015) and is given by

(6)

where is the logarithm with base 10. On the other hand, we arbitrarily choose a quadratic function of for such that it passes at and smoothly joins the other curve at with a slope of zero. The resulting function is given by

(7)

The updated critical curve is shown in Fig. 9.

Moreover, we note that even though particles of can concentrate at a lower solid abundance of , their concentration is one order of magnitude less dense than particles of at . The maximum local solid density reached by particles of is about 300 times the background gas density in the mid-plane , while that by particles of is only about 20–30. By comparing and for particles of (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 can be more than one order of magnitude stronger. At , they reach a peak density of 100–300, on the same order of that by particles of . However, due to the formation of relatively more dense filaments of solids at smaller separations for particles of , the saturated peak density remains comparatively lower than that of particles of .

This observation leads to an interesting scenario for planetesimal formation. Solid particles of smaller sizes might be more predisposed towards gravitational collapse inside self-induced dense filaments of solids than those of larger sizes. To see this, the Roche density is given by {align} ρRρ0 &= 9M⋆4πa3ρ0
&= 270(M⋆M) (a2.5 au)^-3 (ρ010-10g cm-3)^-1, where and are the mass of and the distance to the central star, respectively, and is the gas density in the mid-plane. 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 g cm for au can also be seen in numerical simulations of star-irradiated 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 is well above the Roche density in the terrestrial region of a wide range of disks when the solid abundance , while it may require particularly dense disks or further outward locations for particles of to drive planetesimal formation unless also. Most importantly, we have demonstrated that inside the ice line, strong concentration of mm-sized 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 photo-evaporation of the gaseous disk (Carrera et al. 2017), radial pile-up 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).

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 particle-gas dynamics for a population of particles of various sizes is predominantly driven by the largest ones. They found a critical solid abundance of for particles of seven species ranged from to , and conjectured that it is the total abundance of the particles with that forms the criterion for triggering strong concentration of particles; that is, . However, if the dust coagulation is limited by the bouncing barrier, the largest particles in the distribution may only have (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 particle-gas 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 mid-plane 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 mid-plane (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 , as originally investigated by Carrera et al. (2015). Using the local-shearing-box 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 and , and systematically vary the solid-to-gas column density ratio and the resolution, with two computational domains of and in the radial-vertical plane, where 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 and , the critical solid abundance above which strong concentration of solids occurs is in the range of and , respectively. The timescales required for this process to reach saturation are approximately 400–1000 and 600–2000, respectively, where is the local orbital period. For at a solid abundance of , a maximum local solid density of 20–30 is observed, where is the initial mid-plane density of the gas, while at a solid abundance of , results in a density of 300. We also find a super-linear 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 in Fig. 9. For a solid density of 10, direct gravitational collapse of mm-sized 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.

Acknowledgements.
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 278675-PEBBLE2PLANET. A.J. is grateful for financial support from the Knut and Alice Wallenberg Foundation and from the Swedish Research Council (grant 2010-3710).

References

  • Adachi et al. (1976) Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progress of Theoretical Physics, 56, 1756
  • Bai & Stone (2010a) Bai, X.-N. & Stone, J. M. 2010a, ApJ, 722, 1437
  • Bai & Stone (2010b) Bai, X.-N. & Stone, J. M. 2010b, ApJS, 190, 297
  • Bai & Stone (2014) Bai, X.-N. & Stone, J. M. 2014, ApJ, 796, 31
  • Baillié et al. (2016) Baillié, K., Charnoz, S., & Pantin, E. 2016, A&A, 590, A60
  • Barge & Sommeria (1995) Barge, P. & Sommeria, J. 1995, A&A, 295, L1
  • Baruteau et al. (2014) Baruteau, C., Crida, A., Paardekooper, S.-J., et al. 2014, in Protostars and Planets VI, ed. T. Montmerle, D. Ehrenreich, & A.-M. Lagrange (Univ. Arizona Press, Tucson, AZ), 667–689
  • Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148
  • Bitsch et al. (2015) Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28
  • Bracco et al. (1999) Bracco, A., Chavanis, P. H., Provenzale, A., & Spiegel, E. A. 1999, Physics of Fluids, 11, 2280
  • Brandenburg & Dobler (2002) Brandenburg, A. & Dobler, W. 2002, Computer Physics Communications, 147, 471
  • Brandenburg et al. (1995) Brandenburg, A., Nordlund, A., Stein, R. F., & Torkelsson, U. 1995, ApJ, 446, 741
  • Brauer et al. (2007) Brauer, F., Dullemond, C. P., Johansen, A., et al. 2007, A&A, 469, 1169
  • Carrera et al. (2017) Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16
  • Carrera et al. (2015) Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
  • Dra̧żkowska et al. (2016) Dra̧żkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105
  • Dra̧żkowska & Dullemond (2014) Dra̧żkowska, J. & Dullemond, C. P. 2014, A&A, 572, A78
  • Dutrey et al. (2014) Dutrey, A., Semenov, D., Chapillon, E., et al. 2014, in Protostars and Planets VI, ed. T. Montmerle, D. Ehrenreich, & A.-M. Lagrange (Univ. Arizona Press, Tucson, AZ), 317–338
  • Goldreich & Lynden-Bell (1965) Goldreich, P. & Lynden-Bell, D. 1965, MNRAS, 130, 125
  • Gonzalez et al. (2017) Gonzalez, J.-F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984
  • Haisch et al. (2001) Haisch, Jr., K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153
  • Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
  • Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
  • Hockney & Eastwood (1988) Hockney, R. W. & Eastwood, J. W. 1988, Computer Simulation Using Particles (CRC Press, New York, NY)
  • Ida & Guillot (2016) Ida, S. & Guillot, T. 2016, A&A, 596, L3
  • Jacquet et al. (2011) Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591
  • Johansen et al. (2014) Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, ed. T. Montmerle, D. Ehrenreich, & A.-M. Lagrange (Univ. Arizona Press, Tucson, AZ), 547–570
  • Johansen et al. (2015) Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1500109
  • Johansen & Youdin (2007) Johansen, A. & Youdin, A. 2007, ApJ, 662, 627
  • Johansen et al. (2009a) Johansen, A., Youdin, A., & Klahr, H. 2009a, ApJ, 697, 1269
  • Johansen et al. (2009b) Johansen, A., Youdin, A., & Mac Low, M.-M. 2009b, ApJ, 704, L75
  • Johansen et al. (2012) Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, A&A, 537, A125
  • Kataoka et al. (2013) Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 557, L4
  • Klahr & Hubbard (2014) Klahr, H. & Hubbard, A. 2014, ApJ, 788, 21
  • Lorén-Aguilar & Bate (2015) Lorén-Aguilar, P. & Bate, M. R. 2015, MNRAS, 453, L78
  • Lyra (2014) Lyra, W. 2014, ApJ, 789, 77
  • Mamajek (2009) Mamajek, E. E. 2009, in American Institute of Physics Conference Series, Vol. 1158, American Institute of Physics Conference Series, ed. T. Usuda, M. Tamura, & M. Ishii, 3–10
  • Nakagawa et al. (1986) Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375
  • Okuzumi & Hirose (2011) Okuzumi, S. & Hirose, S. 2011, ApJ, 742, 65
  • Papaloizou & Terquem (2006) Papaloizou, J. C. B. & Terquem, C. 2006, Reports on Progress in Physics, 69, 119
  • Ros & Johansen (2013) Ros, K. & Johansen, A. 2013, A&A, 552, A137
  • Safronov (1969) Safronov, V. S. 1969, Evoliutsiia doplanetnogo oblaka. (Nauka, Moscow)
  • Schäfer et al. (2017) Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69
  • Schoonenberg & Ormel (2017) Schoonenberg, D. & Ormel, C. W. 2017, A&A, in press [arXiv: 1702.02151]
  • Seizinger & Kley (2013) Seizinger, A. & Kley, W. 2013, A&A, 551, A65
  • Simon et al. (2016) Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55
  • Simon et al. (2017) Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, ApJ, submitted [arXiv: 1705.03889]
  • Testi et al. (2014) Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. T. Montmerle, D. Ehrenreich, & A.-M. Lagrange (Univ. Arizona Press, Tucson, AZ), 339–361
  • Turner et al. (2014) Turner, N. J., Fromang, S., Gammie, C., et al. 2014, in Protostars and Planets VI, ed. T. Montmerle, D. Ehrenreich, & A.-M. Lagrange (Univ. Arizona Press, Tucson, AZ), 411–432
  • Wada et al. (2011) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36
  • Weidenschilling (1977a) Weidenschilling, S. J. 1977a, MNRAS, 180, 57
  • Weidenschilling (1977b) Weidenschilling, S. J. 1977b, Ap&SS, 51, 153
  • Whipple (1972) Whipple, F. L. 1972, in Twenty-First Nobel Symposium, From Plasma to Planet, ed. A. Elvius (Wiley Interscience Division, New York, NY), 211
  • Yang & Johansen (2014) Yang, C.-C. & Johansen, A. 2014, ApJ, 792, 86
  • Yang & Johansen (2016) Yang, C.-C. & Johansen, A. 2016, ApJS, 224, 39
  • Yang & Krumholz (2012) Yang, C.-C. & Krumholz, M. 2012, ApJ, 758, 48
  • Youdin & Johansen (2007) Youdin, A. & Johansen, A. 2007, ApJ, 662, 613
  • Youdin (2010) Youdin, A. N. 2010, in EAS Publications Series, Vol. 41, Physics and Astrophysics of Planetary Systems, ed. T. Montmerle, D. Ehrenreich, & A.-M. Lagrange (EDP, Les Ulis), 187–207
  • Youdin & Goodman (2005) Youdin, A. N. & Goodman, J. 2005, ApJ, 620, 459
  • Zhu et al. (2015) Zhu, Z., Stone, J. M., & Bai, X.-N. 2015, ApJ, 801, 81
  • Zsom et al. (2010) Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57

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 so-called super-particles or representative particles is often adopted. Each super-particle models an arbitrarily high number of identical physical particles and assumes the collective dynamical properties of them. The more super-particles are used, the more the system of physical particles is resolved.

On the other hand, under the particle-mesh 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 super-particles 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.

Figure 10: Scale height of the particle layer (top panel) and maximum of the local particle density (bottom panel) as a function of time for two 2D models with different numbers of super-particles. Both models have particles of dimensionless stopping time , a solid abundance of , a computational domain of , and a resolution of 2560. The blue and green lines represent an average number of particles per cell of and 10, respectively.

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 super-particles, that is, with an average number of particles per cell of (see Sect. 2). The model has particles of dimensionless stopping time , a solid abundance of , a computational domain of , and a resolution of 2560. Fig. 10 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 (c.f., 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 ( versus ) and has the same scale height of . On the other hand, the simulation with obtains a saturated maximum local solid density of , roughly a factor of two smaller than that of the simulation with . The radial and vertical velocity dispersions are also slightly lower at and , respectively (c.f., Table 2).

Figure 11: Side-view of the particle density at time from the 2D model of Fig. 10 with an average number of particles per cell of .

The reason for the smaller maximum solid density in the simulation with is that two well-separated axisymmetric filaments of solids are formed in the computational domain, as shown in Fig. 11 (c.f., Fig. 1). We note that the simulation with also exhibits a transient stage of two coexisting filaments around before one final dominant filament emerges (Fig. 2(a)). 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 mid-plane

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 solid-to-gas density ratio . 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 and (Gonzalez et al. 2017).

(a) , , ,
(b) , , ,
(c) , , ,
Figure 12: Final density profiles of gas (blue, left axis) and particles (green, right axis) in the mid-plane for various 2D models with the dimensions of . The background density gradient is added to the gas profile for easy recognition of any pressure bump.

In Figs. 12, we plot the final density profiles of the gas and the particles in the mid-plane for several of our 2D models, in which strong local concentration of solids occurs. The background density gradient (see Sect. 2) is added to the gas profile for easy recognition of any pressure bump in our systems. In all our models, the solid-to-gas density ratio in the mid-plane reaches simply due to the process of sedimentation. For particles of , we find that the perturbation in the gas density is in general small with respect to the background density gradient unless . Even for –100, a pressure bump barely forms near a local density peak of solids (Figs. 11(a) and 11(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 when (Fig. 11(c)). Therefore, it seems that in our systems, and is required to generate dynamically significant perturbation in the gas via the back reaction of the drag force.

Appendix C Non-axisymmetric particle-gas 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 particle-gas dynamics associated with this phenomenon in more detail.

Figure 13: Azimuth-averaged radial profiles of the particle/gas properties in the mid-plane at for the 3D model with particles of and a resolution of 640. 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 mid-plane , while the velocities are weighted by the densities and normalized by the velocity reduction of the gas 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.

Fig. 13 shows the azimuth-averaged radial profiles of the particle/gas properties in the mid-plane at for the 3D model with particles of dimensionless stopping time and a resolution of 640 (c.f., Figs. 5 and 7). The dense filament of solids at this moment is located at . 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 . The azimuthal velocities of the gas and the particles near the filament remain sub-Keplerian. More importantly, the radial velocities of the particles in the filament are predominantly inward, on the order of . For comparison, the outward radial velocity of the filament shown in the left panel of Fig. 7 is about . 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.

(a) particles kinematics
(b) gas kinematics
Figure 14: Streamlines in the mid-plane at for the 3D model with particles of and a resolution of 640 for (a) the particles and (b) the gas. The background shows the densities of the particles and the gas, respectively.

The outward pattern movement of the filament of solids seems to result from the large-scale non-axisymmetric structure of the particle-gas system, which is absent in our 2D models. As shown in Fig. 5, the solid particles concentrate not only radially, but also azimuthally. Fig. 14 shows a slice through the mid-plane at the same time from the same model as in Fig. 13. 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 small-scale vortex motions. The mechanism that leads to this type of non-axisymmetric structure and movement remains unclear and requires future dedicated investigation.

Similar pattern of the particle-gas dynamics described above is also manifested in the 3D models with particles of . 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 particle-gas drag equilibrium solution of Nakagawa et al. (1986). This further supports the idea that the apparent radially outward movement of filaments is exclusively a non-axisymmetric phenomenon.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
226141
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description