Issue 
A&A
Volume 618, October 2018



Article Number  A75  
Number of page(s)  17  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201832783  
Published online  12 October 2018 
Streaming instability of multiple particle species in protoplanetary disks
Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University,
Box 43,
22100
Lund,
Sweden
email: noemi.schaffer@astro.lu.se
Received:
7
February
2018
Accepted:
29
July
2018
The radial drift and diffusion of dust particles in protoplanetary disks affect both the opacity and temperature of such disks, as well as the location and timing of planetesimal formation. In this paper, we present results of numerical simulations of particlegas dynamics in protoplanetary disks that include dust grains with various size distributions. We have considered three scenarios in terms of particle size ranges, one where the Stokes number τ_{s} = 10^{−1}−10^{0}, one where τ_{s} = 10^{−4}−10^{−1}, and finally one where τ_{s} = 10^{−3}−10^{0}. Moreover, we considered both discrete and continuous distributions in particle size. In accordance with previous works we find in our multispecies simulations that different particle sizes interact via the gas and as a result their dynamics changes compared to the singlespecies case. The larger species trigger the streaming instability and create turbulence that drives the diffusion of the solid materials. We measured the radial equilibrium velocity of the system and find that the radial drift velocity of the large particles is reduced in the multispecies simulations and that the small particle species move on average outwards. We also varied the steepness of the size distribution, such that the exponent of the solid number density distribution, dN∕da ∝ a^{−q}, is either q = 3 or q = 4. Overall, we find that the steepness of the size distribution and the discrete versus continuous approach have little impact on the results. The level of diffusion and drift rates are mainly dictated by the range of particle sizes. We measured the scale height of the particles and observe that small grains are stirred up well above the sedimented midplane layer where the large particles reside. Our measured diffusion and drift parameters can be used in coagulation models for planet formation as well as to understand relative mixing of the components of primitive meteorites (matrix, chondrules and CAIs) prior to inclusion in their parent bodies.
Key words: protoplanetary disks / methods: numerical / hydrodynamics / instabilities / turbulence / diffusion
© ESO 2018
1 Introduction
Protoplanets – the building blocks of terrestrial planets, as well as of the cores of gas giants, ice giants and superEarths – form by collisional growth from micrometer dust particles to bodies with sizes of thousands of kilometers. There are, however, several barriers that hinder the formation of planets. One barrier occurs early in the course of planet formation, since pebbles of millimetercentimeter sizes do not stick efficiently when they collide. Their interaction typically results in bouncing, erosion and fragmentation (Güttler et al. 2010).
Even if sticking would be perfect, radial drift limits particle growth. Due to the orbital velocity difference between the solid and gas components, the former feels a headwind and spirals in toward the central star. Meter sized solids have drift times of approximately a hundred years at 1 au (Adachi et al. 1976; Weidenschilling 1977; Youdin 2010). Particles grow maximally to a size where the growth timescale equals the radial drift timescale, resulting typically in centimeter sizes interior of 10 au and millimeter sizes in the outer disk (Birnstiel et al. 2012, 2016; Lambrechts & Johansen 2014).
One way to overcome the radial drift barrier is to concentrate solids into clumps via the streaming instability (Youdin & Goodman 2005; Youdin & Johansen 2007; Johansen & Youdin 2007). Taking into account the backreaction of the solids onto the gas, a small overdensity of particles accelerates the surrounding gas to a higher speed. The radial drift of the given particle clump is thus reduced, since it feels less drag from the surrounding gas. At the same time, growth to even larger clump size is possible by the accumulation of inwards drifting particles. Once the local solid concentration reaches the Roche density, planetesimals form through gravitational collapse (Johansen et al. 2012). The formation of filaments by the streaming instability occurs above a threshold metallicity slightly elevated relative to solar metallicity (Johansen et al. 2009; Bai & Stone 2010b; Carrera et al. 2015; Yang et al. 2017). The metallicity may be increased either by gas photoevaporation (Carrera et al. 2017) or by pileup of drifting solids in the inner regions of the protoplanetary disk (Dra̧żkowska et al. 2016; Ida & Guillot 2016; Schoonenberg & Ormel 2017; Dra̧żkowska & Alibert 2017; Gonzalez et al. 2017).
In recent years, the number of protoplanetary disk observations at different wavelengths has increased dramatically. Observations at millimeter and centimeter wavelengths reveal the presence of pebbles while scattered light observations probe the dust population. Spatially resolved ALMA observations of young disks around Class 0 protostars show evidence of the presence of millimetersized pebbles, which hints to grain growth via coagulation already at an early stage of disk formation (Gerin et al. 2017). ALMA observations of the spectral index in the rings of the disk around HL Tau also suggests local grain growth up to centimeter sizes (ALMA Partnership 2015; Zhang et al. 2015). One of the closest observed systems with a protoplanetary disk, TW Hya, has been studied extensively by several different instruments as well. Menu et al. (2014) reviewed interferometric observations of the disk, from nearinfrared up to centimeter wavelengths, that suggest the presence of both micronsized dust and millimeter–centimetersized pebbles. van Boekel et al. (2017) presented scattered light observations of the TW Hya disk with the SpectroPolarimeter Highcontrast Exoplanet REsearch (SPHERE) instrument on the Very Large Telescope and Rapson et al. (2015) with the Gemini Planet Imager (GPI) on the Gemini South telescope; both these studies showed that the scattered light signal is dominated by dust. This indicates coexistence of solid particles of a range of sizes in young protoplanetary disks.
Even though grain growth is seen in observations, the majority of previous models involving the streaming instability focused on monodisperse particle populations. Bai & Stone (2010a) filled this gap and considered a distribution of particle sizes to study the dynamics of materials in the midplane of protoplanetary disks. They assume a discrete size distribution and saw that the radial drift velocity of all particle species is reduced compared to the one single species would have. Moreover, they show that small grains tend to move outwards in the disk. Dra̧żkowska et al. (2016) also implemented particle size dependent drift velocities into their numerical model and showed that this enhanches particle concentration. Using a fluid formalism, Laibe & Price (2014) provide analytical calculations considering multiple particle species and found the outward migration of small grains in protoplanetary disk settings as well. Johansen et al. (2007) perform numerical simulations of planetesimal formation by the gravitational collapse of locally overdense regions, where they considered multiple solid sizes. They conclude that different particle sizes collapse into the same gravitationally bound system, despite the difference in their aerodynamic properties.
Our goal with this paper is to understand the interaction between the gas and the observed pebble population in protoplanetary disks. We build on the work of Bai & Stone (2010a) and study the dynamics of a wide range of particle sizes with the inclusion of mutual drag forces between the gas and solid species. We performed simulations with various particle size ranges and size distributions. More importantly, we have studied the evolution of both discrete and continuous systems in terms of particle size distributions and compare the results. In Sect. 2 we discuss the numerical model and the dynamical equations. In Sect. 3 we measure the equilibrium drift velocities in our models and discuss the effect of including a range of particle sizes on the radial drift velocities. In Sect. 4 we address the evolution of the particle scale height throughout about 300 orbits and the influence of the particle size distribution on the degree of turbulence in the system. In Sect. 5, we analyze the radial and vertical diffusion of each solid species. In Sect. 6 we discuss potential implications of our results and finally in Sect. 7 we summarize our work.
2 Numerical model
We have modeled a local segment of protoplanetary disks in a shearing box, which corotates with Keplerian velocity, v_{K}, at an arbitrary distance, r, from the central star. Our simulations are twodimensional, in the radialvertical plane, since this is sufficient for capturing the linear and nonlinear evolution of the streaming instability (Youdin & Johansen 2007; Johansen & Youdin 2007; Bai & Stone 2010a). For simplicity, we neglect magnetic effects. Moreover, we focused on the stage before planetesimal formation begins, and hence we neglected selfgravity of solid materials.
2.1 Dynamical equations
In our models, the gas component is described as a fluid with density ρ_{g} and velocity u on a fixed grid. The gas velocity is measured relative to the Keplerian shear motion, which introduces additional advection terms in the dynamical equations (Goldreich & LyndenBell 1965). The equations that govern the motion of the gas are the continuity and momentum equations, respectively (1) (2)
The third term on the lefthand side of Eq. (1) and Eq. (2) represents the background shear flow, where Ω is the local Keplerian angular frequency. The first two terms on the righthand side of Eq. (2) account for the radial component of the stellar gravity and the centrifugal and Coriolis forces, while the third term results from the local pressure gradient. We assumed an adiabatic equation of state with an adiabatic index of 5/3 and with an initial uniform adiabatic sound speed, c_{s}, although in our simulations, c_{s} is close to being constant. The largescale pressure gradient in the disk is represented by the fourth term, where (3)
is the dimensionless parameter that sets the level of radial pressure support in the gas disk (Nakagawa et al. 1986). In Eq. (3) H_{g} is the gas scale height and P is the pressure. The secondtolast term on the righthand side of Eq. (2) represents the vertical component of the stellar gravity, while the last term is responsible for the mutual drag forces between the solid and gas components. The velocity difference between the two components is described by (u − v), where v is the particle velocity. The friction time, τ_{f}, represents the time over which the gas drag changes the velocity of a given particle by a significant amount and ϵ = ρ_{p}∕ρ_{g} is the ratio of solid and gas mass densities.
The solid component is modeled as Lagrangian superparticles, each of which represents a swarm of identical physical particles. The governing equations for each superparticle are (4) (5)
The position of the particle is x_{p} = (x_{p}, y_{p}, z_{p}), while its velocity relative to the Keplerian shear is v = (v_{x}, v_{y}, v_{z}) (Youdin & Johansen 2007).
We considered solids that have constant friction time. Given that in the usual conditions of a protoplanetary disk particle sizes are smaller than the mean free path of gas molecules, we can translate the interaction between the gas and solids to be in the Epstein regime (Weidenschilling 1977). In this case, friction time is (6)
Here, ρ_{∙} and a are the internal density and the radius of the solids. In the rest of the paper, we use the dimensionless Stokes number, τ_{s} = τ_{f}Ω. In terms of the Stokes number and the gas surface density Σ_{g}, Eq. (6) becomes (7)
We write the gas surface density as (8)
where f_{g} is a factor that measures gas depletion, with f_{g} = 1 representinga relatively pristine disk with gas surface density comparable to the minimum mass solar nebula of Hayashi (1981). We note, however, that we use a shallower power law index than the MMSN, motivated by observations of protoplanetary disks (Andrews et al. 2009; Gorti et al. 2011; Bate 2018). Using Eqs. (7) and (8), we can find the physical solid size corresponding to a given Stokes number at different regions of a protoplanetary disk as (9)
Assuming that f_{g} = 1, Eq. (9) shows that with ρ_{∙} = 1 g cm^{−3} at 2.5 au τ_{s} = 10^{−4} corresponds to solids of approximately 160 microns, while τ_{s} = 1 represents solids of about 160 cm.
2.2 Initial conditions and boundary conditions
Both the gas and solid components are initialized following a stratified density profile centered at the disk midplane. The particle scale height of all species is initially set to H_{p} = 0.015H_{g} as in Bai & Stone (2010a).
In terms of initial gas and particle velocity, we have two setups. In the case of discretized particle size distribution, we initialized both the solid and gas drift speed according to the multispecies Nakagawa–Sekiya–Hayashi (NSH) equilibrium solution (Nakagawa et al. 1986; Tanaka et al. 2005; Bai & Stone 2010a). The velocities of N particle species, which are described by and , satisfy the matrix equation (10)
where ηv_{K} represents the global radial pressure support felt by the gas (Eq. (3)). The (sub)matrices A and B are defined as A = Λ^{−1}(1 + Γ)B and . Both are functions of the Stokes number τ_{i} and the dusttogas ratio ϵ_{i} for each species i, through Λ = diag{τ_{1}, τ_{2},..., τ_{N}} and , where . It is useful to nondimensionalize the pressure support in terms of Π = ηv_{K}∕c_{s} (Bai & Stone 2010c). Bai & Stone (2010a) find that Π ≈ 0.05 describes typical pressure support in the inner regions of protoplanetary disks, which is the value we apply in our models (see also Bitsch et al. 2015). The initial gas velocity is then calculated by momentum conservation as(11)
This method works well for discrete particle species, but becomes impractical for the case of a continuous particle size distribution. In that case, the gas component is therefore initialized with a subKeplerian velocity, such that and the particle velocities are set to zero.
The boundary conditions of the simulation domain are reflecting in the vertical direction and periodic in radial direction.
2.3 Numerical method
To simulate the particlegas systems, we used the Pencil Code^{1}. It is a highorder finite difference code used for (magneto–)hydrodynamic flows in astrophysics and it uses thirdorder Runge–Kutta time integration (Brandenburg & Dobler 2002). Eqs. (1), (2), (4), and (5) are solved using explicit finite differencing. The coupling between the particle and gas components is achieved numerically using the triangular shaped cloud (TSC) method, which is a second order particlemesh scheme (Youdin & Johansen 2007).
2.4 Simulation parameters
Table 1 lists the parameters that we used in our model. In the first column, we list the names of our runs. The nomenclature SIMNABp describes some of the parameters of the given run. “SI” stands for the streaming instability, while “M” and “N” describe the minimum and maximum particle sizes such that τ_{s,min} = 10^{−M} and τ_{s,max} = 10^{−N}. Next, “A” describes the exponent, q, of the solid number density distribution,dN∕da ∝ a^{−q}. We refer to a shallow distribution, when q = 3, and a steep distribution when q = 4. The next parameter “B” describes the size of the simulation box, such that L_{x} = L_{z} = 0.BH_{g}. Lastly, “p” describes, whether the given run is continuous (c) or discrete (d) in terms of particle size distribution.
The second and third columns list the three different size ranges we used in our simulations. These Stokes number values correspond to solids ranging from about one hundred micrometers to about one meter in size (see Eq. (9)).
The fourth column in Table 1 lists the slope of the size distribution of solids in our runs. The grain number density follows dN∕da ∝ a^{−q}, where N is the solid number density and a is the solid size. Depending on physical processes considered, q can take updifferent values. In the case of selfsimilar fragmentation cascade in a debris disk, q = 3.5 (Dohnanyi 1969; Williams & Wetherill 1994). This value also agrees with the analytic estimate of size distribution of small asteroids and grains in the interstellar medium (Mathis et al. 1977). In protoplanetary disks gas also plays an important role in the dynamics of solids. As the drag from the gas slows down small particles, coagulation can also take place, which changes the steady state particle distribution. Birnstiel et al. (2011) provided recipes for grain size distribution in the enveloping cases of different collision regimes. In this paper, we study the case of q = 3 (shallow distribution) and q = 4 (steep distribution) in order to stay consistent with various possibilities.
The fifth column in Table 1 encapsulates the different box sizes that we used in our simulations. Our standard dimensions are 0.2 × 0.2 H_{g}, which is sufficient for most of the runs. In some cases, we increased this size in order to study how the box size influences the results.
In the final column we list the number of grid points N_{x} and N_{z} in the x and z directions. As default, we used N_{x} × N_{z} × N_{species} particles, where N_{species} represents the number of particle sizes. Based on Bai & Stone (2010b), this is sufficient to capture the particlegas dynamics with multisized particles. In addition, we used 128 × 128 grid cells (N_{x} × N_{z}) in the radial and vertical directions. We performed convergence tests in terms of grid and particle resolutions. In run SI4142Ad, we increased the number of particles by a factor of four from the default value. In runs SI1042Bd and SI1042Cd the number of cells was increased to 256 × 256 and 512 × 512, respectively.
The solidtogas ratio, Z, is another important factor in the outcome of the dynamics. This parameter is defined as Z = Σ_{p}∕Σ_{g}, where (12)
Here, Σ_{p} and Σ_{g} are the solid and gas column densities, respectively. In our simulations, we used a value of Z = 0.01, which does not result in efficient particle clumping (Bai & Stone 2010a; Carrera et al. 2015; Yang et al. 2017). However, forming planetesimals is not the goal of this paper. Excluding strong solid concentration and thus planetesimal formation allows us to isolatethe effect of the multiple particle sizes on the dynamics of the system.
Simulation parameters.
3 Equilibrium drift velocities and radial diffusion
The relative radial motion of the solids and the gas is driven by their mutual drag force interactions. On one hand, the particles experience radial drift due to the headwind they feel from the subKeplerian gas. On the other hand, the relative motion between particles and gas also triggers the streaming instability and generates disk turbulence. As a consequence of the selfgenerated turbulence, in addition to radial drift, the solid material also experiences radial and vertical diffusion. In this section, we describe our measurements of the radial velocity of particles once the system reaches the equilibrium state, after the vertical settling and diffusion arein balance.
3.1 Discrete particle size distribution
Following Bai & Stone (2010a) we began with a discrete particle size distribution. We first discretized particle sizes into bins of half a dex in Stokes number, τ_{s}, where each bin contained particles of the same τ_{s}. The velocities of both the particles and the gas were initialized according to the NSH solution (Nakagawa et al. 1986; Bai & Stone 2010a) given in Eqs. (10) and (11), respectively. To find the initial velocity profile of the two components at any given height, we needed to find an appropriate value for ϵ (and thus Γ). We calculated the local dusttogas ratio of each size bin as a function of height using the relationship (13)
where the gas density at the midplane is ρ_{g,0}, z is the vertical coordinate, and H_{p} and H_{g} are the particle and gas scale heights, respectively. Here, (14)
is the total mass density of particles of species j in the midplane and the mass fraction f_{j} is (15)
The Stokes number of different particle species is denoted by τ_{s,j}.
3.1.1 Steep discrete size distribution (q = 4)
First, we considered equal particle mass in each logarithmic size bin. The size distribution parameter, q, is related to the mass per logarithmic radius bin through (16)
where M is the total solid mass. This implies that q = 4 corresponds to uniform solid mass per logarithmic size bin.
Figure 1 shows the mean radial equilibrium velocity of runs SI4142d and SI3042d. The measurement was started once the particles settled to a quasisteady scale height, at which point the vertical settling and diffusion are in balance. All particle velocities in the rest of the paper are normalized by ηv_{K}, which is the difference between the gas and Keplerian velocities in the absence of particles.
The vertical lines show the 1σ variation. They represent the level of fluctuations around the mean drift velocity and correspond to the radial diffusion the given particle species experiences. The diffusion of all sizes is relatively uniform in the case of run SI4142d. In the case of run SI3042d the level of diffusion is lower for the largest particles τ_{s} = 10^{−0.5} and τ_{s} = 10^{0} than for the smaller sizes. This may be due to the fact that particles of these sizes are the least coupled to the gas and are not affected by its turbulence to a high degree. We discuss the effect of turbulence in more detail in Sect. 5.
In general, Fig. 1 shows that the radial drift velocity decreases (becomes more negative) with increasing Stokes number. Most interestingly, the smallest particles have positive mean drift velocities, meaning that they move outwards. We compared the measured velocities in Fig. 1 with the theoretical single species solution (red stars) using the NSH solution given by Nakagawa et al. (1986) (17)
We see that the radial drift velocities in our multispecies runs are reduced, especially for the larger Stokes numbers. The deviation of our results compared to the single species NSH solution is due to the fact that different particle species interact with each other through the gas. As large particles drift in, the gas moves out due to momentum conservation. Because small particles are more tightly coupled to the gas, they are entrained in its motion, and contribute to the acceleration of the gas to closer to the Keplerian speed. Thus, larger particles drift in at reduced velocities compared to the single species case, since they feel a weaker headwind from the gas. At the same time, the particles with τ_{s} ≳ 10^{−2} trigger the streaming instability most readily and create turbulence in the surrounding gas, which diffuses the solid materials (Bai & Stone 2010a).
We compared our measurements with the theoretical multispecies solution (see Eq. (10)) using two methods. The first method gives the radial equilibrium velocity near the midplane. Here, we have used the midplane dusttogas ratio of each species calculated under the assumption that the vertical particle distribution is approximately Gaussian. The midplane density for the particles is (18)
where Σ_{p} is the column mass density of the particles. Thus, the local dusttogas ratio for the jth particle bin at the midplane is (19)
where is the average particle scale height in the jth bin and ρ_{g,0} is the gas density at the midplane. We were able to calculate the midplane dusttogas ratio in a given size bin using the global dusttogas ratio and the mass fraction per bin (see Eq. 15), such that (20)
As seen in Fig. 1, the theoretical radial equilibrium velocity near the midplane agrees well in the case of run SI4142d (green circles), but is noticeably more positive than our measured drift velocities in the case of run SI3042d (blue circles). This is likely to be the result of the fact that we used the midplane dusttogas ratio, which is an overestimate of the solid loading, which is further away from the midplane.
In the second method of calculating the analytic multispecies solution, we used the heightaveraged radial equilibrium velocity. In this case, we calculated the drift velocity of each species as a function of height assuming particles of each species are uniformly distributed within its own scale height. Then the dusttogas ratio for species j is given by (21)
where is the timeaveraged scale height of the particle species. Figure 2 shows the equlibrium velocity of the τ_{s} = 10^{−3} species as a function of height in run SI3042d. The drift velocity is largest in two bins around the midplane and is close to the gas velocity in the absence of particles at larger heights. We then averaged the radial drift velocity over height by weighing with the dusttogas ratio in each bin (Bai & Stone 2010a). As seen in Fig. 1 (red and purples squares), the heightaveraged radial equilibrium velocity provides a significantly better match to the results of our simulations.
We conducted convergence tests to ensure that numerical resolution does not influence our results. Figure 3 shows the radial equilibrium velocity of particles in runs SI1042d, SI1042Bd and SI1042Cd. In all three cases, the simulation box size was L_{x} = L_{z} = 0.2H_{g}, but the number of grid cells was successively increased from 128 × 128 to 256 × 256 and finally to 512 × 512. The equilibrium velocity agrees in all three resolutions. Thus, it appears the grid resolution does not have a significant effect on our results.
We also increased the total particle number by a factor of four, to check for convergence. Figure 4 shows that indeed, the equilibrium radial velocity is not appreciably affected by an increase in particle number.
Fig. 1 Mean equilibrium radial drift velocity as a function of Stokes number (τ_{s}) for runs SI4142d and SI3042d measured from the time of saturation. The vertical bars show the 1σ variation of the radial velocity and are measures of the radial diffusion each species experiences. The red stars represent the single species equilibrium solutions. The green and blue circles show the multispecies radial equilibrium velocity near the midplane. The purple and red squares show the height averaged multispecies solution under the assumption that the particle density of each species is distributed uniformly within its scale height. Run SI3042d was shifted horizontally for illustration purposes. 

