© ESO 2018

1 Introduction

Protoplanets – the building blocks of terrestrial planets, as well as of the cores of gas giants, ice giants and super-Earths – 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 millimeter-centimeter 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 back-reaction of the solids onto the gas, a small over-density 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 pile-up 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 millimeter-sized 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 near-infrared up to centimeter wavelengths, that suggest the presence of both micron-sized dust and millimeter–centimeter-sized pebbles. van Boekel et al. (2017) presented scattered light observations of the TW Hya disk with the Spectro-Polari-meter High-contrast 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 two-dimensional, in the radial-vertical plane, since this is sufficient for capturing the linear and non-linear 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 self-gravity 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 & Lynden-Bell 1965). The equations that govern the motion of the gas are the continuity and momentum equations, respectively $$\frac{\partial {\rho }_{\text{g}}}{\partial t}+u\cdot \nabla {\rho }_{\text{g}}-\frac{3}{2}\Omega x\frac{\partial {\rho }_{\text{g}}}{\partial y}=-{\rho }_{\text{g}}\nabla \cdot u,$$(1) $$\begin{array}{ll}\frac{\partial u}{\partial t}+(u\cdot \nabla )u\hfill & -\frac{3}{2}\Omega x\frac{\partial u}{\partial y}=2\Omega u_y\widehat{x}-\frac{1}{2}\Omega u_x\widehat{y}\hfill \\ \hfill & -\frac{\nabla P}{{\rho }_{\text{g}}}+2\eta {\Omega }^{2}r\widehat{x}-{\Omega }^{2}z\widehat{z}-\frac{ϵ}{{\tau }_{\text{f}}}(u-v).\hfill \end{array}$$(2)

