EDP Sciences
Free Access
Issue
A&A
Volume 606, October 2017
Article Number A80
Number of page(s) 16
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201630106
Published online 16 October 2017

© ESO, 2017

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; 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 compact1 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 τs ≡ ΩKts and thus the size of the particles, where ΩK is the local Keplerian frequency and ts 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-4H 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 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 τ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 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 ΩK 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 ts (Whipple 1972; Weidenschilling 1977a) or its dimensionless counterpart τs ≡ ΩKts (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 Code2, 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 (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 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.

For comparison purposes, we investigate both 2D (radial-vertical) 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 ≡ Πcs = 0.05cs, a typical value in the inner region of the MMSN, where cs 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 solid-to-gas 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.

Table 1

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 y-directions, respectively, is the average number of particles per cell, and Nz is the number of cells in z-direction. This implies that the mass density each particle contributes to a cell is on the order of(2)where Δz and Lz are the sizes of each cell and the computational domain in z-direction, 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 Lz = 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 Hp is the scale height of the particles. For , Lz = 0.2H, and Hp = 0.015H, zmax ≃ 1.8Hp, 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 | z | ≲ zmax, the dynamics of the particle-gas system near the mid-plane remains representative. The particles in the high-altitude regions (| z | ≳ zmax) 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 zp: (4)which we use to identify Hp throughout this work. We also note that near the mid-plane, the number of particles per cell is on the order of (5)For , Lz = 0.2H, and Hp = 0.015H, np(0) ≃ 5. In Appendix A, we vary the number of particles in our setup for comparison purposes.

thumbnail 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 mid-plane ρ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. 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 τ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.

thumbnail Fig. 2

Scale height of the particle layer Hp (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 Hp, particle density ρp, and time t are normalized by the vertical scale height of the gas H, the initial mid-plane 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 mid-plane 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 (e-folding) timescale for sedimentation is approximately (2πτs)-1P, 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 tightly-coupled 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 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.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é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., in prep.; see also the discussion on the effect of other aspects of the model below).

thumbnail 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 mid-plane 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 smaller-scale 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 under-resolved, the wavelength of which is on the order of 10-4H (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.

Table 2

Saturation state of the 2D models.

thumbnail 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 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. 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 Zc above which spontaneous strong concentration of particles of τs = 10-2 is triggered lies in 0.01 < Zc < 0.02; this is a somewhat smaller value than the range 0.02 < Zc < 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 mid-plane 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 ~ 100200P, 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 mid-plane 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 mid-plane (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 mid-plane 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 time-averaged 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 Zc needed to trigger strong concentration of particles of τs = 10-3 lies in the range 0.03 <Zc< 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.

thumbnail 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 mid-plane ρ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 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 τs ≪ 1 and the local solid-to-gas density ratio ϵ ≫ 1. 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 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.

thumbnail Fig. 6

Scale height of the particle layer Hp (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 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 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.

thumbnail 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 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 ρ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 lower-resolution 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 higher-resolution 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 50Pt ≲ 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 102ρ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/cs 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 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).

thumbnail 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 Zc for which particles of τs ≪ 1 can spontaneously concentrate themselves, we hereby revise the Zc 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 Zc such that it passes Zc = 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 mid-plane ρ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 self-induced 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 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 10-810-10 g cm-3 for a< 5 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 τ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 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).

thumbnail 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 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 Zc ~ 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 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 τs ≪ 1, 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 τs = 10-2 and 10-3, and systematically vary the solid-to-gas column density ratio Z and the resolution, with two computational domains of 0.2H × 0.2H and 0.4H × 0.4H in the radial-vertical 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 Zc above which strong concentration of solids occurs is in the range of 0.01 <Zc< 0.02 and 0.03 <Zc< 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 mid-plane 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 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 τs ≪ 1 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.


1

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.

2

The Pencil Code is publicly available at http://pencil-code.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 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

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.

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 τ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 Δvp,x = (0.0027 ± 0.0001)cs and Δvp,z = (0.0028 ± 0.0001)cs, respectively (cf., Table 2).

thumbnail Fig. A.1

Scale height of the particle layer Hp (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 super-particles. 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

thumbnail Fig. A.2

Side-view 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 well-separated 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 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 ϵρ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).

thumbnail Fig. B.1

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 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 mid-plane for several of our 2D models, in which strong local concentration of solids occurs. The background density gradient −(2ΠΩK/cs)ρ0x (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 ϵ ~ 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: 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.

thumbnail Fig. C.1

Azimuth-averaged radial profiles of the particle/gas properties in the mid-plane 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 mid-plane ρ0, while the velocities are weighted by the densities and normalized by the velocity reduction of the gas Δv = Πcs 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 azimuth-averaged radial profiles of the particle/gas properties in the mid-plane 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 cs. 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 −0.02Πcs. For comparison, the outward radial velocity of the filament shown in the left panel of Fig. 7 is about + 4 × 10-4Πcs. 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.

thumbnail Fig. C.2

Streamlines in the mid-plane 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 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. Figure C.2 shows a slice through the mid-plane 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 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 τ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 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.

All Tables

Table 1

Model specifications.

Table 2

Saturation state of the 2D models.

All Figures

thumbnail 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 mid-plane ρ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
thumbnail Fig. 2

Scale height of the particle layer Hp (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 Hp, particle density ρp, and time t are normalized by the vertical scale height of the gas H, the initial mid-plane 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
thumbnail 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
thumbnail 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
thumbnail 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 mid-plane ρ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
thumbnail Fig. 6

Scale height of the particle layer Hp (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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. A.1

Scale height of the particle layer Hp (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 super-particles. 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
thumbnail Fig. A.2

Side-view 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
thumbnail Fig. B.1

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 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
thumbnail Fig. C.1

Azimuth-averaged radial profiles of the particle/gas properties in the mid-plane 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 mid-plane ρ0, while the velocities are weighted by the densities and normalized by the velocity reduction of the gas Δv = Πcs 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
thumbnail Fig. C.2

Streamlines in the mid-plane 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.