Open with DEXTER 
Fig. 2 Theoretical equilibrium radial drift velocity of species τ_{s} = 10^{−3} in run SI3042d with respect to height. The drift velocity is highest in the two bins around the midplane and is close to the gas velocity in the absence of large particles further away from the midplane. 

Open with DEXTER 
Fig. 3 Convergence of equilibrium radial velocity of different particle sizes against grid resolution with runs SI1042d, SI1042Bd and SI1042Cd. Both the mean drift velocities (circles) and the level of diffusion (vertical bars) agree between the cases of N_{x} × N_{z} = 128 × 128, 256 × 256, and 512× 512. This impliesthat our results are not dependent on numerical resolution. Runs SI1042Bd and SI1042Cd were shifted horizontally for better visibility. 

Open with DEXTER 
Fig. 4 Convergence of equilibrium radial velocity of different particle sizes against particle number with runs SI4142d and SI4142Ad. The mean drift velocities and the 1σ variation agree well, indicating that the number of particles does not influence our results. In the case of run SI4142Ad, the symbols were shifted horizontally for comparison purposes. 

Open with DEXTER 
3.1.2 Shallow discrete size distribution (q = 3)
We also considered a size distribution with q = 3, with which dM∕dln a ∝ a. Compared to the previous case, we now had more mass in particles with larger Stokes numbers, thus there were more particles available to play an active role in the streaming instability. As before, to initialize the system we used Eqs. (10) and (11).
In Fig. 5, we show the radial equilibrium velocity of runs SI3032d and SI4132d. Similar to the previous case, where all size bins had the same solid mass (runs SI3042d and SI4142d) presented in Fig. 1, all particle species have unique drift velocities. The small particles have positive drift velocities and that of the larger species is reduced relative to the single species solution (marked as red stars). The 1σ bars around the mean drift velocity represent the level of diffusion each species experiences. The radial diffusion is smallest for particles with τ_{s} > 10^{−0.5}, since these are the least coupled to the gas and therefore the least influenced by its turbulence.
In Fig. 5 we have plotted again the theoretical multispecies solution using the two methods presented in the previous section. The green and blue circles correspond to the multispecies radial equilibrium velocity near the midplane. This method produces velocity values that agree well with our measurements in the case of run SI4132d, but deviate in the case of run SI3032d. The reason for the deviation is, again, that using the midplane dusttogas ratio (see Eq. (20)) overestimates the solid loading further away from the midplane thus the drift velocity of the smaller species as well.
We also compare the theoretical multispecies solution using the second method. The red and purple squares in Fig. 5 correspond to the heightaveraged solution, where we assumed that the particle density of each species is distributed uniformly within its own scale height and used Eq. (21) to calculate the dusttogas ratio as a function of height. We then weighed the velocities by the particle density in each height bin. Figure 5 shows that this method provides a significantly better match to our numerical measurements, especially in the case of the smaller species.
Fig. 5 Mean equilibrium radial drift velocity as a function of τ_{s} for runs SI4132d and SI3032d measured from the time of saturation. The color scheme agrees with that of Fig. 1. Run SI3032d was shifted horizontally for illustration purposes. 