The third term on the left-hand 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 right-hand 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 large-scale pressure gradient in the disk is represented by the fourth term, where $$\eta =-\frac{1}{2}{(\frac{H_{\text{g}}}{r})}^2\frac{\partial \text{ln }P}{\partial \text{ln }r},$$(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 second-to-last term on the right-hand 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 super-particles, each of which represents a swarm of identical physical particles. The governing equations for each super-particle are $$\frac{\text{d}{x}_{\text{p}}}{\text{dt}}=-\frac{3}{2}\Omega x_{\text{p}}\widehat{y}+v,$$(4) $$\frac{\text{d}v}{\text{dt}}=(2\Omega v_y\widehat{x}-\frac{1}{2}\Omega v_x\widehat{y}-{\Omega }^2{z}_{\text{p}}\widehat{z})+\frac{u-v}{{\tau }_{\text{f}}}.$$(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 $${\tau }_{\text{f}}=\frac{{\rho }_{•}a}{{\rho }_{\text{g}}{c}_{\text{s}}}.$$(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 $${\tau }_{\text{s}}=\sqrt{2\pi }\frac{{\rho }_{•}a}{{\Sigma }_{\text{g}}}.$$(7)

We write the gas surface density as $${\Sigma }_{\text{g}}(r)=1000{\text{ g cm}}^{-2}{(\frac{r}{\text{au}})}^{-1}{f}_{\text{g}},$$(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 $$a\approx 400\text{ cm }{(\frac{r}{\text{au}})}^{-1}{(\frac{{\rho }_{•}}{1{\text{ g cm}}^{-3}})}^{-1}{f}_{\text{g}}{\tau }_{\text{s}}.$$(9)

Assuming that f_g = 1, Eq. (9) shows that with ρ_∙ = 1 g cm⁻³ at 2.5 au τ_s = 10⁻⁴ 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 ${\gamma }_x={(v_{1x},v_{2x},\mathrm{...},v_{Nx})}^{\text{T}}$ and ${\gamma }_y={(v_{1y},v_{2y},\mathrm{...},v_{Ny})}^{\text{T}}$, satisfy the matrix equation $$\left(\begin{array}{c}{\gamma }_x\\ {\gamma }_y\end{array}\right)=\eta v_{\text{K}}\left(\begin{array}{cc}A& 2B\\ -B/2& A\end{array}\right)\left(\begin{array}{c}0\\ 1\end{array}\right),$$(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 + Γ)B and $B={\{{[{\Lambda }^{-1}(1+\Gamma )]}^2+1\}}^{-1}{\Lambda }^{-1}$. Both are functions of the Stokes number τ_i and the dust-to-gas ratio ϵ_i for each species i, through Λ = diag{τ₁, τ₂,..., τ_N} and $\Gamma ={(ϵ,ϵ,\mathrm{...},ϵ)}^{\text{T}}$, where $ϵ={({ϵ}_1,{ϵ}_2,\mathrm{...},{ϵ}_N)}^{\text{T}}$. It is useful to non-dimensionalize 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$$u=-{\sum }_j{ϵ}_j{v}_j-\eta v_{\text{K}}\$\widehat{y}\$.$$(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 sub-Keplerian velocity, such that $u(t=0)=-\eta v_{\text{K}}\widehat{y}$ 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 particle-gas systems, we used the Pencil Code¹. It is a high-order finite difference code used for (magneto–)hydrodynamic flows in astrophysics and it uses third-order 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 particle-mesh 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 SIMN-A-B-p 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 self-similar 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 particle-gas 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 SI41-4-2A-d, we increased the number of particles by a factor of four from the default value. In runs SI10-4-2B-d and SI10-4-2C-d the number of cells was increased to 256 × 256 and 512 × 512, respectively.

The solid-to-gas ratio, Z, is another important factor in the outcome of the dynamics. This parameter is defined as Z = Σ_p∕Σ_g, where $${\Sigma }_{\text{p},\text{g}}={\displaystyle {\int }_{-\infty }^{\infty }{\rho }_{\text{p},\text{g}}}(r,z)\text{d}z.$$(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.

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 sub-Keplerian 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 self-generated 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 dust-to-gas ratio of each size bin as a function of height using the relationship $${ϵ}_j(z)=\frac{{\rho }_{j,0}\exp (-z^2/2H_{\text{p}}^2)}{{\rho }_{\text{g},\text{0}}\exp (-z^2/2H_{\text{g}}^2)}$$(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, $${\rho }_{j,0}={\rho }_{\text{g},\text{0}}{f}_{j}Z\frac{H_{\text{g}}}{H_{\text{p}}}$$(14)

is the total mass density of particles of species j in the midplane and the mass fraction f_j is $$f_j=\frac{{\tau }_{\text{s},\text{j}}^{4-q}}{{\sum }_k{\tau }_{\text{s},\text{k}}^{4-q}}.$$(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 $$\frac{\text{d}M}{\text{dln }a}\propto a^{4-q},$$(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 SI41-4-2-d and SI30-4-2-d. 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 SI41-4-2-d. In the case of run SI30-4-2-d the level of diffusion is lower for the largest particles τ_s = 10^−0.5 and τ_s = 10⁰ 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) $$v_x=-\frac{2{\tau }_{\text{s}}}{{(1+ϵ)}^2+{\tau }_{\text{s}}^2}\eta v_{\text{K}}.$$(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⁻² 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 dust-to-gas ratio of each species calculated under the assumption that the vertical particle distribution is approximately Gaussian. The midplane density for the particles is $${\rho }_{\text{p}}(z=0)=\frac{{\Sigma }_{\text{p}}}{\sqrt{2\pi }{H}_{\text{p}}},$$(18)

where Σ_p is the column mass density of the particles. Thus, the local dust-to-gas ratio for the jth particle bin at the midplane is $${ϵ}_j=\frac{{\Sigma }_{\text{p},\text{j}}}{\sqrt{2\pi }{\overline{H}}_{\text{p},\text{j}}}\frac{1}{{\rho }_{\text{g},\text{0}}},$$(19)

where ${\overline{H}}_{\text{p},\text{j}}$ 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 dust-to-gas ratio in a given size bin using the global dust-to-gas ratio and the mass fraction per bin (see Eq. 15), such that $${ϵ}_j=\frac{Z{f}_j{H}_{\text{g}}}{{\overline{H}}_{\text{p},\text{j}}}.$$(20)

As seen in Fig. 1, the theoretical radial equilibrium velocity near the midplane agrees well in the case of run SI41-4-2-d (green circles), but is noticeably more positive than our measured drift velocities in the case of run SI30-4-2-d (blue circles). This is likely to be the result of the fact that we used the midplane dust-to-gas 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 height-averaged 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 dust-to-gas ratio for species j is given by $$\begin{array}{l}{ϵ}_j(|z|<{\overline{H}}_{\text{p},\text{j}})=\frac{\sqrt{2\pi }Z{f}_j{H}_{\text{g}}}{2{\overline{H}}_{\text{p},\text{j}}},\\ {ϵ}_j(|z|>{\overline{H}}_{\text{p},\text{j}})=0,\end{array}$$(21)

where ${\overline{H}}_{\text{p},\text{j}}$ is the time-averaged scale height of the particle species. Figure 2 shows the equlibrium velocity of the τ_s = 10⁻³ species as a function of height in run SI30-4-2-d. 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 dust-to-gas ratio in each bin (Bai & Stone 2010a). As seen in Fig. 1 (red and purples squares), the height-averaged 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 SI10-4-2-d, SI10-4-2B-d and SI10-4-2C-d. 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 SI41-4-2-d and SI30-4-2-d 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 SI30-4-2-d was shifted horizontally for illustration purposes.

Fig. 2 Theoretical equilibrium radial drift velocity of species τs = 10−3 in run SI30-4-2-d 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.

Fig. 3 Convergence of equilibrium radial velocity of different particle sizes against grid resolution with runs SI10-4-2-d, SI10-4-2B-d and SI10-4-2C-d. Both the mean drift velocities (circles) and the level of diffusion (vertical bars) agree between the cases of Nx × Nz = 128 × 128, 256 × 256, and 512× 512. This impliesthat our results are not dependent on numerical resolution. Runs SI10-4-2B-d and SI10-4-2C-d were shifted horizontally for better visibility.

Fig. 4 Convergence of equilibrium radial velocity of different particle sizes against particle number with runs SI41-4-2-d and SI41-4-2A-d. 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 SI41-4-2A-d, the symbols were shifted horizontally for comparison purposes.

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 SI30-3-2-d and SI41-3-2-d. Similar to the previous case, where all size bins had the same solid mass (runs SI30-4-2-d and SI41-4-2-d) 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 SI41-3-2-d, but deviate in the case of run SI30-3-2-d. The reason for the deviation is, again, that using the midplane dust-to-gas 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 height-averaged 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 dust-to-gas 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 SI41-3-2-d and SI30-3-2-d measured from the time of saturation. The color scheme agrees with that of Fig. 1. Run SI30-3-2-d was shifted horizontally for illustration purposes.

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 sub-Keplerian 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 $$m_{\text{p},j}=\frac{M_{\text{p},j}}{N_{\text{p},j}}=\frac{Z{\rho }_0\sqrt{2\pi }{H}_{\text{g}}{L}_x{f}_j}{N_{\text{p},j}},$$(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 dust-to-gas ratio in the given size bin calculated under the assumption that the vertical particle distribution is approximately Gaussian. The dust-to-gas 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 SI30-4-2-c.

Figure 6b shows the case of the smaller particle sizes, namely run SI41-4-2-c. 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⁻¹ 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 SI30-4-2-c (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 SI30-3-2-c. 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 SI41-3-2-c. 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 SI30-3-2-c.

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 SI30-4-2-c, SI41-4-2-c, SI30-3-2-c, and SI41-3-2-c 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.

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 SI30-4-2-d and SI30-3-2-d. 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 SI41-4-2-d and SI41-3-2-dand 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 SI30-4-2-d and SI30-4-2-c. 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 SI30-3-2-d and SI30-3-2-c. 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, ${\sigma }_z^2$, thus $$H_{\text{p}}=\sqrt{{\sigma }_z^2}=\sqrt{\langle z_{\text{p}}^2\rangle -{\langle z_{\text{p}}\rangle }^2}.$$(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⁻³ − 10⁰, 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⁻⁴−10⁻¹. 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) $$t_{\text{settle}}=\frac{1}{{\tau }_{\text{s}}\Omega }.$$(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 right-most 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 power-law index of q = 3. In other words, the system now has more mass in larger particle species. The evolution of both SI30-3-2-d and SI41-3-2-d 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 SI30-4-2-d and SI30-3-2-d measured from the time when the saturated state is reached. The level of diffusion is similar in both runs.

Fig. 8 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τs) for runs SI41-4-2-d and SI41-3-2-d 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.

Fig. 9 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τs) for runs SI30-4-2-d and SI30-4-2-c measured from the time when the saturated state is reached. In the case of run SI30-4-2-d, the velocities are averaged to the end of the simulation, while in the case of run SI30-4-2-c, they are averaged over 15 orbits.

Fig. 10 Mean equilibrium radial drift velocity and the 1σ variation of the radial velocity as a function of Stokes number (τs) for runs SI30-3-2-d and SI30-3-2-c measured from the time when the saturated state is reached. In the case of run SI30-3-2-d, the velocities are averaged to the end of the simulation, while in the case of run SI30-3-2-c, they are averaged over 15 orbits.

Fig. 11 Particle scale height evolution of runs SI30-4-2-d and SI30-4-4-d for about 300 orbits. Different particle sizes (marked with different colors) saturate to different scale heights, in such a way that Hp increases with decreasing Stokes number. Left panel: evolution of the run with the larger Stokes numbers (τs = 10−3−100) in a 0.2 × 0.2 Hg 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 scale-height 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.

Fig. 12 Particle scale height evolution of runs SI41-4-2-d, SI41-4-4-d, and SI41-4-2A-d, where the participating particles are τs = 10−4−10−1. Left panel: simulation in a box of 0.2 × 0.2 Hg. Middle panel: evolution of the same system in a 0.4 × 0.4 Hg 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 left-most panel. In runs SI41-4-4-d and SI41-4-2A-d the particle subdisks do not saturate to a constant height within the simulation time.

Fig. 13 Particle scale height evolution of run SI30-3-2-d in a 0.2 × 0.2 Hg simulation box. Compared to the left-hand side of Fig. 11, the particle species saturate to similar scale heights but do not show fluctuations of the same magnitude with time.

Fig. 14 Particle scale height evolution of run SI41-3-2-d in a 0.2 × 0.2 Hg simulation box. The particle species do not reach a constant scale height within the simulation time. Compared to the plot on the left-hand side of Fig. 12, where the power-law index of the distribution q = 4, the particle layers extend to similar scale heights.

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 SI30-4-2-c and SI41-4-2-c. 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 self-generated 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⁻³−10⁰, 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⁻⁴−10⁻¹. 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 SI30-3-2-c, SI30-3-4-c and SI30-3-8-c, 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 SI30-3-4-c 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 shock-like 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 SI30-4-2-c and run SI30-4-4-c. 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.4Hg box in the right panel. The scale heights of the large particles converge while the smallest species saturate to a larger Hp than in the left panel.

Fig. 16 Particle scale height evolution of run SI41-4-2-c and run SI41-4-4-c. The scale height of all species does not saturate within the simulation time in either the 0.2 × 0.2 Hg (left panel), nor the 0.4 × 0.4 Hg box (right panel).

Fig. 17 Particle scale height evolution of run SI30-3-2-c, run SI30-3- 4-c, and run SI30-3-8-c. Panels (from left to right): evolution of the system in simulation boxes of 0.2 × 0.2 Hg, 0.4 × 0.4 Hg, and 0.8 × 0.8 Hg, 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 non-convergent behavior due to the development of shock-like patterns in the radial component of the gas velocity (discussed in Appendix A).

Fig. 18 Particle scale height evolution of runs SI41-3-2-c and SI41-3-4-c. The particle scale height of the small species increases by a factor of two when changing the dimensions from 0.2 × 0.2 Hg in the left panel to 0.4 × 0.4 Hg in the right panel, while that of the largest particles converges to about the same height.

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 self-regulated diffusion driven by the larger species. Compared to single-species simulations, there is also a significant difference in terms of the scale height each particle species saturates to. Carrera et al. (2015) conducted single-species streaming instability simulations of a variety of Stokes numbers. They found that particles of 10⁻³ ≤ τ_s ≤ 10⁰ saturate to H_p ≈ 10⁻² 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 root-mean-square of the vertical gas velocity in terms of height for the runs in three different box sizes, namely run SI30-3-2-c, SI30-3-4-c, and SI30-3-8-c. 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 SI30-3-4-c and SI30-3-8-c 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 SI30-3-2-c, 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 SI41-3-2-c and SI41-3-4-c. 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 SI30-3-8-c, SI30-3-4-c, and SI30-3-2-c 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.

Fig. 20 Horizontal gas velocity fluctuation with respect to height in the case of runs SI30-3-8-c, SI30-3-4-c, and SI30-3-2-c 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 Hg and 0.8 × 0.8 Hg boxes, while in the smallest box δux,rms varies by more than a factor of two.

Fig. 21 Vertical gas velocity fluctuation with respect to height in the case of runs SI41-3-4-c, and SI41-3-2-c 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.

Fig. 22 Horizontal gas velocity fluctuation with respect to height in the case of runs SI41-3-4-c and SI41-3-2-c 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.

5 Radial and vertical diffusion

The mutual drag between the solid and gas components triggers the streaming instability and generates turbulence. This self-regulated 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 self-regulated turbulence. It is calculated following Johansen & Youdin (2007) as $$D_x=\frac{1}{2}\frac{\text{d}{\sigma }_{\text{x}}^{\text{2}}}{\text{dt}}.$$(25)

Here, ${\sigma }_x^2$ 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) $$\frac{H_{\text{p}}}{H_{\text{g}}}=\sqrt{\frac{\delta }{{\tau }_{\text{s}}+\delta }}.$$(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 power-law 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 SI41-4-2-d (blue dashed curve) and SI30-4-2-d (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 SI30-4-2-d and SI41-4-2-d is relatively constant for τ_s ≤ 10⁻², 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 SI41-4-2-d (blue curves) and SI30-4-2-d (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⁻², 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⁻² < τ_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 SI41-4-2-d. 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 SI41-4-2-d and SI30-4-2-d. Here, Dx is approximately constant for small Stokes numbers and increases with increasing particle size as τs ≥ 10−2.5 in the case of run SI41-4-2-d. 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 SI30-4-2-d, 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.

Fig. 24 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI41-4-4-d and SI30-4-4-d. We doubled the box size to 0.4 × 0.4 Hg compared to the simulations in Fig. 23. Compared to the smaller box, both Dx and Dz converge to approximately the same values for each Stokes number. The exception seen in curve Dz,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.

Fig. 25 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI41-4-2-d and SI41-4-2b-d. Both Dx (solid curves) and Dz (dashed curves) show relatively small deviation.

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 SI41-3-2-d and SI30-3-2-d. 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 SI41-3-2-d and SI30-3-2-d. Both Dx and Dz 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 maximum-allowed scale height in the box or do not saturate within the simulation time are indicated with a cross.

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 self-regulated 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 SI30-3-2-c. 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 SI30-3-2-d. Particles with τ_s > 10⁻² drift inwards, while species with τ_s ≤ 10⁻² 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 chondrule-matrix 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 fluid-particle 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⁻⁴−10⁻¹, τ_s = 10⁻³−10⁰, or τ_s = 10⁻¹−10⁰. 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 sub-Keplerian 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 root-mean-square 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 self-generated 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 SI30-3-2-c. 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.

Fig. 28 Evolution of particle radial positions with respect to time in the case of run SI30-3-2-d. 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.

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 724687-PLANETESYS) and the Swedish Research Council (grant 2014-5775). AJ and CCY thank the European Research Council (ERC Starting Grant 278675-PEBBLE2PLANET). 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 SI30-3-8-c 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 shock-like 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 shock-like behavior. This behavior was not seen in the 3D single-species 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 shock-like patterns.

Fig. A.1 Radial gas velocity at t ≈ 1000 Ω−1 corresponding to runs SI30-3-2-c, SI30-3-4-c, and SI30-3-8-c, respectively. In addition to the oscillatory motion driven by the streaming instability, ux develops large-scale shock-like patterns in the largest box, shown in the final panel. Here, we overplot the outline of the two smaller box sizes for comparison.

References

All Tables

All Figures

Fig. 1 Mean equilibrium radial drift velocity as a function of Stokes number (τs) for runs SI41-4-2-d and SI30-4-2-d 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 SI30-4-2-d was shifted horizontally for illustration purposes.

In the text

	Fig. 2 Theoretical equilibrium radial drift velocity of species τs = 10−3 in run SI30-4-2-d 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.
In the text

Fig. 3 Convergence of equilibrium radial velocity of different particle sizes against grid resolution with runs SI10-4-2-d, SI10-4-2B-d and SI10-4-2C-d. Both the mean drift velocities (circles) and the level of diffusion (vertical bars) agree between the cases of Nx × Nz = 128 × 128, 256 × 256, and 512× 512. This impliesthat our results are not dependent on numerical resolution. Runs SI10-4-2B-d and SI10-4-2C-d were shifted horizontally for better visibility.

In the text

	Fig. 4 Convergence of equilibrium radial velocity of different particle sizes against particle number with runs SI41-4-2-d and SI41-4-2A-d. 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 SI41-4-2A-d, the symbols were shifted horizontally for comparison purposes.
In the text

	Fig. 5 Mean equilibrium radial drift velocity as a function of τs for runs SI41-3-2-d and SI30-3-2-d measured from the time of saturation. The color scheme agrees with that of Fig. 1. Run SI30-3-2-d was shifted horizontally for illustration purposes.
In the text

Fig. 6 Radial drift velocity of runs SI30-4-2-c, SI41-4-2-c, SI30-3-2-c, and SI41-3-2-c 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.

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 SI30-4-2-d and SI30-3-2-d measured from the time when the saturated state is reached. The level of diffusion is similar in both runs.
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 SI41-4-2-d and SI41-3-2-d 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.
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 SI30-4-2-d and SI30-4-2-c measured from the time when the saturated state is reached. In the case of run SI30-4-2-d, the velocities are averaged to the end of the simulation, while in the case of run SI30-4-2-c, they are averaged over 15 orbits.
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 SI30-3-2-d and SI30-3-2-c measured from the time when the saturated state is reached. In the case of run SI30-3-2-d, the velocities are averaged to the end of the simulation, while in the case of run SI30-3-2-c, they are averaged over 15 orbits.
In the text

Fig. 11 Particle scale height evolution of runs SI30-4-2-d and SI30-4-4-d for about 300 orbits. Different particle sizes (marked with different colors) saturate to different scale heights, in such a way that Hp increases with decreasing Stokes number. Left panel: evolution of the run with the larger Stokes numbers (τs = 10−3−100) in a 0.2 × 0.2 Hg 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 scale-height 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.

In the text

Fig. 12 Particle scale height evolution of runs SI41-4-2-d, SI41-4-4-d, and SI41-4-2A-d, where the participating particles are τs = 10−4−10−1. Left panel: simulation in a box of 0.2 × 0.2 Hg. Middle panel: evolution of the same system in a 0.4 × 0.4 Hg 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 left-most panel. In runs SI41-4-4-d and SI41-4-2A-d the particle subdisks do not saturate to a constant height within the simulation time.

In the text

	Fig. 13 Particle scale height evolution of run SI30-3-2-d in a 0.2 × 0.2 Hg simulation box. Compared to the left-hand side of Fig. 11, the particle species saturate to similar scale heights but do not show fluctuations of the same magnitude with time.
In the text

	Fig. 14 Particle scale height evolution of run SI41-3-2-d in a 0.2 × 0.2 Hg simulation box. The particle species do not reach a constant scale height within the simulation time. Compared to the plot on the left-hand side of Fig. 12, where the power-law index of the distribution q = 4, the particle layers extend to similar scale heights.
In the text

Fig. 15 Particle scale height evolution of run SI30-4-2-c and run SI30-4-4-c. 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.4Hg box in the right panel. The scale heights of the large particles converge while the smallest species saturate to a larger Hp than in the left panel.

In the text

	Fig. 16 Particle scale height evolution of run SI41-4-2-c and run SI41-4-4-c. The scale height of all species does not saturate within the simulation time in either the 0.2 × 0.2 Hg (left panel), nor the 0.4 × 0.4 Hg box (right panel).
In the text

Fig. 17 Particle scale height evolution of run SI30-3-2-c, run SI30-3- 4-c, and run SI30-3-8-c. Panels (from left to right): evolution of the system in simulation boxes of 0.2 × 0.2 Hg, 0.4 × 0.4 Hg, and 0.8 × 0.8 Hg, 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 non-convergent behavior due to the development of shock-like patterns in the radial component of the gas velocity (discussed in Appendix A).

In the text

	Fig. 18 Particle scale height evolution of runs SI41-3-2-c and SI41-3-4-c. The particle scale height of the small species increases by a factor of two when changing the dimensions from 0.2 × 0.2 Hg in the left panel to 0.4 × 0.4 Hg in the right panel, while that of the largest particles converges to about the same height.
In the text

Fig. 19 Vertical gas velocity fluctuation with respect to height in the case of runs SI30-3-8-c, SI30-3-4-c, and SI30-3-2-c 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.

In the text

Fig. 20 Horizontal gas velocity fluctuation with respect to height in the case of runs SI30-3-8-c, SI30-3-4-c, and SI30-3-2-c 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 Hg and 0.8 × 0.8 Hg boxes, while in the smallest box δux,rms varies by more than a factor of two.

In the text

	Fig. 21 Vertical gas velocity fluctuation with respect to height in the case of runs SI41-3-4-c, and SI41-3-2-c 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.
In the text

	Fig. 22 Horizontal gas velocity fluctuation with respect to height in the case of runs SI41-3-4-c and SI41-3-2-c 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.
In the text

Fig. 23 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI41-4-2-d and SI30-4-2-d. Here, Dx is approximately constant for small Stokes numbers and increases with increasing particle size as τs ≥ 10−2.5 in the case of run SI41-4-2-d. 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 SI30-4-2-d, 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.

In the text

Fig. 24 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI41-4-4-d and SI30-4-4-d. We doubled the box size to 0.4 × 0.4 Hg compared to the simulations in Fig. 23. Compared to the smaller box, both Dx and Dz converge to approximately the same values for each Stokes number. The exception seen in curve Dz,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.

In the text

	Fig. 25 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI41-4-2-d and SI41-4-2b-d. Both Dx (solid curves) and Dz (dashed curves) show relatively small deviation.
In the text

Fig. 26 Radial and vertical diffusion coefficients with respect to Stokes number for runs SI41-3-2-d and SI30-3-2-d. Both Dx and Dz 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 maximum-allowed scale height in the box or do not saturate within the simulation time are indicated with a cross.

In the text

Fig. 27 Evolution of particle radial positions with respect to time in the case of run SI30-3-2-c. 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.

In the text

Fig. 28 Evolution of particle radial positions with respect to time in the case of run SI30-3-2-d. 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.

In the text

	Fig. A.1 Radial gas velocity at t ≈ 1000 Ω−1 corresponding to runs SI30-3-2-c, SI30-3-4-c, and SI30-3-8-c, respectively. In addition to the oscillatory motion driven by the streaming instability, ux develops large-scale shock-like patterns in the largest box, shown in the final panel. Here, we overplot the outline of the two smaller box sizes for comparison.
In the text