Open with DEXTER 
3.2 Continuous particle size distribution
So far we have considered discretized particle size distribution. We now turn to a more realistic case, in which each particle in a given run has a unique size, in other words particle sizes are distributed in a quasicontinuous way. In this case, particle velocities are initialized to be zero and the gas component is set to have an initial subKeplerian velocity of u_{y} = −Πc_{s}, where Π = 0.05.
3.2.1 Steep continuous size distribution (q = 4)
We revisitedthe case of q = 4, such that dM∕dln a is constant, and studied the same systems as in Sect. 3.1.1, but using a quasicontinuous particle size distribution.Here each particle swarm was given a unique mass (per unit azimuthal length) according to the expression (22)
where M_{p,j} is the total mass of a given particle species and N_{p,j} is the number of particles of a given species.
The resulting radial drift velocity of a random subset of all particles is shown in Fig. 6a. In contrast to Fig. 1, now each particle has a unique velocity distributed around the mean shown with red. The level of diffusion is indicated by the scattering of the individual (black) drift velocity points around the mean. Similarly to the discrete case shown with blue in Fig. 1, the level of diffusion is smallest for the largest particles, since these are the least coupled to the gas (see also Youdin & Lithwick 2007).
In Fig. 6a, the blue circles represent the theoretical velocity values of all particles grouped into 16 bins calculated using the multispecies NSH solution (Eq. (10)). For this, we used the midplane dusttogas ratio in the given size bin calculated under the assumption that the vertical particle distribution is approximately Gaussian. The dusttogas ratio in the midplane of each size bin is calculated using Eq. (20). Since we used the midplane equilibrium radial velocity as comparison to our measurements, the theoretical velocities somewhat overestimate the numerical results, especially in the case of run SI3042c.
Figure 6b shows the case of the smaller particle sizes, namely run SI4142c. Compared to the result in Fig. 6a the particles are more concentrated around the mean velocity. This means that the level of diffusion is lower on average, indicating that the particles experience less turbulence. This is due to the fact that now, the largest particles have size τ_{s} ≤ 10^{−1} and more tightly coupled particles generate lower levels of turbulence. Thus the number of particles that are large enough to drive the streaming instability is lower compared to the case of run SI3042c (Bai & Stone 2010a).
3.2.2 Shallow continuous size distribution (q = 3)
We also considered the case of the shallow particle size distribution. Figure 6c shows a randomly selected sample of particles corresponding to run SI3032c. The mean velocity is positive for the smaller particles. Compared to the same run with q = 4 (Fig. 6a) the diffusion of the smaller particles has increased.
We compared this result to the case with smaller Stokes numbers, namely run SI4132c. The drift velocities in Fig. 6d are less dispersed around the mean, contrary to the case of the larger particles in Fig. 6c. Also, both the theoretical and the mean velocities for larger sizes are smaller. As in the case of the steep size distribution, the theoretical midplane velocities overestimate somewhat our measurements, especially in the case of run SI3032c.
Comparing Fig. 6a with Fig. 6c and Fig. 6b with Fig. 6d, it is clear that the use of either q = 3 or q = 4 does not have a large effect on the radial drift velocity. The level of turbulence appears to be higher if the particle size distribution is shallow, as shown by the larger dispersion in Figs. 6c and 6d.
Fig. 6 Radial drift velocity of runs SI3042c, SI4142c, SI3032c, and SI4132c in terms of Stokes number (τ_{s}) at t ≈ 1000 Ω^{−1}. Individual black points show the radial drift velocities of a random subset of approximately 14 000 particles. The particle size distribution is continuous, such that each particle has a unique size and drift velocity. Each red cross shows the mean drift velocity of two hundred particles with similar τ_{s}, which decreases with increasing particle size. The dispersion of points around this mean shows the level of diffusion the particles experience, which goes down with increasing particle size as well. The purple bars show the 1σ variation around the mean velocity. The blue circles show the radial equilibrium velocity near the midplane for each particle species grouped into 16 bins using the multispecies NSH solution in Eq. (10). The theoretical velocity values provide a relatively good fit to the mean drift velocity for small particle sizes but are somewhat off for large Stokes numbers. 

Open with DEXTER 
3.3 Comparison of particle size distributions
In order to understand how parameters such as the exponent of the size distribution (q) and the manner in which the particles are distributed (discrete or continuous) influence the level of turbulence, we compared the 1σ variation around the mean equilibrium radial drift velocities. We do not expect the way the particles are distributed (discretely or continuously), to influence the results significantly. The only difference that stems from the two distribution methods is that in the latter case, each particle has a unique size between τ_{s,min} and τ_{s,max}, whereas in the former case, the particle sizes are binned into discrete sizes.
3.3.1 Comparison of steep and shallow particle size distributions
In Fig. 7, we present the mean equilibrium drift velocity and the 1σ limits in the case of run SI3042d and SI3032d. As the figure shows, there is no significant difference between the fluctuations around the mean velocities whether the particle size distribution is steep (q = 4) or shallow (q = 3). We see a similar tendency in Fig. 8, where we compare runs SI4142d and SI4132dand show that the 1σ variation around the mean equilibrium drift velocities is similar for both particle size distributions.
3.3.2 Comparison of discrete and continuous particle size distributions
Figure 9 shows the mean drift velocities and the 1σ variation around them in the case of run SI3042d and SI3042c. The particle sizes are binned into 16 bins for the continuous run for illustration purposes and the velocities are calculated over 15 orbits. Since the particle sizes in the two runs do not agree completely, the comparison is not entirely straightforward. However, there is no significant difference between the two models in Fig. 9.
In Fig. 10, however, there is some difference between the 1σ limits of runs SI3032d and SI3032c. The run where the particles are distributed continuously shows about a factor of two stronger turbulence than the discrete case.
4 Particle scale height
The mutual drag between the solid and gas components creates turbulence that stirs the particles (Johansen & Youdin 2007). Once the vertical settling and the diffusion are in balance, the system reaches a saturated state, and the particles reach a quasisteadyscale height (Johansen et al. 2009; Bai & Stone 2010a; Yang & Johansen 2014). By measuring the scale height to which each particle species saturates, we can indirectly assess the level of diffusion that they experience and thus the strength of turbulence driven by the streaming instability.
We calculated the scale height of each species as the square root of the particle height variance, , thus (23)
4.1 Discrete size distribution
Figures 11 and 12 show the evolution of the particle scale height in the case of discrete steep particle size distribution. Figure 11 shows the system where the participating particles are τ_{s} = 10^{−3} − 10^{0}, with each particle species differentiated by a unique color. The system is first modeled in a 0.2 × 0.2 H_{g} box, shown on the left. Each particle species saturates to a different scale height in such a way that H_{p} increases with decreasing Stokes number. After saturation, the scale heights show some fluctuations, especially for the smaller species. We ran the same model in a larger box of 0.4 × 0.4 H_{g}, and the result is shown on the right panel in Fig. 11. This was done to check that the turbulent diffusion of the smaller solids is not limited by the size of the simulation box. The saturated scale height of the smallest species is indeed larger for the larger box. This indicates that the vertical box size prevented these solids from reaching the scale heights as they should. On the other hand, the larger species saturate to about the same scale height, indicating convergence.
We compared this to the case with the smaller participating particle species, namely τ_{s} = 10^{−4}−10^{−1}. The result is shown in Fig. 12. Given that the particles are smaller, it takes longer for the system to settle into equilibrium. This is due to the larger settling time (Nakagawa et al. 1986; Dubrulle 1995; Youdin & Lithwick 2007) (24)
As above, we increased the box size to ensure that the vertical diffusion of the particles is not limited by the height of the simulation domain. The panel in the middle of Fig. 12 shows that the scale height of the smallest species increased and that of the larger species decreased compared to the case of the smaller box. In addition to increasing the box size, we conducted another convergence test, in which we increased the total number of particles by a factor of four. The result is shown in the rightmost panel of Fig. 12. None of the species settle to a constant scale height within approximately 300 orbits.
Next, we considered the case of shallow and discrete particle size distribution with a powerlaw index of q = 3. In other words, the system now has more mass in larger particle species. The evolution of both SI3032d and SI4132d is shown in Figs. 13 and 14. Comparing the corresponding panel in Figs. 11 and 12, we see that the particles in the case of both q = 3 and q = 4 reach similar scale heights. The particle scale heights in Fig. 14 and the first panel in Fig. 12 are similar as well. This implies that the level of turbulence is similar both when the particle distribution is steep and shallow, which is consistent with the similar dispersions of the radial drift velocity of particles in the respective models found in Sect. 3.1.
Fig. 7 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τ_{s}) for runs SI3042d and SI3032d measured from the time when the saturated state is reached. The level of diffusion is similar in both runs. 

Open with DEXTER 
Fig. 8 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τ_{s}) for runs SI4142d and SI4132d measured from the time when the saturated state is reached. The level of turbulence does not appear to be affected by the choice of the q parameter. 

Open with DEXTER 
Fig. 9 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τ_{s}) for runs SI3042d and SI3042c measured from the time when the saturated state is reached. In the case of run SI3042d, the velocities are averaged to the end of the simulation, while in the case of run SI3042c, they are averaged over 15 orbits. 

Open with DEXTER 
Fig. 10 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τ_{s}) for runs SI3032d and SI3032c measured from the time when the saturated state is reached. In the case of run SI3032d, the velocities are averaged to the end of the simulation, while in the case of run SI3032c, they are averaged over 15 orbits. 

Open with DEXTER 
Fig. 11 Particle scale height evolution of runs SI3042d and SI3044d for about 300 orbits. Different particle sizes (marked with different colors) saturate to different scale heights, in such a way that H_{p} increases with decreasing Stokes number. Left panel: evolution of the run with the larger Stokes numbers (τ_{s} = 10^{−3}−10^{0}) in a 0.2 × 0.2 H_{g} simulation box. Right panel: convergence of the scale height for the larger species in a box that is twice the size of the original. The scaleheight is relatively similar for the two box sizes, except for the two smallest particle size bins that exhibit an increased scale height when the box size is enlarged. 

Open with DEXTER 
Fig. 12 Particle scale height evolution of runs SI4142d, SI4144d, and SI4142Ad, where the participating particles are τ_{s} = 10^{−4}−10^{−1}. Left panel: simulation in a box of 0.2 × 0.2 H_{g}. Middle panel: evolution of the same system in a 0.4 × 0.4 H_{g} box, where the particle scale heights saturate to different values relative to the smaller box. Right panel: the total particle number was increased by a factor of four from the model in the leftmost panel. In runs SI4144d and SI4142Ad the particle subdisks do not saturate to a constant height within the simulation time. 

Open with DEXTER 
Fig. 13 Particle scale height evolution of run SI3032d in a 0.2 × 0.2 H_{g} simulation box. Compared to the lefthand side of Fig. 11, the particle species saturate to similar scale heights but do not show fluctuations of the same magnitude with time. 

Open with DEXTER 
Fig. 14 Particle scale height evolution of run SI4132d in a 0.2 × 0.2 H_{g} simulation box. The particle species do not reach a constant scale height within the simulation time. Compared to the plot on the lefthand side of Fig. 12, where the powerlaw index of the distribution q = 4, the particle layers extend to similar scale heights. 

Open with DEXTER 
4.2 Continuous size distribution
We turn to studying the evolution of the continuously distributed particle sizes. Since in this case each particle has a unique size, we have grouped them into six bins for illustration purposes in the following discussions.
In Figs. 15 and 16, we show the scale height evolution of runs SI3042c and SI4142c. Similar to the previously presented cases, where the particles are distributed in a discrete way in terms of their size, the large species remain closer to the midplane, while the selfgenerated turbulence diffuses small particles to larger heights. Once again, we conducted convergence tests in terms of simulation box size. As shown in Fig. 15, in the case of size range τ_{s} = 10^{−3}−10^{0}, the smallest species saturate to a larger H_{p} in the 0.4 × 0.4 H_{g} box than in the smaller box, while the larger solids converge to the same scale heights. Figure 16 shows the evolution of the particle scale height, when the particles are τ_{s} = 10^{−4}−10^{−1}. The particle species do not settle to a constant height during the simulation time in either the 0.2 × 0.2 H_{g} box, nor the one double that size, although the smaller species reach large heights in the latter case.
Finally, in Figs. 17 and 18 we show the evolution of the system with shallow continuous particle size distribution. Figure 17 shows the scale height evolution of runs SI3032c, SI3034c and SI3038c, which were performed in boxes of 0.2 × 0.2 H_{g}, 0.4 × 0.4 H_{g}, and 0.8 × 0.8 H_{g}, respectively. We see that in run SI3034c the four smallest particle bins reach heights almost a factor of two larger compared to the run in the smaller box. To test the convergence of the system, we increased the box size again, to 0.8 × 0.8 H_{g}. The scale heights of the three largest species reach a lower scale height than before, while the smaller species are diffused to larger heights. This unusual behavior is accompanied by the formation of shocklike patterns in the radial component of the gas velocity. We discuss this in more detail in Appendix A. Figure 18 shows the evolution of the same model with smaller particles. Since all species settle to a relatively constant height, we ended the run after about 140 orbits. However, when doubling the domain size to 0.4 × 0.4 H_{g} the smallest species bins saturate at heights almost twice as large compared to the smaller box.
Comparing Figs. 17 and 18 with Figs. 15 and 16, respectively, we see that in the latter case where the size distribution is steep (q = 4), all species reach lower scale heights in general. The same trend is present once we increase the box size from our default size. This is likely due tothe difference in the number density of the solids. By increasing the mass in larger particle species, we increase the availability of particles that actively participate in the streaming instability which generates turbulence. Thus, particles in runs with continuous and shallow size distribution experience more diffusion and the particle layers become more puffed up around the midplane.
Fig. 15 Particle scale height evolution of run SI3042c and run SI3044c. Each particle in the simulation has a unique size within the given size range. Particles are grouped into six size bins for illustration purposes. To ensure that the diffusion of the small particles is not limited by the simulation box, we performed the same simulation in a 0.4 × 0.4H_{g} box in the right panel. The scale heights of the large particles converge while the smallest species saturate to a larger H_{p} than in the left panel. 

Open with DEXTER 
Fig. 16 Particle scale height evolution of run SI4142c and run SI4144c. The scale height of all species does not saturate within the simulation time in either the 0.2 × 0.2 H_{g} (left panel), nor the 0.4 × 0.4 H_{g} box (right panel). 

Open with DEXTER 
Fig. 17 Particle scale height evolution of run SI3032c, run SI303 4c, and run SI3038c. Panels (from left to right): evolution of the system in simulation boxes of 0.2 × 0.2 H_{g}, 0.4 × 0.4 H_{g}, and 0.8 × 0.8 H_{g}, respectively. The largest particles converge to the same height in the first two cases, but particles with τ_{s} < 10^{−1} are diffused to heights larger by about a factor of two. The panel on the right, where the simulation box size is the largest, shows unusual nonconvergent behavior due to the development of shocklike patterns in the radial component of the gas velocity (discussed in Appendix A). 

Open with DEXTER 
Fig. 18 Particle scale height evolution of runs SI4132c and SI4134c. The particle scale height of the small species increases by a factor of two when changing the dimensions from 0.2 × 0.2 H_{g} in the left panel to 0.4 × 0.4 H_{g} in the right panel, while that of the largest particles converges to about the same height. 

Open with DEXTER 
4.3 Effect of multiple particle species
The effect of multiple particle sizes in the same system is evident from the change in mean radial drift velocity relative to the single species case and the effect of the selfregulated diffusion driven by the larger species. Compared to singlespecies simulations, there is also a significant difference in terms of the scale height each particle species saturates to. Carrera et al. (2015) conducted singlespecies streaming instability simulations of a variety of Stokes numbers. They found that particles of 10^{−3} ≤ τ_{s} ≤ 10^{0} saturate to H_{p} ≈ 10^{−2} H_{g} (see also Yang & Johansen 2014; Yang et al. 2017). In general, we see that most particle species, especially those with τ_{s} ≤ 10^{−0.5}, are diffused to heights much larger than that. This is the consequence of the strong turbulence driven by the large particles. The vertical velocity fluctuations in the gas are a measure of the level of turbulence driven by the particles. Figure 19 shows the rootmeansquare of the vertical gas velocity in terms of height for the runs in three different box sizes, namely run SI3032c, SI3034c, and SI3038c. The level of turbulence is not only significant near the midplane (z = 0), but also at larger heights. Thus, the smaller particle species that reside at larger heights experience the effect of the stirring from the large particles near the midplane. In Fig. 20 we show the horizontal component of the gas velocity fluctuations as a function of height. The curves in the case of runs SI3034c and SI3038c peak around the midplane and drop to an approximately uniform value beyond ± 0.1 H_{g} until the increase near the walls. In the case of run SI3032c, the fluctuation is at a minimum around the midplane and increases toward the walls of the simulation domain.
In Figs. 21 and 22 we show the vertical and horizontal velocity fluctuations of runs SI4132c and SI4134c. Both δu_{z,rms} and δu_{x,rms} show similarbehavior in terms of box size. The vertical component shows uniformity with respect to height, while the horizontal component peaks near the midplane and decreases toward the walls.
Fig. 19 Vertical gas velocity fluctuation with respect to height in the case of runs SI3038c, SI3034c, and SI3032c averaged over approximately 15 orbits starting at t ≈ 1000 Ω^{−1}. The relatively uniform level of fluctuation further away from the midplane implies that the turbulence is not limited to the region close to the midplane (z = 0), where the large particles reside. The turbulence extends to larger scale heights, thus affecting the small particle species that reside there. 

Open with DEXTER 
Fig. 20 Horizontal gas velocity fluctuation with respect to height in the case of runs SI3038c, SI3034c, and SI3032c averaged over approximately 15 orbits starting at t ≈ 1000 Ω^{−1}. The level of fluctuation peaks around the midplane and is relatively uniform throughout the rest of the domain for the runs in the 0.4 × 0.4 H_{g} and 0.8 × 0.8 H_{g} boxes, while in the smallest box δu_{x,rms} varies by more than a factor of two. 

Open with DEXTER 
Fig. 21 Vertical gas velocity fluctuation with respect to height in the case of runs SI4134c, and SI4132c averaged over approximately 15 orbits starting at t ≈ 500 Ω^{−1}. The level of fluctuation is relatively uniform for the two runs both in terms of height and box size. 

Open with DEXTER 
Fig. 22 Horizontal gas velocity fluctuation with respect to height in the case of runs SI4134c and SI4132c averaged over approximately 15 orbits starting at t ≈ 500 Ω^{−1}. The level offluctuation peaks around the midplane for both box sizes and is relatively uniform throughout the rest of the domain. 

Open with DEXTER 
5 Radial and vertical diffusion
The mutual drag between the solid and gas components triggers the streaming instability and generates turbulence. This selfregulated turbulence in return stirs up the particles which experience diffusion in both the radial and vertical dimensions. To understand the influence of each particle species on the level of turbulence and their response to it, we measured both the radial and vertical diffusion coefficients for each particle size. We did this only for the runs where the particles were distributed in a discrete fashion in terms of their size. The continuous runs cannot be used to measure the turbulent diffusion, since in that case, each particle has a unique size and calculating the diffusion of each species with respect to the rest would be meaningless due to the differential radial drift.
The radial diffusion coefficient, D_{x}, represents the degree of random walk each species experiences due to the selfregulated turbulence. It is calculated following Johansen & Youdin (2007) as (25)
Here, is the varianceof the radial displacement that each species covers at discrete time intervals. We applied linear regression on the resulting curve to find the radial diffusion coefficient of each species.
The vertical diffusion coefficient, D_{z}, is calculatedusing the scale height of the particles, such that (Dubrulle 1995; Youdin & Lithwick 2007) (26)
Here, δ is a dimensionless measure of the diffusioncoefficient, such that D_{z} = δc_{s} H_{g} (Johansen & Klahr 2005; Johansen et al. 2014a).
5.1 Steep size distribution
We first considered the case of steep particle size distribution, with a powerlaw index of q = 4. In Fig. 23 we present the relationship between diffusion coefficient and Stokes number and compare the dependence of D_{x} and D_{z} on particle size. As seen in Fig. 23, D_{z} increases with increasing Stokes number in both run SI4142d (blue dashed curve) and SI3042d (red dashed curve) for τ_{s} ≤ 10^{−1.5}. This trend arises, since as we discussed in Sect. 4, smaller particles saturate to larger scale heights. As a consequence, smaller species can reside in regions where the level of turbulence is somewhat lower than near the midplane, as seen in Fig. 19. The decreasing diffusion coefficients of the smaller species could also be influenced by the box size, such that the random walk of the particles close to the domain walls is restricted. On Fig. 23 we marked these values with crosses. When τ_{s} ≥ 10^{−1.5}, the vertical diffusion coefficient decreases with increasing Stokes number. These species saturate at small scale heights and hence feel the full strength of the turbulent diffusion in the midplane layer.
The radial diffusion coefficient, D_{x}, for run SI3042d and SI4142d is relatively constant for τ_{s} ≤ 10^{−2}, since as Figs. 11 and 12 show, particles of these sizes saturate to similar scale heights. Above τ_{s} ≈ 10^{−2.5}, D_{x} increases with increasing Stokes number, since in contrast with the smaller species, they reside closer to the midplane, where they can experience more turbulence. This behavior can be seen in both Figs. 20 and 22. Comparing the curves that correspond to run SI4142d (blue curves) and SI3042d (red curves), we see that for all particle species both the vertical and radial coefficients are higher in the latter case. This is the consequence of the presence of more particles with τ_{s} ≳ 10^{−2}, which makes it easier for the streaming instability to operate and drive the turbulence (Bai & Stone 2010a). In other words, the larger the abundance of the actively participating particles is, the stronger the streaming turbulence becomes.
As shown in Fig. 12, the vertical diffusion of the smallest particle species is limited by the vertical domain, since their scale height increases as we double the box size. We measured D_{x} and D_{z} in the 0.4 × 0.4 H_{g} box, and show the result in Fig. 24. Compared to the case of the smaller simulation box shown in Fig. 23, there is relatively little change in the radial diffusion coefficients for both runs, which implies convergence. The vertical diffusion coefficient D_{z,41} (blue dashed curves) in both Figs. 23 and 24 take on similar values for each particle species. On the other hand, D_{z,30} (red dashedcurves) takes on slightly different values in both boxes, especially for 10^{−2} < τ_{s} < 10^{−0.5}. This could also be a consequence of the fluctuations of the particle scale height in the smaller box, as see in Fig. 11. We marked the diffusion coefficients that correspond to species which do not saturate within the simulation time with crosses.
In Fig. 25, we present the result of a numerical convergence test performed on run SI4142d. We increased the number of grid cells in the radial and vertical dimensions from our default value of 128 × 128 to 256 × 256, while at the same time leaving the physical box size unchanged, such that L_{x} = L_{z} = 0.2 H_{g}. Figure 25 shows convergence may have been reached against numerical resolution, suggesting that numerical resolution does not appreciably influence the level of turbulence felt by the particles.
Fig. 23 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI4142d and SI3042d. Here, D_{x} is approximately constant for small Stokes numbers and increases with increasing particle size as τ_{s} ≥ 10^{−2.5} in the case of run SI4142d. The vertical diffusion coefficient increases with particle size until τ_{s} ≈ 10^{−1.5}. Above this value it decreases as the particles are less coupled to the gas. Both the radial and vertical diffusion coefficients are higher in the case of run SI3042d, since this model is more abundant in solids with τ_{s} ≳ 10^{−2} which are responsible for triggering the streaming instability. The diffusion coefficients corresponding to species that saturate to scale heights similar to the maximum allowed scale height in the box are marked with a cross. 

Open with DEXTER 
Fig. 24 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI4144d and SI3044d. We doubled the box size to 0.4 × 0.4 H_{g} compared to the simulations in Fig. 23. Compared to the smaller box, both D_{x} and D_{z} converge to approximately the same values for each Stokes number. The exception seen in curve D_{z,30} at intermediate Stokes numbers is likely a consequence of the fluctuations seen in their particle scale height in the smaller box (see Fig. 11). The diffusion coefficients corresponding to species that do not saturate within the simulation time are marked with a cross. 

Open with DEXTER 
Fig. 25 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI4142d and SI4142bd. Both D_{x} (solid curves) and D_{z} (dashed curves) show relatively small deviation. 

Open with DEXTER 
5.2 Shallow size distribution
We also measured the diffusion coefficients in the case of shallow particle size distribution. Figure 26 shows D_{x} and D_{z} of runs SI4132d and SI3032d. Comparing the diffusion coefficients to the ones in Fig. 23, we see that there is not a significant increase in either D_{x}, nor D_{z}. Both Figs. 13 and 14 show the fluctuations of the equilibrium scale heights, which according to Eq. (26), are responsible for setting D_{z}. The diffusion coefficients which correspond to species that either saturate to scale heights similar to the maximum allowed scale height in the box or do not saturate within the simulation time are marked with crosses.
Fig. 26 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI4132d and SI3032d. Both D_{x} and D_{z} take on similar values as in Fig. 23, where the particle size distribution is shallow (q = 4). The diffusion coefficients corresponding to species that saturate to scale heights similar to the maximumallowed scale height in the box or do not saturate within the simulation time are indicated with a cross. 

Open with DEXTER 
6 Implications
Particle scale height is a key parameter in grain growth by coagulation (Zsom et al. 2011; Dra̧żkowska et al. 2013). As shown in Sects. 4 and 5, every particle species saturates to a unique scale height and thus has a response to the selfregulated turbulence. Since coagulation is driven by the relative velocity of the solid materials, as well as their local number density, its efficiency is strongly influenced by the equilibrium particle scale height of the given particle species.
Here we apply the main results presented in the paper to protoplanetary disks by showing the evolution of particle radial position, r, as a function of time, t, taking into account both radial drift and diffusion over approximately 300 orbits.
Figure 27 shows the radial distance traveled by the particles in run SI3032c. In this run, the size distribution of the particles is continuous thus we grouped them into six bins for illustration purposes. The solid lines show the mean position of the particles in each bin, while the dashed lines correspond to the 1σ variation. Each particle size bin is labeled by different colors. We see that, as expected, different particle bins experience varying amounts of radial motion depending on their size. As discussed in Sect. 4, this is partly determined by their vertical position with respect to the midplane. Both radial drift and differential radial drift are strongest for the largest species. The position of the small species increases with time, implying outward transport in the disk. As the particles are distributed continuously in terms of their size, differential radial drift dominates over diffusion in spreading the distribution of particle positions.
Figure 28 shows the evolution of the particle radial positions for run SI3032d. Particles with τ_{s} > 10^{−2} drift inwards, while species with τ_{s} ≤ 10^{−2} move outwards. Compared to Fig. 27, the level of spreading for all species is significantly lower, since differential drift does not play a role.
Overall, our results imply that differential radial drift contributes more to the spreading of particle species than diffusion. However, as is evident in Fig. 27, differential drift does not mix the particle species. In fact, diffusion must be present to achieve mixing between particles of different sizes.
The dynamics of solids driven by the streaming instability, once the realistic case of multiple particle species is considered, may be responsible for the transport of chondrules as well as matrix material in the early solar system. The chondrulematrix complementarity in carbonaceous chondrites implies that both components formed in the same region of the disk and did not drift apart during transport to the chondrite forming region (Hezel & Palme 2010; Jacquet 2014; Johansen et al. 2014b). Our results here show that differential radial drift is generally a more powerful mechanism than diffusion for separating particle species of different sizes.
7 Summary
In this paper we built on the work of Bai & Stone (2010a) and modeled the dynamics of multiple particle species embedded in gas. To do this, we used the Pencil Code and performed 2D fluidparticle simulations that span the radial and vertical plane. We studied three scenarios in terms of particle size, such that the Stokes number of the solid materials in a given model was τ_{s} = 10^{−4}−10^{−1}, τ_{s} = 10^{−3}−10^{0}, or τ_{s} = 10^{−1}−10^{0}. The size distribution of the particles was distributed according to dN∕da ∝ a^{−q}, where N is the total number density and a is the radius of the particles. We considered q = 3 (shallow) and q = 4 (steep) so that our assumed size distributions envelope a variety of distributions predicted by observations and numerical simulations. On one hand, we distributed the particles into discrete size bins, in order to compare with Bai & Stone (2010a). On the other hand, we studied a more realistic case and modeled systems where the particles were distributed continuously in terms of their size.
We show that the dynamics of particles in protoplanetary disks changes once we consider the realistic case of multiple solid sizes embedded in the gas.
Firstly, due to the friction from the subKeplerian gas, the particles experience radial drift. At the same time, the large particles trigger the streaming instability and the generated turbulence drives the radial and vertical diffusion of the solid materials.
Secondly and most interestingly, in Figs. 27 and 28 we see that small particles move outwards as also observed in Bai & Stone (2010a). This behavior is not seen in previous streaming instability simulations that contain particles of a single species for example, Johansen & Youdin (2007). As the larger species stir up the gas, the gas is accelerated locally and moves outwards, taking the strongly coupled (small) dust particles with it. At the same time, larger species move inwards but at drift velocities slower than expected from the single species NSH solution (e.g., Fig. 1).
Moreover, from the comparison of models with different size distribution exponents (q), particle distribution methods (discrete or continuous) and particle sizes, we see that the sizes of the participating particles is the most important factor. The strength of turbulence appears to be independent of the steepness of the size distribution (see Figs. 7 and 8). Whether the particles are discretely or continuously distributed also produces no large differences in the measured turbulence levels (see Figs. 9 and 10).
Finally, as indicated by the relatively uniform level of the vertical component of the gas rootmeansquare velocity away from the midplane (see Figs. 19 and 21), the turbulence generated by the large particles located close to the midplane is not limited to the midplane. The turbulence extends to larger heights and as a result the smaller particle species that reside at these heights also experience stirring generated near the midplane. As a consequence, compared to single species streaming instability simulations (Carrera et al. 2015), the measured particle scale heights in our models are larger by several factors.
The results listed above have important implications for protoplanetary disks. As shown in Fig. 27, in the more realistic case of continuous particle size distribution, both diffusion and differential radial drift contribute to the spreading of particles over time. Contrary to differential drift, diffusion is driven by the selfgenerated turbulence and is present in all our systems, independent of the particle size distribution we consider. The thickness of the dust layer in protoplanetary disks is an important factor that affects the efficiency of coagulation into dust aggregates (Zsom et al. 2011; Dra̧żkowska et al. 2016). In the future, interaction between particle species should be taken into account in planet formation models.
Fig. 27 Evolution of particle radial positions with respect to time in the case of run SI3032c. Since the particles are distributed continuously in terms of their size, we grouped them into six particle size bins. The mean position of the particles in a given size range is marked by the solid curves, while the dashed curves mark the spread of particle positions by differential radial drift and diffusion. Particles with τ_{s} ≥ 10^{−1} drift and diffuse inwards, while particles of τ_{s} ≤ 10^{−1.5} tend to move outwards in the disk. The spreading of each particle species is in fact dominated by the differential radial drift insteadof diffusion. 

Open with DEXTER 
Fig. 28 Evolution of particle radial positions with respect to time in the case of run SI3032d. The smallest species drift outwards, while species with sizes τ_{s} > 10^{−2} drift inwards. Compared to the run with continuous size distribution shown in Fig. 27 the spreading of the positions is smaller because of the lack of differential radial drift in these discrete species. Differential radial drift is therefor completely dominant in spreading the size distribution of particles. 

Open with DEXTER 
Acknowledgements
We thank the anonymous referee for their comments that helped improve the manuscript. NS was funded by the “Bottlenecks for particle growth in turbulent aerosols” grant from the Knut and Alice Wallenberg Foundation (2014.0048). AJ is grateful for support from the KAW Foundation (grant 2012.0150), the European Research Council (ERC Consolidator Grant 724687PLANETESYS) and the Swedish Research Council (grant 20145775). AJ and CCY thank the European Research Council (ERC Starting Grant 278675PEBBLE2PLANET). The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at LUNARC in Lund University.
Appendix A Gas velocity in large simulation box
In Sect. 4 we show the evolution of the particle scale height for run SI3038c in a simulation box of L_{x} = L_{z} = 0.8 H_{g}. Compared to the other panels in Fig. 17 that show the same run in smaller boxes (L_{x} = L_{z} = 0.2 H_{g} and L_{x} = L_{z} = 0.4 H_{g}), the equilibrium scale height of the corresponding species do not converge. The small species are diffused to larger heights than before and the larger ones settle closer to the midplane.
The radial component of the gas velocity may shed some light on the discrepancy. After about 60 orbits, u_{x} develops largescale, diagonal shocklike patterns that persist throughout the rest of the run. In Fig. A.1, we show the radial gas velocity in all three box sizes. The first two panels, which correspond to the smaller boxes, show the expected oscillatory motion due to stirring from the larger particles. The last panel, in addition to oscillatory motion, alsodisplays shocklike behavior. This behavior was not seen in the 3D singlespecies simulations with large boxes in Yang & Johansen (2014) and is likely an artifact of the numerical prescription. As recommended in Li et al. (2018), outflow instead of reflecting boundary conditions in the vertical direction might eliminate the shocklike patterns.
Fig. A.1 Radial gas velocity at t ≈ 1000 Ω^{−1} corresponding to runs SI3032c, SI3034c, and SI3038c, respectively. In addition to the oscillatory motion driven by the streaming instability, u_{x} develops largescale shocklike patterns in the largest box, shown in the final panel. Here, we overplot the outline of the two smaller box sizes for comparison. 

Open with DEXTER 
References
 Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Prog. Theor. Phys., 56, 1756 [NASA ADS] [CrossRef] [Google Scholar]
 ALMA Partnership (Brogan, C. L., et al.) 2015, ApJ, 808, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010a, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010b, ApJs, 190, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010c, ApJ, 722, L220 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R. 2018, MNRAS, 475, 5618 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brandenburg, A., & Dobler, W. 2002, Comput. Phys. Commun., 147, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531 [NASA ADS] [CrossRef] [Google Scholar]
 Dra̧żkowska, J., & Alibert, Y. 2017, A&A, 608, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dra̧żkowska, J., Windmark, F., & Dullemond, C. P. 2013, A&A, 556, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dra̧żkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dubrulle, B. 1995, Icarus, 114, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Gerin, M., Pety, J., Commerçon, B., et al. 2017, A&A, 606, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldreich, P., & LyndenBell, D. 1965, MNRAS, 130, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, J.F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984 [NASA ADS] [Google Scholar]
 Gorti, U., Hollenbach, D., Najita, J., & Pascucci, I. 2011, ApJ, 735, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., et al. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayashi, C. 1981, IAU Symp., 93, 113 [NASA ADS] [Google Scholar]
 Hezel, D. C., & Palme, H. 2010, Earth Planet. Sci. Lett., 294, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Ida, S., & Guillot, T. 2016, A&A, 596, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jacquet, E. 2014, Icarus, 232, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Johansen, A., Youdin, A., & Mac Low M.M. 2009, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, A&A, 537, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014a, in Protostars and Planets VI (Arizona: University of Arizona Press), 547 [Google Scholar]
 Johansen, A., Morbidelli, A., Gounelle, M., Jacquet, E., & Cuzzi, J. 2014b, in Asteroids, Comets, Meteors 2014, eds. K. Muinonen, et al., ACM Conf. Ser. [Google Scholar]
 Laibe, G., & Price, D. J. 2014, MNRAS, 444, 1940 [NASA ADS] [CrossRef] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, R., Youdin, A. N., & Simon, J. B. 2018, ApJ, 862, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Menu, J., van Boekel, R., Henning, T., et al. 2014, A&A, 564, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Rapson, V. A., Kastner, J. H., MillarBlanchaer, M. A., & Dong, R. 2015, ApJ, 815, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Schoonenberg,D., & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tanaka, H., Himeno, Y., & Ida, S. 2005, ApJ, 625, 414 [NASA ADS] [CrossRef] [Google Scholar]
 van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, D., & Wetherill, G. 1994, Icarus, 107, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, C.C., & Johansen, A. 2014, ApJ, 792, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, C.C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Youdin, A. N. 2010, EAS Pub.Ser., 41, 187 [CrossRef] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Publicly available at http://pencilcode.nordita.org/.
All Tables
All Figures
Fig. 1 Mean equilibrium radial drift velocity as a function of Stokes number (τ_{s}) for runs SI4142d and SI3042d measured from the time of saturation. The vertical bars show the 1σ variation of the radial velocity and are measures of the radial diffusion each species experiences. The red stars represent the single species equilibrium solutions. The green and blue circles show the multispecies radial equilibrium velocity near the midplane. The purple and red squares show the height averaged multispecies solution under the assumption that the particle density of each species is distributed uniformly within its scale height. Run SI3042d was shifted horizontally for illustration purposes. 

Open with DEXTER  
In the text 
Fig. 2 Theoretical equilibrium radial drift velocity of species τ_{s} = 10^{−3} in run SI3042d with respect to height. The drift velocity is highest in the two bins around the midplane and is close to the gas velocity in the absence of large particles further away from the midplane. 

Open with DEXTER  
In the text 
Fig. 3 Convergence of equilibrium radial velocity of different particle sizes against grid resolution with runs SI1042d, SI1042Bd and SI1042Cd. Both the mean drift velocities (circles) and the level of diffusion (vertical bars) agree between the cases of N_{x} × N_{z} = 128 × 128, 256 × 256, and 512× 512. This impliesthat our results are not dependent on numerical resolution. Runs SI1042Bd and SI1042Cd were shifted horizontally for better visibility. 

Open with DEXTER  
In the text 
Fig. 4 Convergence of equilibrium radial velocity of different particle sizes against particle number with runs SI4142d and SI4142Ad. The mean drift velocities and the 1σ variation agree well, indicating that the number of particles does not influence our results. In the case of run SI4142Ad, the symbols were shifted horizontally for comparison purposes. 

Open with DEXTER  
In the text 
Fig. 5 Mean equilibrium radial drift velocity as a function of τ_{s} for runs SI4132d and SI3032d measured from the time of saturation. The color scheme agrees with that of Fig. 1. Run SI3032d was shifted horizontally for illustration purposes. 

Open with DEXTER  
In the text 
Fig. 6 Radial drift velocity of runs SI3042c, SI4142c, SI3032c, and SI4132c in terms of Stokes number (τ_{s}) at t ≈ 1000 Ω^{−1}. Individual black points show the radial drift velocities of a random subset of approximately 14 000 particles. The particle size distribution is continuous, such that each particle has a unique size and drift velocity. Each red cross shows the mean drift velocity of two hundred particles with similar τ_{s}, which decreases with increasing particle size. The dispersion of points around this mean shows the level of diffusion the particles experience, which goes down with increasing particle size as well. The purple bars show the 1σ variation around the mean velocity. The blue circles show the radial equilibrium velocity near the midplane for each particle species grouped into 16 bins using the multispecies NSH solution in Eq. (10). The theoretical velocity values provide a relatively good fit to the mean drift velocity for small particle sizes but are somewhat off for large Stokes numbers. 

Open with DEXTER  
In the text 
Fig. 7 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τ_{s}) for runs SI3042d and SI3032d measured from the time when the saturated state is reached. The level of diffusion is similar in both runs. 

Open with DEXTER  
In the text 
Fig. 8 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τ_{s}) for runs SI4142d and SI4132d measured from the time when the saturated state is reached. The level of turbulence does not appear to be affected by the choice of the q parameter. 

Open with DEXTER  
In the text 
Fig. 9 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τ_{s}) for runs SI3042d and SI3042c measured from the time when the saturated state is reached. In the case of run SI3042d, the velocities are averaged to the end of the simulation, while in the case of run SI3042c, they are averaged over 15 orbits. 

Open with DEXTER  
In the text 
Fig. 10 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τ_{s}) for runs SI3032d and SI3032c measured from the time when the saturated state is reached. In the case of run SI3032d, the velocities are averaged to the end of the simulation, while in the case of run SI3032c, they are averaged over 15 orbits. 

Open with DEXTER  
In the text 
Fig. 11 Particle scale height evolution of runs SI3042d and SI3044d for about 300 orbits. Different particle sizes (marked with different colors) saturate to different scale heights, in such a way that H_{p} increases with decreasing Stokes number. Left panel: evolution of the run with the larger Stokes numbers (τ_{s} = 10^{−3}−10^{0}) in a 0.2 × 0.2 H_{g} simulation box. Right panel: convergence of the scale height for the larger species in a box that is twice the size of the original. The scaleheight is relatively similar for the two box sizes, except for the two smallest particle size bins that exhibit an increased scale height when the box size is enlarged. 

Open with DEXTER  
In the text 
Fig. 12 Particle scale height evolution of runs SI4142d, SI4144d, and SI4142Ad, where the participating particles are τ_{s} = 10^{−4}−10^{−1}. Left panel: simulation in a box of 0.2 × 0.2 H_{g}. Middle panel: evolution of the same system in a 0.4 × 0.4 H_{g} box, where the particle scale heights saturate to different values relative to the smaller box. Right panel: the total particle number was increased by a factor of four from the model in the leftmost panel. In runs SI4144d and SI4142Ad the particle subdisks do not saturate to a constant height within the simulation time. 

Open with DEXTER  
In the text 
Fig. 13 Particle scale height evolution of run SI3032d in a 0.2 × 0.2 H_{g} simulation box. Compared to the lefthand side of Fig. 11, the particle species saturate to similar scale heights but do not show fluctuations of the same magnitude with time. 

Open with DEXTER  
In the text 
Fig. 14 Particle scale height evolution of run SI4132d in a 0.2 × 0.2 H_{g} simulation box. The particle species do not reach a constant scale height within the simulation time. Compared to the plot on the lefthand side of Fig. 12, where the powerlaw index of the distribution q = 4, the particle layers extend to similar scale heights. 

Open with DEXTER  
In the text 
Fig. 15 Particle scale height evolution of run SI3042c and run SI3044c. Each particle in the simulation has a unique size within the given size range. Particles are grouped into six size bins for illustration purposes. To ensure that the diffusion of the small particles is not limited by the simulation box, we performed the same simulation in a 0.4 × 0.4H_{g} box in the right panel. The scale heights of the large particles converge while the smallest species saturate to a larger H_{p} than in the left panel. 

Open with DEXTER  
In the text 
Fig. 16 Particle scale height evolution of run SI4142c and run SI4144c. The scale height of all species does not saturate within the simulation time in either the 0.2 × 0.2 H_{g} (left panel), nor the 0.4 × 0.4 H_{g} box (right panel). 

Open with DEXTER  
In the text 
Fig. 17 Particle scale height evolution of run SI3032c, run SI303 4c, and run SI3038c. Panels (from left to right): evolution of the system in simulation boxes of 0.2 × 0.2 H_{g}, 0.4 × 0.4 H_{g}, and 0.8 × 0.8 H_{g}, respectively. The largest particles converge to the same height in the first two cases, but particles with τ_{s} < 10^{−1} are diffused to heights larger by about a factor of two. The panel on the right, where the simulation box size is the largest, shows unusual nonconvergent behavior due to the development of shocklike patterns in the radial component of the gas velocity (discussed in Appendix A). 

Open with DEXTER  
In the text 
Fig. 18 Particle scale height evolution of runs SI4132c and SI4134c. The particle scale height of the small species increases by a factor of two when changing the dimensions from 0.2 × 0.2 H_{g} in the left panel to 0.4 × 0.4 H_{g} in the right panel, while that of the largest particles converges to about the same height. 

Open with DEXTER  
In the text 
Fig. 19 Vertical gas velocity fluctuation with respect to height in the case of runs SI3038c, SI3034c, and SI3032c averaged over approximately 15 orbits starting at t ≈ 1000 Ω^{−1}. The relatively uniform level of fluctuation further away from the midplane implies that the turbulence is not limited to the region close to the midplane (z = 0), where the large particles reside. The turbulence extends to larger scale heights, thus affecting the small particle species that reside there. 

Open with DEXTER  
In the text 
Fig. 20 Horizontal gas velocity fluctuation with respect to height in the case of runs SI3038c, SI3034c, and SI3032c averaged over approximately 15 orbits starting at t ≈ 1000 Ω^{−1}. The level of fluctuation peaks around the midplane and is relatively uniform throughout the rest of the domain for the runs in the 0.4 × 0.4 H_{g} and 0.8 × 0.8 H_{g} boxes, while in the smallest box δu_{x,rms} varies by more than a factor of two. 

Open with DEXTER  
In the text 
Fig. 21 Vertical gas velocity fluctuation with respect to height in the case of runs SI4134c, and SI4132c averaged over approximately 15 orbits starting at t ≈ 500 Ω^{−1}. The level of fluctuation is relatively uniform for the two runs both in terms of height and box size. 

Open with DEXTER  
In the text 
Fig. 22 Horizontal gas velocity fluctuation with respect to height in the case of runs SI4134c and SI4132c averaged over approximately 15 orbits starting at t ≈ 500 Ω^{−1}. The level offluctuation peaks around the midplane for both box sizes and is relatively uniform throughout the rest of the domain. 

Open with DEXTER  
In the text 
Fig. 23 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI4142d and SI3042d. Here, D_{x} is approximately constant for small Stokes numbers and increases with increasing particle size as τ_{s} ≥ 10^{−2.5} in the case of run SI4142d. The vertical diffusion coefficient increases with particle size until τ_{s} ≈ 10^{−1.5}. Above this value it decreases as the particles are less coupled to the gas. Both the radial and vertical diffusion coefficients are higher in the case of run SI3042d, since this model is more abundant in solids with τ_{s} ≳ 10^{−2} which are responsible for triggering the streaming instability. The diffusion coefficients corresponding to species that saturate to scale heights similar to the maximum allowed scale height in the box are marked with a cross. 

Open with DEXTER  
In the text 
Fig. 24 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI4144d and SI3044d. We doubled the box size to 0.4 × 0.4 H_{g} compared to the simulations in Fig. 23. Compared to the smaller box, both D_{x} and D_{z} converge to approximately the same values for each Stokes number. The exception seen in curve D_{z,30} at intermediate Stokes numbers is likely a consequence of the fluctuations seen in their particle scale height in the smaller box (see Fig. 11). The diffusion coefficients corresponding to species that do not saturate within the simulation time are marked with a cross. 

Open with DEXTER  
In the text 
Fig. 25 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI4142d and SI4142bd. Both D_{x} (solid curves) and D_{z} (dashed curves) show relatively small deviation. 

Open with DEXTER  
In the text 
Fig. 26 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI4132d and SI3032d. Both D_{x} and D_{z} take on similar values as in Fig. 23, where the particle size distribution is shallow (q = 4). The diffusion coefficients corresponding to species that saturate to scale heights similar to the maximumallowed scale height in the box or do not saturate within the simulation time are indicated with a cross. 

Open with DEXTER  
In the text 
Fig. 27 Evolution of particle radial positions with respect to time in the case of run SI3032c. Since the particles are distributed continuously in terms of their size, we grouped them into six particle size bins. The mean position of the particles in a given size range is marked by the solid curves, while the dashed curves mark the spread of particle positions by differential radial drift and diffusion. Particles with τ_{s} ≥ 10^{−1} drift and diffuse inwards, while particles of τ_{s} ≤ 10^{−1.5} tend to move outwards in the disk. The spreading of each particle species is in fact dominated by the differential radial drift insteadof diffusion. 

Open with DEXTER  
In the text 
Fig. 28 Evolution of particle radial positions with respect to time in the case of run SI3032d. The smallest species drift outwards, while species with sizes τ_{s} > 10^{−2} drift inwards. Compared to the run with continuous size distribution shown in Fig. 27 the spreading of the positions is smaller because of the lack of differential radial drift in these discrete species. Differential radial drift is therefor completely dominant in spreading the size distribution of particles. 

Open with DEXTER  
In the text 
Fig. A.1 Radial gas velocity at t ≈ 1000 Ω^{−1} corresponding to runs SI3032c, SI3034c, and SI3038c, respectively. In addition to the oscillatory motion driven by the streaming instability, u_{x} develops largescale shocklike patterns in the largest box, shown in the final panel. Here, we overplot the outline of the two smaller box sizes for comparison. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.